Temperature Structures Associated with Different Components of the Atmospheric Circulation on Tidally Locked Exoplanets

Observations of time-resolved thermal emission from tidally locked exoplanets can tell us about their atmospheric temperature structure. Telescopes such as JWST and ARIEL will improve the quality and availability of these measurements. This motivates an improved understanding of the processes that determine atmospheric temperature structure, particularly atmospheric circulation. The circulation is important in determining atmospheric temperatures, not only through its ability to transport heat, but also because any circulation pattern needs to be balanced by horizontal pressure contrasts, therefore implying a particular temperature structure. In this work, we show how the global temperature field on a tidally locked planet can be decomposed into contributions that are balanced by different components of the atmospheric circulation. These are the superrotating jet, stationary Rossby waves, and the divergent circulation. To achieve this, we partition the geopotential field into components balanced by the divergent circulation and the rotational circulation, with the latter comprising the jet and Rossby waves. The partitioned geopotential then implies a corresponding partitioning of the temperature via the hydrostatic relation. We apply these diagnostics to idealised general circulation model simulations, to show how the separate rotational and divergent circulations together make up the total three-dimensional atmospheric temperature structure. We also show how each component contributes distinct signatures to the thermal phase curve of a tidally locked planet. We conclude that this decomposition is a physically meaningful separation of the temperature field that explains its global structure, and can be used to fit observations of thermal emission.


INTRODUCTION
Recently, the collection of known planets has grown significantly. This is due to the discovery of numerous planets orbiting stars other than the Sun, called extrasolar planets (or exoplanets for short).
The most observationally accessible exoplanets tend to be on close-in, short-period orbits. Many are sufficiently close to their host star that the gravitational interaction between star and planet synchronises their orbital and rotational periods (Dole 1964;Guillot et al. 1996). A planet in this orbital state is referred to as tidally locked. As tidally locked planets always present the same face to their host star, they have a permanent 'day side' and 'night side'.
These planets can be characterised by observing their thermal emission (Burrows 2014;Crossfield 2015). Timeresolved observations of thermal emission are known as thermal phase curves, and show the disk-integrated infrared flux emitted from a planet as it orbits its host star (Cowan & Agol 2008;Demory et al. 2016), revealing the longitudinal structure of its brightness temperature. Eclipse maps show the two-dimensional structure of brightness temperature on a planet's day side (Majeau et al. 2012).
Atmospheric temperatures and atmospheric circulation are closely related. This is not only due to the ability of an atmosphere to transport heat (Koll & Abbot 2016), and thus maintain a temperature structure against radiative heating and cooling , but also through the pressure gradient forces in the horizontal momentum equations. These pressure gradient forces mean that any distribution of winds must be balanced by a particular temperature distribution. The most familiar example of such a balance is the geostrophic thermal wind relation (Holton 2004, Chapter 3) (1) which can be derived by assuming the flow is in hydrostatic and geostrophic balance. Above, p is pressure, a is the planetary radius, R is the specific gas constant, and f = 2Ω sin ϑ is the Coriolis parameter, with Ω the planetary rotation rate. λ is longitude and ϑ is latitude. The thermal wind relation states that vertical variation in the horizontal wind u = (u, v) must be accompanied by horizontal variation in the temperature T . The utility of geostrophic balance itself is limited to scenarios where the Coriolis term is dominant in the momentum equations, but the requirement that accelerations in the horizontal momentum equations are balanced by horizontal pressure contrasts is a generic one. We will show how an understanding of how different circulation features are balanced by atmospheric temperature structure aids interpretation of observations of thermal emission from extrasolar planets. Tidally locked planets are heated on their day side and cool from their night side, and this pattern of heating and cooling strongly influences their atmospheric circulation (see reviews by Showman et al. 2013;Pierrehumbert & Hammond 2019;Zhang 2020). Heating on the day side drives divergent motion , which in turn generates stationary Rossby waves (Showman & Polvani 2010; see also Sardeshmukh & Hoskins 1988). Interaction between these two circulation components can lead to momentum convergence towards the equator, and the acceleration of a westerly superrotating jet (Showman & Polvani 2011).
In Hammond & Lewis (2021), we showed that each of these circulation components can be isolated from one another by separating the horizontal velocity u into a rotational ('divergence free') component u r and a divergent ('vorticity free') component u d (called a Helmholtz decomposition) The zonal-mean zonal jet and stationary Rossby waves are contained within the zonal-mean and eddy components of u r = u r +u r respectively, and the divergent circulation (comprising, e.g., thermally-direct overturning; cf. Hammond & Lewis 2021, and Kelvin waves; see Appendix A) is contained within u d . Above, χ is the velocity potential and ψ is a streamfunction, defined in terms of the divergence δ = ∇ p · u and vorticity ζ = k × ∇ p u via the relations ∇ 2 p χ = δ and ∇ 2 p ψ = ζ. ∇ p is the horizontal gradient operator acting on isobaric surfaces, and k is the unit vector in the vertical direction.
In this paper, our objective is to use the Helmholtz decomposition to understand how different components of atmospheric circulation contribute to the distribution of atmospheric temperature, and by extension thermal emission, on tidally locked planets. This is achieved by decomposing the geopotential φ = gz into a component that is balanced by purely rotational winds ('the rotational geopotential'), and a component that is due to divergent winds and rotationaldivergent interactions ('the divergent geopotential') (Section 2; see also Trenberth & Chen 1988). We then use the hydrostatic relation to split the global temperature field into rotational and divergent components (Section 3). Finally, we show how each circulation component contributes separately to simulated thermal phase curves by linearly expanding F = σT 4 around the horizontally averaged temperature T 0 (Section 4). A summary of our results is included at the end of the manuscript (Section 5).
The rotational and divergent geopotential components are defined by assuming different balances of terms in the divergence equation (Equation 6; introduced in Section 2). Appendix A illustrates the relationship between balances of terms in the divergence equations, and balances of terms in decomposed momentum equations for the rotational and divergent winds, using the simple example of the linearised shallow water equations. By deriving momentum equations for the rotational and divergent velocities, we show that a steady linear Kelvin wave can exist on a sphere in the absence of drag, in contrast to the equatorial beta-plane where drag is a condition for a linear Kelvin wave (Showman & Polvani 2011).
In each section, the diagnostics are applied to idealised general circulation model (GCM) simulations of atmospheric circulation on dry, temperate terrestrial tidally locked planets. The model is constructed using the Isca modeling framework (Vallis et al. 2018), and is described in detail in Appendix B. We use results from four simulations with different rotation periods, P = 4, 8, 16, and ∞ days (i.e., zero rotation). The zero-rotator is included to illustrate the effect of rotation in the other simulations. We also analysed output from a simulation with P = 32 days. The results from this analysis were nearly identical to those for the P = 16 days simulation and so are not shown.
In this work, we only present results from simulations of terrestrial planets. However, there is no reason why the techniques described herein cannot be applied to gaseous planets (e.g., 'hot Jupiters' or 'mini Neptunes') and this is something we aim to do in the near future.

Definition of rotational geopotential
The horizontal momentum equations can be written in the form (Holton 2004, Chapter 4) where |u| 2 = u 2 + v 2 , ω = Dp/Dt is the pressure vertical velocity, and f = 2Ω sin ϑ is the Coriolis parameter. k f is a linear drag coefficient that is non-zero close to the surface. Equations for the vorticity ζ and divergence δ can then be obtained by taking k · ∇ p × (4) and ∇ p · (4), respectively, The geopotential φ does not appear in the vorticity equation (Equation 5), and thus all of the dynamics that determines φ is described by the divergence equation (Equation 6).  Figure 1. Geopotential height (contours) and horizontal velocity (quivers) at p = 480 hPa for each GCM simulation. For each simulation, full fields are shown in the top row. The next three rows show the zonal-mean rotational circulation, eddy rotational circulation, and divergent circulation. Note that the circulation in the P = ∞ days simulation has no rotational component.
In order to separate the geopotential into components associated with the rotational and divergent winds, we substitute u = u r + u d into Equation (6), which yields where J(ψ, χ) = k · (∇ p ψ × ∇ p χ) is the Jacobian, and β = a −1 ∂f /∂ϑ, with a the planetary radius. A derivation of Equation (7) is given in Appendix C.
To decompose the geopotential into components associated with the rotational and divergent circulation we expand it as φ = φ 0 + φ r + φ d , where φ 0 (t, p) is the horizontal-mean geopotential, which does not contribute to Equation (7) as ∇ 2 p φ 0 = 0. We then define the 'rotational geopotential' as that which is balanced by purely rotational circulation (with no vertical motion), whence we obtain This equation is called the non-linear balance equation (Lorenz 1960) and takes the place of the divergence equation for purely non-divergent horizontal circulation (see, e.g., Holton 2004, Chapter 11). For stationary circularlysymmetric flow, it is equivalent to the gradient wind approximation. If ψ (i.e., the rotational wind distribution) is known, then Equation (8) can be inverted to obtain φ r (e.g., by expressing Equation 8 in terms of spherical harmonics).

Divergent geopotential
Given the rotational geopotential as defined by Equation (8), we define the 'divergent geopotential' to be that obtained as a residual from The remaining terms in Equation (7) which determine φ d involve terms that are purely divergent, and terms due to interaction between the rotational and divergent winds. . Geopotential height (contours) and velocity (u, ω; quivers) in the longitude-pressure plane, averaged between ±40 • latitude, shown for each GCM simulation. The top row shows the full circulation. Contributions to the full circulation from the eddy-rotational and divergent components are shown in the next two rows. Data is averaged over the tropics only in order to show the baroclinic structure of the eddy rotational wind.

Note on the relationship between balances in the divergence equation and the momentum equations
In the next subsection we will use Equations (7), (8), and (9) to partition the geopotential in our GCM simulations into rotational and divergent components, and will discuss the terms in Equations (7) and (8) that contribute the most to the decomposed height fields. In anticipation of this, we note that a large term in, e.g., Equation (8), which defines φ r , does not necessarily correspond to a dominant balance in the momentum equation for u r . This is because an additional term is introduced if Equation (8) is integrated to obtain an equation for u r ; the same applies if Equation (7) (minus the solely rotational terms) is integrated to obtain an equation for u d . This is illustrated for the linearised shallow water equations on the sphere in Appendix A. Figure 1 shows the horizontal wind u = (u, v) and geopotential height z in the mid-troposphere (p = 480 hPa) for each GCM simulation. The full circulation is shown in the top row, and subsequent rows show the rotational circulation -split into contributions from the zonal-mean (jet) and eddies (stationary waves) -and the divergent circulation. Figure 2 shows the vertical structure of the full circulation, eddy rotational circulation, and divergent circulation averaged between ±40 • latitude. In the P = 4 days simulation, the circulation is mostly contained within the rotational component, whereas in the P = 8 and P = 16 days simulations it is more evenly divided into rotational and divergent con-tributions. Finally, the circulation in the P = ∞ days simulation is purely divergent. In this subsection, we describe the important features of each component of the circulation that can be identified in Figures 1 and 2. We refer to Appendix D throughout, which contains additional analysis of the contributions of individual terms in the expanded divergence equation (Equation 7) to each of the rotational and divergent height fields.

Application to simulations
The zonal-mean rotational height field φ r in each of the P = 4, 8, and 16 days simulations is in approximate geostrophic balance with the zonal-mean rotational wind u r : which gives rise to a local height maximum centred on the equator, associated with the superrotating equatorial jet (Hammond & Pierrehumbert 2018). In the P = 4 days simulation there are additional retrograde jets poleward of ±50 • in each hemisphere which lead to a secondary maximum in the zonal-mean rotational height field at each pole. The eddy-rotational circulation in each simulation has a horizontal wind and geopotential structure characteristic of a stationary equatorial Rossby wave. The longitude of the eddy rotational height maximum is located eastwards of the substellar point due to a Doppler-shift by the zonal-mean zonal wind (Tsai et al. 2014). For the rapidly rotating P = 4 days case, the structure and phase offset of the Rossby waves is very similar to that obtained in a linearised shallow water model with a prescribed background jet (e.g., Figure 10 in Hammond & Pierrehumbert 2018). As the rotation period is reduced, the Rossby lobes with a positive geopotential anomaly become elongated with respect to those with a negative geopotential anomaly. Similar phenomenology is also exhibited in the GCM simulations presented in Edson et al.
(2011, see their Figure 5). In Appendix D we show that this is due to non-linear effects by partitioning the rotational height into contributions from the linear and non-linear terms in the non-linear balance equation (Equation 8). In the vertical, the stationary Rossby waves have a wavenumber-1 modulation (Tsai et al. 2014), as can be identified in the second row of Figure 2 (most clearly evident in the P = 4 and 8 days simulations). The divergent circulation generally displays the most complex structure of all of the circulation components in Figures 1 and 2. The exception to this is the P = ∞ days simulation, for which the circulation is purely divergent and takes the form of a single isotropic, thermally-direct, overturning cell that features air rising on the day side, near the substellar point, and sinking on the night side. The geopotential structure associated with the overturning circulation is predominantly associated with the ∇ p · (ω∂u d /∂p) term in Equation (7), which is dominant in regions of ascent and descent. In the boundary layer the surface friction term is important, and in the outflow regions aloft the horizontal nonlinear −∇ 2 p [|∇ p χ| 2 /2] term is important (see Appendix D). The balance of terms in Equation (7) for the P = 16 days is similar to that in the zero-rotation simulation. In particular, the ∇ p · (ω∂u d /∂p) term is still important in these simulations in the region of strong ascent near the substellar point. However, additional terms in Equation (7) that describe interactions between the rotational (mostly u r ) and divergent circulations are also important (Appendix D), and our interpretation is that in the P = 16 days case the overturning circulation is 'tipped over' by the eastward superrotating jet.
In the P = 4 days simulation, the structure of the divergent height field is very different from that in the P = 16 or P = ∞ days simulations, and is harder to interpret. Both its horizontal and vertical structure are somewhat wave-like, although very different to the structure of a linear Kelvin wave on the equatorial beta plane (see, e.g., Matsuno 1966). In Appendix D we show that the total divergent height structure shown in Figures 1 and 2 for the P = 4 days simulation is mostly due to the linear u d β term, which takes the form of a Kelvin wave in a linearised shallow water model (Appendix A), and terms describing interaction between rotational and divergent horizontal winds. Finally, in the P = 8 days simulation, the divergent height field has features that are characteristic of both the slowly rotating simulations, with thermally-direct overturning, and the P = 4 days simulation with non-linear and wave-like contributions to the height, suggesting that it resides in an intermediate regime.

RELATIONSHIP BETWEEN GEOPOTENTIAL COMPONENTS AND TEMPERATURE
In this section, we use the hydrostatic relation to obtain components of the atmospheric temperature that are associated with the rotational and divergent geopotential components. This is achieved by requiring that each component of the geopotential (φ 0 , φ r and φ d ) is in hydrostatic balance with the corresponding component of the temperature (T 0 , T r and T d ), i.e., (11) Understanding how each component of the circulation contributes to atmospheric temperature structure is important for interpreting observations of thermal emission from tidally locked planets, such as eclipse maps (Majeau et al. 2012) and phase curves (Cowan & Agol 2008;Demory et al. 2016;Morris et al. 2021). Figure 3 shows the temperature structure in the lower troposphere (p = 750 hPa) for each GCM simulation. As with the geopotential height shown in Figure 1, the temperature is decomposed into contributions from the zonal-mean and eddy rotational circulations, and the divergent circulation. Figure 4 shows the vertical structure of the full circulation, eddy rotational circulation, and divergent circulation, averaged meridionally over all latitudes for consistency with the disk-integrated phase curves shown later.
Generally, the horizontal and vertical structures of the temperature components shown in Figures 3 and 4 are similar to the structure of the height components shown in Figures 1 and 2. The main different between them is a phase offset in their vertical structures, which appears as the temperature is obtained from the vertical derivative of the geopotential (this is why the temperature in Figure 3 is plotted deeper in the atmosphere than the height field in Figure 1).
In each of the simulations with rotation, there is an eastward 'hot spot' shift from the substellar point. In the P = 8 and 16 days simulations, this is almost exclusively due to the rotational standing Rossby waves which are Doppler-shifted eastwards from the substellar point by the zonal-mean jet (Tsai et al. 2014;Hammond & Pierrehumbert 2018). The zonal-mean rotational temperature that balances the zonalmean rotational wind cannot contribute a longitudinal shift to the temperature structure by definition, although it does contribute towards 'brightening' the equator with respect to the poles. The divergent temperature structure in the P = 16 days simulation is similar to that of the P = ∞ days simulation which only features overturning circulation. The overturning circulation is a direct response to heating at the substellar point, and it appears that this strong coupling to the surface prevents it from being Doppler-shifted much by the jet.
By contrast, in the P = 4 days simulation, the divergent temperature structure has a greater eastward shift than the eddy rotational circulation (which is now only shifted by a few degrees longitude from the substellar point), although its amplitude is comparatively lower. The eastward shift of the divergent temperature in this simulation is consistent with the wave-like nature of its divergent circulation and height field,  Figure 5. Phase curves calculated using the broadband thermal emission from a radiating level at prad = 750 hPa, for each GCM simulation. The secondary and primary eclipses would be at phase 90 • and −90 • . The solid black curve shows emission calculated using the full temperature field, and the dashed grey curve shows emission using the linear approximation given by Equation (13). The dashed red and blue curves show the contributions to the linear approximation from the rotational and divergent temperature fields, respectively. discussed in the previous section. The temperature structure associated with the zonal mean rotational circulation is also different to that in the P = 8 and 16 days simulations, due to the retrograde jets that exist in this simulation in midlatitudes, and now contributes a brightening to both the equator and the poles, relative to the mid-latitudes.
4. THERMAL PHASE CURVES Thermal phase curves show the infrared flux emitted from a planet as it orbits its star. As a planet progresses through its orbit we see emission from different angles, so these observations characterise the longitudinal distribution of temperature on extrasolar planets (Cowan & Agol 2008;Koll & Abbot 2015). Thermal phase curves are already available for hot Jupiters and high-temperature rocky planets (Demory et al. 2016;Kreidberg et al. 2019), and should increase in quality in the coming years (Bean et al. 2018) following the successful launch of the James Webb Space Telescope. With JWST it should also be possible to produce thermal phase curves for smaller and cooler terrestrial planets (Deming et al. 2009).
In this section, we compute thermal phase curves from the decomposed temperature structures to show how each component contributes to their observable features. We approximate that the broadband thermal emission I ↑ originates from a specific radiating level p rad , so that We do this to simulate observing at a wavelength at which the atmosphere is optically thick, instead of using the model's actual semi-grey outgoing longwave radiation (OLR), because in our simulations the OLR is dominated by the surface emission and so is not an instructive demonstration of the potential contribution to the phase curve from the atmosphere. This also corresponds to the case of gaseous planets where the thermal emission is entirely due to the atmosphere and its circulation. As Equation (12) is non-linear in T , we cannot easily separate out contributions to I ↑ from each of T r and T d . We therefore linearise Equation (12) about the horizontal mean temperature T 0 , which separates contributions to I ↑ from T r and T d . Normalised thermal phase curves are then computed according to where ξ is the phase angle of the planet in its orbit (Cowan & Agol 2008). Figure 5 shows phase curves computed for each GCM simulation using Equation (14), for a radiating level located at p = 750 hPa, which corresponds to the level for which the horizontal temperature structure was shown in Figure 3. Solid black curves show the phase curve computed using the non-linear expression given by Equation (12), and dashed grey curves show the phase curve computed using the linear approximation given by Equation (13). There is generally good agreement between the phase curves computed using the non-linear and linear expressions, and the difference between the two is reduced with increasing rotation period as deviations from the horizontal-mean temperature become smaller. This shows that the decomposition of the phase curve into rotational and divergent components by the linear approximation is useful.
The contributions from the rotational and divergent components of the temperature to the linear approximation are shown in Figure 5 as dashed red and blue lines, respectively. In each simulation with rotation, the amplitude of the rotational contribution to the phase curve is larger than that of the divergent contribution, and leads to a negative offset in the phase curve peak from ξ = 90 • (corresponding to an eastward hot-spot shift in the temperature maps shown in Figure 3). As discussed in Section 3, this is entirely due to the Doppler-shifted Rossby waves, as the temperature structure that balances the jet cannot itself contribute a direct hot-spot shift by definition (although it is the velocity of the jet that Doppler-shifts the waves).  The negative phase offset associated with the stationary Rossby waves increases monotonically with increasing period. However, this trend is not reflected in the full phase curves, due to contributions from the divergent circulation. In the P = 4 days simulation, the divergent circulation, which in this simulation resembles a non-linear wave that is Doppler-shifted eastwards, also contributes an eastward phase shift, that is in fact larger than that associated with the rotational circulation. This means that the peak of the full phase curve has a greater negative phase offset from ξ = 90 • than it would if only the rotational circulation (stationary Rossby waves) contributed. By contrast, in the P = 16 days simulation, the divergent circulation now takes the form of overturning circulation, and is restricted to be close to the substellar point, due to strong coupling between the overturning circulation and the instellation. This means that the divergent contribution to the phase curve is not offset from ξ = 90 • , and the full phase curve has a lesser negative phase offset than it would if due to the stationary Rossby waves alone.
Our analysis of the vertical temperature structure in the previous section shows that the phase curves in Figure 5 will depend strongly on our choice of radiating level. To illustrate this we show an additional phase curve in Figure 6 for the P = 4 days simulation, this time for p rad = 200 hPa. In this figure, the sign of the normalised flux for each curve (full, rotational, and divergent) is essentially flipped with respect to that shown in Figure 5, and the peak of the full phase curve has a positive offset from ξ = 0 • , which would correspond to a westward hot-spot shift in temperature. In our simulation, this is due to the wavenumber-1 structure of the stationary Rossby waves in the vertical (see Figure 4), which means that they can contribute either a negative or a positive shift to the phase curve peak depending on the exact pressure level from which they contribute to thermal emission. This is a key result of understanding thermal phase shifts as due to stationary wave shifts rather than heat advection by a jet, and is a potential explanation for unexpected observations of westward phase shifts (Dang et al. 2018). 5. SUMMARY Hammond & Lewis (2021) showed that the Helmholtz decomposition u = u r + u d is a useful tool for separating out different components of the atmospheric circulation on tidally locked planets. The jet and stationary Rossby waves are contained within u r and divergent circulation is contained within u d , which takes the form of thermally-direct overturning for slowly rotating planets (e.g., our P = 8 and P = 16 days simulations; see also Hammond & Lewis 2021), but has a complex linear Kelvin / non-linear wave structure for more rapidly rotating planets (cf. our P = 4 days GCM simulation, and the shallow water results in Appendix A).
The objective of this work was to build upon Hammond & Lewis (2021) by relating the atmospheric temperature structure to each of the circulation components contained within u r and u d . This was achieved by defining rotational and divergent components of the geopotential using balance relations in the divergence equation, which were then used to define a temperature structure via the hydrostatic relation. To illustrate the utility of decomposing the temperature in this way, we applied it to output from idealised GCM simulations of atmospheric circulation on terrestrial tidally locked planets (run using Isca; Vallis et al. 2018), considering four rotation periods: P = 4, 8, 16, and ∞ days (i.e., zero rotation).
The temperature maps we obtain (Figure 3) show that both the rotational and divergent circulations make non-negligible contributions to the horizontal temperature structure. Understanding what determines the relative strength of these circulations, and how they depend on properties of a particular exoplanet, will be crucial for correctly interpreting observations of thermal emission such as phase curves (Cowan & Agol 2008) and eclipse maps (Majeau et al. 2012).
Our analysis shows that both circulation components can contribute to the hot spot shifts that have been inferred from observations of thermal emission (e.g., Demory et al. 2016) (Figures 3 and 5). In our simulations, both the rotational circulation and the divergent circulation contribute an eastward hot spot shift if thermal emission is assumed to come from the lower troposphere (p = 750 hPa). For our P = 8 and P = 16 days simulations, the hot spot shift is almost entirely due to the rotational stationary Rossby waves (which are Doppler-shifted eastwards by the superrotating jet). In the P = 8 days case, this is because the amplitude of the rotational temperature component is much larger than the divergent temperature component. In the P = 16 days case, both components have a similar amplitude, but the divergent component does not contribute a hot spot shift as it is dominated by thermally-driven overturning, which means it is strongly coupled to the pattern of instellation which is centered on the substellar point. By contrast, in the P = 4 days simulation, both the rotational component and the divergent component (which is now more wave-like) contribute an eastward hot spot shift if emission is from the lower troposphere. However, we show that the rotational component can also contribute a westward hot spot shift if emission comes from higher up in the atmosphere (p = 200 hPa in our simulations; Figure 6), which is due to its wavenumber-1 structure in the vertical (Figure 4).
From this analysis, we suggest that the temperature structure of the atmosphere of a tidally locked planet can be divided into four physically meaningful components: 1. The overturning divergent component, providing a temperature peak at the substellar point and uniform temperature elsewhere.
2. A wave-like divergent component, providing a zonalwavenumber-1 sinusoidal temperature modulation on the equator which can be Doppler-shifted by a zonal jet.
3. The Rossby-wave-like eddy rotational component, providing a zonal-wavenumber-1 sinusoidal temperature modulation with maxima off the equator, which can be Doppler-shifted by a zonal jet.
4. The zonal-mean rotational component, i.e. the zonal jet, providing a zonally uniform temperature field with a meridional gradient determined by the gradient of the jet itself.
These components could be used to fit thermal phase curves or eclipse maps with physically meaningful functions. We hope that the techniques and analysis presented herein will be of use for interpreting observations of thermal emission on tidally locked planets.

NTL was supported by Science and Technology Facilities
Council grant ST/S505638/1. MH was supported by a Junior Research Fellowship at Christ Church, Oxford. The authors are grateful to Dr. Man-Suen Chan for IT support. In this appendix, we use the example of the linearised shallow water equations on the sphere to show how partitioning the height field into rotational and divergent components can be used to obtain momentum equations for the rotational and divergent wind. We use these equations to interpret the existence of a steady linear Kelvin wave in the spherical geometry, in the absence of drag, which does not exist on the equatorial beta plane.
In the absence of friction, the linearised shallow water equations may be written (Vallis 2017, Chapter 3) ∂u ∂t Above, η is the free surface displacement (or height), H is the mean fluid depth, and Q is a mass source or sink term. As with the primitive equations (Section 2), equations for the vorticity and divergence can be obtained by taking k·∇× (A1) and ∇ · (A1), respectively, Assuming a steady state, we can partition the divergence equation into two equations that define the rotational and divergent height, η r and η d : which can then be integrated to yield momentum equations for the rotational and divergent circulations in a steady state where γ is defined by Equation (A8) can be obtained by taking the curl of either Equation (A6) or (A7), and the second equality comes from the vorticity equation in steady state. Equations (A6) and (A7) each contain an additional term involving the potential function γ, when compared with the momentum equation for the full horizontal velocity u. These terms arise as an integration constant and are due to the variability of f with latitude, and their proper representation is contingent on working in a spherical geometry (Trenberth & Chen 1988). Figure 7 shows horizontal velocity and height fields from a linear shallow water simulation using the GFDL 1 spectral shallow water model. We ran the model at T85 resolution for 500 days with a timestep of 300 seconds, and plot output averaged over the final 100 days of the run. The mass source term is configured for a tidally locked planet as following Perez-Becker & Showman (2013), where η eq is given by η eq = H + ∆η eq cos λ cos ϑ on the day side,  ∆η eq is the difference in the radiative equilibrium height between the substellar point and the terminator. In Equation (A9), τ rad is the radiative relaxation timescale. The GFDL model solves the full non-linear shallow water equations, and the solution was kept in a linear regime by using a small forcing amplitude, η eq /H 1, following Perez-Becker & Showman (2013). An important feature of the model is that it does not include a linear drag (i.e., τ drag = ∞ in Perez-Becker & Showman 2013). The planetary parameters we used for this simulation were: a = 6.371 × 10 6 m, P = 8 days, and g = 9.81 m s −2 . The equilibrium height was H = 10000 m, the magnitude of the forcing was ∆η eq = 10 m, and the radiative timescale was τ rad = 0.1 days. We chose these parameters for demonstrative purposes as they give rotational and divergent components of similar magnitude, but our conclusions are not sensitive to the chosen parameters as long as the forcing is linear. In Figure 7 the velocity and height is split into eddy rotational and divergent components. The zonal mean rotational component is not shown as it is small.
The decomposed rotational and divergent circulations shown in Figure 7 have structures that resemble the linear solutions for equatorial Rossby and Kelvin waves, respectively, on the equatorial beta plane (Matsuno 1966), although they are modified by the spherical geometry. This partitioning of the Rossby waves into the rotational component and the Kelvin waves into the divergent component is consistent with Vanneste & Vial (1994), who show that in the limit of small Lamb parameter, Λ ≡ 4Ω 2 a 2 /(gH) 1 (in our simulation Λ ≈ 0.1), Rossby waves are purely rotational and Kelvin waves are purely divergent. The Kelvin wave obtained here has non-zero meridional velocity, which derives from the spherical geometry (Yamamoto 2019). The Rossby component is very similar to that obtained by Showman & Polvani (2011) on the equatorial beta plane in the absence of friction (as is the case here; cf. their Figure 14). Notably, Showman & Polvani (2011, their Appendix C) show that in the absence of drag, an equatorial Kelvin wave cannot exist on the equatorial beta plane.
In a steady state, the zonal momentum equation is As f = 0 at the equator (and also v = 0 due to the hemispheric symmetry of the problem), Equation (A11) requires that ∂η/∂λ = 0 at the equator. This is satisfied in our linear shallow water simulation (see the left-hand panel in Figure  7). However, Showman & Polvani (2011) additionally show that on the equatorial beta plane, equatorial Rossby wave solutions have no longitudinal η variation at the equator, which by Equation (A11), means that an equatorial Kelvin component cannot exist (as it would introduce a non-zero ∂η/∂λ). By contrast, both the equatorial Rossby and Kelvin components in Figure 7 have a non-zero height modulations that cancel one another so that the total ∂η/∂λ at the equator is still zero. This is due to the effect of the spherical geometry on the waves, which now are governed by a zonal momentum equation of the form (e.g., for u r , from Equation A6) If ∂γ/∂θ is non-zero at the equator (which will be the case if v is symmetric about the equator; γ obtained in our simulation is shown in Figure 8 for reference), then neither will ∂η r /∂λ| ϑ=0 = −∂η d /∂λ| ϑ=0 . As a consequence, on the sphere and in the limit of small Lamb parameter, the existence of a rotational stationary Rossby wave requires the existence of a divergent stationary Kelvin wave, and vice versa. Showman & Polvani (2011) introduced a linear drag to their beta plane model to generate the stationary Kelvin wave necessary for the acceleration of a zonal jet, but we have shown here that this linear drag is not necessary for the formation of a Kelvin (-like) wave on a sphere rather than on a beta plane.

B. DESCRIPTION OF GENERAL CIRCULATION MODEL
The GCM simulations analysed in the main text were run using Isca, which is a a framework for building idealised general circulation models of varying complexity (Vallis et al. 2018). Isca is built on the GFDL spectral dynamical core, which is used to integrate the primitive equations of motion forwards in time. In the present study we use a dry GCM configuration, i.e., water vapour and moist processes are not included in the model atmosphere. The model is configured to simulate the atmospheric circulation of terrestrial planets with a lower boundary modeled as a 'slab ocean'. All planetary parameters aside from the rotation rate and properties of the stellar irradiation (e.g., radius, surface pressure, gravitational acceleration) are set to be those of the Earth. Four rotation periods are used, P = 4, 8, 16, and ∞ days (as described in the main text), and the distribution of stellar irradiation is set to be that of a tidally locked planet (defined below). The main components of the model are described below.

B.1. Dynamical core
The dynamical core integrates the primitive equations forwards in time, using a semi-implicit leapfrog scheme with a Robert-Asselin time filter. The equations are solved in vorticity-divergence form on a thin spherical shell (similar to Equations 5 and 6) using a pseudospectral method in the horizontal (prognostic fields are represented by a triangular truncation of spherical harmonics), and a finite difference method in the vertical (see, e.g., Satoh 2014, Chapters 21-23). The vertical coordinate used is a terrain-following 'sigma' coordinate, defined σ = p/p s , where p is pressure and p s is surface pressure. The simulations presented here use a T42 triangular truncation, which corresponds to approximately 2.8 • latitude-longitude resolution. We use 25 unevenly spaced levels in the vertical, distributed according to σ = exp[−5(0.05z + 0.95z 3 )] wherez is equally spaced in the unit interval. A ∇ 8 hyperviscosity term is applied to the momentum and thermodynamic equations operating on a timescale of 0.1 days at the grid scale. The data analysed in the main text is interpolated from the model's sigma coordinate to a pure pressure coordinate.

B.2. Radiative transfer
We use a semi-grey radiative transfer code to model radiative heating and cooling. Radiative fluxes are split into two wavelength bands: one covering longwave (thermal emission) wavelengths, and the other covering shortwave (stellar heating) wavelengths.
In the shortwave band, the atmosphere is assumed to be transparent, and all incoming radiation is absorbed at the surface (i.e., the surface albedo has been folded into the solar constant). The instellation at the top-of-atmosphere is imposed with the following distribution, appropriate for a tidally locked planet: where λ 0 = 0 • is the substellar longitude, and S 0 = 1000 W m −2 .
In the longwave band, upward and downward radiative fluxes are computed according to where τ = τ 0 (p/p 0 ) is the longwave optical depth, and we set τ 0 = 1 and p 0 = p s . The longwave radiative heating in the model is then The boundary conditions for the longwave fluxes are F ↓ (p = 0) = 0 and F ↑ (p = p s ) = σT 4 s where T s is the surface temperature.

B.3. Surface energy budget, boundary layer sensible heat transport, and convection
The surface is modeled as a static slab ocean that can exchange energy with the atmosphere in the vertical but cannot transport heat horizontally. The surface temperature evolves according to where C is a specified heat capacity that corresponds to an ocean mixed-layer depth of 2.5 m, F ↓ sfc is the downward flux of longwave radiation incident on the surface, and H is the surface sensible heat flux. H is computed according to where ρ is the density, the subscript 'a' indicates quantities that are evaluated at the lowest model level, and C = 0.001 is a dimensionless bulk transfer coefficient. Turbulent vertical transport of heat within the boundary layer is parametrised following Thatcher & Jablonowski (2016) as where z is height, w = Dz/Dt is the vertical velocity, and θ is the potential temperature. The diffusion coefficient K is calculated according to p pbl = 850 hPa is the top of the boundary layer, and p strat = 100 hPa controls the rate of decrease of boundary layer diffusion with height. The contribution to the model temperature tendency from the boundary layer diffusion is then Finally, the model includes a dry convection scheme, which instantaneously restores the atmospheric temperature structure to the dry adiabat whenever it is convectively unstable.

B.4. Surface friction and sponge layer
Surface friction is parametrised following Held & Suarez (1994) as a linear Rayleigh drag where k f is defined as k f,s = 1 day −1 is the drag timescale at the surface, and the top of the frictional boundary layer is located at σ b = 0.7.
D. SUBDIVISION OF HEIGHT FIELD This appendix includes some additional analysis of the terms in the divergence equation (Equation 7) for our simulations, to identify those that contribute most strongly to the rotational and divergent height components. In particular, we analyse the divergent circulation by highlighting the term that corresponds to overturning circulation, and comparing it to the terms that correspond to wave-like divergent circulation and non-linear rotational-divergent interactions In the P = 4 and P = 8 days simulations, the eddy rotational height (Figure 1) is dominated by the linear ∇ p ·[f ∇ p ψ] term (not shown). It has the structure of a stationary equatorial Rossby wave, that is Doppler-shifted eastwards by the zonal-mean jet (Tsai et al. 2014;Hammond & Pierrehumbert 2018). In the P = 16 days simulation, the eddy rotational height also resembles a stationary Rossby wave, but the lobes with a positive geopotential anomaly are noticeably elongated with respect to those with a negative geopotential anomaly (Figure 1). Figure 9 show the eddy rotational height field for the P = 16 days simulation, divided into contributions from the linear term, and the non-linear ∇ p · [∇ 2 p ∇ p ψ] − ∇ 2 p [|∇ p ψ| 2 /2] terms. This subdivision of the height field shows that the elongation of the stationary Rossby is mostly due to non-linear rotational-rotational interactions.
The divergent height field has a complex structure in each of the simulations presented in Figures 1 and 2. In order to interpret it, we show the divergent height field decomposed into contributing terms in Figures 10 and 11. The P = ∞ days simulation represents a limiting case where the divergent circulation takes the form of an isotropic, overturning cell (see Figures 1 and 2). In this scenario, Figure 11 shows that the divergent height structure is due to the ∇ p · (ω∂u d ∂p) term in regions of ascent near the substellar point, and descent on the night side. Meanwhile, the k f δ (surface friction) contribution is important in the boundary layer, and non-linear horizontal divergent-divergent interactions are important in the outflow area surrounding the region of ascent near the substellar point. The balance of terms in the P = 16 days simulation is similar to that in the P = ∞ days case, although there is a contribution to the divergent height field from the horizontal non-linear rotational-divergent interaction terms, −∇ 2 p [J(ψ, χ)]−J(χ, ∇ 2 p ψ) (mostly through u r ). Our interpretation of the divergent height structure in the P = 16 days simulation is that it is consistent with a thermally-direct overturning circulation that is 'tipped over' by the superrotating jet.
By contrast, the most important contributions to the divergent height in the P = 4 days simulation are from the linear u d β term, and non-linear horizontal rotational-divergent interactions (note that surface friction is still important in the boundary layer). Figure 10 shows the horizontal structure of the linear u d β term and horizontal non-linear rotationaldivergent interaction terms for the P = 4 days simulation. The height field associated with the u d β term is qualitatively similar to that of the linear Kelvin wave obtained in the linear shallow water simulation presented in Appendix A. In the vertical, the height associated with this term has a wavenumber-1 structure ( Figure 11). However, unlike in the  Figure 11. Subdivision of divergent geopotential height into contributions from different terms in the divergence equation (Equation 7 in the main text). The top row shows the total divergent height field for each simulation, and subsequent rows show contributions to the divergent height from different terms (labelled to the left of each row). Contributions from terms in the full divergence equation that are due to solely rotational circulation are not shown, as they do not contribute to the divergent geopotential by definition.
linear shallow water simulation, the divergent height field is modified by the contribution from the non-linear rotationaldivergent interaction term (largest in the upper troposphere), which yields the complex structure shown in Figure 1. The divergent height in the P = 8 days simulation has contribu-tions from all of the terms discussed for both the P = 4 days and P = 16 days simulations, and it appears that this simulation is an intermediate case between the Kelvin/nonlinear wave regime of the P = 4 days simulation, and the thermally-direct overturning regime of the P = 16 and P = ∞ days simulations.