Long‐Term Variability of Mars' Exosphere Density Based on Precise Orbital Analysis of Mars Reconnaissance Orbiter and Mars Odyssey

The variability of Mars exosphere over monthly to solar‐cycle scales at 251 and 412 km altitude is quantified by analysis of 41‐Ls mean densities derived from precise orbit determination of the Mars Reconnaissance Orbiter (MRO) and Mars Odyssey (MO) satellites, respectively. The data encompass 2006–2020 (MRO) and 2002–2020 (MO). At both altitudes, most of the variance is captured by cos(Ls–ϕ), where ϕ ≈ 258°. This term represents the effects of solar heating changes due to the eccentricity of Mars orbit around the Sun, and climatological changes in heating due to lower‐atmosphere dust loading, which does not play a significant role. The remaining variability is connected with the “irregular” variability of solar flux over monthly time scales. For MO, the presence of Helium disrupts a clean correlation with these sources.


Introduction
A process fundamental to the behavior of planetary upper atmospheres is the absorption of solar radiation from the Sun and the variability that it produces.In addition to photoionzation and the creation of ionospheres, solar radiation absorption also heats the atmosphere and drives the variability of neutral densities over time scales ranging from minutes to hours (e.g., solar flares), to the ∼27-day rotation of the Sun, to the seasonal and solar-cycle (∼11year) scales that encompass the focus of the present paper.In Earth's thermosphere, it is the variability of extreme ultraviolet (EUV) radiation originating in the solar corona that is most relevant to density variability.EUV/UV radiation absorption is also important at Mars, but the origins of its solar-driven density variability also reside in the eccentricity of Mars orbit around the Sun.This introduces an annual dependence in the solar radiation received at Mars (the so-called ``orbital effect''), which is maximum at Mars' closest approach to the Sun (perihelion, at solar longitude Ls = 251°).
Infrared radiation (IR) from the Sun does not change significantly, and therefore does not play a role in the variability of Earth's thermosphere, although IR cooling to space by CO 2 is relevant to its long-term behavior.However, IR radiation and in fact the whole spectrum of radiation received from the Sun at Mars are subject to the orbital effect.The absorption of IR radiation by CO 2 in Mars' middle atmosphere (Bougher & Dickinson, 1988) can lead to thermal expansion of the overlying atmosphere, thus contributing to the annual variation of density in the thermosphere and exosphere (Fang et al., 2022).Similarly, radiation absorption by dust in the lower and middle atmosphere, which climatologically has a seasonal variation close to that of the orbital effect, can also contribute to solar-driven density variability in the thermosphere (e.g., Liu et al., 2018).
Efforts to quantify the above long-term variability have involved analyses of densities obtained from precise orbit determination (POD) of Mars Global Surveyor (MGS) at 390 km during 1999-2004 (Forbes et al., 2008), and of • Dust-induced inter-annual variability is small, but the analysis method does suppress potentially important regional-scale events

Supporting Information:
Supporting Information may be found in the online version of this article.
Mars Odyssey (MO) at 405 km during 2004-2010 (Bruinsma et al., 2014).Bruinsma et al. (2014)  A similar POD analysis was not considered for MAVEN, since periapsis has not remained fixed (most recently raised to 180-200 km in 2020), similar analyses can be performed at fixed heights using NGIMS data (e.g., Fang et al., 2022), and the range of solar activity is not yet comparable to those experienced by MRO and MO.

Data
An overview of the data sets considered in the present study is presented in Figure 1, along with the MGS densities analyzed in Forbes et al. (2008).The MRO densities correspond to its periapsis altitude of 251 km and extend from 2006 to 2020.The MO and MGS densities correspond to mean altitudes of 412 and 394 km, and 2002-2020 and 1999-2006, respectively.All three satellites are in low-eccentricity near-polar day/night Sun-synchronous orbits with maximum drag weighted toward middle to high latitudes in the Southern Hemisphere (SH).The 81-day mean 10.7 cm solar radio flux adjusted to Mars is also shown in Figure 1, and shows the relation of these density data sets to solar cycle.The quasi-2-year oscillations reflect the relationship between solar flux and exosphere density imposed by the eccentricity of Mars' orbit around the sun, which has a period of 1.88 Earth years.
The derivation of densities from MRO POD follows the same basic procedure described for MO in Bruinsma et al. (2014) and for MGS in Bruinsma and Lemoine (2002), and the reader is referred there for details.Briefly, a scale factor is applied to a Mars thermosphere empirical model (DTM Mars, Bruinsma & Lemoine, 2002) to minimize the differences between the calculated and observed orbits of these satellites.The positions of these satellites are provided by tracking data from the Deep Space Network (DSN), which are corrected for effects introduced by Earth's atmosphere, solid-Earth and ocean tides.The densities are derived over 96 hr arcs for MO, and thus represent approximately 4-sol day-night and longitudinal means.Due to the lower altitude of MRO and the larger atmospheric drag, much shorter arcs are possible, even as short as 1-orbit cadence.In the present study, the mean densities are confined to 60°S-90°S latitude with 24 hr (∼1-sol) duration.Therefore, the MRO data also represent day-night and longitudinal means, but over this restricted latitude range.Information on sampling and distribution of the densities, which shows the bias toward the SH of the three data sets but notably for MRO, and the three solar local time planes for MO, is given in Table S1, Figures S1, and S2 in Supporting Information S1.The derivation of densities requires accurate modeling of forces on the satellites other than atmospheric drag.These include Mars' gravity field, including perturbations due to the Sun, planets, and Martian moons; solar and planetary radiation pressures; thruster firings; and momentum exchange between the satellite surface and the colliding particles, which is modeled through a drag coefficient.The drag coefficient involves modeling the satellite profile projected into the direction of the velocity vector, which also requires knowledge of satellite shape (including satellite bus, solar array, high-gain antenna, etc.) and orientation, the latter being provided through telemetered quaternions.Drag coefficients are stable for MO and MRO, as can be seen in Figure S3 of the Supporting Information S1.
Compared to the MO results presented in Bruinsma et al. (2014), the MO data are much less noisy due to the improvement of force models used, especially the gravity field of Mars, the more precise knowledge of the attitude of the satellite (including solar panels and high gain antenna) and a new strategy for adjustment of residual accelerations during desaturation of reaction wheels.Uncertainties for single-arc density values for all 3 satellites are now less than 15% (RMS), and are mainly due to tracking errors and solar radiation pressure modeling.DSN observations become very noisy near conjunction (Earth-Sun-Mars angle = 0), and more data tend to be rejected and reflect greater variability around these periods.These effects are mitigated by taking 41-Ls running means, reducing the random errors from 15% to 15/ ̅̅̅̅̅ 10 √ = 4.7% for MO, and 15/ ̅̅̅̅̅ 40 √ = 2.4% for MRO.A constant bias of up to 10% is also possible due to uncertainties in the assumed drag coefficient; however, this does not affect analyses and conclusions with regard to density variability.

Analysis Method
We use the method of multiple linear regression to quantify the density variabilities displayed in Figure 1.The method consists of fitting a function of the form: in a least-squares sense to either ρ or log 10 ρ where ρ is the daily 41-Ls mean density; MF10.7′ is the "residual" (defined below) 41-Ls mean 10.7-cm solar radio flux adjusted to Mars' orbit; f(CDOD) represents the influence of lower-atmosphere dust where CDOD is the column dust optical depth (https://www-mars.lmd.jussieu.fr/mars/dust_climatology/) from TES and THEMIS nadir observations integrated over all longitudes and between 0°and 90°S latitude; and ϕ is the phase of an annual variation.Note that ϕ = 251°(perihelion, Mars' closest approach to the Sun) corresponds to the so-called "orbital effect" on the solar flux; that is, the modulation of solar irradiance through the inverse square of Mars' distance to the Sun.As shown by Fang et al. (2022), the orbital effect includes not only density changes due to in-situ absorption of extreme ultraviolet (EUV) and UV radiation in the thermosphere, but also the absorption of infrared (IR) radiation between 50 and 120 km altitude (e.g., Bougher & Dickinson, 1988).However, while the former varies with solar activity and Mars' distance from the Sun, the latter is only subject to the orbital effect.MF10.7, which is commonly used as a proxy for EUV/UV solar radiation variability, contains variability due to both solar activity and the orbital effect.However, the cos(Ls-ϕ) term in (1) also represents variability due to IR heating and heating associated with lower atmosphere dust (CDOD also varies annually with maximum around perihelion; see below).Therefore MF10.7 is to some degree correlated with these other sources of variability, which is problematic when applying least-squares regression, in particular when allocating variances to the different sources of variability.To address this problem, we first fit cos(Ls-ϕ) into MF10.7 to obtain the MF10.7 residuals (MF10.7′)from the cos (Ls-ϕ) fit; then, MF10.7′ is used in the fitting instead of MF10.7.The same approach is applied to CDOD to obtain CDOD′ as the fitting parameter.This leaves the responses to the EUV/UV and IR flux variability due to the orbital effect, and the purely annual variation due to dust heating, all contained in the cos(Ls-ϕ) term.
In addition to redefining variables in Equation 1 to de-correlate the fitting variables, there are a few other aspects of the present work that differ from those applied to MGS in Forbes et al. (2008) and to MO in Bruinsma et al. (2014).Our prior analysis of MO data revealed low densities during the 2006-2010 solar minimum period that could not be well-captured by MF10.7.In this paper, we also analyze the MO data using the S10.7 index (Tobiska et al., 2008)  The reason is that ρ and CDOD are not close to the sinusoidal variation expressed in the cos(Ls-ϕ) term in (1), but log 10 ρ and log 10 CDOD are.The sinusoidal nature of log 10 CDOD is demonstrated in Figure S5 in Supporting Information S1.The sinusoidal nature of the density variability is obvious for MRO in Figure 1, but not as much for MO.In fact, as shown below, for MO the best correlations are found with respect to ρ and CDOD.

Results
Figure 2 illustrates results for log 10 ρ of MRO.Table 1 contains numerical values of coefficients in Equation 1 and information on the variance captured by the fit and sub-elements of Equation 1.A stepwise approach is used to properly include the influences of CDOD on the fit, as explained below.The top panel (a) of Figure 2 illustrates the fit to log 10 ρ including the cos(Ls-ϕ) and MF10.7′ terms, which achieve a correlation coefficient of R = 0.99.
From Table 1 (MRO-(b)), the variance of the data is 0.24, and the fit captures 97.% of the variance.The cos(Ls-ϕ) term with ϕ = 259.ºaccounts for 91.% of the total fit variance (MRO-(d)), while the MF10.7′term accounts for 6.0% of the total variance; the sum of these two terms thus accounts for all 97.% of the total variance associated with the fit.
Our first attempt at adding the influence of CDOD to the fit took the form f(CDOD) = log 10 CDOD′, where log 10 CDOD′ was the residual of a fit to log 10 CDOD of the form cos(Ls-ϕ).This was done to ensure that the first and last terms of Equation 1were uncorrelated with each other.log 10 CDOD′ represents the CDOD variability that is not represented by the annual term.The fit variance in this case was 0.23 (98.%) (MRO-(g) in Table 1), and variances associated with the cos(Ls-ϕ), MF10.7′ and log 10 CDOD′ terms were 67.%, 6.3%, and 2.2% (not shown in Table 1), totaling 75.% which is significantly short of 98.%.These differences result from the fact that log 10 CDOD and MF10.7′ were fortuitously to some degree correlated.Therefore, f(CDOD) = log 10 CDOD′ in Equation 1 was replaced by f (CDOD) = log 10 CDOD′′, where log 10 CDOD′′ is the residual from the fit to both MF10.7′ and cos(Ls-ϕ).The result is that all 3 terms in Equation 1 are uncorrelated with each other.
The results of including log 10 CDOD′′ in the fit are summarized in Figure 2b and Table 1, MRO(f)-(j).The correlation coefficient is now also R = 0.99, and the variances captured by the annual and solar flux terms remain at 91.% and 6.0%, respectively, and where ϕ remains unchanged at 259. o .The variance captured by the log 10 CDOD′′ term, the aggegate influence of interannual variability of lower-atmosphere dust heating, is only 0.27%.
To further explore the dust effect, the residuals between the fit and the data in Figure 2a were examined against log 10 CDOD (see Figure S6 in Supporting Information S1).That figure shows that the two peaks in the log 10 ρ residuals from the fit during MY 28 and MY 34 are fully captured by the log 10 CDOD′′ term in Equation 1 and in Table 1 (MRO-(f)), and to first order log 10 ρ ≈ 0.1⋅log 10 CDOD.If we assume that the annual variations of log 10 ρ and log 10 CDOD follow in proportion to their ratios during these dust storm events, then the 0.72 amplitude of log 10 CDOD (cf. Figure S5 in Supporting Information S1) implies an annual amplitude of 0.072, or 10% of the C1 coefficient in MRO-(f) in Table 1.It follows that solar heating associated with dust contributes an estimated 9.1% of the total 91.% attributable to the annual term, the remaining 82.% being accounted for by the orbital effect on solar heating due to absorption of solar UV/EUV and IR radiation throughout the atmospheric column.
Figure 3 depicts the same information as Figure 2, except for MO over MY25-35, and the corresponding numerical data are provided in Table 1.For MO the better correlation was found with respect to ρ (R = 0.84) as opposed to log 10 ρ (R = 0.78); it is the results for ρ that are presented here.Note also that the solar flux term is represented by MF10.7′, as it was found that MS10.7′ actually resulted in a slightly worse correlation.From Table 1, 70.% of the variance was captured, partitioned between the annual term with ϕ = 257.o (42.%) and the solar flux term (27.%) (i.e., MO(d)-(e)).When the CDOD′ term is added, which is uncorrelated with MF10.7′ (see Figure S5 in Supporting Information S1), a mere 0.31% out of 70.% of the fit variance is captured, and virtually nothing else changes.When fit with respect to log 10 ρ instead of ρ, a total of only 60% of the variance was captured by the fit (not shown).
For MO, the residuals between the fits and the data were examined against CDOD, similar to the procedure followed for MRO.No perceptible relationships were found between CDOD and the density residuals, leading us to conclude that the effects of dust heating are not detectable at the 412 km altitude of MO.

Summary and Conclusions
Through least-squares multiple linear regression, we have identified the major contributors to long-term density variability in the exosphere of Mars.To obtain the best fits, the dependent variable is 41-Ls mean log 10 ρ for MRO and 41-Ls mean ρ for MO.The independent variables include an annual term with maximum near perihelion, that captures the variability due to heating by the absorption of UV/EUV radiation in the thermosphere and IR radiation in the middle atmosphere due to Mars orbit around the Sun, and the annual variability associated with heating associated with solar radiation absorption by dust.A second term represents the density response to heating by the "irregular" component of UV/EUV solat flux absorbed in the thermosphere; that is, the part that is not captured by the annual term associated with the "orbital effect."The 10.7 cm solar radio flux adjusted to Mars  • 98.% of the observed variance is captured by the fit.The main contributors are the annual term with maximum at Ls = 259.º(91.%), and the irregular part of the solar UV/EUV flux variability (6.0%).• According to this fit, the irregular component of dust heating contributes only 0.27% to the total variance captured by the fit.However, if we assume that the proportionality between log 10 ρ and log 10 CDOD in the annual variation is the same as that during the planet-encircling events of MY28 and MY34, then an estimated 9% of the above 91.%can be attributed to dust heating.
With respect to ρ at 412 km at in the SH from MO POD during 2002-2020: • 70.% of the observed variance is captured by the fit.The main contributors are the annual term with maximum at Ls = 257.º(42.%), and the irregular part of the solar UV/EUV flux variability (27.%).• The irregular component of dust heating contributes only 0.44% to the total variance captured by the fit.No additional contributions of dust heating, for example, to the annual variation in density, were detectable as is the case for MRO.• Given the less robust fit for MO as compared with MRO, especially around solar minimum, the S10.7 parameter at Mars (Tobiska et al., 2008) was also used instead of F10.7 for MO.However, use of S10.7 did not result in an improved fit.• Based on the earlier MO data set for 2002-2010, Bruinsma et al. (2014) ) speculated that residuals between the fit and the data may be due to the presence of He.Based on subsequent measurements and modeling (Elrod et al., 2017(Elrod et al., , 2023;;Gupta et al., 2021), the presence of He in MO-measured densities is likely.However, in addition to affecting total mass density, the influences of He on satellite drag coefficients are open to considerably uncertainty.Therefore, these two effects combined are not quantifiable, and how much of the remaining 32% of the variance not captured by the total fit to MO densities may be ascribable to He remains undetermined.Finally, it is important to recognize that our use of 41-Ls running means throughout this paper suppresses the potential influences of solar variability and lower atmosphere dust at shorter time scales.For instance, while 27day mean MF10.7 correlate very highly (R = 0.96) with those of 17-22 nm EUV fluxes measured from MAVEN, visibly more scatter exists with respect to daily values (Fang et al., 2022, Supporting Information).This means that solar-induced variability in densities likely exist that are omitted in our analysis, such as 27-day solar rotation (e.g., Forbes et al., 2006;Hughes et al., 2022).Similarly, use of 41-Ls mean CDOD averages over the SH suppresses any potential relationships between density and shorter-duration regional-scale dust events that are known to exist during many MY (e.g., Kass et al., 2016).

•
Densities from precise orbit determination of Mars Reconnaissance Orbiter and Mars Odyssey are used to quantify long-term variability of Mars' exosphere • The annual term is the greatest contributor to total variance, followed by the irregular component of EUV/ UV flux intercalibrated these two data sets to create a unified data series of densities at 405 km extending from 1999 to 2010.More recently, Fang et al. (2022) analyzed CO 2 , Ar, N 2 , and O densities measured by the NGIMS instrument on MAVEN between 180 and 275 km during 2015-2020 to deduce the origins of their long-term variability.In this paper, we analyze a new POD data set consisting of Mars Reconnaissance Orbiter (MRO) densities around periapsis at 251 km and between 60°S-90°S during 2006-2020.In addition, we present results for a new improved and extended (2002-2020) MO data set.
adjusted to Mars (MS10.7),which does account for the low solar fluxes that occurred during this unusually deep solar minimum.The S10.7 index is based on the integrated 26-34 nm solar irradiance measured by the Solar Extreme-ultraviolet Monitor (SEM) instrument on the NASA/ESA Solar and Heliospheric Observatory (SOHO) research satellite, and is designed for use in the JB2008 empirical satellite drag model(Bowman, Tobiska, Marcos, et al., 2008;Bowman, Tobiska, & Kendra, 2008) for Earth.FigureS4in Supporting Information S1 compares MF10.7 with MS10.7 for MY25-35, and also includes CDOD.The noteworthy points are the deeper solar minima for MS10.7 versus MF10.7, and the large dust enhancements during the planetencircling dust storms during MY 28 and MY 34.In addition, for MRO we have found that improved correlations are obtained by fitting log 10 ρ instead of ρ, and similarly in this case log 10 CDOD is used instead of CDOD.
C0 = 9.03, C1 = 3.66, ϕ = 257., C2 = 0.172, C3 = 3as a proxy for the UV/EUV flux.A third term provides an aggregate measure of the variability associated with the irregular component of height-integrated dust content in Mars atmosphere.These latter two terms are also represented by their 41-Ls means.With respect to log 10 ρ at 251 km at high latitudes in the SH from MRO precise orbit determination (POD) during 2006-2020:

Table 1
Coefficients in Equation 1 Obtained by Fitting MRO and MO Data With ((f)-(j)) and Without ((a)-(e)) the f(CDOD) Term, the Associated Variances, and Relative Contributions (R 2 ) to the Total Fit Variance Associated With Each Term