Longitudinally asymmetric stratospheric oscillation on a tidally locked exoplanet

Using a three-dimensional general circulation model, we show that the atmospheric dynamics on a tidally locked Earth-like exoplanet, simulated with the planetary and orbital parameters of Proxima Centauri b, support a longitudinally asymmetric stratospheric wind oscillation (LASO), analogous to Earth's quasi-biennial oscillation (QBO). In our simulations, the LASO has a vertical extent of 35--55 km, a period of 5--6.5 months, and a peak-to-peak wind speed amplitude of -70 to +130 m/s with a maximum at an altitude of 41 km. Unlike the QBO, the LASO displays longitudinal asymmetries related to the asymmetric thermal forcing of the planet and to interactions with the resulting stationary Rossby waves. The equatorial gravity wave sources driving the LASO are localised in the deep convection region at the substellar point and in a jet exit region near the western terminator, unlike the QBO, for which these sources are distributed uniformly around the planet. Longitudinally, the western terminator experiences the highest wind speeds and undergoes reversals earlier than other longitudes. The antistellar point only experiences a weak oscillation with a very brief, low-speed westward phase. The QBO on Earth is associated with fluctuations in the abundances of water vapour and trace gases such as ozone which are also likely to occur on exoplanets if these gases are present. Strong fluctuations in temperature and the abundances of atmospheric species at the terminators will need to be considered when interpreting atmospheric observations of tidally locked exoplanets.


INTRODUCTION
Rocky planets orbiting M-class stars have become a prime target for exoplanet research due to their prevalence and potential habitability (Yang et al. 2020;Schwieterman et al. 2018;Kopparapu 2013;Heath et al. 1999;Scalo et al. 2007;Wandel 2018). As the habitable zones (HZ) of cool M stars are found close to the host star, planets within the HZ have a high chance of being tidally locked (Barnes 2017;Pierrehumbert & Hammond 2019). These planets receive stellar irradiation on only one hemisphere, an orbital configuration not found for planets in the solar system, which has a substantial impact on the planet's atmospheric dynamics and chemistry.
An accurate understanding of atmospheric dynamics is important both to assessing the theoretical habitability of tidally locked rocky planets and interpreting observations. Theoretical studies of tidally locked rocky exoplanets have used 3-D atmospheric general circulation models (GCMs) to explore their possible atmospheric dynamics and climate states (Koll & Abbot 2016;Boutle et al. 2017;Wolf 2017;Komacek & Abbot 2019). A key finding of GCM studies is that tidally locked planets within the HZ of M-class stars can retain temperate climates in spite of their uneven thermal forcing as a result of winds redistributing heat from their hot daysides to their cold nightsides (Merlis & Schneider 2010;Pierrehumbert & Hammond 2019;Wordsworth 2015;Joshi et al. 1997). This day-night temperature gradient of a tidally locked exoplanet has already been observed in some cases, even for a rocky planet (Demory et al. 2016). However, most past studies have examined mean climate states, masking time-dependent phenomena which may impact observations by upcoming instruments such as the James Webb Space Telescope's Near Infrared Spectrograph (NIRSpec) (Mollière et al. 2017;Morley et al. 2017) and the ESA Atmospheric Remote-sensing Infrared Exoplanet Large-survey mission (Venot et al. 2018;Tinetti et al. 2018). As observations are performed in a narrow time window, modes of variability with long periods such as analogues of the Earth's quasi-biennial oscillation (QBO) may result in observed temperatures or abundances that differ by an order of magnitude or more from predictions based on mean GCM states. Other studies have used GCMs, coupled with synthetic planetary spectrum generators, to simulate observations of transmission spectroscopy for transiting exoplanets (Suissa et al. 2020;Arney et al. 2017;May et al. 2021;Boutle et al. 2020;Lines et al. 2018), enhancing the value of these models as an aid to interpreting observational data (Brogi et al. 2016;Louden & Wheatley 2015). However, gaps in our understanding of the physics driving atmospheric phenomena on tidally locked planets, partly driven by differences in model structure and parameterisations (Sergeev et al. 2020;Fauchez et al. 2021;Turbet et al. 2021;Sergeev et al. 2021), may cause misinterpretations of atmospheric observations. Global 3-D models of tidally locked planets exhibit an area of deep convection around the substellar point (Hammond & Lewis 2021;Yang et al. 2013;Kopparapu et al. 2016), leading to permanent cloud cover and heavy precipitation in this region (Boutle et al. 2017;Sergeev et al. 2020;Labonté & Merlis 2020;Yang et al. 2013). The tropospheric and stratospheric circulation is characterised by stationary Rossby waves of zonal wavenumber 1 induced by the day-night thermal forcing (Merlis & Schneider 2010;Showman & Polvani 2011;Hammond & Lewis 2021). These waves have been identified as a driving source of the superrotating equatorial jets that develop in tidally locked simulations of both Earth-like planets and hot Jupiters (Showman & Polvani 2010, 2011Tsai et al. 2014;Debras et al. 2019Debras et al. , 2020. The jets are associated with two pairs of cyclones-anticyclones: one pair forms in the northern hemisphere, rotating clockwise on the dayside and anticlockwise on the nightside, while the other pair forms in the southern hemisphere and rotates in the opposite directions. The cyclones contribute a zonal wind vector component near the equator that points westward in the substellar region and eastward in the antistellar region. This results in a background flow which is more eastward on the nightside than on the dayside, a pattern that holds true in both the troposphere and the stratosphere. The day-night thermal forcing additionally creates a region of rising air on the dayside and subsiding air on the nightside (Hammond & Lewis 2021;Joshi et al. 1997).
On Earth, the background flow determines the zonal structure of the QBO. In the QBO, the stratospheric zonal wind forms a pattern of vertically stacked jets flowing in opposite directions. The pattern propagates downwards from roughly 20 − 40 km at a speed of 1 km per month such that the wind direction at a given height reverses over time with a mean periodicity of 26 − 28 months, with substantial variation from cycle to cycle (Dunkerton et al. 2015;Braesicke 2015). The oscillation of the zonal wind is accompanied by oscillations in air temperature and trace gases. The basic mechanism driving the QBO was described by Lindzen & Holton (1968) and Plumb (1977) and reproduced in a laboratory setting by Plumb & McEwan (1978). The minimum requirements for the QBO to occur are a background flow that is a function of altitude and at least two gravity waves propagating in opposite directions horizontally, as well as upwards vertically (Plumb 1977). An eastward propagating gravity wave will travel upwards until it encounters an eastward perturbation in the background flow, where it is absorbed, accelerating the flow in the eastward direction. Westward propagating waves pass through this area of eastward flow until they reach a westward perturbation and are absorbed in turn. The resulting pattern of stacked jets propagates downwards over time as the bottom of each shear zone is accelerated and the lowest jet dissipates due to viscous diffusion (Plumb 1977;Dunkerton et al. 2015). Although equatorially trapped Rossby and Kelvin waves contribute some momentum flux, it is believed that up to 80% of the momentum flux driving the QBO is provided by convectively generated gravity waves (Dunkerton 1997;Lane 2015). Gravity waves are generated when the forces of gravity or buoyancy react against vertical perturbations in a fluid medium. The waves travel internally within the fluid until they are absorbed at critical levels where the velocity of the background flow is comparable to the wave's horizontal phase speed (Booker & Bretherton 1967), making the planet's mean flow an important element in the phenomenon.
Here, we describe the first reported simulation of a stratospheric oscillation, analogous to the QBO on Earth, on a tidally locked planet with an Earth-like atmospheric composition. This new stratospheric oscillation consists of periodic reversals in the direction of the zonal wind, an associated cooling and warming of the equatorial atmosphere, and fluctuations in the stratospheric water vapour content. While analogous to the QBO in its causal mechanism, the oscillation displays longitudinal asymmetries arising directly from the planet's tidally locked state. Accordingly, we refer to the oscillation on a tidally locked planet as a longitudinally asymmetric stratospheric oscillation (LASO) and reserve the term QBO for the phenomenon on Earth. To put our results into context, throughout the paper we compare and contrast key features of the tidally locked LASO with the QBO on Earth (Baldwin et al. 2001).
Section 2 describes the global 3-D model we use to describe atmospheric dynamics on a tidally locked planet, the metrics we use to characterize the LASO, and the radiative transfer model we use to simulate JWST observations of the fluctuations in atmospheric water vapour content caused by the LASO. Section 3 describes the mechanism underpinning the LASO and characteristics of the LASO that differ from those of the QBO. We also consider the effect of LASO-related water vapour fluctuations on transit observations. In Section 4 we discuss the implications of the LASO for future spectral atmospheric observations and the sensitivity of the LASO to differences in sub-grid scale model parameterisations. We conclude the paper in Section 5.

Model description
We use the Global Atmosphere 7.0 configuration of the Met Office Unified Model to simulate the atmosphere of a tidally locked Earth-like exoplanet, nominally Proxima Centauri b (Walters et al. 2019). For brevity, we provide a basic description of the model and refer the reader to Mayne et al. (2014), Boutle et al. (2017), andSergeev et al. (2020) for more details. We run the model at a horizontal resolution of 2 • latitude by 2.5 • longitude, and with 60 vertical levels from the surface to up to 85 km ( Table 1). The atmospheric levels are identical to those in Boutle et al. (2017) up to 38 km. The remaining 47 km are accounted for by an additional 22 levels. Table 2 shows the orbital and planetary parameters used in this study, following Boutle et al. (2017), with data taken from Anglada-Escudé et al. (2016) and Turbet et al. (2016). Proxima Centauri b is simulated as an aquaplanet with a slab ocean (Frierson et al. 2006), with a heat capacity of 10 7 J K −1 m −2 representing a mixing layer of 2.4 m. We specify an atmospheric composition of 100 % nitrogen with a trace fixed CO 2 abundance of 5.941 × 10 −04 kgkg −1 . The model includes interactive water vapour and a full water cycle with unlimited evaporation from the slab ocean, precipitation in the liquid and ice phases, and a cloud scheme. The stellar spectrum for Proxima Centauri was taken from BT-Settl (Rajpurohit et al. 2013) with T ef f = 3000 K, g= 1000ms −2 , and metallicity = 0.3 dex. The star is treated as quiescent. We spin up the model for 900 Earth days from an equilibrium state taken from a previous run of the Proxima Centauri b model. The simulation was determined to be in equilibrium when the incoming and outgoing radiation at top-of-atmosphere were balanced, evaporation balanced precipitation in the global mean, and the surface temperature no longer showed a long-term upward or downward trend. The experiments we report were run for 1800 Earth days in total, and we analysed a sampling period of 900 days after equilibrium, which is long enough to include multiple oscillations. We define the substellar point to be 0 • longitude and latitude. The antistellar point is located at 180 • longitude and the eastern and western terminators at 90 • E and 90 • W, respectively. Hereafter, days refer to Earth days.
Earth GCMs typically require both resolved and parameterised waves to generate a QBO (Garfinkel et al. 2021). Past studies have shown that the Unified Model can reproduce a realistic oscillation (Scaife et al. 2000(Scaife et al. , 2002. Our model can sustain resolved long wavelength gravity waves and inertio-gravity waves generated by physical and numerical imbalances in the flow (as discussed in Section 3). The model also includes a sub-grid scale non-orographic gravity wave parametrisation (Warner & McIntyre 1996, 2001Bushell et al. 2015). Within this parameterisation, an unsaturated spectrum of gravity waves is launched near the surface and carries a vertical flux of horizontal wave pseudo-momentum upwards in equal amounts towards the north, south, east, and west. Upward-propagating flux is diminished by Doppler modification and density amplification. Wave pseudo-momentum is conserved by an opaque model lid which sets the outgoing flux to zero (Bushell 2021). The scheme includes a scale factor for the vertical flux and a characteristic 11.186 Rotation speed (rad s −1 ) 6.501×10 −6 Eccentricity (·) 0 Obliquity (·) 0 Radius (km) 7160 Acceleration due to gravity (m s −2 ) 10.9 Table 2. Simulation parameters used to describe an Earth-like tidally locked Proxima Centauri b.
(spectrum peak) wavelength which must be tuned to the vertical resolution of the model to achieve a QBO. We use a factor enhancement for the vertical flux of horizontal pseudo-momentum of 1.2 and a characteristic wavelength of 4.3 km with model level heights shown in Table 1. These are the standard values used in the GA 7.0 for Earth with the model levels in Table 1 and we retain them in order to simulate Proxima Centauri b as an Earth-like planet. The gravity wave scheme can be run with two different gravity wave source locations: 1) a standard scheme in which the gravity wave launch flux is globally invariant, distributed uniformly around the model; and 2) an experimental scheme in which the gravity wave launch flux is determined by local precipitation, developed based on an empirical correlation found between gravity waves and precipitation on Earth (Bushell et al. 2015). Our study is based on the standard globally invariant scheme. To test whether our results are robust against the choice of scheme, we ran the experimental precipitation-dependent scheme and found that the results are comparable with the exception of the amplitude of the longitudinal asymmetry of the mean flow, which we discuss in Section 3.
In addition to the control experiment described below, we performed a sensitivity run to understand the influence of horizontal resolution by increasing it to 2 • latitude by 2 • longitude. We performed a separate sensitivity test by decreasing the timestep from 20 to 2 minutes. We found that these sensitivity runs did not result in major changes from the control experiment.

GRAVITY WAVE-INDUCED ACCELERATION
To confirm that the reversal in the direction of the stratospheric jets is caused by absorption of gravity waves, we calculate the acceleration of the atmosphere in each model gridbox due to gravity waves, F , as: where u ′ is the high-pass filtered zonal anomaly of the zonal wind, w ′ is the time anomaly of the vertical wind, and z is the height. This definition is adapted from the formulation of the wave momentum flux described by Plumb (1977) and used in its original form by Showman et al. (2019) to identify shear transition zones between alternating jets. We modified u ′ to better describe the characteristic flow on a tidally locked planet. The zonal anomaly is usually defined as the deviation of that quantity from its zonal mean and, in a zonally symmetric baseline flow, identifies transient phenomena. However, tidally locked planets lack zonal symmetry. The planet's stationary Rossby waves skew the zonal anomaly so that it represents the Rossby waves rather than transient flows. To adapt u ′ to this asymmetrical background flow, we use a high-pass filter to remove the influence of waves with zonal wavenumber less than 5 (i.e., long Rossby waves), following Sugimoto et al. (2021). The remaining waves are resolved long wavelength gravity and inertio-gravity waves. Additionally, the zonal anomaly of the vertical wind w would be impractical for identifying waves because rising air on the dayside and subsiding air on the nightside cancel each other out when a zonal mean is calculated. We instead use the time anomaly to distinguish transient, wave-like fluctuations in w from the background. The resulting metric identifies wave-induced acceleration caused by resolved gravity and inertio-gravity waves only.

LATITUDINAL EXTENT
We expect that a QBO-like oscillation on an exoplanet may extend to higher latitudes than Earth's QBO, covering a larger portion of the terminators, which may affect transit spectroscopy. To characterise the latitudinal extent L of the LASO, we use a formula derived by Haynes (1998) for the latitudinal extent of the QBO: where σ is the frequency of time variation of the applied force, α is the radiative damping rate, N is the buoyancy frequency, D is the depth scale of the force, and β 0 is the meridional gradient of the Coriolis parameter evaluated at the equator. We assume that all quantities except β 0 are identical to those for Earth, and focus our analysis on the impact of differences in the Coriolis factor.

PERIOD
The period of the LASO is of interest because the duration of the temperature and atmospheric composition oscillations may affect observations and in particular could cause discrepancies between repeated observations of the same target exoplanet. We determine the period T of the LASO using a wavelet analysis of the wind velocity time series at the altitude where the LASO has its greatest amplitude. Following Plumb (1977), the period of the oscillation is inversely proportional to the momentum flux: wherek,ĉ,N ,μ, andF are scaling factors for the wavenumber and phase speed of the wave driving the oscillation, the buoyancy frequency, the thermal dissipation rate, and the momentum flux, respectively. The parameter η is a dimensionless number that depends on the specific details of the problem. In Plumb's simple numerical model, the wavenumber and phase speed of the waves can be specified, while in the real atmosphere the oscillation is driven by a broad spectrum of waves whose constituents vary with location and time. However, the inverse relationship between the period and the momentum flux at the lower boundary of the QBO shear zone is also accepted to exist for the real QBO (Dunkerton et al. 2015). We assume all other values are identical to those for Earth and focus on the impact of the wave momentum fluxF on the period.

LAGRANGIAN ROSSBY NUMBER
To understand how a tidally locked orbital configuration affects where gravity waves are generated on the planet we use the Lagrangian Rossby number, Ro (L) , to identify resolved gravity wave sources. Ro (L) is not an exclusive measure of gravity wave sources, but is often used to detect unbalanced flows which give rise to gravity waves (Zhang et al. 2000;Lin 2007;Sugimoto et al. 2021). It is defined as the ratio of the acceleration of a parcel of air to the Coriolis acceleration and essentially measures the local departure from geostrophy: where v h is the horizontal wind vector and f is the Coriolis parameter. We modify Ro (L) in several respects to adapt the metric to a slowly rotating planet. Following Sugimoto et al. (2021), we disregard the local wind tendency term, ∂v h ∂t , and fix the value of f near the equator (≤ ±10 • ) equal to f at 10 • to extend the metric to cover the equatorial region. As the slowly rotating Proxima Centauri b is not in geostrophic balance, the resulting values of Ro (L) are very large. We normalise Ro (L) by dividing by the minimum value, which we consider the baseline ageostrophy for our simulation. A larger value indicates a relative local increase in imbalances in the flow. We also exclude the top five model levels and latitudes above ±40 • to avoid interference from instability at the model top and poles.

NASA Planetary Spectrum Generator
To describe the spectrum of the atmospheric water vapour content that might be observed by the James Webb Space Telescope's NIRSpec instrument, we use the NASA Planetary Spectrum Generator (Villanueva et al. 2018), publicly available at https://psg.gsfc.nasa.gov/. The NIRSpec instrument has a range of 0.6 to 5.3 µm, covering water vapour features in the infrared at 1.4, 1.8, and 2.7 µm, with spectral resolutions (R) of 100, 1000, and 2700. As inputs to the PSG, we use the basic orbital, planetary, and stellar parameters for the Proxima Centauri system, as well as pressure, temperature, altitude, H 2 O, N 2 , and CO 2 data from our UM simulations. We generate a transit spectrum for each latitude at both the eastern and western terminators (i.e., every 2 • of latitude, for a total of 180 individual spectra) and average the outputs to get a final spectrum, following Suissa et al. (2020). Figure 1 shows the wind speed and temperature oscillations in our simulation of the atmosphere of the tidally locked Proxima Centauri b. The period of the LASO is much shorter (5 − 6.5 months) than that of the QBO on Earth (26 − 28 months). There is also an asymmetry between the dayside and nightside. Nightside zonal winds are more positive (eastward) than dayside winds and have shorter westward phases. Consistent with the theory that describes the QBO, the dayside-nightside asymmetries in the mean planetary-scale circulation, in particular the dayside-nightside differences in wind speeds, are the primary factor responsible for the asymmetry of the LASO. We find that the location of gravity wave sources also affects the local flow. Figure 1c shows temperature fluctuations with a peak-topeak amplitude of ∼ 60 − 70 K, more than an order of magnitude larger than the variation of ∼ 4 K associated with the QBO (Dunkerton et al. 2015).

Longitudinally asymmetric stratospheric oscillation
To confirm that the LASO is generated by the same physical mechanism as the QBO, we calculated the waveinduced acceleration due to gravity waves and compared it to the change in zonal wind direction over time. Where the oscillation is driven by gravity waves, the direction of the acceleration caused by gravity waves should match the direction of the change in the zonal wind throughout the vertical extent of the LASO. Figure 2a shows the filtered u ′ field with zonal wavenumbers smaller than 5 removed (Section 2.2). We find positive and negative fluctuations in wind speeds due to resolved long wavelength gravity waves, with shear zones at 30 km, 41 km, and 50 km that correspond to regions where the flow changes direction. Waves with a leftward tilt correspond to westward-propagating modes and those with a a rightward tilt correspond to eastward-propagating modes. The band of reduced u ′ intensity between 15 km and 30 km corresponds to the permanent superrotating equatorial jet in the troposphere. Figure 2b shows the zonal mean gravity wave-induced acceleration calculated using Equation 1 and averaged over 90 days, corresponding to approximately half of a LASO cycle. The contour lines show the change in zonal wind speed over the same 90 days. Areas of positive (negative) wind direction change match well with areas of gravity wave-induced eastward (westward) acceleration. Regions below 25−30 km show little change in zonal wind speeds and negligible gravity wave-induced acceleration.

Latitudinal extent of LASO
The LASO extends to higher latitudes than the QBO on Earth. The latitudinal extent of the QBO, defined here by the full width at half maximum of a Gaussian fit to the horizontal cross section of the zonal wind at the height of the QBO maximum (Schenzinger et al. 2017), is roughly 15 • S-15 • N (Braesicke 2015). Figure 3 shows that the LASO extends from about 60 • S-60 • N by this definition. As the LASO on a slowly rotating tidally locked planet extends to higher latitudes, the oscillations in air temperature and in abundances of atmospheric species will cover more of the terminators and potentially have a larger impact on observations. Previous work suggested that at higher latitudes the Coriolis force would counteract the acceleration due to gravity waves, limiting the QBO to equatorial regions on the Earth (Lindzen & Holton 1968). Equation 2 gives a value for Earth of about L Earth =10 • (Baldwin et al. 2001), which underestimates the actual extent of 15 • (Dunkerton & Delisi 1985) by ≃ 33%. The Coriolis force on Proxima Centauri b is an order of magnitude weaker than on Earth due to the planet's slower rotational period of 11.2 days (Ω=6.5×10 −6 rad/s). If we assume the values of the other quantities in Eqn. 2 for Proxima Centauri b are identical to those for Earth, this gives a latitudinal extent of √ 10L Earth or ≃47 • . Scaling this value to account for the 33% underestimation gives the latitudinal range simulated by the model. The broader extent of the LASO compared to the equatorial QBO agrees with the description of slowly rotating tidally locked planets as 'all tropics' (Showman et al. 2013). A LASO on a slow rotator may therefore have a larger impact on transit spectroscopy than one on a more rapidly rotating planet.  Another major difference between the LASO and the QBO is the period of oscillation. The period of the QBO has a large variability, with a mean value of 26-28 months (Braesicke 2015). Our calculations show that the LASO in our simulation of Proxima Centauri b's atmosphere has a period of 5-6.5 months. The period of the QBO is not directly related to the planet's rotation rate (Plumb & McEwan (1978) and Eqn. 3). Because a tidally locked planet does not rotate with respect to its star, the oscillation is also not related to seasonal and annual variations. The shorter period of LASO instead implies stronger wave momentum flux and wave-induced acceleration.

Period of LASO
On Earth, up to 80% of the QBO forcing originates from convectively generated gravity waves (Lane 2015). Figure  4 shows that deep convection in our simulation of Proxima Centauri b occurs within an area of 30 • of the substellar point. This area has the shape of a pointed oval and is asymmetrical, with convection consistently more intense at the western edge of the region. This region of deep convection is collocated with a gravity wave source (as shown in Figure 6b). The accepted dependence of the QBO's period on the gravity wave momentum flux as expressed in Equation 3 supports the conclusion that greater wave momentum flux, driven by the intense convection, accounts for the shorter period of the LASO compared to the QBO on Earth. Accordingly, we find higher wind speeds in the LASO than the QBO. The QBO has maximum eastward winds of 15 ms −1 and maximum westward winds of -30 ms −1 (Dunkerton et al. 2015), while the LASO's eastward wind speeds reach 100 ms −1 and westward winds reach -75 ms −1 .

Longitudinal asymmetry of the LASO
The LASO exhibits longitudinal asymmetries that do not have an analogy on Earth. We show below that this is related to underlying asymmetries in the steady-state atmospheric dynamics of tidally locked planets and ultimately to the asymmetrical thermal forcing.
As described above, the emergence of the QBO is associated with gravity wave absorption at critical layers in the atmosphere where the speed of the background flow is similar to the phase speed of the wave, leading to wave breaking and the deposition of momentum. Longitudes and heights where the background flow is faster will absorb faster gravity waves and experience stronger acceleration, perpetuating the asymmetry. This accounts for the stronger nightside eastward winds seen in Figure 1b, as the nightside experiences higher eastward winds due to the Rossby wave contribution to the zonal wind, as well as for the height at which the LASO has its maximum amplitude (41 km), which is also where the Rossby waves are most intense.  In the Unified Model, convection is diagnosed by the surface buoyancy flux and an undilute parcel ascent from grid points where this flux is positive. It is then categorised as deep if the parcel reaches neutral buoyancy above 2.5km or the freezing level (whichever is higher), and if the vertical velocity is greater than 0.02ms −1 within a layer at least 1.5km thick above the level of neutral buoyancy (Walters et al. 2019).
The zonal wind does not change direction simultaneously at every longitude. Figure 5 shows the temporal variations of the 41 km zonal wind along the equator at two phase transitions, one from eastward to westward and the other from westward to eastward flow. In both transitions, longitudes close to the western terminator transition first, in addition to exhibiting higher wind speeds. Figure 5a shows that the westward winds begin to appear around 90 • W on days 125-130, while other longitudes transition around days 135-140. The converse is also true (Figure 5b). Wind speeds are consistently highest at 90 • -120 • W. Close inspection also shows that the antistellar point has a very weak westward phase, with wind speeds of only 5-10 ms −1 in some places, compared to 100-120 ms −1 during the eastward phase. Periods of westward flow at the antistellar point are comparatively brief. The weakness of the westward phase here implies lower westward acceleration, which may be due to the absorption of westward-propagating gravity waves before they travel around the planet to the nightside, a lack of critical levels available to absorb westward-propagating modes, an asymmetry in the underlying resolved gravity wave spectrum itself, or a combination of these mechanisms.
The amplitude of the LASO is largest at the western terminator and smallest at the antistellar point ( Figure 5). This pattern cannot be explained directly by the presence of Rossby waves because the equatorial wind vector component contributed by the Rossby waves is larger at the antistellar point than at the western terminator. A possible explanation for the greater amplitude and early phase switching of the oscillation at the western terminator is the presence of a second gravity wave source in this region. The western terminator is a jet exit region where both stratospheric and tropospheric jets decrease rapidly in speed. The steep gradient in wind speeds is caused by the convergence of the two opposing wind vectors contributed by the Rossby wave cyclone-anticyclone pair. Jet exit regions are well-known to be sources of gravity waves (Plougonven & Zhang 2014). To investigate further, we calculated the Lagrangian Rossby number (Eqn. 4). Figure 6a shows that at 41 km, where the LASO reaches its maximum amplitude, a strong gravity wave source is present at the western terminator. In the troposphere, unbalanced flows are found in the region of deep convection (Figure 6b and Figure 4). The equatorial altitude-longitude cross-section (Figure 6c) confirms that areas of instability that can generate gravity waves occur preferentially westward of the substellar point in the troposphere and around the western terminator in the lower stratosphere. To link wave acceleration to the jet exit region, Figure 6d shows the gravity wave-induced acceleration at the equator (as calculated in Figure 2), overlaid with the zonal wind. The strongest acceleration is found at the jet exit region for both the tropospheric and the lowest stratospheric jets between 20-45 km. The western terminator both generates gravity waves due to the slowing of the jet speeds in this region and also absorbs high-energy gravity waves because this is where the critical levels exhibit the highest wind speeds.   The result is an unusually unstable region where zonal winds reverse earlier than elsewhere and may oscillate back and forth in shorter time periods, as seen in Figure 5a between 125-145 days and 120 • -90 • W .
We also ran the model using the experimental precipitation-dependent gravity wave scheme to investigate the impact of moving the gravity wave source on the characteristics of the LASO, as described in Section 2. In the precipitationdependent scheme, parameterised gravity waves are launched preferentially from the convection zone around the substellar point on the dayside, creating an asymmetry in the parameterised gravity wave source that does not exist in the standard scheme. While the period, latitudinal extent, amplitude, resolved gravity waves, and water vapour oscillation did not differ between the precipitation-dependent and standard schemes, the longitudinal asymmetry in  wind speeds and phase transitions decreased and the LASO became more zonally symmetric. In the standard scheme, the wind velocity maximum at the western terminator and the minimum at approximately 30 • E differ by ∼100 ms −1 at 41 km, while in the precipitation-dependent scheme they differ by only ∼20 ms −1 . The stronger gravity wave source at the substellar point appears to play a similar role as the source at the jet exit region in modifying the local flow at higher altitudes, counteracting the westward component contributed by the dayside cyclone. Although the overall mean circulation remains longitudinally asymmetrical, this experiment demonstrates the non-linearity in outcomes resulting from interactions between Rossby waves and gravity waves in the simulation. In summary, as discussed above, the LASO is shaped by interactions between multiple components of the circulation typical of tidally locked planets. Figure 7 shows a schematic illustration of the interplay between Rossby waves, the superrotating tropospheric jet, and the gravity wave sources in the substellar region and at the western terminators. The diagram represents the mean circulation generated with the standard, globally invariant gravity wave source.

Atmospheric water vapour oscillations
On Earth, the oscillation of the zonal wind direction is associated with variations in the amount of atmospheric water vapour of the order of ±10% (Randel et al. 1998;Wang et al. 2020). An analogous oscillation on exoplanets is of interest because large fluctuations in atmospheric water vapour may be observable. In our simulations, specific humidity at 41 km varies about ±50% on average, with peak oscillations up to 200% of the mean. This results in oscillations with minimum and maximum values that can differ by an order of magnitude.
To test whether a fluctuation of this magnitude would be detectable, we generated synthetic planetary spectra as observed by the NIRSpec instrument aboard JWST for the above maximum and minimum water vapour concentrations. We found the fluctuations are too small to be observable for an Earth-sized planet like Proxima Centauri b, but as QBOlike phenomena have also been detected on Jupiter and Saturn in the solar system (Showman et al. 2019;Orton et al. 1991;Leovy et al. 1991), a LASO could potentially occur on a larger planet with greater potential for observation. This could be particularly important for short-period, potentially tidally locked, super-Earth or mini-Neptune-type planets detected by missions such as the Transiting Exoplanet Survey Satellite (see Bluhm et al. (2021) for a recent example) Using a global 3-D climate model we show that Earth-like tidally locked planets can exhibit a longitudinally asymmetric stratospheric wind oscillation (LASO), analogous to the quasi-biennial oscillation (QBO) on Earth. As QBO-like oscillations are also known to exist on Jupiter and Saturn (Showman et al. 2019;Orton et al. 1991;Leovy et al. 1991), our finding suggests a similar conclusion for gas giants. The LASO reported in our study has implications for the climatology and observations of tidally locked exoplanets, as well as for our understanding of GCM results.

Oscillations in meteorological and atmospheric composition
It is well established that the QBO on Earth is responsible for oscillations in the abundance of many atmospheric species, notably ozone and methane (Bowman 1989;Huang et al. 2008;Randel et al. 1998;Li et al. 2009). In this initial study, we have focused on the properties and acceleration of the LASO. However, by including atmospheric chemical kinetics, future work could also explore LASO-induced variations in species important for observations. For example, Earth's total ozone column fluctuates by ±1 ppm (∼ 4%) on the 26-28 month timescale of the QBO. A variation of this magnitude is unlikely to be detectable with existing instrumentation, but an ozone column oscillation may be larger in amplitude on a tidally locked planet than on Earth, in line with our simulation results for water vapour. Studies of ozone photochemistry on tidally locked planets (Yates et al. 2020) have also shown that the planet's nightside accumulates ozone in its cold traps, resulting in a zonally asymmetrical ozone layer without analogy on Earth and high levels of ozone in one hemisphere. Future work could explore how this tidally locked ozone chemistry would interact with oscillations in dayside production of ozone through stratospheric photochemistry. In addition, LASOrelated temperature anomalies in our simulation reach 40-50 K (∼ 25% deviation from the mean), while the QBO of temperature on Earth is limited to 1-2 K (Tegtmeier et al. 2020). Large variations in air temperature will also potentially affect measurements of thermal emission spectra from the atmosphere.

Sensitivity to model parameterisations
Gravity wave parametrisation schemes and differences in ability to generate a QBO-like phenomenon may account for some of the differences found in modelling studies of the same exoplanet. For example, predicted nighttime surface temperatures for Proxima Centauri b diverge by 50 K in different studies (e.g. Turbet et al. (2016); Boutle et al. (2017)), a discrepancy that can only be partly explained by differences in stellar irradiation used by these studies. Boutle et al. (2017) instead proposed that these temperature discrepancies could be explained by differences in convection parameterizations and model resolutions, which affect the water vapour profile and consequently the dayside-nightside heat transport. Gravity wave parameterisations likewise vary between models and Watkins & Cho (2010) have shown that gravity waves can have a significant impact on the development of the background flow, temperature, and heat redistribution from dayside to nightside on tidally locked planets. Gravity waves can become trapped within jets and travel longitudinally from the gravity wave source on the dayside until the jet speed develops a critical level where the wave can be absorbed (Pitteway & Hines 1965;Nappo 2012). Differing gravity wave schemes may therefore result in different predictions for heat redistribution and dayside-nightside temperature contrast (Watkins & Cho 2010).

The importance of the western terminator for atmospheric dynamics
Our results identify the western terminator as a region of particular complexity. At the equatorial western terminator, the zonal wind components contributed by the standing Rossby waves on the dayside and nightside meet and oppose each other. We have also shown in Figure 6 that the western terminator is a jet exit region where the flow is unbalanced and internal gravity waves are generated and locally reabsorbed. These factors lead to significant turbulence and variability due to the interplay between these elements. As transit spectroscopy is expected to be a key tool in constraining the composition of exoplanet atmospheres, understanding the dynamics of the terminator regions of tidally locked planets is important to interpreting the results of observational studies. If the time scale of an oscillation is longer than the observing window, observers may view the planet in a low-temperature/high-temperature or lowabundance/high-abundance phase of an oscillation. As we noted above, a LASO on a hot Jupiter could be intense enough to be detectable. The high amplitude of the oscillation at the western terminator, a critical region for transit spectroscopy, suggests that the LASO should be studied on tidally locked hot Jupiters to determine its possible impact on observations.

CONCLUSION
Using a global 3-D general circulation model, we reported the discovery of a phenomenon analogous to Earth's quasibiennial oscillation (QBO) in simulations of the atmosphere of an Earth-like tidally locked planet, nominally Proxima Centauri b. As the phenomenon exhibits asymmetries related to the tidal locking, we refer to it as the longitudinally asymmetric stratospheric oscillation (LASO). The LASO begins above the tropospheric equatorial superrotating jet and extends vertically from roughly 35-55 km. It consists of a periodic reversal in the direction of the zonal winds within this region, accompanied by oscillations in temperature. We found the LASO in our simulation to have a latitudinal extent of ± 60 • and a period of 5-6.5 Earth months. The shorter period compared to its Earth counterpart is explained to a first order by stronger wave driving from gravity wave generation in the deep convection zone at the substellar point. We also identified a secondary gravity wave source in the jet exit region at the western terminator. Our analysis found that wind speeds in the LASO reach a maximum at the western terminator and that reversals of direction at the western terminator precede those at other longitudes. We further found that the antistellar point experiences a very weak oscillation: zonal winds are preferentially eastward-flowing, with only brief and low-speed westward phases.
Our results highlight the significance of atmospheric variability in global 3-D simulations of exoplanets and point towards a need for further studies to investigate this largely uncharted territory. In this initial study, we focused on the characteristics and mechanism of the LASO and did not explore the potential effect of oscillations in abundances of atmospheric species other than water vapour on exoplanet observations. Future work in this area could use a 3-D model with online atmospheric chemistry such as presented in Yates et al. (2020) to explore such effects. While we performed sensitivity tests by varying the horizontal and temporal resolution, we expect that the LASO's fundamental dependence on the planet's mean circulation will make it sensitive to factors that alter the planet's atmospheric dynamics, such as different atmospheric compositions, the presence of land, and the planet's rotation rate. On the Earth, the QBO is also known to interact with other dynamical phenomena such as the semi-annual oscillation, Madden-Julian oscillation (Martin et al. 2021), and polar vortex. We leave such potential interactions between the LASO and similar sources of variability in the atmosphere for future work.
As the body of observational data from exoplanets grows in the coming decades, our understanding of the parameter space for planetary climate states must grow with it. This parameter space includes time-dependent atmospheric phenomena such as the one we describe in this paper which, even if they are familiar features on Earth, will inevitably take new forms on tidally locked planets. Characterisation of the sources and properties of atmospheric variability on tidally locked exoplanets promises to refine our interpretations of observations and, in the case of terrestrial planets, deepen our understanding of the Earth in the context of the larger population of planets inside and outside the solar system.