Planetary‐Scale Waves Seen in Thermal Infrared Images of Venusian Cloud Top

Planetary‐scale waves are thought to play a crucial role in the maintenance of the atmospheric circulation of Venus. They have been mainly studied using ultraviolet images. In this study, cloud top temperature oscillations associated with planetary‐scale waves were extracted from thermal images taken over 142 Earth days by the Longwave Infrared Camera onboard JAXA's Akatsuki spacecraft. A spectral analysis showed that multiple planetary‐scale waves overlap. The latitudinal structures of four distinct waves with periods of 3.6, 5.0, 5.4, and 6.1 days were studied. The 3.6‐day wave is considered to be a Kelvin wave and the other waves are considered to be hemispherically symmetric Rossby waves. Although the Rossby waves have similar periods, their latitudinal structures are different. The Rossby waves are considered to originate from different sources based on a comparison with analytical solutions.

• Planetary-scale waves were observed at the Venusian cloud top in thermal infrared images taken by the Venus orbiter Akatsuki • Four distinct modes were identified in the spectra; one is attributed to a Kelvin wave and the others are attributed to Rossby waves • The Rossby waves are thought to originate in different regions based on a comparison with analytical solutions

Supporting Information:
Supporting Information may be found in the online version of this article.
The brightness temperature images of the cloud top taken by LIR potentially provide complementary information on the structures of planetary-scale waves. However, the images taken by LIR are relatively featureless, and thus the extraction of planetary-scale waves other than the thermal tides (Fukuya et al., 2021;Kouyama et al., 2019) and stationary topographic waves Kouyama et al., 2017) has been unsuccessful. The Pioneer Venus infrared spectrometer data show the existence of temperature oscillations with periods of 2.9-5.8 d around the cloud top and above in the northern hemisphere (Apt & Leung, 1982;Taylor et al., 1980). Precise measurements of the wave periods and the extraction of the waves' latitudinal structures from the temperature field are desired. In this study, we extract temperature oscillations associated with planetary-scale waves from LIR data using a newly developed method. The complete local time coverage of LIR combined with longterm observation enables unprecedented continuous sampling, which is advantageous for precise spectral analysis. The latitudinal structures of the detected waves are compared with analytical solutions of equatorial waves to identify the wave sources.

Data Set
We used global images of Venus taken by LIR. The details of LIR are described in Fukuhara, Taguchi, et al. (2017) and Taguchi et al. (2007). LIR measures infrared radiation in the wavelength range of 8.0-12.0 μm, where the influence of sunlight scattered by Venus is negligible, and thus all local time regions can be observed. The observed radiance is expressed by the brightness temperature, which is interpreted as the temperature of the cloud top. The contribution function of the radiation measured by LIR has a maximum at an altitude of around 65 km with a full width at half maximum of ∼10 km. LIR has a field of view of 16.4° × 12.4° and captures images with a detector format of 328 × 248 pixels. The pixel resolution is 0.6°-3.0° in both latitude and longitude for the data used in this study. The systematic error, which appears as an offset between images, is ∼3 K. The relative error among the pixels in a given image is ∼0.3 K. The data used are LIR Level 3x products in the Akatsuki Data Archive Ogohara et al., 2017).
The Akatsuki spacecraft orbits the equatorial region with a period of ∼10.5 d. Because of its extended elliptical orbit, the spacecraft spends most of the time near the apoapsis. The direction of the major axis of the orbit is almost unchanged in the inertial frame. From this orbit, LIR acquires images every 2 hr. Observations are suspended for 8 hr each day during telecommunication with the tracking station. For about one month or more every Venus year (∼225 d), when the apoapsis is located above the midnight region, stray light from the Sun prevents observations. Because data should be as regular and continuous as possible for a spectral analysis, 272 images taken at roughly 12-hr intervals (UTC 05:30 ± 2 hr and 17:30 ± 2 hr) were extracted from the 142-d period from May 18 to 6 October 2017. This period had 12 data gaps and five images taken outside the range above were included. The data are listed in Table S1. To continuously observe the same side of the planet in the inertial frame of reference, the right ascension of the spacecraft direction as seen from the center of Venus was limited to the range 190°-270°, which includes the apoapsis (Figure 1). The standard deviation of the right ascension is ∼12°, which corresponds to a phase difference of 0.034 cycles for a wave with a zonal wavenumber of unity. The influence of the uneven sampling and the movement of the spacecraft position is examined later. The frequency resolution in the later spectral analysis is given by the reciprocal of the length of the time series, which is 1/(142 d) ∼ 0.007 d −1 .

Method
Planetary-scale waves are expected to be observed as periodic variations in the brightness temperature of the cloud top. However, the change in brightness temperature from one image to another is difficult to measure because of the systematic error of the measurement (∼3 K), which is comparable to or larger than the wave amplitudes in a model that reproduces the observed velocity amplitudes (Kouyama et al., 2015) and the amplitudes suggested by infrared spectroscopy (Taylor et al., 1980). Furthermore, the statistical properties of the systematic noise are unclear. To deal with this difficulty, we focus on the brightness temperature gradient in the east-west direction in each image, which is affected only by the relative error between pixels (∼0.3 K). The zonal propagation of a planetary-scale wave causes an oscillation of the longitudinal temperature gradient, which can be detected by finding a peak in the spectrum of the temperature gradient time series. Note that the observed altitude region depends on the emission angle; a larger emission angle leads to a higher peak altitude of the contribution function Sato et al., 2014). Therefore, pixels belonging to a given emission angle range were used in the calculation of the longitudinal temperature gradient. The procedure used to obtain spectra is as follows.
1. In each image, 11 latitudinal bins were set from 50°S to 50°N in intervals of 10°. Along each latitudinal circle, the average of the brightness temperature in the emission angle range of 40-60° was individually calculated for the eastern and western sides of the Venus disk ( Figure 2). Typically, 10-20 data points were averaged in the emission angle range. The longitudinal gradient of the temperature was obtained by dividing the difference in the mean temperature between the eastern and western sides by the longitudinal angular distance. The longitudinal angular distance ranges from around 50° to 100°. The above procedure was repeated for all images to obtain the temperature gradient time series at the 11 latitudinal bins (Figure 3). Although temperature gradient data exist for every ∼12-hr period, out of 284 data points, 12 data points are missing at each of 50°N, 40°N, 30°N, 20°N, 10°N, 0°, and 10°S, 14 data points are missing at 20°S, 13 data points are missing at 30°S, 19 data points are missing at 40°S, and 46 data points are missing at 40°S. 2. Each time series included a prominent long-time-scale component with a time scale of 50-100 d, which is considered to have been caused by the combination of the local-time-fixed structure (thermal tides) and the change in the local time of the apoapsis in the course of Venus's revolution. In order to remove artificial spectral features that appear when Fourier transform is applied to a finite-length time series, this component  was suppressed by fitting a third-order polynomial to the data and subtracting it from the original time series. It was confirmed that the wave properties did not change when the order of the polynomial was changed to 4 or 5. 3. The time series at all latitudinal bins were spectrally analyzed using the discrete Fourier transform and the Lomb-Scargle periodogram (Scargle, 1982). The purpose of Fourier transform is to apply narrow-band filtering for the estimation of the amplitude and phase (see Section 5).
In the Fourier transform, the sampling interval was regarded as 12 hr and the data gaps were filled by linear interpolation in advance. In the Lomb-Scargle periodogram, the original data with uneven sampling were used. The effect of data gaps on the Lomb-Scargle periodogram was examined by analyzing sinusoidal functions with an actual sampling window; it was confirmed that no spurious peaks appear. The Python functions numpy.fft.fft and signal.lombscargle were used for the Fourier transform and the Lomb-Scargle periodogram, respectively.
The influence of the uneven sampling and the orbital motion of the spacecraft around Venus was found to be negligible according to a simulation. We let a sinusoidal wave with a zonal wavenumber of unity propagate westward along the equator with a period of 5.071 d in the inertial frame, which coincides with one of the frequency bins in the Fourier spectrum. Then, we calculated the temperature gradient at each actual observation time based on the difference between the two longitudinal regions with an emission angle of 50° as seen from the spacecraft position (Figure 4a), although the spacecraft position was projected onto the equatorial plane for simplicity. The Fourier transform of this time series showed a distinct peak at a frequency that corresponds to the assumed wave period (Figure 4b). The spectral amplitudes for neighboring bins separated by a spectral resolution of 0.0070 d −1 dropped to almost zero. Simulations with other wave periods gave similar results. The apparent absence of the influence of the orbital motion can be explained by the fact that almost the same side of the planet continues to be observed due to the extended elliptical orbit. Although the observed phase is disturbed by up to ∼0.10 cycle when the spacecraft approaches Venus, the phase development is quite regular in most of the observation period. The planetary rotation relative to the major axis of the spacecraft's orbit, which is almost fixed to the inertial frame, causes the overestimation of the observed frequency by the reciprocal of the orbital period of 243.0 d, that is, 0.0041 d −1 . The observed frequencies were corrected by this small amount.

Spectra of Longitudinal Temperature Gradient
The latitudinal dependence of the Fourier spectrum is shown in Figure 5a and that of the Lomb-Scargle periodogram is shown in Figure 5b. The frequency interval of the Lomb-Scargle periodogram can be chosen arbitrarily. Since the effective frequency resolution is considered to be similar to that of the Fourier spectrum, setting the interval much smaller than this might cause artificial features. On the other hand, if the frequency interval is set equal to that of the Fourier spectrum, the heights of spectral peaks can be underestimated if the peaks are located right in the middle of the frequency bins. Considering these, we set the interval to half that of the Fourier spectrum as a compromise. The recurrence period of the background wind at the cloud top, which is also plotted in Figure 5, was calculated with the following model zonal velocity ( Figure A1): where 0 = 95 m s −1 , 1 = 40 m s −1 , and is the latitude. This model profile is roughly aligned with the mean zonal velocity profile obtained by cloud tracking from 365-nm images taken by Akatsuki UVI during March 2017 to January 2019  and reproduces the observed equatorial velocity of ∼100 m s −1 and broad midlatitude jets at around 45° latitude (e.g., Hueso et al., 2015;Khatuntsev et al., 2013;Machado et al., 2017;Peralta et al., 2007). The 95% significance interval is plotted based on the 2 -distribution of the red noise spectrum (Gilman et al., 1963). The spectra plotted in a wider spectral range ( Figure B1) show, in addition to peaks with periods of several days under study, a broad peak around the orbital period of ∼10.5 d, which might be created by the combination of the local-time-fixed structure (thermal tides) and the spacecraft motion. Figures 5a and 5b show similar spectral features. Several narrow spectral peaks extend over a wide latitudinal range, which are indicative of planetary-scale waves. Oscillations with periods shorter than that of the background wind tend to be confined to low latitudes. Assuming that the oscillations have a zonal wavenumber of unity based on the previous studies (Section 1), this implies that the oscillations are manifestations of Kelvin waves. Oscillations with periods longer than that of the background wind show large amplitudes at mid to high latitudes. This implies that the oscillations are manifestations of Rossby waves. The widths of the spectral peaks are comparable to the resolution of the discrete Fourier transform, which is ∼0.11 d for a 4-d wave and ∼0.25 d for a 6-d wave. This suggests that the frequencies of the waves were stable over tens of days.
The similarity of the two spectra suggests that the effect of the uneven sampling is negligible; the Lomb-Scargle periodogram provides slightly more resolved spectra. Based on this result, in the next section we obtain the amplitude and phase structures for the distinct modes indicated by the arrows in Figure 5 by applying the inverse Fourier transform to the narrow spectral intervals of the Fourier spectra.

Latitudinal Structures of Observed Waves
Here, we focus on the four distinct spectral peaks indicated by arrows in Figure 5. The corresponding waves, which have periods of 3.60 ± 0.09, 4.99 ± 0.18, 5.37 ± 0.20, and 6.06 ± 0.25 d, are hereafter referred to as the 3.6-d wave, 5.0-d wave, 5.4-d wave, and 6.1-d wave, respectively. The uncertainty in the period was taken to be the interval between the frequency bins in the discrete Fourier transform. The 5.0-d and 5.4-d waves are considered as different spectral modes since their frequencies are separated by twice the spectral resolution and the waves have different latitudinal structures. To extract each quasi-monochromatic wave from the time series, we applied the inverse Fourier transform to the corresponding spectral components. Because each spectral peak typically includes three frequency bins, the inverse Fourier transform was applied to the complex Fourier amplitudes at the peak frequency bin and its two neighboring bins. Examples of the time series obtained by this procedure are shown in Figures C1 and C2; the extracted dominant mode reproduces to some extent the observed notable oscillations in the period where the mode is enhanced relative to the other periods.
The amplitude at each latitude was obtained from the root-mean-square of the quasi-monochromatic component multiplied by √ 2 . The phase at each latitude was calculated by dividing the time lag between the quasi-monochromatic components at the relevant latitude and at the equator by the wave period. The time lag was estimated by finding the maximum time-lagged correlation coefficient between the time series. In the plots shown below, the phase at the equator was taken to be . The time interval used for computing the amplitude and phase was taken to be the interval where each wave has a relatively large amplitude, namely 50-120 d for the 3.6-d wave, 50-140 d for the 5.0-d wave, 0-142 d for the 5.4-d wave, and 0-80 d for the 6.1-d wave.
The obtained latitudinal structures are shown in Figure 6. We present the amplitude in the temperature, ̂ , instead of that in the temperature gradient, ̂ . ̂ was calculated using the relation ̂= (2 ∕ )̂ , where is the zonal wavelength of a wavenumber-1 wave at each latitude. The sign of the phase is such that an increase (decrease) of the phase represents a westward (eastward) phase shift. The amplitude and phase are not plotted at 50° latitude for the 3.6-d and 6.1-d waves because the spectral peak was low compared to the background component at high latitudes and thus not statistically significant. The phase is nearly constant with respect to the latitude for all waves, suggesting that zonally propagating modes were acquired. The amplitudes are around 1 K or smaller, in accordance with a model that reproduces the observed velocity amplitudes (Kouyama et al., 2015). The error bars 10.1029/2021JE007047 7 of 13 in the amplitude were calculated from the root-mean-square of the amplitudes of the background noise outside the spectral peaks. The error in the phase was calculated as where is the amplitude of the wave (signal), is the amplitude of the noise, and is the instantaneous phase difference between the signal and the noise. The instantaneous phase error is given by ∕ sin , where is considered to be distributed uniformly in the range 0-2 .
The 3.6-d wave shows a slight increase in amplitude near the equator. This, combined with faster propagation compared to the background wind, suggests that the wave is a Kelvin wave. Similar low-latitude modes with periods of 3.5-4.0 d have been found in the UV brightness variation and the zonal velocity variation deduced from UV images (e.g., Imai et al., 2019;Kouyama et al., 2012;Nara et al., 2019;Rossow et al., 1990).
The amplitude of the 5.0-d wave increases with latitude up to 50°S and 50°N and seems to reach its maximum value at higher latitudes in both hemispheres. The 5.4-d wave shows similar latitudinal variation, although the amplitude may reach its maximum value at 40° in the southern hemisphere. The 6.1-d wave seems to have maximum amplitudes at around 20-30°. These characteristics, combined with slower propagation compared to the background wind, imply that the waves are hemispherically symmetric Rossby waves. Similar midlatitude modes with periods of 5.0-5.5 d have been found in the UV brightness variation and the meridional velocity variation deduced from UV images (e.g., Imai et al., 2019;Kouyama et al., 2015Kouyama et al., , 2013Rossow et al., 1990). Imai et al. (2019) identified Rossby-wave-like vortices centered at ∼40° latitude. The 5.0-d and 5.4-d waves found in our study might correspond to such midlatitude modes. The 5.0-d and 5.4-d waves might also correspond to the temperature oscillation with a period of 5.3 ± 0.4 d that peaks at 60-70° observed in the northern hemisphere by the Pioneer Venus infrared spectrometer (Apt & Leung, 1982).

Discussion and Conclusions
We found for the first time that at least four planetary-scale waves exist simultaneously at the Venusian cloud top around 65 km altitude. As described above, the 5.0-d, 5.4-d, and 6.1-d waves are considered to be Rossby waves based on their periods and latitudinal structures. The 5.0-d and 5.4-d waves have maximum amplitudes at a latitude of around 50° or higher, and the 6.1-d wave has a maximum amplitude at around 20-30° latitude. In this section, the latitudinal structures are compared with analytical solutions of equatorial waves. The energy of any free mode in spherical geometry becomes concentrated toward the equator as the Lamb parameter increases (Longuet-Higgins, 1968). Equatorial wave theory can thus be a useful approximation for modes with small equatorial deformation radii. The latitudinal gradient of the Coriolis parameter, , that appears in the analytical solution is replaced by the background vorticity gradient.
Because the analytical solution for an equatorial wave (Matsuno, 1966) cannot take into account the latitudinal variation of the angular velocity of the background wind, we selected two representative latitudes for the background rotation, namely 20° and 50°. The former corresponds to the lowest latitude of the maximum amplitude and the latter is the highest latitude where the amplitude was retrieved. Waves that are consistent with the background wind in the midlatitude should be reproduced in this parameter range. The recurrence period and the vorticity gradient of the background wind were calculated for these latitudes using the velocity profile (Equation 1) ( Figure A1). The intrinsic frequency of the wave was calculated from the observed wave's period and the recurrence period of the background wind. We first determined the vertical wavenumber m corresponding to the wave's period, namely 5.0, 5.4, or 6.1 d, for the gravest equatorial Rossby wave (n = 1) by solving the following dispersion relation (Andrews et al., 1987): where is the intrinsic frequency, N is the buoyancy frequency, and k is the horizontal wavenumber. The equatorial deformation radius l E was calculated from m using the relation 10.1029/2021JE007047 8 of 13 where H is the scale height. The latitudinal profile of the temperature amplitude was obtained using the relation where y is the northward distance from the equator.
The latitudinal structures of waves with periods of 5.0, 5.4, and 6.1 d for the representative latitudes of 20° and 50° are plotted in Figure 7. The adopted parameters and the calculated vertical wavelengths are summarized in Table 1. We assumed H ∼ 5.4 km and N ∼ 0.018 s −2 based on the Venus International Reference Atmosphere (VIRA) at an altitude of ∼65 km (Seiff et al., 1985). In all solutions, the maximum amplitude occurs around latitudes 20-35° and the variation of the amplitude with latitude is less than ∼20% at <30° latitude. The observed 6.1-d wave is consistent with these characteristics, whereas the observed 5.0-d and 5.4-d waves have maximum amplitudes at latitudes of ∼50° or higher, and have amplitudes at high latitudes that are more than twice as large as those near the equator.
It is unlikely that Rossby waves of similar frequencies originating in the same region could have such different latitudinal structures. A possible explanation is that the 6.1-d wave was generated locally around the cloud top and that the 5.0-d and 5.4-d waves were generated deep in the clouds and propagated upward to the cloud top. Imamura (2006) showed, using a numerical model in a spherical geometry, that an equatorial Rossby wave forced at an altitude of 60 km expands poleward during upward propagation in the presence of midlatitude jets, and has maximum amplitudes at around latitude 60° where the vorticity gradient peaks ( Figure A1) at the cloud top level. Imamura (2006) also showed that the latitude of the maximum amplitude shifts to ∼45° in the absence of the midlatitude jets. It is probable that the 5.0-d and 5.4-d waves originated at lower altitudes and were trapped around the vorticity gradient maxima at high latitudes. The 5.0-d and 5.4-d waves might also have been generated by planetary-scale instability (Iga & Matsuda, 2005;Wang & Mitchell, 2014) or by unstable midlatitude jets (Piccialli et al., 2012) and then trapped by the jets.
In the present study, multiple global-scale waves were simultaneously observed in the spectra and at least four distinct modes were identified. The dominance of planetary-scale waves is considered to be a characteristic of the Venusian atmosphere, where the Rossby radius of deformation is large because of the slow rotation and the background atmosphere is relatively homogeneous due to the small effects of topography and the seasonal cycle. Venusian general circulation models exhibit such a superposition of planetary-scale waves (Lebonnois et al., 2016;Yamamoto & Takahashi, 2003), although the sources of the waves have not been well identified. Coupled Rossby-Kelvin modes that emerge from planetary-scale instability (e.g., Iga & Matsuda, 2005;Wang & Mitchell, 2014) might play an important role. Planetary-scale streak structures seen in middle-and lower-cloud images at near-infrared wavelengths taken by the Akatsuki 2-μm camera (Limaye et al., 2018;Peralta et al., 2018;Figure 7. Latitudinal structures of equatorial Rossby waves with periods of 5.0 d (black), 5.4 d (red), and 6.1 d (blue) for the background wind conditions at latitudes 20° (solid) and 50° (dashed). The amplitudes were normalized by the maximum amplitude for each wave. Note. Background wind conditions at latitudes 20° and 50° are considered.  Figure C1, and the sum of the 5.0-d, 5.4-d, and 6.1-d waves is shown in Figure C2.

Data Availability Statement
The Akatsuki LIR Level 3x data used in this study are available in the Data Archives and Transmission System (DARTS) of JAXA  at https://doi.org/10.17597/ISAS.DARTS/VCO-00019. An internal version of the data (v20181001) was used, but it is essentially the same as the publicly available version. The data presented in the figures and the time series of the longitudinal gradient of the brightness temperature derived from the LIR data are given in the online repository Zenodo .

Acknowledgments
The data used in this study were provided by the Akatsuki project of JAXA. This work was supported by JSPS KAKENHI grant numbers 20H01958, 19H05605, and 19K14789. The authors would like to thank one anonymous reviewer and Pedro Machado for their constructive reviews and Editor Anni Määttänen for editorial handling.