Observational Verification of High‐Order Solar Tidal Harmonics in the Earth's Atmosphere

This study combines 8 years of middle atmospheric wind data observed at 52°N latitude from two radars in different longitudinal sectors to investigate solar tides. The power spectral density of horizontal winds exhibits a −3 power law within the frequency range 2.0 < f < 7.0 cpd (equivalent to periods 3.6 − 12.0 hr). Particularly noteworthy are the 4.8‐ and 4‐hr tides, exhibiting signal‐to‐noise ratios ranging between 13 and 16 dB, surpassing the 0.01 significance level. This challenges their previous oversight in literature, possibly due to inadequacies in prevailing noise models. Cross‐spectra between longitudinal sectors emphasize the dominance of sun‐synchronous components in the six lowest‐frequency tides. Composite spectra indicate that tidal enhancements during SSWs resemble regular seasonal variations. Intriguingly, year‐to‐year spectral variations suggest that these enhancements are more influenced by seasonal dynamics than by SSW, contrasting with established literature. These findings underscore the need to reevaluate tidal harmonics and consider appropriate noise models in future studies.


Introduction
In the universe, celestial objects often engage in orbital motions, such as galaxies orbiting supermassive centers and binary star systems in mutual embrace.These celestial ballets generate periodic gravitational gradients and radiation heating, resulting in tides, a phenomenon seen as fluid movements and geological deformations.Within our solar system, tides play vital roles in regulating heat and mass circulation.They drive the dust transport on Mars (Wu et al., 2020), intensify volcanic on Jupiter's moon Io (Lainey et al., 2009), and asymmetrically heat Saturn's moon Enceladus (Nimmo & Pappalardo, 2006).Earth experiences the influence of both solar and lunar tides, stemming from either gravitational or thermal excitation and impacting diverse processes, such as carbon storage (Wang et al., 2019), glaciation and climate changes (Green et al., 2020), thermal circulation and ice shelf melting (Richter et al., 2022), orbital evolution (Mitchell & Kirscher, 2023), geophysical activities and hazards (Kasahara, 2002), meteorological patterns (Kohyama & Wallace, 2016), Earth's angular velocity and biological evolution (Mitchell & Kirscher, 2023).Despite their pivotal roles, existing literature has predominantly focused on a few lowest-frequency tidal harmonics, with higher-frequency tides receiving less attention, partly due to challenges in distinguishing high-frequency tides from regional buoyancy fluctuations.
However, beyond these four harmonics, the existence of higher-order harmonics remains uncertain.While higherorder harmonic signatures have tantalizingly emerged in just two instances (He, Forbes, et al., 2020;Hedlin et al., 2018), numerous unresolved questions remain about these high-order signals.It is imperative to ascertain whether the reported signals genuinely originate from tidal forces rather than sporadic buoyancy waves (namely, gravity waves).Moreover, questions arise concerning the regularity of these tidal occurrences and the total count of regularly occurring solar tidal harmonics.To tackle these inquiries, this study analyzes jointly 8 years of middle atmospheric wind observations sourced from two distinct longitudinal sectors.Compelling evidence reveals the statistically existence of the first six orders of sun-synchronous tidal harmonics.

Data and Methods
In this study, we analyzed zonal and meridional wind data (u and v) covering the years 2012-2019.These observations were obtained from two meteor radar stations, one situated in Mohe (M) at coordinates 122°E, 53.5°N, as detailed in Yu et al. (2013), and the other located in Collm (C) with coordinates 13.0°E, 51.3°N, as documented in Jacobi (2012) and Stober et al. (2021).
Hourly wind estimates were derived within a vertical altitude range (h) spanning from 80.5 to 95.5 km, following the methodology outlined in Hocking et al. (2001).These wind estimates were obtained from meteor observations collected across an all-sky range, approximately a 400-km diameter circle centered on the radar antenna.It is noteworthy that the spatial resolution of our data exceeds the latitudinal separation of 2.2°between the two radar stations, making the latitudinal separation inconsequential for this study.
We emphasize that our analysis is specifically focused on altitudes below 96 km.Beyond this altitude, tides may undergo distinct processes, including thermospheric absorption of solar radiation within the Schumann-Runge bands and Schumann-Runge continuum, which are phenomena prevalent around the mesopause region.Additionally, dissipation associated with ion drag becomes increasingly significant at these higher altitudes (Becker, 2017), further justifying our limitation to altitudes below 96 km.

Power Spectral Density Function (PSD) of Meteor Winds
The power spectral density function (PSD) of a time series x(t) is defined as: Here f σ ≔ 1 Δt measures the spectral frequency resolution, where Δt ≔ t max − t min represents the duration or temporal span of the time series; F{x} ≔ ∫ tmax t min x(t)e − i2πf t dt Δt denotes the Fourier transform of x.The "2" in the numerator accounts for both the real and imaginary parts contributing equally.(Alternatively the PSD can be expressed using the sampling frequency f s and data length L as Using Equation 1, we calculate the PSD for both wind components u and v collected from 2012 to 2019 at various altitude indexed as j = 1, 2, …, J, and each longitude indexed as k = 1, 2, resulting in S{u j,k } and S{v j,k }.The sum of these two PSDs is further averaged across all altitude levels and both longitudinal sectors, as described by: (2)

Power Law and Signal-To-Noise Ratio
The obtained P( f ) in the frequency range of 2 < f < 7 cpd undergoes a robust least-square regression analysis utilizing the model P( f ) = Cf β , facilitating the estimation of Ĉ and β.The robust regression employs Bisquare weights to minimize the impact of outliers.Subsequently, P( f ) ≔ Ĉf β is presented as the solid red line in Figure 1a.
Our algorithms for estimating P( f ) and β are expected to be accurate, as they have been validated by detecting a f − 5/3 power law in a different scenario, as shown in Figure 6 in He, Vogt, et al. (2017).
Referring to P( f ) as the measure of background noise power, we define a signal-to-noise ratio as SNR ≔ P P. The 99th percentile of the SNR distribution, denoted as the dashed line in Figure 1, serves as the α = 0.01 significance level.

Cross-Spectral Analysis
We calculate the cross spectra of horizontal winds between the two longitudes following the same approach as used in plotting Figure 3 in He, Chau, et al. (2020) and Figure 4 in He et al. (2021).To achieve this, cross-wavelet spectra are calculated for each altitude and for both zonal and meridional wind components.These spectra are then summed between the two components and sampled at a period of T = 1 day.Combining the samples from all altitudes allows for the compositing of a cross-spectrum for the harmonic in the date-altitude representation, covering the altitude range of 80.5-95.5 km and spanning the years 2012-2019.This procedure is repeated to composite spectra at other harmonic periods (T = 1 n day, where n = 2, 3, …, 7).The resulting seven spectra are presented in Figures 2a-2g.The grayscale intensity of the plot denotes the spectral magnitude while the color denotes the phase difference between our Mohe and Collm stations, which is a function of the zonal wavenumber m of the underlying wave and longitudinal separation between our Mohe and Collm stations, Note that the spectra in Figure 2 have been customized individually so that in each panel, green represents the argument arg{e i ⋅ nλ Δ } which equals the remainder of nλ Δ 2π and therefore corresponds to a sun-synchronous tide, see Section 2.5.

Composite Analysis
Each of 8-year spectra, illustrated in Figure 2, is averaged according to the day of year.The resulting averaged spectra, displayed in the left panels of Figure 3, serve to represent the seasonal variation.
Subsequently, the spectra from Figure 2 are averaged for another round, this time with respect to the days leading up to the nearest polar vortex weakening event (PVW), as specified in Table S1.The PVW serves as an indicator of the stratosphere sudden warming event (SSW), and its determination is based on the analysis of westward wind patterns at an altitude of 48 km and at 70°N latitude, utilizing MERRA reanalysis data, as documented in Siddiqui et al. (2015).The PVW average, presented in the right panels of Figure 3, is interpreted as the response to the SSWs within the scope of this study.

Phase Difference Technique
As demonstrated in a series of studies (e.g., He, Chau, et al., 2020;He & Forbes, 2022), the cross spectrum allows us to estimate the zonal wavenumber (m) of the underlying wave using the phase difference technique (PDT) as follows: Here, C represents the cross spectra displayed in Figures 2 and 3; Z ∈ Z is an integer indicating the whole-cycle ambiguity corresponding to the aliasing associated with Nyquist's theorem in space.The PDT entails both singleand long-wave assumptions (e.g., He & Forbes, 2022).The former implies the existence of a dominant wave within the defined time and frequency window, while the latter necessitates the wavelength of this dominant wave to be longer than twice λ Δ to ensure Z = 0.However, prior knowledge, such as m ∈ Z, can relax the long-wave assumption (He & Forbes, 2022).Readers are referred to He (2023) for a plain-language introduction to the principles and assumptions of the PDT, as well as a practical demonstration using synthetic data.In the current work, we relax the long-wave assumption to the constraint m ∈ {n, n ± 1, n ± 2} in which m = n corresponds to the sun-synchronous component and m ∈ {n ± 1, n ± 2} embraces all sun-asynchronous (i.e., non-migrating) components potentially originating from nonlinear interactions between sun-synchronous tides and stationary planetary-scale patterns with zonal wavenumber |m| = 0, 1, or 2.
Under this constraint, and considering the green customization of the argument arg{e i ⋅ nλ Δ } as defined in Section 2.3, the green spectra in Figures 2 and 3 correspond to estimations of m = n because Displayed are the sum of cross-spectra of u and v.The grayscale intensity and color in each panel denote the magnitude and phase of the spectrum.The phase is a function of zonal wavenumber m as specified in the color code.The color hue is adjusted so that in all panels, green denotes the corresponding sun-synchronous component.The vertical dash lines indicate the SSW central day referring to PVWs as specified in Table S1.In each panel, the solid lines represent a single isoline, as specified by the horizontal line in the corresponding colormap.
Figure 3. (a-g) Seasonal variation of 24-, 12-, 8-, 6-, 4.8-, 4-, and 3.4-hr cross-spectra composited from the spectra in Figures 2h-2n A similar composite as (a-g) but with respect to the SSW central day.In (a-g), the dashed and dotted lines display the average and standard deviation of the day specified in the third column in Table S1.
Note that in (a-g), the composites run from July to June so that the winter maxima in (b-f) are located in the center.In each panel, the solid lines represent a single isoline, as specified by the horizontal line in the corresponding colormap.
The tidal components with m = n share the same zonal phase velocity v p = 1 mT = n m * 1cpd = n m * Ω where Ω ≔ 1cpd is the apparent motion of the sun.These components are inherently synchronized with the sun's motion and are collectively referred to as migrating tides.

Power-Law of the PSD
The power spectral density (PSD) of middle atmospheric horizontal winds from 2012 to 2019, as depicted in Figure 1a and detailed in Section 2.1, exhibits a power-law frequency pattern with PSD ∝ f β within the frequency range of 2 < f < 7 cycles per day (cpd).This range falls between the local inertial frequency f c = 1.6 cpd and the local Brunt-Väisälä frequency f B ≈ 300 cpd (Wüst et al., 2017).In the range from f c to f B , atmospheric fluctuations display universal power spectra in terms of frequency and horizontal wavenumber (e.g., Tulloch & Smith, 2006;VanZandt, 1982), characterized by a canonical spectral slope of − 5/3.However, in contrast to this expected − 5/3 slope, our PSD reveals a power-law exponent of β = − 2.98 ± 0.02 at 2 < f < 7 cpd.This estimation is derived through a least-squares (LS) regression analysis, as outlined in Section 2.2.We attribute this discrepancy between our estimation and the canonical slope to the inherent limitations of the horizontal resolution of our wind observations.
The distribution of kinetic energy as a function of horizontal wavelength has been extensively investigated, showing a slope of − 3 at synoptic scales and − 5/3 at mesoscales (e.g., Callies et al., 2014;H. L. Liu, 2019;Nastrom et al., 1984;Tulloch & Smith, 2006).While the − 3 power law aligns with classical two-dimensional turbulence theory (Kraichnan, 1967), the − 5/3 sub-range on the shorter scale side of the − 3 sub-range remains a subject of debate (e.g., Koshyk & Hamilton, 2001;Tung & Orlando, 2003).Simulations with a general circulation model at a horizontal resolution of 35 km (km) revealed that the transition between these two slopes occurs at wavelengths ranging from 500 to 4,000 km, with variations depending on altitude (Koshyk & Hamilton, 2001).The horizontal resolution of our observations, approximately 400 km (see Section 2), limits our ability to fully capture mesoscale fluctuations, emphasizing the presence of the synoptic − 3 power law.
The PSD in Figure 1a displays distinct peaks at frequencies f = 1, 2, …, 6cpd, with signal-to-noise ratio (SNR) exceeding 14 dB above the background PSD estimated through LS regression (represented by the red solid line in Figure 1b).These six peaks significantly exceed the α = 0.01 significance level, as defined in Section 2.2.Given that the spectral cascade should not favor any specific frequency, we suspect that these subharmonics correspond to tidal phenomena.It is worth noting that existing literature has extensively reported the first four harmonics while largely overlooking the 5-th and 6-th (4.8-and 4-hr) harmonics.We attribute this oversight to the common use of frequency-independent white noise models for assessing spectral significance (e.g., He, Forbes, et al., 2020).In accordance with the white noise model, represented by the dashed blue line in Figure 1a, the 5-and 6-cpd peaks do not surpass the α = 0.01 significance level and are thus often overlooked.We therefore recommend employing a frequency-dependent noise model when assessing the significance of tidal spectral components.

Tides Versus Regional Buoyancy Waves
The cross-spectra corresponding to the six harmonics, as depicted in Figures 2a-2f, are dominated by green peaks.The persistent presence of this color within these cross-spectra indicates coherent oscillations between the two longitudinal sectors housing our radars.The coherence implies that the underlying wavelength is approximately of the same order of magnitude as the longitudinal separation which is 109°.Consequently, it can be deduced that these coherent oscillations are not associated with regional buoyancy waves.Specifically, the prevalence of the green color suggests that the coherence is linked to sun-synchronous tides, as explained in Sections 2.3 and 2.5.Regarding the 7-cpd harmonic, even though it marginally surpasses the significance level, as illustrated in Figure 1, its cross-spectra, depicted in Figures 2g and 3g, do not exhibit any predominant color, and therefore, further discussion on this harmonic will not be pursued.

Seasonal Incidence Versus SSW Association
The 8-year cross-spectra presented exhibit annual recurrences, which are notably evident in the seasonal variations depicted in Figures 3a-3f.Among these seasonal variations, the most pronounced patterns include winter enhancements in December-January observed in the 2-, 3-, 4-, 5-, and 6-cpd harmonics (Figures 3b-3f) and equinox enhancements observed in the 1-, 2-, 3-, and 5-cpd harmonics (Figures 3a, 3c, and 3e).The equinox enhancements can be attributed to the prominent semiannual variation in tropospheric latent heat release, with a peak during the equinox months, which has previously been used to explain the seasonal variability of diurnal sun-synchronous tide (Hagan & Forbes, 2002).
As for the winter enhancements, two potential factors come into play.First, there is the enhanced forcing of ozone heating in the winter hemisphere (e.g., Hagan et al., 1999;Schneider et al., 2005;Xu et al., 2012).Second, the influence of mean winds on dissipation plays a role.The prevailing eastward wind in the winter atmosphere reduces susceptibility to dissipation of the sun-synchronous tides by increasing their vertical wavelengths through Doppler shifts (e.g., Becker, 2017;Forbes et al., 2020;Forbes & Vincent, 1989).The reversal of the prevailing eastward wind to westward can account for the termination of the winter enhancements.This Doppler shift effect of zonal wind has also been utilized to elucidate the modulations of the quasi-biennial oscillation and El Ni ǹo-Southern oscillation on tidal amplitudes (e.g., Kogure & Liu, 2021;H. Liu et al., 2017;Pancheva et al., 2021;Xu et al., 2009).
While the first four tidal harmonics have been reported to be active during Sudden Stratospheric Warming events (SSWs) (e.g., Fuller-Rowell et al., 2011;Goncharenko et al., 2012;He & Chau, 2019;He, Forbes, et al., 2020;Pedatella et al., 2014;Smith, 2012;van Caspel et al., 2023), our results reveal that the winter enhancements of the 2-6-th harmonics appear to be independent of the occurrence and timing of SSWs.For instance, in the winters of 2013/2014, 2014/2015, and 2015/2016 (Figure 2f (Figure 2f), the enhancements of the 6-th harmonic are observed in the absence of, after, and before the Polar Vortex Weakening (PVW), respectively.(The PVW is a reference day of the SSW.See Section 2.4.).Although the composite responses of the harmonics during SSWs (Figures 3h-3n) display green peaks around the PVW, these responses statistically resemble the winter enhancements, as SSWs themselves are predominantly a winter phenomenon occurring mostly in January-February.It is challenging for observational studies to effectively disentangle the influence of SSWs from seasonal variations.

On the High-Order Harmonics
Numerous experimental and modeling endeavors have been dedicated to the exploration of 1-4-cpd tidal components (e.g., Azeem et al., 2016;Du & Ward, 2010;Jacobi et al., 2017Jacobi et al., , 2018;;Lilienthal & Jacobi, 2019;G. Liu et al., 2020;M. Liu et al., 2015;Moudden & Forbes, 2013;Pancheva et al., 2013Pancheva et al., , 2021;;Xu et al., 2012;Yue et al., 2013).Tidal models successfully capture the distinctive features of these sun-synchronous components (Hagan et al., 1995(Hagan et al., , 1999(Hagan et al., , 2001)), when forced by solar heating alone.This statement pertains to observations made below about 100 km.It might be tempting to consider the study of higher-order tides as an esoteric pursuit with little practical significance, given the conventional belief that the diurnal variation of solar heating is adequately described by Fourier series coefficients in the range of 1-6-cpd, with descending amplitudes.However, Smith and Ortland (2001) and Smith et al. (2004) argue that nonlinear tide-tide interactions also contribute to the forcing of 3-and 4-cpd tides under specific conditions.These tides' activities are explained in terms of spatial and temporal variations of heat sources, zonal-mean winds, and nonlinear interactions with gravity waves and planetary waves (e.g., Forbes, 1982;Hagan, 1996;He, Chau, et al., 2017;McLandress, 2002;Mlynczak & Solomon, 1993;Xu et al., 2010).
Recent literature increasingly supports the notion that nonlinear tide-tide interactions play a crucial role in generating high-order solar tides, attaining significant amplitudes in the E-region (100-150 km altitude).This implies the potential generation of dynamo electric fields that map to the F-region and drive variability in E × B plasma drifts and densities at approximately 4-6 hr time scales.A key observation in this context is presented by Azeem et al. (2016), who reports 4-cpd tidal temperature amplitudes reaching 10-30 K at 120-130 km within ±30°latitude.This contrasts with 2-and 3-cpd amplitudes of approximately 20-40 K and 20-30 K, respectively. Angelats I Coll and Forbes (2002) provide a theoretical basis for this observation, who report 4-cpd tidal meridional wind amplitudes up to 27 ms − 1 at 130 km and 20°S latitude due to 2-cpd tidal nonlinear selfinteraction according to their simulation.Support for tidal self-interaction as a viable mechanism exists in simulations by Huang et al. (2007), who demonstrate that an equatorial 2-cpd tidal temperature amplitude of 18 K at 106 km results from 1-cpd tidal self-interaction.The 3-cpd tide is generated, in part, by the 1-cpd and 2-cpd tidal interactions (see the summary in Moudden & Forbes, 2013).Therefore it is reasonable to assume that 4-, 5-, and 6-cpd tides may also arise partially through tidal interactions.

Conclusions
In summary, this study presents six tidal harmonics derived from 8-year middle atmospheric wind observations obtained by two mid-latitude meteor radars separated by 109°in longitude.The observed coherence between these radars strongly suggests that these harmonics are predominantly associated with sun-synchronous tides, rather than sporadic regional buoyancy waves.Notably, the 5-and 6-th (4.8-and 4-hr) tidal harmonics have been largely overlooked in previous investigations, primarily due to the prevalent use of inappropriate frequencyindependent noise models.In contrast, our approach utilizes a noise model characterized by a power law with a frequency dependence proportional to f − 3 , as determined from the observational PSD diagram.
Atmospheric thermal tides play a crucial role in shaping celestial bodies, influencing planetary habitability, climate fluctuations, geophysical activities, and heat and mass circulation.Hence, our study initiates exploration into the broader significance of high-order solar thermal tides, underscoring the potential for further research to unveil their comprehensive implications for planetary systems.

•
Wind spectrum reveals 6 tidal harmonics significantly higher than background noise with − 3 frequency power law • Coherence between two longitudinal sectors reveals that the harmonics are synchronized with the Sun • Winter tidal enhancements seem to be influenced by seasonal factors rather than SSW, presenting a contrast to existing literature Supporting Information: Supporting Information may be found in the online version of this article.

Figure 1 .
Figure 1.(a) Power spectral density function (PSD) of horizontal wind within 2012-2019 averaged geometrically across all altitude levels and two longitude sectors, according to Equation 2. (b) Signal-to-noise ratio (SNR) with reference to the noise level defined by the solid red line in (a) that is least-squares fitted according to the model P ∝ f β in the frequency range 2 < f < 7 cpd.The blue dashed and red dashed lines indicate the α = 0.01 significance level referring to a frequency-independent and -dependent noise model, respectively.(See Sections 2.1 and 2.2).In (a), the label S n (n = 1,2,…, 7) denotes the period T = 1 solar day n .Note that the spectra exhibit crowding as a consequence of the heightened frequency resolution associated with the 8-year sampling windows: f σ ≔ 1 Δt = 1 8 years = 3.4 × 10 − 4 cpd.