The Equatorial Jet Speed on Tidally Locked Planets: I -- Terrestrial Planets

The atmospheric circulation of tidally locked planets is dominated by a superrotating eastward equatorial jet. We develop a predictive theory for the formation of this jet, proposing a mechanism in which the three-dimensional stationary waves induced by the day-night forcing gradient produce an equatorial acceleration. This is balanced in equilibrium by an interaction between the resulting jet and the vertical motion of the atmosphere. The three-dimensional structure of the zonal acceleration is vital to this mechanism. We demonstrate this mechanism in a hierarchy of models. We calculate the three-dimensional stationary waves induced by the forcing on these planets, and show the vertical structure of the zonal acceleration produced by these waves, which we use to suggest a mechanism for how the jet forms. GCM simulations are used to confirm the equilibrium state predicted by this mechanism, where the acceleration from these waves is balanced by an interaction between the zonal-mean vertical velocity and the jet. We derive a simple model of this using the"Weak Temperature Gradient"approximation, which gives an estimate of the jet speed on a terrestrial tidally locked planet. We conclude that the proposed mechanism is a good description of the formation of an equatorial jet on a terrestrial tidally locked planet, and should be useful for interpreting observations and simulations of these planets. The mechanism requires assumptions such as a large equatorial Rossby radius and weak acceleration due to transient waves, and a different mechanism may produce the equatorial jets on gaseous tidally locked planets.


INTRODUCTION
Tidally locked planets always present the same face to the star that they orbit. Their atmospheric circulation is dominated by an equatorial jet, the strength of which determines directly observable features like the hot-spot shift and day-night contrast. Showman & Polvani (2011) showed that the jet is produced by the day-night instellation gradient, which induces stationary equatorial waves that transport prograde momentum towards the equator.
No studies so far have used this process to predict the equilibrium jet speed on these planets, as the process that balances this acceleration has not been identified. In this study we propose a mechanism by which this jet forms on terrestrial tidally locked planets, which does not rely on frictional drag. This provides a estimate of the jet speed that only depends on the basic atmospheric and planetary parameters. Our primary aim is to demonstrate the mechanism by which the jet forms, and to derive how the jet speed with the planetary parameters. Our estimate of an exact jet speed only applies to the idealised atmospheres we consider in this study, and will not apply to planets with thick atmospheres, significantly different heating profiles, or strong moisture effects.
The structure of this paper is as follows. In Section 2 we review previous work on the atmospheric circulation of tidally locked planets, and show their typical global circulation in a GCM simulation. In Section 3 we introduce an idealised model of the threedimensional stationary waves induced in the atmosphere of a tidally locked planets by its day-night instellation gradient. The mechanism we propose for the formation of the jet relies entirely on the zonal acceleration caused by these stationary waves, so we aim to isolate them in this idealised model. We solve the primitive equations on a beta-plane using the Dedalus software package, and show the structure of these waves and the zonal acceleration that they produce. The vertical profile of the zonal acceleration is then used in Section 4 to propose a mechanism for the formation and equilibration of the jet. Leovy (1987) and Zhu (2006) proposed similar mechanisms for the formations of zonal jets on Venus and Titan.
Section 6 uses a suite of GCM simulations of terrestrial tidally locked planets to test the proposed theory. We show that the equilibrium zonal momentum budget matches the expected balance from the proposed mechanism. We also show that the scaling of equatorial jet speed with instellation approximately matches the predicted speed from the stationary wave calculation. In Section 5 we derive a simple estimate of the maximum equatorial jet speed on a terrestrial tidally locked planet: for planetary radius a, acceleration due to gravity g, Stefan-Boltzmann constant σ, specific gas constant R, Brunt-Väisälä frequency N * , surface pressure p 0 , specific heat capacity c p , and instellation F 0 . This estimate corresponds to the jet speed at a height z = H, where H is the atmospheric scale height. The constant of proportionality depends on the heating profile at the substellar point, so the predicted jet speed is a scaling relation in general, and only a numerical prediction for planets with zero albedo and a forcing profile with a vertical wavelength ≈ 2H.
In Section 7 we discuss how this mechanism relies on several assumptions and simplifications, and suggest how other sources of acceleration such as transient waves could affect the jet speed. We also show how the different properties of gaseous "hot Jupiter" exoplanets could complicate the formation of a jet via this mechanism, which we will investigate in a forthcoming study.
We conclude that this mechanism is a good description of the formation of the equatorial jet on a terrestrial tidally locked planet with a dry, cloud-free atmosphere, and can predict the approximate jet speed for these planets. The key assumptions required by this mechanism are that the zonal acceleration is initially dominated by the contribution from stationary waves, and that once the jet forms it does not strongly affect the magnitude of this zonal acceleration. This mechanism could be extended to describe the jet formation and speed on planets with thicker atmospheres, clouds with strong radiative effects, or significant moisture content.

THE GLOBAL CIRCULATION OF TERRESTRIAL TIDALLY LOCKED PLANETS
Many of the best planetary candidates for atmospheric characterisation or potential habitability are expected to be tidally locked, but their atmospheric circulation is not fully understood. This section reviews previous work related to the equatorial jets of terrestrial tidally locked planets, and shows the typical features of a simulation of the global circulation on such a planet.

Review of Previous Work
The atmospheric circulation of tidally locked planets is measurable through observations such as thermal phase curves (Parmentier & Crossfield 2017). This study is motivated by the need to understand the formation of their equatorial jets, which are the dominant dynamical feature of this circulation (Pierrehumbert & Hammond 2019). Numerical atmospheric modelling has been invaluable to understanding the composition and climates of tidally locked gaseous planets such as hot Jupiters (Mayne et al. 2014;Showman et al. 2015;Drummond et al. 2018;Mendonça et al. 2018;Debras et al. 2020) and terrestrial planets (Joshi et al. 1997;Hammond & Pierrehumbert 2017;Boutle et al. 2017). The lack of observational data makes verifying the accuracy of these models difficult, compared to models of the Earth and the other planets of the Solar System. This has led to studies using different modelling approaches (Cho et al. 2015) and different approximations (Mayne et al. 2019) as it is not clear what techniques give the most realistic results. This study aims to provide a theoretical basis for the equatorial jet on these planets, which may help to guide modelling choices -for example, ensuring that an imposed surface drag or top-of-atmosphere sponge layer does not interfere with the momentum balance of the jet, or using a high enough upper boundary to fully resolve the zonal momentum fluxes that produce the jet.
Previous studies have demonstrated different circulation regimes on these planets, varying properties such as instellation and rotation rate to show their effect on the global circulation and its observable features (Kataria et al. 2014;Showman et al. 2015;Carone et al. 2015). Other studies have compared suites of simulations to observations such as phase curves of terrestrial planets (Demory et al. 2016;Hammond & Pierrehumbert 2017) and hot Jupiters (Arcangeli et al. 2019). More detailed measurements of circulation are becoming possible, such as measuring wind speeds via Doppler spectroscopy (Louden & Wheatley 2015;Brogi et al. 2016;Flowers et al. 2019), and measuring multi-dimensional temperature maps by eclipse mapping and spectral phase curves (Majeau et al. 2012;Stevenson et al. 2014).
An eastward superrotating equatorial jet is a common feature of almost all simulations of tidally locked planets. Read & Lebonnois (2018) define local superrotation as a state with a local excess of eastward atmospheric angular momentum relative to solid-body rotation at the equator. Any eastward flow at the equator is therefore superrotating and must be produced by up-gradient angular momentum transport towards the equator, a requirement know as Hide's Theorem (Hide 1969). This process has been investigated for Solar System planets such as Venus (e.g. (Fels & Lindzen 1974)) and the Earth (e.g. (Shell & Held 2004)). Showman & Polvani (2010) showed that on tidally locked planets this transport is provided by planetary-scale stationary waves, similar to the equatorial waves present in the tropics of the Earth (Matsuno 1966). In a pioneering study, Showman & Polvani (2011) developed this concept further and demonstrated it in a 2D linear shallowwater model, a 2D nonlinear shallow-water model and a 3D GCM. Their linear model relied on a linear drag to produce the appropriate stationary waves for an eastward equatorial acceleration (specifically, the equatorial Kelvin wave), and did not produce a zonal equatorial acceleration without this linear drag. Their non-linear model did produce a zonal equatorial acceleration without linear drag, which we will explore further in Section 3.2. Heng & Workman (2014) and Perez-Becker & Showman (2013) further explored two-dimensional models of this system, showing how non-linear balance or a more realistic forcing field gives different solutions, while still preserving the formation of the equatorial jet. This study focuses on the three-dimensional structure of the zonal acceleration, following authors such as Mendonça (2020) and Debras et al. (2020) who analysed the vertical structure of the zonal acceleration on tidally locked planets, identifying the key role of the vertical transport of zonal momentum. Tsai et al. (2014) followed the approach of Wu et al. (2000) to construct a three-dimensional linear model where the stationary wave response to stellar forcing is composed of separable vertical modes coupled to twodimensional shallow-water systems. This showed that a uniform eastward zonal flow shifts the stationary equatorial waves eastward, producing the equatorial hotspot shift seen in observations (Parmentier & Crossfield 2017). Tsai et al. (2014) also showed that on hot Jupiters the eastward shift of these equatorial waves can reduce the acceleration they produce, allowing the equatorial jet to reach a steady state. Hammond & Pierrehumbert (2018) used a two-dimensional linear shallowwater model to show how a zonal flow with meridional shear, and an associated geostrophically balanced geopotential perturbation, produces the characteristic shape of the atmospheric circulation and hot-spot shift on tidally locked planets.
The global circulation on a slowly rotating tidally locked planet has similarities to the tropical circulation on the Earth due to the large Rossby number in both cases (Pierrehumbert & Hammond 2019). This leads to behaviour that can be approximated by the "Weak Temperature Gradient" (WTG) regime (Pierrehumbert 2010a;Koll & Abbot 2016), where a non-linear balance in the zonal momentum equation leads to weak horizontal geopotential gradients (and therefore weak temperature gradients). This will be a key simplifying assumption later in this study. The tropics of the Earth also host similar equatorial stationary waves to those on tidally locked planets (Matsuno 1966;Gill 1980), which can lead to similar equatorial superrotation (Norton 2006). Lutsko (2018) considered a very similar system to this study, with a localised tropical heat source on the equator of the Earth rather than a planetary-scale day-night instellation gradient.
This study also builds on other work on the magnitude and scaling behaviour of the velocity and temperature fields on these planets. Komacek & Showman (2016) and Komacek et al. (2017) introduced a predictive theory for the scaling of temperature and velocity perturbations in the atmospheres of tidally locked planets, and successfully applied it to explain the observed scaling of day-night temperature differences on hot Jupiters. Koll & Abbot (2015) produced a similar theory for terrestrial planets based on the WTG approximation, which Kreidberg et al. (2019) used to interpret observations of the thermal phase curve of a terrestrial planet. Zhang & Showman (2017) derived scaling relations for properties of the atmospheric circulation of tidally locked planets, based on the relations of Komacek & Showman (2016). Koll & Komacek (2018) treated the global circulation of a tidally locked planet as a heat engine, to predict the wind speeds on hot Jupiters. Many of these theoretical predictions of observable quantities depend on the equatorial jet speed, which previously needed to be diagnosed from GCM simulations. This study aims to provide a predictive theory for this jet speed on terrestrial tidally locked planets, to enable the prediction of many other observable quantities. Figure 1 shows the typical global circulation on a terrestrial tidally locked planet, similar to the planets of the Trappist-1 system (Gillon et al. 2017). The simulation was run in the GCM Exo-FMS Hammond & Pierrehumbert 2017, 2018Pierrehumbert & Hammond 2019). Section 6.1 describes the details of the numerical modelling in this study in more depth; we introduce a single simulation here to show the key features of its circulation. Its main parameters are a radius of 6 × 10 6 km, a rotation period of 10 days, and instellation at the substellar point of 300 Wm −2 . Figure 1a shows the geopotential and velocity fields at the 500 mbar pressure level of this simulation, at the peak of the equatorial jet. Tsai et al. (2014) and Hammond & Pierrehumbert (2018) explained that the geopotential and velocity fields represent a "Matsuno-Gill" pattern of stationary equatorial waves (Matsuno 1966;Gill 1980), which are shifted eastwards by the equatorial jet. Hammond & Pierrehumbert (2018) showed how the velocity field is primarily a combination of a meridionally sheared (but zonally uniform) equatorial jet, plus a stationary wave response with zonal wavenumber 1. This results in a strong eastward velocity at −90 • where these two components combine, and a region of weak flow at 90 • where these two components cancel.

Typical Global Circulation
The eastward shift of the peak of the geopotential in Figure 1a corresponds to a shift in the peak of the temperature field eastwards from the substellar point. The maximum shift in the geopotential is at the level of the peak of the equatorial jet; the maximum shift in the temperature field is at a different pressure level as it is out of vertical phase with the geopotential field due to the hydrostatic relation. The shift of the hot-spot and the difference in temperature between the day-side and the night-side are observable quantities (Parmentier & Crossfield 2017;Komacek & Showman 2016), which should depend on the speed of this jet (Zhang & Showman 2017). The next section introduces the idealised model that we use to calculate the three-dimensional stationary wave response to a day-night instellation gradient.

THREE-DIMENSIONAL STATIONARY WAVE RESPONSE TO FORCING
In this section we calculate the three-dimensional stationary waves produced in the primitive equations on an equatorial beta-plane by an idealised forcing, representing the atmosphere of a tidally locked planet. Our aim is to show the structure of the zonal acceleration produced by these stationary waves, which we will use in the next section to suggest a mechanism for the formation and equilibration of the equatorial jet. Unlike previous studies, we find the response to forcing without imposing a linear drag on the horizontal velocities, which we will show to be vital to matching the magnitude of the velocity perturbations and zonal acceleration in our GCM simulations.

Idealised Beta-Plane Model
The adiabatic, inviscid primitive equations in height (log-pressure) coordinates (x, y, z) are (Vallis 2006): where the advective derivative is: and the horizontal derivative at constant z is: The variables in these equations are the temperature T , the horizontal velocity u = (u, v), the vertical velocity w, and the geopotential Φ. The log-pressure z coordinate can be transformed to the pressure coordinate p = p s e −z/H , and the vertical velocity w can be transformed to the vertical pressure velocity ω = −(p/H)w.
The Coriolis parameter is f = 2Ω sin φ, where Ω is the planetary rotation rate and φ is the latitude (which we will represent with the y coordinate on the betaplane, so that f = 2Ω sin φ = βy). The density is ρ = ρ 0 exp(−z/H), for a surface density ρ 0 (determined from the surface pressure, for an ideal gas) and a scale height that we approximate as H = RT 0 /g. R is the specific gas constant, g is the acceleration due to gravity, and T 0 is the planetary equilibrium temperature for the instellation F 0 . The Brunt-Väisälä frequency N * is  Figure 1. The equilibrium circulation of a simulation of a terrestrial tidally locked planet with instellation 300 W m −2 , discussed in Section 2.2. It is plotted at the 390 mbar level, for consistency with later figures. The substellar point is located at 0 • longitude. This is a typical circulation pattern for a tidally locked planet, with an eastward equatorial jet producing an eastward hot-spot shift. The eastward equatorial jet is centered at about 500 mbar. defined by N 2 * = R H (κT /H + ∂T /∂Z), where the dry adiabatic exponent κ = R/c p for the specific heat capacity c p . We approximate N 2 * to have constant value of 5 × 10 −4 s −2 , to approximate the value used in the tropics of the Earth by Wu et al. (2000), and to match the magnitude of the velocity perturbations in our GCM simulations in Section 6.
We recast these equations onto the beta-plane, add a forcing Q to represent that on a tidally locked planet, and introduce a linear radiative cooling term and a Rayleigh drag term (which we set to zero for the "nonlinear" solutions later). We impose a second-order divergence damping, second-order hyperdiffusion, and a sponge layer to stabilise the horizontal velocity fields (Jablonowski & Williamson 2011), which are all represented by the terms D u and D T . This results in the system: where α rad = 1/τ rad , for a constant radiative timescale τ rad : for an equilibrium temperature T 0 defined by F 0 = σT 4 0 . The factor e −2/3 p s corresponds to the pressure level at which the longwave optical depth is 2/3 in our GCM simulations, i.e. the radiating level of the outgoing longwave radiation (the structure of the solutions is not directly sensitive to this parameter). The dynamical damping rate α dyn is an arbitrary parameter corresponding to the rate of a linear Rayleigh drag, which we set to zero for most of the calculations. The second-order diffusive damping and second-order divergence damping applied to the velocity field has the form: and the second-order diffusive damping applied to the temperature field is: where the vertical distribution of the second-order diffusion terms is: to give a sponge layer that prevents the reflection of waves in the vertical direction from the upper boundary of the model. The sponge layer height is z s = 3H and the length scale over which the sponge layer increases to its maximum value is z w = H/3. We choose the coefficients to be K 2 = 10 7 and ν 2 = 3 × 10 7 to stabilise the calculations, without playing a major role in the momentum balance (Jablonowski & Williamson 2011).
The forcing Q is a three-dimensional field representing the heating applied to the atmosphere on a terrestrial tidally locked planet. We use an idealised vertical heating profile similar to that used by Wu et al. (2000) to represent heating by convective plumes in the tropics of the Earth, as the atmosphere in the simulations in this study is forced by absorption of longwave radiation from the surface and by dry convective adjustment: Gill & Philips (1986) uses the same profile to represent localised tropical heating on the Earth that excites waves with a vertical wavelength 2H. We found that the qualitative structure of the resulting stationary waves did not depend strongly on the exact structure of the vertical heating profile, as long as the forcing primarily excited waves with this vertical wavelength. The discontinuity in the first derivative of this vertical profile did excite high-order modes than smoother forcing profiles; we still used this profile as we found that the abrupt end to the forcing at a height H produced zonal acceleration profiles that matched those in the later GCM simulations well. In Tsai et al. (2014) a similar system of stationary waves is forced by a well-defined heating profile due to absorption of stellar radiation in the atmosphere of a hot Jupiter.
The forcing has the same horizontal distribution at all vertical levels, following the cosine of the longitude θ to represent the instellation on a tidally locked planet, and following a Gaussian envelope in the meridional direction: where the latitudinal scale of the forcing y 0 is set to be √ 2a, where a is the radius of the planet, in order to generate planetary-scale stationary waves similar to those in the GCM simulations and the shallow-water model of Showman & Polvani (2011).
We normalise the three-dimensional forcing field Q = Q x,y,z = Q 0 Q z (z)Q xy (x, y) such that the columnintegrated energy absorbed at the substellar point is equal to the instellation at the top of the atmosphere; i.e. that there is zero albedo and all of the instellation acts to heat the atmosphere. For an instellation F 0 , this requires: For a surface density ρ 0 = p s M/(RT 0 ) and a scale height H = RT /M g, this gives Q 0 ≈ 2.53gF0 pscp . The constant of proportionality depends on the vertical profile of heating at the substellar point, and the albedo of the planet, so in general this is a scaling relation Q 0 ∼ gF0 pscp rather than an exact equality.
We solve Equation 5 in a three-dimensional spectral basis using the Dedalus package (Burns et al. 2016(Burns et al. , 2020. This is an open-source Python package for the numerical solution of partial differential equations with spectral methods. We find the stationary solution of this forced system using different basis functions in each of its three dimensions. The x direction is represented with a Fourier basis, from which we exclude all zonally uniform modes with zonal wavenumber n = 0. This prevents the formation of a zonally uniform flow that affects the stationary wave response, as we are trying to isolate the acceleration due to stationary waves, without any background zonal flow. We will investigate the feedback of a zonal-mean flow on the stationary wave response in Section 4.3. The y direction is represented with a combination of sine and cosine functions (the "Sin/Cos" basis in Burns et al. (2020)), which allows for even or odd parity to be imposed on the variables in this direction. The u, w, Φ, and T fields are required to be symmetric about the equator, and the v field is required to be antisymmetric. The z direction is represented with a Chebyshev basis, where we impose the boundary conditions w = 0 at the top and bottom of the model. This system of equations is similar to the primitive equations solved in a GCM; however, the solution is conceptually closer to the results of shallow-water models such as those in Matsuno (1966), Showman & Polvani (2011), and to the results of the three-dimensional stationary wave model of Tsai et al. (2014). This system is on a beta-plane rather than a sphere, has strict parity conditions on all of its variables about the equator, is restricted to non-zero zonal wavenumbers, and is forced by a pre-defined stationary sinusoidal function. It therefore recovers only the three-dimensional stationary waves due to this forcing, rather than producing a full planetary circulation as a GCM would. The solution should be considered as "pseudo-equilibrium", as it satisfies the constrained set of equations we are solving, but in reality would produce a zonal-mean zonal flow which would modify the solution.

Stationary Wave Structure
Previous studies have found the stationary wave response and resulting zonal acceleration in systems with a linear drag applied to the velocity fields (Tsai et al. 2014). Showman & Polvani (2011) found that a lin-  . This level corresponds to the peak of the "Stationary Horizontal" acceleration in Figure 3. The solution in the first plot has a linear drag α dyn u for a drag timescale τ dyn = 1/α dyn = 1d, and the second plot has no linear drag. Away from the equator, where the momentum balance is governed by the Coriolis force, the two solutions are similar. Near the equator, the solutions are governed by linear drag and nonlinear balance respectively. This distinction is vital to accurately calculating the equatorial velocities and acceleration.
ear shallow-water system with no linear drag produced no zonal acceleration at the equator, as it could not produce the required Kelvin wave response there without drag. In this study, we calculate the solution with no linear drag in order to match our GCM simulations, and to derive an acceleration and equatorial jet speed that does not depend on an arbitrarily chosen linear drag timescale. We will show that our non-linear model can produce a zonal acceleration at the equator without linear dynamical drag. Showman & Polvani (2011) showed that this was possible in two-dimensional nonlinear shallow-water simulations, but did not examine the resulting stationary wave structure. Figure 2 shows the temperature and velocity fields of the equilibrated solutions to Equation 5 in the Dedalus software package, with and without a linear drag α dyn . Both of the solutions have a planetary radius 6 × 10 6 m, acceleration due to gravity 10 m s −2 , day (and year) length 10 d, instellation 10 3 W m −2 , specific heat capacity 10 3 J kg −1 K −1 , molar mass 28 g mol −1 , and surface pressure 10 5 Pa. The square of the Brunt-Väisälä frequency is set to a constant value N 2 * = 5 × 10 −4 s −2 . The magnitude of the velocity field in the solution is sensitive to this parameter, but it is difficult to estimate accurately for a atmosphere in general (unless the atmosphere is isothermal, which is not the case here). We suggest that this value is appropriate for these "Earth-like" atmospheres, as it gives velocities that approximately match the GCM simulations. It is important to note that the overall magnitude of the acceleration fluxes is sensitive to this parameter, but if the actual mechanism is not.
These parameters correspond approximately to the GCM test in Section 6 with instellation 10 3 W m −2 . Figure 2a shows the "linear" solution with a strong linear dynamical damping with a timescale of 1 d. While gaseous planets such as hot Jupiters may have be affected by magnetic drag that can be represented by this term (Komacek & Showman 2016), there is no reason to expect that the equatorial jet speed of terrestrial planets at the temperatures we consider will be affected be a uniform Rayleigh drag (apart from the drag near their surface, which will not affect the momentum budget at the level of the jet). We include the case with linear drag for comparison with the two-dimensional shallow-water solutions of Matsuno (1966) and Showman & Polvani (2011), as well as the three-dimensional stationary wave calculations of Tsai et al. (2014). The geopotential and velocity fields in Figure 2a are similar to those in Matsuno (1966) and Showman & Polvani (2011), with a stronger response on the equator than the solution in Matsuno (1966) due to the higher radiative and dynamical damping rates. Figure 2b shows the "non-linear" solution with the same parameters as the solution in Figure 2a, but with no dynamical damping. This system is governed by the same balances in the momentum and thermodynamic equations in Equation 5 as the GCM simulations later. Komacek & Showman (2016) discusses this non-linear (or "advective") momentum balance in detail, and compares it to other momentum balances via the Coriolis term and linear damping. The 'linear" and "non-linear" solutions are similar far away from the equator, where   Figure 2, which are the basis of the mechanism in Section 4. The qualitative forms of the profiles are the same, but their magnitudes are different as the zonal momentum equation is governed by different balances in each case. The local peaks of the acceleration profiles at about 0.4 bar correspond to the acceleration in two-dimensional shallow-water models such as in Showman & Polvani (2011). The peak of the "Stationary Vertical" term at about 0.6 bar is the main contribution to the unbalanced initial acceleration, which produces the equatorial jet.
the momentum balance is governed by the Coriolis term f × u.
In the linear case, the geopotential gradient ∇ z Φ is balanced by the linear drag α dyn u. This leads to a linear scaling between the velocity and the geopotential perturbation (and, via the hydrostatic balance equation and the thermodynamic equation, between the velocity and the forcing): This gives the linear relation between the velocity perturbations and forcing in Matsuno (1966) and Showman & Polvani (2011), which means that the zonal acceleration caused by these perturbations scales quadratically with the forcing. In the "non-linear" case this Rayleigh drag term is not present and the balance must be different. We discuss this non-linear balance in more detail in Section 5, and show how it leads to a much weaker dependence of the velocity perturbations and zonal acceleration on the magnitude of the stellar forcing. This approximation will be invalid when the WTG approximation does not apply globally on sufficiently rapidly rotating planets, or for sufficiently strong forcing that leads to non-negligible non-linear terms in the thermodynamic equation in Equation 5 (Pierrehumbert & Hammond 2019).
An important difference between the "linear" and "non-linear" cases is that in the linear case the magnitude of the velocity field is entirely dependent on the arbitrarily chosen α dyn damping parameter. This means that for a given forcing Q (and constant planetary parameters), the magnitude of the equatorial acceleration and the resulting jet speed will be entirely determined by the choice of the dynamical damping rate α dyn . Conversely, in the "non-linear" case, the magnitude of the equatorial velocities is determined only by the planetary parameters. The non-linear term plays the role of the linear dynamical damping, balancing the gradient of the geopotential Φ. The non-linear solution therefore has more predictive power than the linear solution. The linear case is still useful for emulating the behaviour of the non-linear case, and providing analytically soluble systems, as the damping parameter can represent the non-linear term if an appropriate value is chosen (for example, if α dyn u ∼ u∂ x u).
Both solutions are sensitive to the Brunt-Väisälä frequency N * . If the frequency is large enough that the w term is the dominant balance in the thermodynamic equation in Equation 5 (which is the case in the solutions plotted in this study), the qualitative form of the solutions should not depend on the exact magnitude of N * . The absolute magnitude of the vertical velocity perturbation (and by extension, the horizontal velocity perturbations) will depend on its magnitude. If the Brunt-Väisälä frequency is small enough that the w term is not the main term balancing the forcing in the thermodynamic equation (i.e. the WTG approximation does not apply), the stationary wave response to forcing may be very different. The Brunt-Väisälä frequency depends on the temperature structure of the atmosphere; it is simple to calculate for an isothermal atmosphere but we found that assuming an isothermal atmosphere did not produce stationary wave solutions that matched the GCM simulations, where the convective atmosphere (on the day-side) is much closer to neutral stability. The Brunt-Väisälä frequency is zero for a neutrally stable atmosphere, and it is difficult to estimate an accurate value to represent the entire atmosphere of a tidally locked planet.
We chose a constant value of N 2 * = 5 × 10 −4 s −2 to give solutions that approximately matched the magnitude of the stationary waves in the GCM, and to be consistent with the value used to represent the tropics of the Earth by Wu et al. (2000). The absolute value of the jet speed we predict is only as accurate as the value chosen for this frequency, although the qualitative mechanism is the same whatever value is chosen. An improved representation of this frequency could be a way to improve the accuracy of this model, highlighting the importance of understanding the generation of static stability in the atmospheres of tidally locked planets.

JET FORMATION MECHANISM
The three-dimensional stationary waves induced by the day-night instellation gradient produce a zonal acceleration at the equator (which would accelerate an equatorial jet, if we had not suppressed this in our idealised calculations). In this section, we calculate the vertical profiles of the different terms contributing to this acceleration. We will then propose a mechanism in which the jet forms due to this acceleration, then interacts with the vertical velocity to produce a deceleration that balances this acceleration and produces equilibrium. This mechanism is similar to that used to describe the formation of superrotating flows on Venus and Titan by Leovy (1987) and Zhu (2006).

Acceleration Profile Structure
The velocity fields in the solution to Equation 5 produce a zonal-mean zonal acceleration in spherical pressure coordinates (θ, φ, p) (Lutsko 2018): for the acceleration terms "Mean Horizontal" (MH), "Mean Vertical" (MV), "Stationary Horizontal" (SH), "Stationary Vertical" (SV), "Transient Horizontal" (TH), and "Transient Vertical" (TV). On the beta-plane (x, y, p) this is: where overbars denote zonal-mean values, and asterisks denote deviations from the zonal-mean value ("stationary" terms). All quantities such as u are timemeans, apart from the "Transient" terms where the primes denote deviations from the time-mean ("eddy" terms), and the square brackets denote a time-mean taken after the two eddy terms are multiplied together. Note that these terms are simplified compared to the equivalent terms in other studies such as Mayne et al. (2017); the "Mean" terms have been expanded and partially cancelled by combination with the zonallyaveraged continuity equation.
For forcing that is symmetric about the equator, v will be zero on the equator. We assume that the stationary forcing will produce stationary waves that are much larger than the transient waves, so the acceleration terms associated with the transient waves can be neglected (although we will still include these terms when diagnosing the momentum budget in the GCM simulations). These assumptions lead to the following simplified expression for the zonal-mean zonal acceleration on the equator: Figure 3 shows the terms in Equation 16, for the solution to Equation 5 in the linear limit and the nonlinear limit (which are the solutions shown in Figure 2). The "Mean Vertical" term is zero in both cases as the zonal-mean zonal velocity is zero. The linear and nonlinear cases have qualitatively similar vertical acceleration profiles, although the non-linear case has a larger acceleration, as the linear case must have smaller velocity perturbations from its additional dynamical damping. The "Stationary Horizontal" term corresponds to  a transport of eastward momentum towards the equator at about 0.4 bar, as shown by Showman & Polvani (2011) in a two-dimensional model corresponding to this level. The "Stationary Vertical" term corresponds to a transport of eastward momentum down from this level to about 0.6 bar, where it will accelerate the initial jet at the initialisation of the atmosphere from rest in the GCM. The structure of this acceleration is similar to that shown by Debras et al. (2020) for hot Jupiters, who identified the role of the "Stationary Vertical" term in accelerating the jet, rather than only decelerating it as in Showman & Polvani (2011). In the next section, we will describe a mechanism where this initial jet produces a new "Mean Vertical" acceleration term, which moves it up back towards the 0.4 bar level, and eventually balances the "Stationary" terms in equilibrium.

Predicted Jet Speed
The acceleration profiles caused by momentum transport from stationary waves in Figure 3 correspond to the zonal acceleration when there is no zonal-mean zonal flow. We suppose in this section that the jet does not strongly feed back on these stationary waves and affect the resulting acceleration, an assumption that we test in Section 4.3. Instead, the primary effect of the zonalmean jet u is to increase the "Mean Vertical" term in Equation 16, which is zero when the atmosphere is at rest. If this is the only feedback from the jet on the zonal momentum budget, the jet will accelerate until the "Mean Vertical" term exactly balances the sum of the "Stationary" terms, shown as the dashed black line in Figure 3.
Setting ∂u ∂t = 0 in Equation 16 shows that the zonalmean jet u required for the "Mean Vertical" term to balance the "Stationary" terms is: (17) Figure 4 is a schematic of how this mechanism works in practice. The first panel shows the evolution of the zonal-mean zonal flow to equilibrium, where line 3 is the equilibrium jet predicted by Equation 17, and lines 1 and 2 are example jet profiles chosen to represent the spin-up of the jet. The second panel shows the resulting "Mean Vertical" acceleration from each of these velocity profiles. As the jet speed increases over time, this term increases until it balances the sum of the "Stationary Horizontal" and "Stationary Vertical" terms in the zonal-mean momentum equation and equilibrium is reached. Line 3 and the dashed line do not cancel exactly above z = H due to the way we estimate ω, which we will discuss later.
The mechanism can also be understood qualitatively by considering the direction of zonal momentum transport. The stationary waves form in response to the dayside instellation and night-side cooling. These produce a "Stationary Horizontal" transport of eastward angular momentum towards the equator, with a peak around one scale height above the surface. This is opposed by the "Stationary Vertical" transport of this eastward momentum towards the surface. This accelerates a jet closer to the surface than the peak of the horizontal momentum transport. As this jet forms it interacts with the zonal-  Figure 5. Equilibrium zonal-mean zonal velocity profiles predicted by the "non-linear" (zero Rayleigh drag) calculation using the Dedalus software described in Section 3, according to the mechanism in Section 4. The first panel shows the profiles for different values of instellation, and the second panel shows the maximum value of each profile versus the instellation. The jet speed only depends weakly on instellation, due to the non-linear balance governing the magnitude of the velocity perturbations at the equator, its inverse dependence on the zonal-mean vertical velocity. mean vertical velocity associated with a Hadley-like circulation (which is not zero, as there is a zonal-mean component to the instellation, unlike in the idealised sinusoidal forcing in Showman & Polvani (2011)). This interaction produces a "Mean Vertical" accleration, as the jet is moved upwards by the zonal-mean vertical velocity. As the jet increases, this acceleration term balances the "Stationary" terms, until they are exactly balanced in equilibrium. If the "Stationary" terms exactly cancel at a certain pressure level, the peak of the equilibrium jet will form there, as the "Mean Vertical" term will be zero there as the jet has no vertical shear at that point.
The terms denoted by asterisks in Equation 17 are determined entirely by the calculation of the stationary wave response in Section 3. This calculation specifically excludes the zonal-mean quantities with zonal wavenumber 0, so does not determine the zonal-mean vertical velocity w. We instead estimate w to be the zonal mean of the day-side vertical velocity, as we expect there to be a Hadley-like overturning on the day-side due to instellation there. The sinusoidal variation of forcing with longitude in Section 3 implies an equal and opposite overturning on the night-side, giving no zonal-mean vertical velocity. A more realistic forcing field like that in Perez-Becker & Showman (2013) would be uniform on the night-side, and would only support overturning on the day-side, giving a non-zero zonal-mean vertical velocity. Our approximation of w in this stationary wave calculation is therefore: ω(y, p) = 1 π π/2 −π/2 ω(x, y, p)dx.
We only calculate u up to z = H using this approximation, as above this level ω can cross zero, giving an undefined solution to Equation 17. Figure 4b shows that neglecting u above z = H is a reasonable approximation, as the majority of the "Stationary" acceleration is below this level, so the portion of u above z = H is not important to the zonal-mean momentum budget. Figure 5 shows the zonal velocity profiles that are predicted by this mechanism using the stationary wave calculations in Dedalus, as discussed in Section 3. We vary the instellation (the magnitude of the forcing) but keep all other parameters the same. It is notable that the predicted zonal-mean jet only depends weakly on instellation. This is due to the nonlinear dependence of the velocity perturbations on the instellation, as well as the inverse relation between the jet speed and the zonalmean vertical velocity (which increases with increasing instellation) given by Equation 17. In reality, and in the GCM simulations, the evolution of the global circulation will be more complicated as the jet will affect the "Stationary" terms and be affected by the "Transient" terms, among other non-linear processes. In the next section, we investigate how strongly the jet affects the "Stationary" terms.

Jet Feedback on Stationary Acceleration Terms
This study proposes that the jet on terrestrial tidally locked planets reaches equilibrium when it becomes strong enough that the "Mean Vertical" acceleration term cancels the "Stationary" acceleration terms in Maximum SH Acceleration Fraction 10 2 Wm 2 10 3 Wm 2 10 4 Wm 2 (b) Fractional "Stationary Horizontal" term versus imposed jet Figure 6. The first panel shows the stationary wave response to forcing given an instellation 10 3 W m −2 and a background jet u = u0e −(y/a) 2 , where u0 = 30 m s −1 (similar to the equivalent test in the GCM). The geopotential and velocity fields appear very different to the solution in Figure 2b, which is the same apart from the imposed background jet. This case with a background jet has almost exactly the same zonal acceleration at the equator, as shown by the second panel, which plots the fractional change in the maximum "Stationary Horizontal" acceleration relative to the acceleration with zero zonal-mean background flow. The vertical dashed lines show the maximum jet speed in the corresponding GCM simulation.
waves produced by the instellation far enough eastward, that the "Stationary" acceleration terms decrease to zero. These acceleration terms require a phase difference between the on-equator stationary Kelvin waves and the off-equator stationary Rossby waves; if they are shifted towards the same position the acceleration will decrease (Tsai et al. 2014;Hammond & Pierrehumbert 2018).
So far, we have neglected this effect and assumed that the jet does not affect the stationary waves and the resulting acceleration strongly enough to affect the zonal momentum budget. In this section, we find the effect of an imposed background flow on the non-linear solutions to Equation 5, to see if this feedback is an important effect. In Section 3 the zonal mean of all variables was set to zero, but we now impose a non-zero zonal-mean background flow with a Gaussian profile in the meridional direction, similar to the background flow in Hammond & Pierrehumbert (2018): where a is the planetary radius. The exact meridional profile of this flow is not critical, as we focus on its effect at the equator only. The background flow must satisfy Equation 5, so we also impose a geopotential perturbation: We impose the same zonal flow at all vertical levels; in reality the jet would vary in the vertical direction. Figure 6a shows the geopotential and velocity fields for the non-linear solution with this background jet now imposed, with otherwise the same parameters as the solution with instellation 1000 W m −2 in Figure 2b. The imposition of the jet gives a visually different solution (with a structure explained in Hammond & Pierrehumbert (2018), but the magnitude of the zonal acceleration is almost the same in both cases. This shows that it can be difficult to estimate the zonal acceleration from the visual appearance of the stationary wave field. Figure 6b shows the effect of the zonal-mean zonal velocity on the maximum value of the "Stationary Horizontal" acceleration term, for a range of values of instellation. For the calculation with instellation 10 3 W m −2 , the "Stationary Horizontal" term is not strongly affected by the imposed background flow until the flow is greater than about 50 m s −1 . As the jet speed in the corresponding GCM simulation is approximately 38 m s −1 (as shown in Figure 6b), this feedback should not play a large role in the zonal momentum budget. The feedback from the jet is also small in the case with instellation 10 2 W m −2 , where the Dedalus calculations imply that the jet speed in the relevant GCM simulation is smaller than required to significantly affect the acceleration.
However, Figure 6b suggests that this feedback is relevant to the case with instellation 10 4 W m −2 . The calculation implies that the jet speed in this test in the GCM is strong enough to reduce the acceleration to approximately half its strength in the absence of a jet. The jet should still reach equilibrium via the proposed mechanism, but may have a lower maximum speed than predicted. Although this implies that the jet speed could be half that predicted by our estimate (where we neglect this feedback effect), in Section 6 we find that our estimate still matches the GCM simulation results reasonably well. This may mean that we have overestimated the effect of this feedback in this section, possibly due to our use of a vertically uniform background flow at the "jet speed". In reality, the flow only has this speed at one pressure level, and is significantly less at the peaks of the "Stationary" acceleration terms" (as we explain later, the peak jet speed does not necessarily coincide with either of the maxima of these acceleration terms). In addition, in these Dedalus calculations the imposed jet does not significantly affect the vertical structure of the acceleration terms beyond decreasing their magnitude. This is consistent with Tsai et al. (2014), where the stationary waves Doppler-shifted by an imposed jet retained the same vertical structure.
In summary, this section shows that for a terrestrial tidally locked planet with instellation 10 3 W m −2 the feedback of the jet on the stationary waves and the resulting "Stationary" acceleration terms is negligible for our simulations with instellation of 10 3 W m −2 and below. It may be significant for the simulations with instellation higher than this, but we find later that the effect may not be as great as predicted by our idealised calculations in this section, which apply the "peak" jet speed at all vertical levels and so may overestimate the effect of this feedback.

PREDICTING JET SPEED
In this section, we estimate the magnitude of the equatorial acceleration for given planetary parameters, and use this to derive the jet speed that will balance this acceleration. In Section 3, we predicted the equilibrium jet speed in Equation 17 by estimating the zonal flow that would produce a "Mean Vertical" acceleration that would balance the positive "Stationary Vertical" acceleration peak shown in Figure 3. In this section, we approximate Equation 17 as: This required us to approximate that the numerator (the "Stationary" terms) and the denominator (the zonal-mean vertical velocity) of the integrand have a similar vertical profile. This is supported by Figure 3b, where the sum of the "Stationary" terms has a similar vertical profile to the forcing Q used in Section 3. This means that the vertically varying parts of the numerator and denominator approximately cancel, and the integral in Equation 17 becomes a multiplication by H.
To estimate the speed of the jet, we find the magnitude of the velocity perturbations u * , v * , w * and the zonal-mean vertical velocity w, for a given forcing Q. The behaviour of the system is partly governed by the dominant balance in the thermodynamic equation in Equation 5. In this study, the dominant balance is between the w term and the forcing Q (Holton 2004;Vallis 2006), which is a state known as the "Weak Temperature Gradient" (WTG) regime (Pierrehumbert 2010a;Koll & Abbot 2016;Pierrehumbert & Hammond 2019). Imposing the WTG approximation and equating these two terms gives Approximating the derivatives in the continuity equation in Equation 5 as ∂/∂x ∼ 1/a and ∂/∂z ∼ 1/H, gives an estimate of the horizontal velocity perturbation (also known as the "eddy" component, the perturbation to the zonal-mean value): This is the same as the "advective balance" in Komacek & Showman (2016), and allows us to estimate the magnitude of the "Stationary Vertical" acceleration term in Equation 21: This is the same as the magnitude of the "Stationary Horizontal" term if it is approximated as u * u * /a. This means that the magnitude of the total acceleration in Equation 21 scales in the same way with forcing as the "Stationary Vertical" term alone.
We estimate the magnitude of the zonal-mean vertical velocity to predict the jet speed from Equation 17. We assume that the zonal-mean form of the thermodynamic equation is governed by the same WTG balance: We also need to estimate the magnitude of the zonalmean vertical velocity w. This depends on the zonal mean of the forcing field Q. The idealised forcing field used in the calculation in Section 3 had no zonal mean, but we can estimate the magnitude of an equivalent "realistic" field where the night-side forcing is uniform (or zero), as in Perez-Becker & Showman (2013). If we assume that in reality only the day-side forcing contributes to the zonal mean of the forcing (and therefore to the zonal-mean vertical velocity), we can estimate the zonal mean of the forcing to be: This gives the following expression for w: Equation 21 then gives the expected jet speed, using the expressions for the magnitude of the acceleration and vertical velocity we have derived: We use the normalisation of the forcing Q from Section 3 to write this in terms of the instellation: This corresponds to a maximum jet speed at z = H, as by approximating the vertical derivative ∂u/∂z as u/H we have assumed that the jet speed increases from 0 at the surface to u over the distance H. The predicted equatorial jet speed is comparable to the magnitude of the velocity perturbations ∼ u * on the equator in a state of nonlinear balance. This arises from the cancellation of the vertical velocity term in the "Stationary Vertical" term with the vertical velocity in the denominator of Equation 17, as well as the cancellation of the vertical derivative in the "Stationary Vertical" term with the vertical integral in Equation 17. It is important to note that u * is not the gravity wave speed √ gH, as might be expected. Instead, it is the magnitude of the horizontal velocity providing nonlinear zonal momentum balance on the equator. As discussed in Section 3.1, the constant of proportionality in Equation 29 depends on the vertical profile of the heating at the substellar point, so will vary for planets without simple heating profiles (and zero albedo) like those in our GCM simulations. We expect that the overall mechanism and scaling relations will remain the same, unless the heating profile has a scale very different to H or the atmosphere is very thick.
In the next section we will show how this prediction approximately matches the jet speeds predicted by the calculation in Dedalus in Section 3, and also matches the GCM results that we will present later.
This prediction of the jet speed can be arrived at even more simply, by interpreting the proposed mechanism as requiring that the magnitude of the acceleration due to "Stationary Vertical" transport, ∂ ∂p (u * ω * ), to be equal to the magnitude of the gradient due to the "Mean Vertical" interaction between the zonal-mean zonal flow and the zonal-mean vertical velocity, ω ∂u ∂p . Setting the "Stationary Vertical" and "Mean Vertical" terms to be equal, and approximating the derivatives as before gives This leads to Equation 29 by the path outlined above. The jet speeds predicted by the calculation in Dedalus and by Equation 29 depend on the value of the Brunt-Väisälä frequency N * . These values are therefore only as accurate as the approximation of a constant N * , and the accuracy of the value of the chosen N * .
The jet speed is strongly sensitive to the acceleration due to gravity g, scaling with g 3 . Each of these factors of g enters through the scale height H, which affects the magnitude of the vertical velocity, the horizontal velocity, and determines the vertical scale of the forcing. This cubic dependence on g is misleading, as the Brunt-Väisälä frequency for an (isothermal) atmosphere depends on g 2 , which means that in reality the jet speed only depends linearly on g. We have not replaced N * with the exact expression for an isothermal atmosphere in Equation 29 as this is not a good estimate for a realistic terrestrial atmosphere, which is far from isothermal. This again shows the importance of an accurate estimate of the Brunt-Väisälä frequency in estimating the stationary wave strength and the resulting jet speed.
The jet speed depends inversely on p s , implying that a very high surface pressure gives a very weak jet. This comes about because we assume that the jet forms at a pressure level comparable to p s on terrestrial planets. For a gaseous planet without a set surface pressure, it would instead be more appropriate to use a pressure level at which the atmosphere is heated strongly by shortwave radiation. In summary, this estimate of jet speed relies on several assumptions about the planet and its atmosphere, and will not apply to planets with different properties. We will investigate the jet speed on gaseous tidally locked planets in a study to follow.

GCM SIMULATIONS OF JET FORMATION
In this section, we examine how well the preceding reasoning fares in explaining the zonal momentum budget of GCM simulations. It is to be hoped that at some stage in the future, the hierarchy of understanding can  Figure 7. The jet profile and acceleration terms averaged over the first 10 days of the GCM simulation with instellation 10 3 W m −2 . The shape of the zonal-mean velocity profile matches the shape of the sum of the "Stationary" acceleration terms, supporting our proposed mechanism, where the the initial acceleration is only due to these terms.
be completed by comparing the GCM simulations to observations of terrestrial tide-locked planets. We will show that the jet speed in the GCM simulations can be approximately predicted by our previous calculations of the stationary wave response.

Numerical simulations
We ran simulations of the atmospheres of terrestrial tidally locked planets in the general circulation model (GCM) ExoFMS , built on the GFDL FMS (Flexible Modelling System) and the associated cubed-sphere dynamical core (Lin 2004). This solves the hydrostatic primitive equations on a cubed-sphere grid with a hybrid "sigma-pressure" coordinate system in the vertical, and uses a semi-grey radiative transfer scheme and a dry convective adjustment scheme (Pierrehumbert 2010b;Hammond & Pierrehumbert 2018). To generate plots, the simulation results are regridded from the cubed-sphere grid to a latitude-longitude grid, and interpolated to a pressure grid in the vertical.
We use 100 vertical levels, to better resolve the vertical profiles of the momentum fluxes shown in Figure 8. The mechanism should still apply to simulations with a lower vertical resolution, and this high vertical resolution would be unnecessary for most other studies. We use a cubed-sphere grid with 48 grid cells on each of the six faces, which corresponds approximately to a spectral resolution of T63. We apply a fourth-order horizontal divergence damping to the velocity fields in the model, with an Earth-like non-dimensional damping coefficient of 0.16 (Lin 2004).
The planetary atmospheres simulated all have a radius of 6 × 10 6 m, a rotation period of 10 d, acceleration due to gravity of 10 m s −2 , specific heat capac-ity 10 3 J kg −1 K −1 , molar mass 28 g mol −1 , and surface pressure 10 5 Pa, in order to match the stationary wave calculations in Dedalus in Section 3. The simulations have different values of instellation at their substellar point, from 10 2 W m −2 to 10 4 W m −2 . This range is comparable to the range of instellation values for the Trappist-1 system (Gillon et al. 2017), which could provide an observational test of the dependence of the strength of global circulation on instellation. The simulations were spun up for 1000 days to ensure the outgoing longwave radiation and zonal-mean zonal velocity reached equilibrium, then data was recorded every 2 days for 1000 days, to ensure a long enough measurement of the time-mean quantities, with enough time resolution to capture the transient quantities. Figure 7 shows the state of the model in the first 10 days of the spin-up of the test with instellation 10 3 W m −2 . The "Mean Vertical" term is very small, as the equatorial jet is very weak at this early stage. The shape of the zonal-mean zonal velocity profile in Figure 7a almost exactly matches the shape of the sum of the two "Stationary" terms in Figure 7b. This is a key result that supports our proposed mechanism, where the initial jet is only due to these "Stationary" terms. Note that these terms have not yet reached their full strength in equilibrium; this is because the stationary velocity perturbations do not reach their maximum amplitude immediately and take some time to evolve. This means that rather than the "Stationary" acceleration terms forming completely, and then the jet forming in response, the jet instead evolves in tandem with the "Stationary" terms. The final balance should still be the same as if the jet only evolved to balance the to-  Figure 8. The equilibrium jet profile and acceleration terms for the simulation with instellation 10 3 W m −2 in the GCM. The jet peaks at about p/ps = 0.4. The zonal momentum budget is dominated by the balance described in Section 4, where the "Stationary" terms produce a jet that provides a "Mean Vertical" acceleration that balances them.

Equilibrium Momentum Budget
tal "Stationary" terms after they formed completely, as long as it does not affect these terms strongly (as shown in Section 4.3). Figure 8 shows the state of the simulation in equilibrium. Figure 8a plots the zonal-mean zonal velocity at the equator of this test in equilibrium, showing that its maximum is at about 400 mbar. In the mechanism we propose, this equatorial jet forms to balance the unbalanced "Stationary" acceleration terms. Figure 8b shows the terms in the zonal-mean momentum budget (Equation 16) for this simulation. Note that the "Mean Vertical" term is zero at the peak of the jet, where the vertical gradient of the jet is zero. The peak of the jet therefore corresponds to the pressure level where the "Stationary Horizontal" and "Stationary Vertical" terms originally cancel, matching the two-dimensional shallow-water model of Showman & Polvani (2010) and Showman & Polvani (2011). The amplitude of these two-dimenensional models corresponds to the amplitude of the first baroclinic mode with vertical wavelength 2H, which is the dominant mode excited by the forcing with wavelength 2H that we introduced in Section 3.
The result that the peak of the jet corresponds to the pressure level where the peaks of the two "Stationary" terms always cancel (in Figures 3 and 8) is notable. If the jet were equilibrated by a linear Rayleigh drag, the peak of the jet would be at the pressure level that had the highest initial acceleration. Instead, it is at a level that has zero zonal-mean zonal acceleration at the initialisation of the model.
In these GCM simulations the "Mean Vertical" term does not totally cancel the peak of the "Stationary Vertical" term, unlike in the idealised calculations in Section 3. The remainder of the momentum balance in the GCM is mostly due to the Rayleigh drag near the surface (not shown explicitly in the plot, but included in the "Sum" line). The "Transient" terms in Equation 14 also contribute to the momentum balance (and are included in the "Sum" term), but are small compared to the other terms. In addition, the "Stationary Horizontal" term in equilibrium in Figure 8 is negative at around 700 mbar, and plays a role in equilibrating the jet.
This does not match the idealised solutions in Section 3, where this term is always positive. This difference appears to arise in the GCM in equilibrium at the pressure level where the vertical shear of the zonal-mean zonal jet is strongest; this could be modifying the stationary wave response by a mechanism similar to that shown by Kato & Matsuda (1992).This process could be similar to the equilibration process suggested by Tsai et al. (2014) that we discussed in Section 4.3, where the jet increases until it decreases the "Stationary" terms enough to reach equilibrium.
The momentum budgets shown in this section have shown the importance of considering the threedimensional structure of the stationary waves induced by forcing in these atmospheres, and the resulting vertical structure of the zonal acceleration. The next section shows how the equatorial jet speed scales with the magnitude of the instellation in these simulations, and compares it to the speed predicted by our calculations in Section 3. Figure 9 plots the zonal-mean zonal velocity at the equator of each test in the GCM, showing how the jet speed increases weakly with increasing forcing, and how the peak moves to lower pressures. This increase in the height of the jet is likely related to an increase in the atmospheric scale height at higher equilibrium temperatures.

Jet Speed Scaling with Forcing
The atmospheric scale height will increase at higher temperatures, increasing the height of the jet according to the theory in Section 5. However, the idealised scale height H should always correspond to the same pressure p s e −1 , for surface pressure p s . This means that the decrease in the pressure of the core of the jet at higher temperatures seen in Figure 9 cannot be explained by an increase in the scale height. It implies that the height of the jet is increasing slightly faster than the scale height H increases with temperature, which cannot be explained by the idealised theory in Section 5. The pressure level of the center of the jet only changes by a factor of 2 over the three orders of magnitude of forcing we model, so we suggest that our approximation that its height is H is fairly accurate. Modelling the exact height of the jet would require a more sophisticated theory, but would be very useful for observable planets where the effect of the jet is measured in phase curves corresponding to particular pressure levels (Parmentier & Crossfield 2017).
The jet speed only depends weakly on forcing -the simulation with the highest instellation has 10 3 times more energy deposited in the atmosphere than the simulation with the lowest instellation, but its jet is only a few times faster. It might be expected from the linear model in Showman & Polvani (2011) (the case when their forcing is weak) that a forcing 10 3 times stronger would produce velocity perturbations 10 3 times stronger as they are linearly related in that model. This would give an acceleration 10 6 times stronger, giving a jet 10 6 times faster. This is not the case in our GCM simulations, as the nonlinear balance at the equator means that the magnitude of the equatorial acceleration depends more weakly on the forcing. In the mechanism we propose, the equilibrium jet speed also depends inversely on the zonal-mean vertical velocity, which increases with increasing instellation, further reducing the dependence of the jet speed on the instellation. Figure 9b compares the results of our hierarchy of models. It shows the maximum jet speed of each test in the GCM, and compares them to the jet speed of the equivalent tests using Dedalus in Section 3, as well as the jet speed predicted by Equation 29. The Dedalus calculations approximately match the magnitude and scaling of the jet speeds in the GCM simulations, suggesting that the mechanism we propose is a good description of the process for jet formation in the GCM. The speed predicted by Equation 29 also approximately matches the speeds in the GCM, suggesting that this is a rea-sonable estimate of the jet speed on these planets (given the approximations discussed in Section 5).
The jet speed in the GCM scales less strongly than the prediction of Equation 29 at high instellations. This may be due to the feedback discussed in Section 4.3, where the faster jet at higher instellations reduces the magnitude of the "Stationary Horizontal" acceleration term that transports momentum towards the equator. A higher instellation and temperature will also increase the relative importance of the nonlinear terms in the thermodynamic equation which we neglected in our derivation of Equation 29. This will reduce the resulting temperature perturbation and stationary wave strength, reducing the jet speed in the GCM relative to our simple estimate. At the lower end, the jet speed increases more rapidly than predicted by Equation 29. We could not find a convincing explanation for this, but suggest that it could be due to a linear balance, rather than nonlinear balance, controlling the magnitude of the stationary waves. At lower instellation, the stationary waves and the jet are closer to the ground and more affected by the linear Rayleigh drag applied there. This rapid increase could also be due to the onset of the resonance identified by Tsai et al. (2014) that we discussed in Section 4.3.
In summary, the hierarchy of models approximately agree, but the exact jet speed deviates from the prediction at very high and low forcing values. Given the model complexity required in other studies to accurately predict a single jet speed in specific planets (Held & Hou 1980;Leovy 1987;Zhu 2006), we suggest that this approximate match over three orders of magnitude in forcing is a reasonable confirmation of the proposed mechanism.

DISCUSSION
The mechanism we propose relies on several approximations and assumptions about which processes are dominant in the zonal momentum budget of the atmosphere, and which are negligible. This section discusses the validity of these approximations and considers other processes which could accelerate or decelerate the equatorial jet.

Complications to this mechanism
The mechanism proposed in this paper is an idealised representation of a complex, time-dependent process. We assume that the "Stationary" terms in the zonal momentum equations form to their full strength instantly, then the jet evolves over time to balance them. In reality, the temperature and velocity perturbations of the stationary waves produced by the day-night instellation  Figure 9. The scaling of zonal-mean zonal equatorial velocity versus instellation in our hierarchy of models. The first panel shows how the jet speed depends weakly on instellation in the GCM simulations, scaling approximately with F 1/2 0 . The peak of the jet moves to slightly lower pressures at higher values of instellation. This may be due to a larger scale height at higher temperatures, giving a longer vertical scale for the forcing and therefore a longer vertical wavelength for the dominant vertical modes in the stationary waves induced by the forcing. The second panel shows how the maximum jet speed in the GCM simulations scales with instellation, and compares this to the calculation with Dedalus and to the prediction from Equation 29. gradient evolve on a timescale of tens of days in the GCM simulations. This is a similar timescale to the evolution of the jet itself, so the "Stationary" acceleration terms and the jet evolve in tandem rather than sequentially. The equilibrium zonal momentum balance in Figure 8 still behaves as if the jet evolved to balance the fully-formed "Stationary" terms, so this approximation is reasonable. We also approximated that the jet does not strongly affect these "Stationary" terms once it has formed; we showed this to be true for the terrestrial planets in this study in Section 4.3, but it may not be true for higher temperature planets, or gaseous planets, which might have stronger jets.
Some of the terms in Equation 14 that we neglected could affect the jet speed. For example, we did not consider the effect of the "Transient" terms, as they were small for the planetary parameters we used. If the travelling waves were comparable in magnitude to the stationary waves we considered, they could accelerate the equatorial jet (Read & Lebonnois 2018) or decelerate it (Showman et al. 2015). Other processes imposed in GCM simulations may also affect the speed, such as the surface drag. This drag extends to 700 mbar in our model; if it were to extend higher or if the jet were to form lower, the drag would affect the zonal momentum budget and the speed of the jet.
We have also neglected the effect of the planetary rotation rate, which may affect the jet speed (Showman et al. 2015;Pierrehumbert & Hammond 2019). The rotation rate determines the horizontal scale of the stationary waves, and therefore the magnitude of the gradients that lead to the "Stationary" acceleration terms. It will also set the horizontal length scale of the acceleration terms, determining the meridional width of the equatorial jet. The calculations in this study assumed that the horizontal scale of the velocity and temperature perturbations was the planetary radius, which is only true in general for atmospheres with sufficiently low Rossby numbers.
At higher rotation rates (for example, for terrestrial Earth-sized planets with periods of less than one day), the (equatorial) Rossby radius is significantly smaller than the planetary radius, and the meridional scale of the equatorial waves will therefore be smaller as well. This should produce a stronger gradient in the velocity fields, giving a higher acceleration and a faster jet, as shown in Pierrehumbert & Hammond (2019). The calculations in this study therefore only apply to planetary atmospheres with stationary waves on the scale of the entire planet. The rotation rate could be included in the estimate of jet speed given by Equation 29, by replacing the planetary radius a with an estimate of the actual meridional scale of the stationary waves induced by the forcing. The equatorial Rossby radius would be a possibility for this scale, but would give a much stronger dependence on rotation rate than was seen in Showman et al. (2015) and Pierrehumbert & Hammond (2019).

Gaseous Planets
The atmospheric circulation of gaseous tidally locked exoplanets such as hot Jupiters is generally more easily observed that that of terrestrial planets. This makes predicting the jet speed of these planets more observationally relevant than for terrestrial planets. In a prelim-inary investigation, we found that not all of the simplifying assumptions we made in this study apply to gaseous planets. Firstly, there appeared to be a non-negligible feedback from the zonal jet on the "Stationary" terms in the zonal-mean momentum budget, as investigated in Section 4.3 and demonstrated in Tsai et al. (2014). This effect appeared strongest where the vertical shear of the jet was strongest, which could be due to an interaction with the vertical external mode (Kato & Matsuda 1992).
Second, we found that the assumption of uniformly upward zonal-mean vertical velocity did not apply to the simulation of a hot Jupiters; as in Mayne et al. (2017), there were regions of zonal-mean downwelling on the equator. This affects the "Mean Vertical" term that equilibrates the "Stationary" terms in our mechanism, which relies on an upward zonal-mean vertical velocity at all vertical levels at the equator. In addition, hot Jupiters may be strongly affected by magnetic drag, which is often represented by a linear Rayleigh drag (Perna et al. 2010). This would modify the mechanism further (in fact, this could make the mechanism simpler as the scaling of the velocities and the equilibration of the jet would be due to this linear drag only).
We will investigate these issues in a study to follow, and are aiming to produce a simple estimate of the magnitude and scaling of the equatorial jet speed on hot Jupiters. This would provide a basis for estimates of observable quantities, such as the predictions of hotspot shift and day-night contrast in Zhang & Showman (2017).

CONCLUSION
This study aimed to explain the formation of the equatorial jet on terrestrial tidally locked planets, and to predict its speed. We proposed that the jet forms by a mechanism in which the day-night forcing induces stationary waves which accelerate a jet, which then interacts with the zonal-mean vertical velocity to produce a deceleration that balances the acceleration from stationary waves in equilibrium. We derived the structure of the zonal acceleration by calculating the threedimensional stationary wave response to the forcing on a tidally locked planet using the Dedalus software. This calculation allowed us to predict the equatorial jet speed for given planetary parameters.
This mechanism was verified by GCM simulations in the ExoFMS model, where the dominant zonal-mean momentum balance was the same as in the proposed mechanism. We ran a suite of simulations with different values of instellation, and showed that the zonal momentum balance and resulting jet speed was approximately the same as predicted by the idealised calculations us-ing Dedalus. With this confirmation of the proposed mechanism, we derived a simplified expression for the jet speed using the WTG approximation: u ∼ 2.53πag 3 σ 1/2 RN 2 p 0 c p F 1/2 0 .
The exact jet speed predicted by this expression depends on parameters such as the Brunt-Väisälä frequency N * , which is difficult to estimate accurately in general. The constant of proportionality also depends on the vertical heating profile at the substellar point, and the planetary albedo. However, we expect that the relation u ∼ F 1/2 0 should hold as long as the mechanism we propose in this study holds.
We discussed why this mechanism may not apply to gaseous tidally locked planets such as hot Jupiters, and suggested how it could be modified to describe their behaviour. The feedback of the jet on the "Stationary" acceleration terms may need to be taken into account for other types of planet, particularly those with high instellation. The mechanism could be modified in other ways to describe other types of tidally locked planetary atmospheres, such as those with cloud layers or significant moisture. These may have different vertical heating profiles and different thermodynamic balances to those assumed in this study of dry, cloud-free atmospheres.
We conclude that this mechanism describes the formation of the equatorial superrotating jet on terrestrial tidally locked planets. It provides a simple prediction of the approximate jet speed, and explains why the jet speed only depends weakly on instellation. We hope that the prediction for jet speeds will be useful for interpreting observations and directing modelling, and will investigate the equivalent mechanism for the formation of equatorial jets on gaseous tidally locked planets in a study to follow.