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

We use a linear shallow-water model to investigate the global circulation of the atmospheres of tidally locked planets. Simulations, observations, and simple models show that if these planets are sufficiently rapidly rotating, their atmospheres have an eastward equatorial jet and a hot-spot east of the substellar point. We linearize the shallow-water model about this eastward flow and its associated height perturbation. The forced solutions of this system show that the shear flow explains the form of the global circulation, particularly the hot-spot shift and the positions of the cold standing waves on the night-side. We suggest that the eastward hot-spot shift seen in observations and 3D simulations of these atmospheres is caused by the zonal flow Doppler-shifting the stationary wave response eastwards, summed with the height perturbation from the flow itself. This differs from other studies which explained the hot-spot shift as pure advection of heat from air flowing eastward from the substellar point, or as equatorial waves travelling eastwards. We compare our solutions to simulations in our climate model Exo-FMS, and show that the height fields and wind patterns match. We discuss how planetary properties affect the global circulation, and how they change observables such as the hot-spot shift or day-night contrast. We conclude that the wave-mean flow interaction between the stationary planetary waves and the equatorial jet is a vital part of the equilibrium circulation on tidally locked planets.


INTRODUCTION
Tidally locked planets always present the same face to the star they orbit, so only receive starlight on one side. If such a planet has an atmosphere, its circulation transports heat from the permanent day-side to the permanent night-side. Simulations of such atmospheres show they have an eastward jet on the equator and a hot-spot to the east of their substellar point, rather than at the substellar point as might be expected. They also show a pair of cold low pressure lobes on the night-side, peaked in the mid-latitudes and symmetrically disposed about the equator. In this paper, we use a shallow-water model linearized about the equatorial jet to investigate this circulation, which has previously been interpreted to be composed of stationary planetary-scale waves (Showman & Polvani 2011) (Tsai et al. 2014). We show that the We built on these previous studies with a shallowwater model linearized about both a non-uniform equatorial jetŪ (y) and its associated height perturbation H(y). We used a pseudo-spectral method to solve these equations on an equatorial beta-plane and on a sphere, to find the stationary wave response to a day-night forcing. Our solutions agreed with the interpretation of Tsai et al. (2014) that the hot-spot shift is caused by the eastward jet Doppler-shifting the stationary waves eastwards.
We introduce the linearized shallow-water equations in Section 3. In Section 4 we use a pseudo-spectral method to solve the shallow-water equations linearized aboutŪ (y) andH(y) on an equatorial beta-plane (Boyd 2000) (Boyd 2017). We find the free modes of the system and its response to stationary forcing. This lets the response pattern be understood in terms of the spectrum of waves excited by the forcing. We explain how the jet shifts the wave response eastwards and changes the pattern of planetary standing waves. Our results match our simulations much better than previous linearized shallow-water calculations in zero or uniform background flow.
In Section 5 we solve the shallow-water equations linearized aboutŪ (φ) andH(φ) in a spherical coordinate system, which exactly represents the latitudinal direction, unlike the beta-plane. We compare this to our solutions on the beta-plane, and show that the solutions are similar despite the approximations of the beta-plane Section 6 shows how observables such as the hot-spot shift and day-night contrast depend on the properties of the planet and atmosphere. We derive simple onedimensional scaling relations on the equator, and compare them to the effects of changing the parameters of the full two-dimensional shallow-water model. Section 7 compares our linear shallow-water solutions to simulations in our General Circulation Model (GCM) Exo-FMS. We decompose its output to isolate the Rossby and Kelvin components for comparison with the linear shallow-water model. Our results show how the standing waves are Doppler-shifted as the zonal jet forms during the spin-up from rest, and match the linear shallow-water model when the GCM reaches equilibrium.
We conclude that linearizing about the eastward jet is vital to the shallow-water model of a tidally locked atmosphere. Our model gives a new interpretation of the main processes controlling the circulation on such planets.

THE ATMOSPHERIC CIRCULATION OF TIDALLY LOCKED PLANETS
Planets such as gas giants in short-period orbits ("Hot Jupiters"), close-in terrestrial planets around G-stars ("lava planets"), and a range of terrestrial planets around M-dwarfs are all expected to be tidally locked. The atmospheres of these planets are very different to those in the Solar System. Their permanent dayside receives much more radiation than the permanent night-side, which drives a strong global circulation. Articles such as Showman et al. (2012), Heng & Showman (2015), and Pierrehumbert & Hammond (2018) review the atmospheric circulation of tidally locked planets.
Shallow-water models (Showman & Polvani 2011) Tsai et al. (2014, GCM simulations (Showman et al. 2012), and observations (Louden & Wheatley 2015) (Crossfield 2015) (Parmentier & Crossfield 2017) suggest that tidally locked atmospheres have a superrotating eastward equatorial jet and an eastward hot-spot shift. Showman & Polvani (2011) set out a theory of the global circulation of atmospheres of tidally locked planets. They focused on the initial forcing and planetary waves which formed an equatorial jet, and did not consider the effect of the jet on the waves after it had formed. Tsai et al. (2014) discussed the effect of a uniform mean flowŪ 0 on the equatorial waves, and identified how the jet Doppler-shifts the waves eastwards (Arnold et al. 2012). They suggested that the hot-spot shift is caused by this Doppler-shift, rather than free wave propagation or advection of air from the substellar point. We build on their work by linearizing the shallow-water equations about a shear flowŪ (y) and a height pertur-bationH(y), which are in geostrophic balance or gradient wind balance. Boyd (1978) and Wang & Xie (1996) solved the linear shallow-water equations for planetary waves in shear flows, and showed how the shear affects the latitudinal extent of the waves. The shear flow in our linear solutions affects both the latitudinal and longitudinal structure of the forced response. We reach the same conclusions as Tsai et al. (2014) about the mechanism of the hot-spot shift, but our linear model matches our GCM results more closely, as the sheared flow and the jet height perturbation have a large effect on the response to stationary forcing. Figure 1 summarises the problem that we aim to solve in this study. The first plot shows typical GCM results of the mean height (analogous to temperature here) for a tidally locked Earth-sized planet. The global circulation is similar to that seen in Hot Jupiter simulations, e.g. . The second plot reproduces the linear shallow-water model of Showman & Polvani (2011). This linear model did not match equilibrated GCM results similar to those in Figure 1, but did match their GCM a short time after it was started from rest. Specifically, the model did not match the positions and wind directions of the hot (high height) and cold (low height) parts of the GCM results. The third plot reproduces the linear shallow-water model in Tsai et al. (2014) which showed how a uniform background flow U 0 could shift the maxima and minima of the forced response. This Doppler-shifted response matches the position of the height maximum (i.e. the hot-spot in the GCM). Our shallow-water model linearized about a shear flow U (y) and associated height perturbationH(y) matches the overall form of our equilibrated GCM results. We use the new model to predict how the global circulation depends on the planetary parameters like rotation rate or damping rate.
Understanding the mechanism behind the circulation on these planets is crucial to interpreting observations of them. Komacek & Showman (2016), Komacek et al. (2017), and Zhang & Showman (2017) developed a simplified model of the circulation on a tidally locked planet using a one-dimensional balance of advection and damping, to predict how the day-night temperature difference and the hot-spot shift scale with planetary parameters. They tested this against GCM results and observations, and found that their scaling worked well. The advective model of these studies has the potential to explain the equatorial hot-spot shift, but does not contain the wave dynamics required for the off-equator response. This offequator response can significantly affect the phase curve -for example, the coldest parts of the circulation at the level of the jet are often the night-side off-equator Rossby lobes. In this paper, we suggest that the wave dynamics are important, as in our linear model, the hotspot shift is caused by a Doppler-shift of the the wave response (see Section 4). In Section 6 we discuss how our wave-based approach leads to some of the same scaling relations identified in Zhang & Showman (2017), but with a different physical basis.
The global circulation on tidally locked planets has been shown to affect habitability and climate stability (Kite et al. 2011) (Yang et al. 2013) (Kopparapu et al. 2017). We hope to show that simple atmospheric models are key to understanding the processes driving the global circulation on these planets, and to interpreting observations of them.

LINEARIZED SHALLOW-WATER EQUATIONS IN ZERO FLOW ON THE BETA-PLANE
We used the linear shallow-water equations on a onelayer equatorial beta-plane to model the atmosphere of a tidally locked planet. These equations describe the motion of a single layer of fluid of constant density where the horizontal scale of its flow is much greater than the depth of the fluid. The linear form of these equations describe small perturbations to this layer (Vallis 2006). We model the atmosphere of a tidally locked planet with a similar shallow-water model to Showman & Polvani (2011). The model corresponds to an active upper layer following the single-layer shallow water equations, above a quiescent layer which can transport mass and momentum to and from the upper layer. The forcing due to stellar heating is represented by Q, a relaxation to the radiative equilibrium height field.
In this section, we derive the wave response to stationary forcing on the beta-plane (Matsuno 1966). The beta-plane system approximates the Coriolis parameter as linear, which is only accurate at low latitudes but leads to more intuitive and useful solutions than the full spherical geometry. We solve the equations in a spherical geometry in Section 5, and show that the betaplane approximation leads to very similar solutions, as in other studies of the atmospheres of tidally locked planets (Showman & Polvani 2011) (Heng & Workman 2014).
All the quantities are linearized as the sum of a zonally mean background value F (y) and a perturbation with the form f (y)e i(kxx−ωt) (unlike Matsuno (1966), who uses the less conventional form f (y)e i(kxx+ωt) ). Throughout this paper, we will use capital letters for mean zonal quantities such asŪ andH, and lower-case letters for perturbations to this background, such as u and h (unless otherwise specified, such as the forcing Q). The shallow-water equations for these perturbations with zero background flow are: where h is the height, c = √ gH is the gravity wave speed (Matsuno 1966), and there is no friction or damping. Non-dimensionalizing with time scale 1/cβ and length scale c/β (the equatorial Rossby radius of deformation L R ), and assuming all quantities have the form f (y)e i(kx−ωt) , the free equations are: 45°S   0°4   5°N   90°N   2250  2350  2450  2550  2650  2750  2850 Height (m) GCM simulation of the atmosphere of a tidally locked planet, showing the hot-spot shift in the mid-atmosphere height field.
The height field in the analytic linear model in Showman & Polvani (2011), which explained the eastward equatorial acceleration.
The height field in the linear model in Tsai et al. (2014) which showed how a uniform flowŪ 0 = 1.0 shifts the forced response. Figure 1. Summary of previous work on this topic, showing how previous linear shallow-water models (the second and third plots) did not fully explain the equilibrium circulation in simulations of such planets. Areas of high height are red, and areas of low height are blue. In the first plot, the substellar point is the white cross at 0 • latitude and 0 • longitude. In the second and third plots, the forcing Q is overlaid in white, with solid contours for positive forcing and dashed contours for negative forcing.

Free Solutions in Zero Flow
The dispersion relation for the free modes is: and the free modes ξ are: (4) where ψ n = e − 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:

Forced Solutions in Zero Flow
for constant non-dimensional forcing σ = (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: The boundary conditions are: Matsuno (1966) shows that a forced solution χ can be expressed as a sum of the free solutions ξ: 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 forcing magnitude Q 0 = 1 to match 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.

Formation and Effect of 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: where overbars are zonal means, primes are deviations from zonal means, and the "thickness-weighted zonal average of any quantity A" is A * ≡ hA/h (Showman & Polvani 2011). R u is the zonal component of the vertical momentum transport R: Showman & Polvani (2011) show why the tilted wave pattern (such as the first plot in our Figure 6) transports eastward zonal momentum towards the equator. The tilts result in velocities "such that u v < 0 in the northern hemisphere and u v > 0 in the southern hemisphere", so the term in u v 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 ∂Ū (y)/∂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 time scales (Matsuno 1966), the equatorial acceleration of 0.1 implies a dimensional acceleration of 0.01 ms −2 on an Earth-like planet or 0.1 ms −2 on a Hot Jupiter. GCM simulations show that the jets of order 100 ms −1 on Earthlike planets (see Section 7 and jets of order 1000 ms −1 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.
The net equatorial acceleration in Figure 2 is eastward, 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 . 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 eastward to a maximum of π/2 ahead of the forcing. The flow shifts the frequency ω of the free modes by k xŪ , which then affects the imaginary part of the coefficients projecting the free modes onto the forced solution: A positive eastward flowŪ shifts the wave response eastwards. The third plot in Figure 1 shows how a uniform flow shifts the maximum of the forced response towards +90 • longitude. For a significant shift, the frequency shift k xŪ must be on the order of ω m , i.e. the jet speedŪ must be on the order of ω m /k x . The frequencies are the solutions of ω 2 − k 2 x − k x /ω = 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Ū ≈ 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 ω m ∼ k xŪ and the Rossby mode is significantly shifted isŪ crit = (gH) 1/2 /4.
For an Earth-like planet, g ∼ 10 ms −2 and h e ∼ 5 km, soŪ crit = 50 ms −1 . GCM simulations show equilibrium jet speeds of order 100 ms −1 , so the Doppler shift should be significant. A Hot Jupiter has g ∼ 20 ms −2 and h e ∼ 200 km (Showman & Polvani 2011), soŪ crit = 500 ms −1 . GCMs and observations show equilibrium jet speeds of order 1000 ms −1 , so the Doppler shift should again be significant. In Section 4, we will linearize the shallow-water equations about a non-uniform jet.

LINEARIZED SHALLOW-WATER EQUATIONS IN SHEAR FLOW ON THE BETA-PLANE
In this section, we discuss the main results of this paper -our solutions to the shallow-water equations linearized around a zonally uniform shear flowŪ (y) and heightH(y) (Boyd 2017).H(y) andŪ (y) satisfy the second line of Equation 1, so are geostroph-ically balanced (H y (y) = −yŪ (y)). For our Gaussian jetŪ (y) = U 0 e −y 2 /2 , the height perturbation is H(y) = U 0 e −y 2 /2 . As before, we use a forcing magnitude Q 0 = 1 and equal radiative and dynamical damping rates α rad = α dyn = 0.2 (Matsuno 1966). We will discuss the effect of unequal radiative and dynamical damping rates in Section 4.4. As explained in Section 3.3, the forced response is significantly shifted eastwards by zonal flows of order unity, so we vary our non-dimensional velocity U 0 from 0 to 1.
The forcing magnitude Q 0 must be comparable to the velocity and height produced by the jet, which are both of order U 0 = 1, so we set Q 0 = 1. This large forcing results in a zonal jet acceleration of order 0.1 in Figure 2, which implies that the jet will accelerate about ten times faster than in our GCM. This is due to the large forcing required on the beta-plane discussed previously. In Section 5 we will show that a more realistic forcing magnitude in a spherical geometry produces a more physically reasonable acceleration.
The new background ofŪ (y) andH(y) changes the shallow-water equations in Section 3 to: The solutions have the form A(y)e i(kxx−ωt) . To consider the free (unforced) modes of the system, we set Q(y) = 0 and ∂/∂t = −iω, giving the free linear system: To find the stationary response to steady forcing (in this case, a day-night heating difference) we set Q(y) = Q 0 e −y 2 /2 (Matsuno 1966) and ∂/∂t = 0, giving the forced linear system: Appendix A explains how we solved both of these sets of equations using a pseudo-spectral method, by expanding the solutions in terms of parabolic cylinder functions. Appendix A.3 shows the accuracy of this method, demonstrating that the pseudo-spectral method identifies the exact solution in the case with no background jet, and that the solution with a background flow changes by less than 1 part in 10,000 for any modes past n = 30. The pseudo-spectral method finds N m solutions (the number of modes used in the calculation) for the eigenvalue equation governing the free modes. Many of these are spurious, but we can distinguish them from the physical modes by inspecting their eigenvalues. The pseudo-spectral method only produces a single solution for the forced linear system, which is simpler to interpret.
The mean shallow-water height H 0 (non-dimensionalized to unity above) is determined by the vertical modes excited by the vertical heating profile of the atmosphere. Tsai et al. (2014) showed that in tidally locked atmospheres, the vertical heating profile excites a continuum of vertical modes, leading to a continuum of horizontal modes. Most of the energy is contained in one of these vertical modes, so it is sufficient to only consider this vertical mode and associated horizontal mode with height H 0 .

Free Solutions in Shear Flow
In this section, we will 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δ(y)x . The phase shift δ(y) in a shear flow is equivalent to the phase shift (ω m − k xŪ ) 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Ū 0 , the free modes are linearly Doppler-shifted as discussed in Section 3.3 (Tsai et al. 2014). For a shear background flowŪ (y), we found the free modes of Equation 14 using the method described in Appendix A, 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Ū (y) = U 0 e −y 2 /2 with a variable magnitude U 0 . 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 towards a maximum shift of +90 • and contributing to the hot-spot shift. 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 hot-spot shift. Figure 3 shows that in this shear flowŪ (y), the n = 1 Rossby mode eigenvalue increases but does not become positive for U 0 = 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 hot-spot shift. So, the forced response and the hot-spot 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 hot-spot 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 substel-lar 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. Non-zero 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 suffi-ciently 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 Figure 5 and 4, resemble the free solutions with zero shear flow (Matsuno 1966), with their latitudinal structure slightly changed by a weak shear flow U = 0.1e −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 A.2 in Appendix A). Boyd (1978) discusses how a shear flow affects the meridional structure of these modes in more detail.

Forced Solutions in Shear Flow
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Ū (y) andH(y). Both plots were calculated with the method in Appendix A, 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.
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 hot-spot.
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 hot-spot 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 0 = 1, so does not shift past the substellar point. This means that in Figure 6 the hot-spot shift is smaller at high latitudes. The higher-order Rossby modes are shifted further east, contributing more to the hot-spot 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 hot-spot 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 nondimensionalized 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.

Jet Acceleration and Equilibrium
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 which accelerates the flow initially (Showman & Polvani 2011).
As the mean eastward zonal flowŪ (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 accelera-Forced solution without a background jet, calculated with the pseudo-spectral method but exactly the same as the analytic solution in the second plot in Figure 1.
Forced solution with a background jet -now, the high height (hot) region matches our GCM results, as do the the low height (cold) Rossby lobes. tion 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 non-zero. 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 towards +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 −ū * /τ 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).
With the length and time scales on the beta-plane discussed previously, a non-dimensional velocity of U 0 = 1 corresponds to a velocity of order 100ms −1 on an Earth-

Effect of Damping
In the previous section we assumed that the radiative damping α 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, MHD damping, or extra nonlinear terms (Heng & Workman 2014). In this section, we suggest that al-though 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 .

Effect of Shear Flow on Hot-Spot Location and Global Circulation
Understanding the mechanism behind the global circulation and eastward hot-spot shift observed on many tidally locked planets (Parmentier & Crossfield 2017) is vital to interpreting observations of such planets, particularly why the magnitude of the hot-spot shift varies. Showman & Polvani (2011) suggest that the hot-spot 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 nonperiodic 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 hot-spot shift on tidally locked planets. This model predicted a hot-spot shift on the equator, but did not represent wave dynamics directly so did not include their effect on the hot-spot 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Ū 0 , and suggested that the hot-spot 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 non-uniform flowŪ (y) and and a height perturbation 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 perturbationH(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 hot spot pattern. The hot-spot 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 hot-spot shift is caused by an eastward Doppler-shift of the stationary Kelvin and Rossby waves towards +90 • , rather than eastward propagating Kelvin waves or pure heat transport from air in the eastward flowing jet.
The analytic solutions in Figure 9 are approximately the first-order terms in the pseudo-spectral expansion in Appendix A. The pattern in Figure 6 is more complex than the sum of patterns in Figure 9 as the flowŪ (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 hot-spot shift. In their case (a), the maximum zonal speed is approximately 100 ms −1 , 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 −1 , 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 hot-spot in the top right panel of their Figure 8. If the cases are Forced response in zero flow with α dyn = 0.
Forced response in flowŪ (y) = e −y 2 /2 with α dyn = 0. Figure 8. The same plots as Figure 6, but with α dyn = 0. The first plot shows the forced wave response in zero background flow, where the Rossby mode now dominates. The second plot shows the forced response on a background jet withŪ (y) = e −y 2 /2 -despite the lack of dynamical damping and weak Kelvin response, the general form of the solution is the same as in Figure 6.
sufficiently strongly damped, then a zonal flow may still not produce a hot-spot shift. Equation 12 shows that the projection coefficient will be mostly real if α is much larger than k xŪ , so the maximum of the wave response will be in phase with the forcing (at the substellar point, with no hot-spot shift).

LINEARIZED SHALLOW-WATER EQUATIONS IN SHEAR FLOW ON A SPHERE
Linearizing the shallow-water equations about a jet on an equatorial beta-plane preserves the link to the intuitive, exact solutions of Matsuno (1966). However, the beta-plane approximation is very inaccurate at high latitudes, and is not well suited to represent the effect of parameters like rotation rate on the latitudinal extent of the waves and circulation. The y coordinate in the previous sections is non-dimensionalized to the equatorial Rossby radius of deformation.
Converting the y coordinate on the beta-plane to a fraction of planetary radius shows how the solution varies with latitude, which works well for planets where the equatorial Rossby radius is smaller than the planetary radius. This is why the previous beta-plane so-lutions resemble both the solutions in spherical coordinates in this section, and the GCM simulations below. However, this conversion leads to serious inaccuracies for planets with a larger equatorial Rossby radius than planetary radius, as the Matsuno (1966) beta-plane solution assumes that the scale of the heating (the size of the tidally locked planet in our case) is the same as the equatorial Rossby radius.
In this section, we solve the shallow-water equations linearized about an equatorial jet in spherical coordinates. This represents the latitudinal dependence of our solutions better, and lets us compare them directly to GCM simulations. We use the same pseudo-spectral method as before, with the Legendre polynomials as our basis set. We emphasize that the beta-plane solutions are still the most useful for understanding the global circulation, as they expand the forced solutions in term of the free modes of the shallow-water system in zero background flow (Matsuno 1966). The spherical solutions lose the link to these intuitive solutions, but do represent the Coriolis force properly, and let us relate our results to real planets more closely.
The non-dimensional linearized spherical shallowwater equations are: where θ is the longitude, φ is the latitude, m is the zonal wavenumber, F is the forcing, Ω is the angu-lar speed of the planet, and G = gH 0 /a 2 Ω 2 where (a) Height field of the exact forced solution (Matsuno 1966).  As before, the zonal flowŪ (φ) is balanced by a height fieldH(φ). SubstitutingŪ (φ) into the second line of Equation 16 gives a height field in gradient wind balance with the zonal flow: The background flow is set to be: The cos(φ) factor is required as the flow must be zero at the poles for Equation 17 to be valid. We assume the perturbations have the form f (φ) exp[i(mφ − ωt)]. We set ∂/∂t = 0 and insert the forcing Q, to find the response to stationary forcing: We use the same method as in Appendix A, with the Legendre polynomials rather than the parabolic cylinder functions (Wang et al. 2016). We solve the equations with y-coordinate µ = sin φ, in the domain −1 < µ < 1. This matches the domain of the Legendre polynomials and reduces most of the trigonometric terms and derivatives in the calculation to powers of µ. We rescale the height variable h in the calculation to avoid singularities at the poles, which we explain in Appendix A.4. Figure 10 shows the non-dimensional spherical solutions for the forced problem with no shear, the Dopplershifting of this solution, and the total forced solution in a shear flow. They have the same form as Figure  6, showing that the beta-plane approximation produces much the same results as the spherical geometry. It appears that the n = 1 Rossby mode is shifted past the substellar point in this case, unlike the beta-plane system discussed in the previous section. The beta-plane system could be made to match the spherical solution more closely by choosing different values for parameters like zonal flow speed or damping rates.
In the figures in this section, we used the same damping rates as before, α rad = α dyn = 0.2, which are also used in Matsuno (1966) and Showman & Polvani (2011). The time and length scales specified above require that the radius is a = 1 and the angular frequency is Ω = 1. We set the jet width φ 0 = π/3, but the exact value of this is not vital to the form of the solution (unless it is unphysically narrow, or so broad as to be uniform). An eastward shift of approximately 90 • requires U 0 ≈ 0.5, so we set this as our maximum flow speed.
We set the forcing Q 0 = 0.5 so that the background jet height and forced wave response had comparable magnitude, which is more illuminating than either one dominating the other. This value of Q 0 is comparable to those used in Showman & Polvani (2011), and we will show the effect of varying it later.
We will consider Earth-like tidally locked planets with rotation rates in the range 1 to 10 days, which have values of G = gH/a 2 Ω 2 varying from 0.3 to 30. For comparison, Hot Jupiters with rotation rates of 1 to 2 days have a value of G ∼ 0.1. In Figure 10 we set G = 1. Later we will show how varying all of these parameters changes the forced solution and often makes either the jet or the wave response dominate.

SCALING RELATIONS
In this section, we will simplify our linear shallowwater model to a 1D model of the most dominant terms on the equator. In Section 6.1, we use this 1D model to predict simple scaling relations between the planetary parameters (damping rate, jet speed), and the phase and amplitude of the equatorial height variation (i.e. the hot-spot shift and day-night contrast). In Section 6.2 we will show that these simple relations qualitatively describe the behaviour of our pseudo-spectral solutions to the full 2D linear equations.

Scaling Relations in 1D Model
The pseudo-spectral solutions to Equation 15 such as 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 hot-spot. 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 Then as before, we set ∂/∂t = 0 and ∂/∂x = ik x : So then the real part of the full solution h(x, y = 0) = h(y = 0)e ikxx is: (23) This corresponds to a sinusoidal height perturbation on the equator. The magnitude of this height perturbation is approximately: The position of the hot-spot x 0 is the maximum of this height perturbation, where ∂h/∂x = 0: Note that this is the same as the approximation to the hot-spot shift from Zhang & Showman (2017), λ s = tan −1 ( τ rad τ adv ), if we set τ rad = 1/α and τ adv = k x /U 0 . This is because their prediction for the hot-spot shift is calculated from Equation 41 in Zhang & Showman (2017): This is equivalent to our Equation 22, if we set T ≡ h, τ adv = k x /U 0 , λ ≡ x, and (T eq (λ) − T )/τ 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 hot-spot shift only matches the approximate solution to their equation.
This 1D model predicts a hot-spot shift varying from 0 • to 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 hot-spot 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 0 the 1D model predicts a hot-spot 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 hot-spot 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 hot-spot 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 Height field of the forced solution with U (φ) = 0.
Forced wave solution with a background flow, minus the flow itself.
considering the full disk-integrated phase curve, where the off-equator response is very different to the 1D oneequator model.

Scaling Relations in 2D Model
In the previous section, we used a 1D model to calculate simplified predictions of how observables such as hot-spot 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 0 and y 0 , the forcing strength Q 0 , 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 0 and Q 0 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 0 and the jet U 0 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.
Section 6.1 predicts that increasing the damping rate α should decrease the hot-spot 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.
The linear shallow-water equations in a spherical geometry described in Section 5 also have the additional parameter G = (c/aΩ) 2 = gH/a 2 Ω 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 WTG parameter, which Pierrehumbert & Ding (2016) used to characterise 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 >> 1 (a low rotation rate Ω) waves dominate the forced response, giving a large day-night contrast. For G << 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 (Koll & Abbot 2015)  ) (Carone et al. 2015). Q 0 = 2.5, a high forcing strength gives a strong response to forcing and a wave-dominated response.
Q 0 = 0.1, a low forcing strength gives a weak response to forcing and a jet-dominated response. Figure 11. Height fields of the response to forcing in spherical coordinates for different values of forcing strength Q0, which controls the relative magnitude of the response to forcing and the background jet height and velocity.
α rad = α dyn = 1, a high damping gives a weaker response to forcing which is more in phase with the forcing Q.
α rad = α dyn = 0.04, a low damping gives a stronger response to forcing which is shifted further eastwards. This 2D linear model does not give such a straightforward prediction of the hot-spot shift and day-night contrast as the 1D model, but does capture the important off-equator behaviour and the effect of the planetary waves. It would be ideal for comparison with observational methods which retrieve 2D maps of the planets (Rauscher et al. 2018).

COMPARISON OF SHALLOW-WATER MODEL AND GCM RESULTS
In this section, we discuss the results of tests comparing 3D simulations in our GCM Exo-FMS (Hammond & Pierrehumbert 2017) to the 2D pseudo-spectral solutions described in this paper. Exo-FMS is based on the GFDL Flexible Modelling System. Our tests had a 96 by 144 by 40 grid, a timestep of 100 seconds, and greygas radiation. The tests took at most 50 Earth days to reach their equilibrium circulation pattern and at most 300 days to reach a steady energy balance. Our "equilibrium" results in this section were an average from 500 to 1000 days. Our "spin-up" results were taken over the first 100 days.
We simulated Earth-sized planets with 1 bar atmospheres with the thermodynamical properties of nitrogen, longwave optical depth 1, shortwave optical depth 0, and the same incoming stellar radiation as on Earth and a 10 day orbital period unless otherwise specified. We chose a 10 day orbital period with these parameters for the default test as this appears to give a pattern where the jet height and velocity (the wavenumber-zero component of the mean GCM output) is comparable to the wavenumber-one wave response, similar to Figures  6 and 10. Very high or very low rotation rates lead to a less illuminating global circulation dominated either by a zonally homogeneous jet or a zonally inhomogeneous G = 1/30, for an orbital period of 8 hours (or a low gravity or scale height, or a high radius). G = 30, for an orbital period of 10 days (or a high gravity or scale height, or a low radius). Figure 13. Height fields for low and high G. A low value of G corresponds to a stronger jet height perturbation and a more zonally homogeneous global circulation with a smaller latitudinal extent, as discussed in Section 5. Temperature (K) Figure 14. Mean temperature and velocity fields from the GCM simulation discussed in Section 7, corresponding to the half-surface pressure level. The substellar point is the white cross at 0 • latitude and 0 • longitude. This simulation has similar parameters to the linear solution in Figure 13.
wave reponse. We will show the effect of varying the 10 day orbital period later. Figure 14 shows the mean temperature and wind pattern in our simulations of a planet with the parameters given above, and a 10 day orbital period. The temperature distribution and wind direction qualitatively match our linear shallow-water solutions. The maximum zonal wind speed is about 100ms −1 The GCM results differ from the shallow-water model in the higher temperature at the substellar point, which may be due to rising hot air from the substellar point, which our shallow-water model does not account for.
Some GCM tests with high rotation rates have mid-latitude instabilities resembling travelling planetary waves, or baroclinic or barotropic instabilities (Polichtchouk & Cho 2012)  ) (Pierrehumbert & Hammond 2018). The time-averaged circulation of most of these planets still resemble the stationary wave pattern in Figure 14. We suggest that the instabilities lead to equilibrated travelling waves which move through the stationary wave pattern without greatly changing it (as discussed in Section 4.1), which is why our linear shallow-water model matches the mean global circulation pattern.
The shallow-water solutions and our plots of the midatmosphere in the GCM do not represent the temperature of every level of the atmosphere. GCM simulations of terrestrial planets have a surface temperature which is closely coupled to the incoming stellar radiation. Hot Jupiters have a varying temperature distribution with height, with an upper atmosphere dominated by radiative balance, sometimes with a temperature inversion. The mid-atmosphere circulation, where the wave patterns show up most clearly, is still very important to understanding observations and climate. Phase curve observations of Hot Jupiters normally find a hot-spot shift (Crossfield 2015), showing that the observations are of a level of the atmosphere where the jet and waves are important. On terrestrial planets, the strong midatmosphere circulation does affect the surface temperature, cooling the day-side and warming the night-side. The global circulation may have an even stronger effect on the climate if clouds are present (Kopparapu et al. 2017).

Equilibrium Global Circulation
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 sec-ond plot in Figure 1, but is in phase with the third plot, which has been Doppler-shifted eastwards by a background flow.
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 (2018) 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).

Spin-up of Global Circulation
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 day 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 90ms −1 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.
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 towards +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), Haqq-Misra et al. (2018), and Pierrehumbert & Hammond (2018) varied the parameters of simulated tidally locked planets and identified changes in the global circulation.

Global Circulation Regimes
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 hot-spot shift will be larger. This gives a smaller day-night contrast and larger hot-spot 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 daynight 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 hot-spot 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 hot-spot 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 hot-spot 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 hot-spot shift, which makes sense as the hot air transported in their jet model will cool quickly on the night-side.   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 hot-spot 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 hot-spot 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 with-out an eastward jet. They may be limited by nonlinear constraints on geostrophic adjustment (Pierrehumbert & Hammond 2018), and the linear theory in this paper will not apply.

CONCLUSIONS
We modelled the global circulation of a tidally locked planet using a shallow-water model linearized about an eastward equatorial jet and its associated height perturbation. This built on previous studies which linearized about a state of rest, or a uniform flow with uniform height. We found that the non-uniform flow and height were crucial to forming the global circulation pattern seen in GCM simulations.  We solved these shallow-water equations using a pseudo-spectral method, first on the beta-plane for simplicity, and then on a sphere to properly represent the Coriolis parameter and the latitude coordinate. Both models showed that the shear flow shifted the wave response eastwards, explaining the hot-spot shift seen in simulations. They also showed that the non-uniform height perturbation of the jet was crucial to the overall form of the global circulation. We varied the parameters of the model such as stellar forcing and planetary rotation rate, to show how they affect the scaling of observables such as hot-spot shift and day-night contrast. We showed how these scaling relations match both a simpler one-dimensional model, and the results of simulations in a GCM.
In this study, we only considered planets that orbit their host star closely and are exactly tidally locked. Thermal tides could prevent such a planet from reaching an exactly tidally locked state. Penn & Vallis (2017) used a shallow-water model to investigate planets which are rotating slightly faster or slower than if they were tidally locked, and showed that this affects the position of the hot-spot. It will be important to consider other possible orbital and climate states for any observed exoplanets, rather than assuming tidal locking.
Further work could build on this linear model by adding the effect of vertical structure (Tsai et al. 2014) or vertical shear (Boyd 1978). Instabilities seen in the GCM simulations could be investigated further by recreating them in a shallow-water model and finding how they scale with planetary properties.
In summary, this paper described a linear shallowwater model of the atmospheric circulation of a tidally locked planet. It included the shear and height perturbation of the eastward equatorial jet which forms on such planets. It matches the results from 3D general circulation models better than previous shallow-water models, especially the position of their shifted hot-spot and night-side cold spots.
We thank the anonymous reviewer for their time and valuable comments, which greatly improved the paper. M.H. was supported by an S.T.F.C. studentship. This work also received support from European Research Council project EXOCONDENSE.
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. f is a vector made up of N subvectors f i , which are the forcing terms in each equation evaluated at each collocation point. H is the same as the matrix in Equation A4 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: u(y) = N n=0 a n φ n ; v(y) = This gives a linear matrix equation with one solution corresponding to the coefficient vectors a n , b n , c n of the forced solution.
Without forcing, the shallow-water equations define an eigensystem where the eigenvalue is the frequency ω.
The pseudo-spectral equation is then: R is an M × N square matrix with elements: 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.
This gives an eigenvalue matrix equation, with N eigenvalues and eigenvectors, corresponding to the frequencies and coefficient vectors a n , b n , 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 .

A.3. Beta-plane solutions
We use the parabolic cylinder functions ψ n (y) (Showman & Polvani 2011) as defined in Equation A15 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) (Boyd 2000).
Their collocation points are at their zeros (which are just the zeros of the Hermite polynomials H n ). Figure A.2 shows the first few parabolic cylinder functions. ψ n (y) = e −y 2 /2 H n (y) (A15) Figure 20 shows the magnitude of the coefficients (Equation A10) 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 non-zero. 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.

���� ������
Coefficients forŪ (y) = exp(−y 2 /2). Figure 20. Coefficients of the pseudo-spectral solution on the beta-plane coordinates with and without a background jet (the plots in Figure 6). The method identifies the exact solution in the first case, and converges rapidly to an accurate solution in the second case.
For non-zero 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.

A.4. Spherical solutions
We use the Legendre polynomials as a basis set for the pseudo-spectral method in a spherical geometry (Equation 16). Figure A.2 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 γ = h/ cos φ (Iga & Matsuda 2005). We replaced h with γ cos φ in Equation 16, solved as normal, then multiplied the solution for γ by cos φ 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.

���� ������
Coefficients calculated with rescaled height γ = h/ cos φ. Figure 21. Coefficients of the pseudo-spectral solution in spherical coordinates (the first plot in Figure 10), with the height variable h and the rescaled height γ. Rescaling the height makes the method converge to a smooth solution at the poles.