Density Fluctuations in the Lower Thermosphere of Mars Retrieved From the ExoMars Trace Gas Orbiter (TGO) Aerobraking

The upper atmosphere of Mars is constantly perturbed by small-scale gravity waves propagating from below. As gravity waves strongly affect the large-scale dynamics and thermal state, constraining their statistical characteristics is of great importance for modeling the atmospheric circulation. We present a new data set of density perturbation amplitudes derived from accelerometer measurements during aerobraking of the European Space Agency’s Trace Gas Orbiter. The obtained data set presents features found by three previous orbiters: the lower thermosphere polar warming in the winter hemisphere, and the lack of links between gravity wave activity and topography. In addition, the orbits allowed for demonstrating a very weak diurnal variability in wave activity at high latitudes of the southern winter hemisphere for the first time. The estimated vertical damping rates of gravity waves agree well with theoretical predictions. No clear anticorrelation between perturbation amplitudes and the background temperature has been found. This indicates differences in dissipation mechanisms of gravity waves in the lower and upper thermosphere.


Introduction
The Trace Gas Orbiter (TGO), launched on 14 March 2016, is the first mission of the European Space Agency's (ESA's) ExoMars program, with the primary objective to search for the evidence of methane and other trace gases. On 19 October 2016, TGO was inserted into an elliptical orbit. TGO performed an 11-month-long aerobraking campaign between 14 March 2017 and 20 February 2018 to gradually reduce its orbital period from approximately 24 h to 2 h, by passing through the Martian atmosphere 952 times at the pericenters of the initially highly elliptical orbit, thereby accumulating a total drag ∆V of about 1017 m s −1 with a limited propellant usage. An extensive overview of the TGO aerobraking design and operations is given in the paper of Denis et al. [1], while Castellini et al. [2] elaborated on this challenging campaign from a Flight Dynamics perspective. Our work exploits the fact that TGO reached the denser layers in the thermosphere and the on-board accelerometers registered variations of the atmospheric drag with high temporal and spatial resolution. This allowed for determining small-scale thermospheric density fluctuations associated with gravity waves (GWs).
GWs are ubiquitous in the Martian and all other planetary atmospheres studied so far [3]. They have been detected at various altitudes during all seasons: in the lower atmosphere using remote sensing techniques such as radio occultations [4], dayglow imaging [5], and infrared radiometer [6], and in the upper atmosphere by in situ measurements with mass spectrometers, e.g., in [7][8][9], and satellite aerobraking [10,11]. GWs are generated at all heights and their amplitudes grow exponentially with altitude upon vertical propagation to less denser layers in the middle and upper atmosphere. Waves originated in the lower atmosphere are of particular interest, because their amplitudes increase by several orders of magnitude. Observations indicate that magnitudes of GW-induced density perturbations with respect to the background density ρ /ρ reach tens of percent in the thermosphere [10,12] and, occasionally, exceed 100% as shown in Figure 2 in [7].
GWs not only continuously disturb the atmosphere, but significantly affect its large-scale dynamics and thermal state. The essential role of these waves has originally been recognized in the terrestrial atmosphere (e.g., see reviews [13,14]), and their role in the general circulation of other planetary atmospheres is increasingly appreciated [3]. The first estimates of momentum forcing created by dissipating GWs in the Martian thermosphere have been inferred from Mars Global Surveyor (MGS) and Mars Odyssey (ODY) aerobraking densities [10]. However, much of the knowledge concerning the impact of GWs on the atmosphere of the Red Planet comes from numerical simulations. Atmospheric models provide the advantage that different atmospheric layers can be modeled continuously, and interactions of various processes can be accounted for. Global simulations of GW effects in Earth's thermosphere [15], have provided motivations to search for similar GW signatures and effects in the upper atmosphere of Mars. Indeed, GWs strongly decelerate atmospheric jets and even reverse them in the Martian middle and upper atmosphere ("GW drag") [16][17][18]. They contribute to the upper-atmosphere polar warmings [19] discovered by ODY [20] and provide stability of the thermosphere by cooling it down in addition to molecular thermal conduction [19,21]. Cooling and heating rates produced by GWs alternate with height, which was initially shown in the context of Earth's thermosphere [22]. GWs have significant magnitudes, as was established from both modeling [23] and measurements [12], and therefore cannot be ignored in the energy balance. GWs also facilitate CO 2 ice cloud formation by providing a net cooling and forming pockets of cold air in the Martian mesosphere and lower atmosphere [24][25][26].
Typical horizontal scales of GWs are between several tens and several hundreds of kilometers. Sources of GWs in the lower atmosphere are highly irregular and include flow over topography, atmospheric convection, and instabilities of weather systems. Larger-scale background winds strongly influence the vertically propagating harmonics via wave refraction and filtering, thus modulating the spatiotemporal and spectral characteristics of the GW field. Some of these phenomena have been explored in the lower and middle atmosphere using a high-resolution (GW-resolving) Martian general circulation model (MGCM) [27][28][29]. Thus, GW activity in the thermosphere is a complex result of the local wave excitation and processes of wave generation and propagation from below. As vertical damping determines the rate of dynamical and thermal forcing imposed by the waves on the large-scale flow, observations can help quantify the mechanisms of their dissipation. Knowledge of the GW field is required for increasing the predictive capabilities of atmospheric models (by constraining and validating GW parameterizations), which is essential for planning and executing future aerobraking campaigns. This paper presents a new set of aerobraking density profiles and outlines the method of extracting information on small spatial-scale disturbances interpreted as GWs. The organization of the article as follows. Section 2 describes the methods of the observations used in this study; Section 3 presents the observational analysis of GW activity; and discussion and conclusions are given in Section 4.

Coverage in Parameter Space
TGO is the fourth spacecraft that was put through an extensive aerobraking campaign in the Martian atmosphere following Mars Global Surveyor (MGS) [30][31][32], Mars Odyssey (ODY) [33][34][35][36], and Mars Reconnaissance Orbiter (MRO) [37][38][39][40][41]. Mars Atmosphere and Volatile Evolution Mission (MAVEN) reached its science orbit with the pericenter altitude of~155 km without aerobraking, but also collected valuable accelerometer data about the lower altitudes during the "deep dip" campaigns when the pericenter dropped below 135 km [42][43][44]. Tolson et al. [45,46] compared the particular circumstances of these missions and their results focusing on the accelerometer data. As both smalland large-scale variations in the atmosphere show dependencies on parameters such as the local solar time (LST), solar zenith angle (SZA), solar longitude (L s ), latitude, etc., all observations that expand the scarce coverage of the whole parameter range are rather valuable. Figure 1 illustrates the evolution of some of these parameters during the TGO aerobraking operations, as also listed below.  Panel A shows that the aerobraking was interrupted after the first 125 passes during the solar conjunction period that divided the campaign into two distinct phases. Most of the braking took place in the second phase, when the orbital period initially decreased rapidly, then at a slower pace, as the orbit was becoming shorter.

•
Except for the walk-in phases when the altitude was lowered gradually, the pericenter altitude (defined with respect to the Mars reference ellipsoid [47]) was between 100 km and 112 km. In the analysis to follow, only passes with pericenter altitudes within this interval are considered.

•
The spacecraft passed through the atmosphere mostly in the southern winter hemisphere (Panel B and C). The pericenter moved into the polar region for a significant number of passes, where the atmospheric behavior markedly differs from that at lower latitudes.

•
Although the LST covered the entire day, the SZA changed within a more constrained interval, because the Sun remained below the horizon irrespective of the LST in the winter polar region (Panel D).

•
The second phase of the aerobraking can further be divided into two parts depending on whether the pericenter was moving southward or northward (Panel B). As can be seen in Panel D, these two parts are symmetric with respect to the evolution of the SZA. However, because of the polar winter, the LST evolved differently.

Acceleration Measurements
Accelerometers are mandatory sensors for aerobraking operations serving multiple purposes: monitoring the aerobraking passes, allowing for the on-board estimation of the time of passage as well as the prediction of the time of future passes, and the on-board detection of an excessive deceleration required to automatically trigger an altitude raising manoeuvre when the aerodynamic constraints are violated [1]. TGO is equipped with two inertial measurement units (IMU), each containing three ring laser gyroscopes (RLG) and three accelerometers.
The raw accelerometer data recorded by the on-board software at the 10 Hz clock frequency were fully dumped to ground as part of the telemetry information. The 10 Hz data were processed in multiple steps for operational purposes, among which the first two are relevant for this analysis [2]: 1.
The effect of angular rates (extracted from the RLG data) on the linear accelerations at the IMU location was removed, to obtain the components of the acceleration of the center of mass.

2.
The accelerometer biases averaged before and after the aerobraking passes were also removed.
The aerobraking mode of TGO was designed to avoid undesired pitch/yaw actuations during the passes using aerodynamic stability to control its attitude. Once the denser atmosphere was left in each pass, the aerodynamic torques were no longer sufficient for bringing the spacecraft back to its target attitude. The residual rates led to the increase of the off-pointing, which eventually activated the thrusters. However, our analysis concerns only the denser regions of the atmosphere and, except for five passes having taken place at higher altitudes at the beginning of the second phase of aerobraking, the thrusters were not activated in the investigated regimes in the vicinity of the pericenters. These five passes have been promptly removed from the dataset.

Processing the Accelerometer Data
The procedure described in the previous subsection yielded a semi-raw 10 Hz accelerometer data set, which was quantized due to the finite resolution of the IMU and, therefore, needed to be smoothed. The Gaussian filters have been selected for this purpose due to their superior qualities in the frequency space compared to those of a simple moving average. In the following, the 2σ values of the respective Gaussian filters are provided for easier interpretation, as the frequency response of a Gaussian filter with a given σ closely approximates the central lobe of the response of a window representing a rolling average with a width of 2σ. When this filter is applied to a Q quantity, Q 2σ stands for the respective filtered time series.
The semi-raw data is smoothed with a 2σ = 0.6 s filter (a 0.6 ). So as to analyze the perturbations in this signal, a mean background value has to be estimated. Many previous studies working with similar observations utilized an approximately 40 s-long rolling mean, e.g., [45,48]. However, as many acceleration peaks are narrow in the current data set, applying the 2σ = 40 s filter often introduces a non-negligible bias in the vicinity of the maximum acceleration, as smaller values are inherently included in the filtering window ( Figure 2). To circumvent this issue, the mean is defined with a seventh-order polynomial fit following the work of Yigit et al. [7]. The relative perturbations are then estimated by removing this background (a) from the a 0.6 signal and then normalizing the result with respect to the background (a' 0.6 = (a 0.6 − a)/a). Only those parts of the passes have been considered, in which a 0.6 is larger than 20% of the maximum acceleration so as to avoid higher numerical and fitting errors at the lower acceleration values. This restriction also constrains the altitude range to about 10 km above the pericenter. The latter makes the assumption of a constant scale height plausible within this height interval.  To evaluate atmospheric disturbances, the obtained a' 0.6 signals need further processing, as it has been found that several profiles exhibit a peculiar regular oscillation pattern (e.g., see pass number 446 in Figure 3). It is unlikely that this behavior is caused by the atmosphere. These oscillations could be removed by increasing the width of the filtering window, however, through this process information of interest might be lost. It would be difficult to argue that a particular σ value is the appropriate choice without understanding the cause of the observed oscillation.
The frequency spectra of all a' 2σ signals have been analyzed using the fast Fourier transformation (FFT). As it is seen in Figure 4, there is a pronounced peak at~0.2 Hz. This finding is in agreement with the spacecraft documentation provided by Thales Alenia Space, according to which the frequency of the first flexible mode of the solar arrays is estimated to be 0.16 Hz with an uncertainty of ±15%. Therefore, it is more appropriate to filter out the peaks in this frequency range, rather than simply increase the width of the Gaussian filter, which affects other frequencies as well ( Figure 4B). This approach was applied to all processed data presented in this paper.  The densities have been calculated from the accelerometer data using the aerodynamic parameters of TGO [1,2]. The above described procedure for acceleration was then applied for evaluating density fluctuations. For each orbit, the standard deviation of the relative densities was calculated to quantify the amplitude of GWs on a single orbit. This measure is referred to as "GW intensity" or "GW activity" in the remaining of this paper.

Estimation of the Scale Height
For all passes, the part that is relevant to our analysis occurred within a narrow altitude range around the pericenter. Therefore, following the recommendation provided in previous publications [45], we derived the density scale height H for each individual pass by fitting a vertical exponential function with constant H to the filtered density profile ρ 0.6 .

Gravity Wave Activity
Next, we present some results based on the analysis of the entire set of aerobraking passes. For convenience, we subdivided it into three phases: (1) before the conjunction (passes 1 to 125, see Figure 1); (2) after the conjunction with the pericenter moving southward (passes 134 to 629); and (3) after the conjunction with the pericenter moving northward (passes 630 to 952). In the figures to follow, the corresponding data are indicated in red, blue, and green. Figure 5 shows the distribution of relative density variations as a function of longitude. It is seen that the amplitudes of fluctuations vary between a few and~20% from pass to pass, occasionally exceeding 25%. To better visualize the cloud of individual passes, we introduced the longitudinally running mean amplitudes and plotted them with solid lines. They vary mainly between 5 and 10%, with the exception of the local maximum up to 13% around 200 • and 250 • longitude during Phase 3. It is not clear whether this maximum is associated with topographical features or caused by large amplitudes detected during a few passes. It is seen from the lower panel that the variability of density fluctuations increases when amplitudes are large. We attribute this local maximum to very large amplitudes detected over a few passes, which are shown with dots in the upper panel. If these passes are neglected, the GW intensity weakly depends on longitude.  Next, the latitudinal dependence of density fluctuations is presented in Figure 6. If the subdivision into phases is neglected, the cloud of dots is distributed very similarly to that obtained from MRO aerobraking depicted in the lower left panel of Figure 12 in [48]. This similarity can be appreciated as the coverage in the orbital parameter space is also quite similar. If the distinction among the phases is introduced, then a clear difference emerges. It is seen that the running averages corresponding to Phases 2 and 3 coincide southward of 60 • S. This is likely because these measurements have been taken over a narrow interval of times corresponding to the southern hemisphere winter (L s~1 00 • to 120 • , see Figure 1A). Northward of 60 • S, the averaged amplitudes for Phases 2 and 3 depart from each other. A closer inspection shows an enhanced variability of amplitudes with respect to the mean during Phase 3 (green line). Essentially, the increase of the mean amplitude at the northward part of the track is created by the same passes with enhanced amplitudes that manifested in the local maximum in Figure 5. The enhanced GW activity at high latitudes just over the winter polar jet has been predicted with Martian GCMs accounting for GWs e.g., Figure 3d in [26]. Thus, this finding seems to be the first observational confirmation for the prediction.

GW intensity [%]
P hase 1 P hase 2 P hase 3 Figure 6. The same as in Figure 5, except for the dependence on latitude. The running means have been calculated by averaging over passes with pericenter latitudes positioned within the ±2 • interval around each pass. No more than 40 passes closest in time have been considered in order to avoid averaging over hundreds of passes at high latitudes.
More insight into the distribution of GW activity in the parameter space can be gained by considering its dependence on the local solar time in Figure 7. It is seen that the magnitudes of fluctuations during the pericenter southward of 60 • S (denoted by the dashed line) remained flat over a broad range of LST-between 5 and 16 h. This indicates a very weak diurnal variation of wave activity at high latitudes of the southern hemisphere in winter. GW propagation is tightly controlled by the background atmospheric fields and, especially, by the mean zonal wind. Thus, the lack of diurnal variations of GW intensities can be linked to the weak solar tide at the winter polar latitudes [49]. The behavior changes northward of 60 • S and differs for the two phases. The activity is stronger from the nighttime to early morning hours during Phase 3 (green line in Figure 7A) at mid-latitudes of the southern hemisphere. These are the passes with large but highly variable amplitudes of disturbances ( Figure 7B) discussed above. The latter result is consistent with MAVEN-NGIMS measurements, which showed that the nighttime GW activity is larger than the daytime one [7]. Note that the MAVEN-NGIMS study focused on higher altitudes than those of the TGO aerobraking campaign. During Phase 2 (blue lines), the activity drops in the late evening hours.  It is seen that these measurements at pericenters at lower latitudes were taken much earlier in the season (starting from L s ≈ 60 • ). The amplitude of density fluctuations steadily grows from ∼5 to 10% (in the average sense) during the southern winter. This is consistent with the recent simulations using a GW-resolving Martian GCM [29], which also demonstrate the increase of GW activity at the southern high-latitudes during the aphelion season. Note that the results of simulations are shown at~80 km, whereas the majority of aerobraking data come from passes with pericenter altitudes between 105 km and 111 km (see red lines in Figure 1A). More information on the altitude dependence of density fluctuation amplitudes is given in Figure 9. The plot clearly demonstrates that ρ /ρ decreases with height, which is consistent with the dissipation of vertically propagating GWs. The estimated vertical damping rate of approximately (1/20) km −1 agrees well with the theoretical predictions for the dissipation due to nonlinear and molecular diffusion in the lower thermosphere obtained with the GW parameterization operational in terrestrial and Martian GCMs e.g., see Figure 4 in [50]. The pass-to-pass variations of fluctuation amplitudes also decay with increasing height by more than a factor of two, as is seen in Figure 9b. Note that all passes with pericenters above 111 km took place at the beginning of the Phase 2 of aerobraking, at latitudes of 28 • S .

Mean Temperature
The following relation can be used to calculate the mean temperature around the pericenter altitude (referred to as temperature hereafter) based on the fitted scale height H: where g = 3.711 m s −2 is the acceleration of gravity, M = 43.3 g mol −1 is the molar mass of air, and R = 8.314 J (mol K) −1 is the universal gas constant. In the analyses of temperatures to follow, we have used the same averaging methods as were applied for GW activity. Figure 10 presents the latitudinal distribution of thus derived background temperature for each pass (cloud of dots) and the running mean along with the RMS variations. The southern hemisphere is covered very well by the observations, except for a narrow band of latitudes (15-25 • S). Whereas temperature decreases rapidly (up to -15 K) at the low-latitudes away from the equator, as demonstrated by Phase 1, a different type of temperature variation is revealed by Phases 2 and 3. Overall, there is a clear temperature increase from the mid-, to high latitudes in the southern hemisphere in Phase 2. This thermospheric polar warming is consistent with that found from density measurements by ODY [20] (see also Figure 12 in [48]), which performed aerobraking at similar heights. The Phase 3 observations suggest a much weaker temperature variation between the high-and mid latitudes. It is remarkable that the standard deviation of temperature is the largest at mid-latitudes, with a clear increasing tendency from Phase 2 to 3. This suggests the temperature variability is larger in the early morning hours. Thus, early in the morning, dynamical processes, such as GW perturbations, can play a more dominant role than during the daytime or late in the evening (see Figure 7B). The vertical distribution of the mean temperature is given in Figure 11. It shows that the average temperature remained constant between 105 km and 110 km. In other words, the pericenter altitudes of the majority of passes coincided with the mesopause between 105 km and 110 km, although the RMS variation about the mean is sufficiently large: ±15 K. Running mean background temperature f or all phases Figure 11. The derived mean temperature as the function of the pericenter altitude for all passes. Dots represent individual orbits, and the solid line is the running average within the ±0.5 km interval.
Using the data collected by Neutral Gas Ion Mass Spectrometer (NGIMS) on-board MAVEN, Terada et al. [8] found a clear anticorrelation between the magnitudes of relative density disturbances and background temperature. They explained it as a manifestation of GW damping/saturation due to convective instability. In Figure 12A, we also plotted the amplitudes of density fluctuations for all passes as a function of the background temperature. No correlation is seen in it (more precisely, the correlation coefficient −0.003 is very small). To explore the issue further and to exclude possible variability, in Figure 12B, we plotted the same quantities, but for the running mean magnitudes. The latter have been computed in the same manner as above by averaging over passes within ±2 • latitude. The plot shows different behaviors over different phases: positive correlation during Phases 1 and 2 (0.689 and 0.729, correspondingly) and anticorrelation (−0.26) during Phase 3. This result is in line with that of Vals et al. [48], who analyzed the aerobraking data from MGS, ODY, and MRO. They discussed possible reasons for the difference with the results based on the MAVEN-NGIMS data. One of them was that the arguments of Terada et al. [8] heavily rely on temperature being isothermal with height. This assumption can well be applied to MAVEN-NGIMS measurements, which were performed close to the exobase (above 150 km). On the other hand, our results reveal no anticorrelation despite the mean temperature being isothermal as well. This finding challenges the explanation based on MAVEN-NGIMS measurements. One of the possible reasons worthy of further exploration is that the physics of GW damping is more complex. In the upper thermosphere, wave damping by molecular diffusion and thermal conduction is more dominant than convective instability. In the altitude range covered by our study, the nonlinear damping (instability) competes with damping by molecular diffusion for different harmonics [21].

Discussion and Conclusions
This paper has presented an analysis of GW-like disturbances of density fluctuations in aerobraking measurements by ESA's TGO mission. The 11-month-long aerobraking campaign provided an extensive coverage of the Martian southern hemisphere and has presented further evidence that the Martian thermosphere is continuously populated by small-scale GWs. The presented data set complements the measurements in the lower thermosphere during aerobraking of three earlier Mars orbiters: MGS, ODY, and MRO [48]. The results are in many senses consistent with the previous ones. This includes the magnitudes of relative density fluctuations associated with GWs and their lack of dependence on longitude, which demonstrates the absence of a correlation with topographic features and thus also with sources of GWs in the lower atmosphere. The TGO data also show the increase of temperature at high latitudes of the winter hemisphere (the so-called thermospheric polar warming). In addition, the peculiarity of the TGO orbits, namely, a large number of passes at high southern latitudes during Phases 2 and 3, has allowed to find a very weak diurnal variation of GW amplitudes. This finding along with the small variation of temperature with LST (not shown here) imply that the solar tide has a small amplitude at high latitudes during the southern winter.
The presented TGO measurements in conjunction with Martian GCMs can be further used for studying the mechanisms of GW propagation in the lower thermosphere and the dynamical and thermal forcing of the large-scale circulation produced by them. The finding that the inferred vertical damping rate is consistent with theoretical predictions provides optimism that the current description of GW processes is not far from reality. An intriguing issue is to understand the difference between the wave dissipation processes in the lower and upper thermosphere in general, and the lack of correlation between GW amplitudes and the background temperature in particular. The latter seems to be clearly demonstrated by the presented data set.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: