Energetic Neutral Atoms near Mars: Predicted Distributions Based on MAVEN Measurements

Determining ENA fluxes from CEX rates requires estimated ion and neutral distributions along a path where CEX occurs and ions flow, determined by some LOS. The density of the Martian exosphere is significant out to large distances from the planet, yet too low to be measured by the Neutral Gas and Ion Mass Spectrometer on MAVEN (Mahaffy et al. 2015). Instead, we use models for the cold and hot components of the H and O coronae, as well as thermospheric CO₂ densities, by Valeille et al. (2010). Specifically, we adopt the one-dimensional parameterization by Modolo et al. (2016), shown in Figure 1(a).

Zoom In Zoom Out Reset image size

The SupraThermal And Thermal Ion Composition (STATIC) instrument on MAVEN (McFadden et al. 2015) measures the local instantaneous ion flux distribution with 4 s cadence. STATIC's ion-optical assembly consists of an electrostatic deflector array, which sweeps the entrance elevation angle ± 45° for a top-hat electrostatic energy analyzer capable of selecting the ion energy per charge (E/q) between 0.1 eV/q–30 keV/q with an energy bandpass resolution ΔE/E = 16%, although the energy table used varies with time and location. A time-of-flight (TOF) system with a 15 kV preacceleration gap separates ions of differing mass per charge (m/q), with sufficient resolution to resolve major ion species such as H⁺, O⁺, and ${{\rm{O}}}_{2}^{+}$. The ions are recorded by microchannel plate (MCP) detectors with 16 azimuthally sectored anodes, resolving the instrument's field of view (FOV) in 225 bins. The full FOV is 360° × 90° for energies <5 keV, decreasing gradually to 360° × 15° at 30 keV. STATIC natively resolves each ion distribution in 16 azimuths, 16 elevations, 64 energies, and 1024 TOF bins, though binned to lower-resolution data products in order to reduce the required data downlink volume. Here, we use the STATIC d1 data product, in which each distribution is binned in 4 elevations, 16 azimuths, 32 energies, and 8 m/q bins, which can be unambiguously defined in the Mars‒Sun‒Orbit (MSO) reference frame. In Cartesian MSO coordinates, X_MSO is defined by the direction to the Sun, Z_MSO is parallel to the normal vector of the Martian heliocentric orbit, and Y_MSO completes the right-handed system.

To accurately determine the average ENA flux, the instantaneous (4 s) distributions need to be systematically organized according to the average global ion flow, which at high altitudes is mainly driven by the Lorentz force distribution in the induced magnetosphere. The orientation of the induced magnetosphere is controlled by the azimuthal clock angle of the upstream interplanetary magnetic field (IMF) (Riedler et al. 1989) and well represented by the corotating Mars‒Sun‒Electric field (MSE) reference frame (Dubinin et al. 1996). In Cartesian MSE coordinates X_MSE is antiparallel to the solar wind velocity vector, v _sw, Z_MSE is defined by the direction of the solar wind motional electric field, E _mot = − v _sw × B _IMF, and Y_MSE completes the right-handed coordinate system. Here, B _IMF is the IMF vector, which changes systematically and stochastically on a wide variety of timescales (Marquette et al. 2018), often shorter than the response time of the system. Lacking a dedicated solar wind monitor at Mars, we estimate the time-varying upstream IMF using measurements from the magnetometer investigation (MAG) on MAVEN (Connerney et al. 2015). Specifically, vector magnetic field measurements on orbit segments in the magnetosheath and undisturbed solar wind are low-pass filtered using a 30 min scale-length inverse exponential window function as described by Ramstad et al. (2020), providing interpolated estimates of the IMF vector in MSO coordinates also for orbit segments inside the induced magnetosphere and ionosphere. We estimate the IMF clock angle using the interpolated IMF vectors as ${\phi }_{\mathrm{IMF}}={\tan }^{-1}({B}_{\mathrm{IMF},{\rm{z}}}/{B}_{\mathrm{IMF},{\rm{y}}})$, though only include measurements for which the IMF is stable with a running standard deviation of ${\sigma }_{{\phi }_{\mathrm{IMF}}}\lt 45^\circ$, derived using the same window function used for ϕ_IMF.

In order to map the instantaneous ion distributions measured by STATIC in the MSE reference frame, we first correct the measured distributions for spacecraft potential and spacecraft velocity such that each measurement represents a sample of the local phase-space distribution in the Martian inertial frame. The ion distributions are subsequently transformed to the MSE reference frame by correcting for the average 4° solar wind aberration angle and rotating each measured distribution by −ϕ_IMF around the X-axis. We discretize the six-dimensional MSE phase-space domain by defining the average differential flux, $\bar{j}$, as

$$\begin{eqnarray}&&{\bar{j}}_{{S}^{+},i,j,k,l,m,n}={\bar{j}}_{{S}^{+}}({X}_{i},{Y}_{j},{Z}_{k},{\theta }_{l},{\phi }_{m},{E}_{n}).\end{eqnarray} \tag{ 1 }$$

Here, S⁺ is the ion particle species (H⁺, O⁺, or ${{\rm{O}}}_{2}^{+}$), while the subscripts i, j, k, l, m, and n are indices that refer to elements in the discretized phase-space coordinate system. X_i , Y_j , and Z_k represent MSE coordinates in the three Cartesian spatial dimensions, with each consecutive element separated from the next by a step size ΔX = ΔY = ΔZ = 0.1 R_M. θ_l , ϕ_m , and E_n are spherical coordinates in the energy-equivalent MSE velocity space, where each angular element is separated by Δθ = Δϕ = π/7 rad, and E_n is logarithmically distributed from 0.01 eV to 30 keV in 33 steps. We go through all available MAVEN orbits from 2014 November to 2020 October and average element-wise all ion distributions measured by STATIC, excluding angular elements blocked by spacecraft surfaces and all measurements from sector 11, which is blocked by the harness for the instrument's attenuator (McFadden et al. 2015). The average differential flux is determined as

$$\begin{eqnarray}&&{\bar{j}}_{{S}^{+},i,j,k,l,m,n}=\displaystyle \frac{1}{{N}_{{S}^{+},i,j,k,l,m,n}}\sum _{\iota =1}^{{N}_{{S}^{+},i,j,k,l,m,n}}{j}_{{S}^{+},\iota },\end{eqnarray} \tag{ 2 }$$

where ι = 1, 2,..., N_{S,i,j,k,l,m,n} are indices of the instantaneous differential flux measurements of ion species S⁺, made at times t_ι (which are disregarded), that sample the phase-space element [X_i , Y_j , Z_k , θ_l , ϕ_m , E_n ]. Phase-space elements of the H⁺ distribution upstream of MAVEN's orbit apoapsis are extrapolated by averaging all covered elements of the same velocity-space coordinates at X_MSE > 2 R_M. To demonstrate the resulting average ion distributions, we also integrate the average distributions over solid angle and energy in all spatial locations and display the resulting net H⁺ and O⁺ fluxes in Figure 2. The maps reproduce large-scale features expected to appear in the MSE reference frame. For example, Figures 2(d) and (h) clearly show the "plume" of O⁺ ions with energies of a few keV (Dubinin et al. 2006; Dong et al. 2015), scavenged from the ionosphere and picked up by the motional electric field of the solar wind flow, thus flowing largely in the +Z_MSE ∣∣ E _mot direction near the planet. Panels (a) and (b) show the deflection of the solar wind H⁺ bulk flow on the dayside with the divergence point offset in the + Z_MSE direction due to mass loading of the solar wind near the planet (Dubinin et al. 2018). The deceleration and heating of solar wind H⁺ population are clearly apparent when comparing the velocity-space distributions between panels (e) and (f) as the distribution widens in the latter.

Zoom In Zoom Out Reset image size

ENA instruments measure the differential flux, j, of a neutral particle species, S, with energy, E, along an LOS. To produce comparable quantities, we estimate the relevant differential ENA production rates, P, and loss rates, L, integrated along the LOS (C) from an arbitrary vantage point, r, near the planet, yielding the total differential flux along the LOS as

$$\begin{eqnarray}&&{j}_{S}({\boldsymbol{r}},\theta ,\phi ,E)={\int }_{C}({P}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )-{L}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )){ds}.\end{eqnarray} \tag{ 3 }$$

Here, ds is an infinitesimal part of path C at spatial location ${\boldsymbol{r}}^{\prime}$, whereas θ, ϕ, and E are spherical coordinates in energy-equivalent velocity space. At low energies, the planet's gravity has a significant effect on the particles' path in phase space from the point of production to the point of detection. We account for the effects of gravity and backtrace the particles' position in phase space along C, such that $\theta ^{\prime}$, $\phi ^{\prime}$, and E' are the velocity-space coordinates (${\boldsymbol{v}}^{\prime}$) of particles traversing along C at ${\boldsymbol{r}}^{\prime}$ corresponding to particles that would be detected at [ r , θ, ϕ, E].

The CEX ENA production rate of species S can be calculated from the differential flux of its corresponding ion species S⁺, and the densities of secondary neutral species, ${n}_{{S}_{s}}$, as

$$\begin{eqnarray}&&{P}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )={\bar{j}}_{{S}^{+},i,j,k,l,m,n}{{\rm{\Sigma }}}_{{S}_{s}}{n}_{{S}_{s}}({\boldsymbol{r}}^{\prime} ){\sigma }_{\mathrm{CEX},{S}^{+},{S}_{s}}(E^{\prime} ),\end{eqnarray} \tag{ 4 }$$

where S_s can be any of the neutral species included here (CO₂, O, or H) and the indices i, j, k, l, m, and n correspond to the phase-space element containing $[{\boldsymbol{r}}^{\prime} ,\theta ^{\prime} ,\phi ^{\prime} ,E^{\prime} ]$. The average differential flux ${\bar{j}}_{{S}^{+}}$ is empirically derived from MAVEN/STATIC measurements per Equation (2). The cross sections for CEX between S⁺ and S_s , ${\sigma }_{\mathrm{CEX},{S}^{+},{S}_{s}}(E^{\prime} )$, are energy dependent and experimentally determined by Lindsay & Stebbings (2005), shown in Figure 1(b). Scattering effects are only included as kinetic-collisional loss rates, L_S,col, with cross sections σ_col, contributing to the total loss rate,

$$\begin{eqnarray}&&{L}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )={L}_{S,\mathrm{CEX}}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )+{L}_{S,\mathrm{col}}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} ).\end{eqnarray} \tag{ 5 }$$

For the CEX neutralization rate, L_S,CEX, we cannot approximate the ion species as stationary compared to the ENAs and have to solve the Boltzmann collision integral. The loss rate depends on the phase-space densities associated with the ENA flux and the ion distribution,

$$\begin{eqnarray}&&{f}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )={j}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )/{v}_{S}^{{\prime} },\end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray}&&{\bar{f}}_{{S}^{+},i,j,k,l,m,n}={\bar{j}}_{{S}^{+},i,j,k,l,m,n}/{v}_{{S}^{+},n},\end{eqnarray} \tag{ 6b }$$

here defined in units of [cm⁻³ sr⁻¹ eV⁻¹]. The factors ${v}_{S}^{{\prime} }=\sqrt{2E^{\prime} /{m}_{S}}$ and ${v}_{{S}^{+},n}=\sqrt{2{E}_{n}/{m}_{{S}^{+}}}$ are the speeds of the ENA and ion species, respectively. The relative ENA velocity with respect to each phase-space element is ${\rm{\Delta }}{{\boldsymbol{v}}}_{S,{S}^{+},l,m,n}\,={{\boldsymbol{v}}}_{{S}^{+},l,m,n}-{{\boldsymbol{v}}}_{S}^{{\prime} }$ with corresponding energy ${E}_{{\rm{\Delta }}v,S,{S}^{+},l,m,n}\,={m}_{{S}^{+}}| {\rm{\Delta }}{{\boldsymbol{v}}}_{S,{S}^{+},l,m,n}{| }^{2}/2$. The CEX neutralization rate can be determined by integrating the differential rates over velocity space as

$$\begin{eqnarray}&&\begin{array}{l}{L}_{S,\mathrm{CEX}}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )=\sum _{{S}^{+},l,m,n}| {\rm{\Delta }}{{\boldsymbol{v}}}_{S,{S}^{+},l,m,n}| {f}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )\times \\ {\bar{f}}_{{S}^{+},i,j,k,l,m,n}{\sigma }_{\mathrm{CEX},S,{S}^{+}}({E}_{{\rm{\Delta }}v,S,{S}^{+},l,m,n}){\rm{\Delta }}{{\rm{\Omega }}}_{l,m}{\rm{\Delta }}{E}_{n},\end{array}\end{eqnarray} \tag{ 7 }$$

where ${\rm{\Delta }}{{\rm{\Omega }}}_{l,m}={\rm{\Delta }}{\theta }_{l}{\rm{\Delta }}{\phi }_{m}\cos ({\phi }_{m})$ is the space angle of the phase-space angular element [θ_l , ϕ_m ]. Lacking experimentally determined CEX cross sections for reactions neutralizing ${{\rm{O}}}_{2}^{+}$, we assume these are similar to the equivalent cross sections for O⁺. For the collisional loss term, L_S,col, we only include collisions between the ENAs and other (secondary) neutrals, S_s , as these are by far more numerous at nearly all altitudes. The collisional cross sections are determined by the kinetic radii, r_S , of each particle pair as ${\sigma }_{\mathrm{col},{\rm{S}},{{\rm{S}}}_{s}}=\pi {({r}_{S}+{r}_{{S}_{s}})}^{2}$. The collisional loss rate can be determined by adding up the collision rates with all secondary species as

$$\begin{eqnarray}&&{L}_{S,\mathrm{col}}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} )={j}_{S}({\boldsymbol{r}}^{\prime} ,{\boldsymbol{v}}^{\prime} ){{\rm{\Sigma }}}_{{S}_{s}}{n}_{{S}_{s}}({\boldsymbol{r}}^{\prime} ){\sigma }_{\mathrm{col},S,{S}_{s}}(E^{\prime} ).\end{eqnarray} \tag{ 8 }$$

With production and loss rates determined along the LOS, the ENA distribution at any spatial location r can be determined by numerically evaluating the line integral in Equation (3) for a discretized domain of directions and energies. The step size (ds) cannot be infinitesimal in numerical integration. Instead, Equation (3) is iteratively evaluated forward along C as a Riemann sum using a finite step size (Δs) selected to be smaller than the typical length scales of the plasma and neutral distributions at any corresponding altitude, h. For computational efficiency, the values for Δs depend on altitude range as Δs(h < 300 km) = R_M/ 100, Δs(300 km < h < 500 km) = R_M/40, for Δs(500 km < h <6800 km) = R_M/20 and for Δs(h > 6800 km) = R_M/5. The loss rates at each step (Equations (8) and (5)) are calculated based on the differential ENA flux derived from the previous step with the initial value starting at zero, thus the derived ENA flux is unable to continuously grow in optically thick (i.e., collisional) regions.

Examples of ray-traced ENA fluxes are shown in Figure 3, here providing parts of the full distributions at the two shown convergence points. The rays correspond to the spatial location of C in Equation (3) and bend slightly near the planet due to gravity affecting the path in phase space. The nadir-looking rays require the smallest number of iterations (25 steps to 2 R_M from the ray start at 100 km altitude), while more iterations are required for all other rays. The H-ENA rays are traced out to 30 R_M upstream of the planet by extrapolating the upstream solar wind H⁺ distributions, and O-ENA rays are limited to the domain covered by MAVEN, i.e., within the ∼2.8 R_M apoapsis of its original science orbit.

Zoom In Zoom Out Reset image size

The accuracies of the synthetic ENA distributions derived here are inherently limited by the accuracy of the ion distributions, neutral profiles, and cross sections used to derive them. For instance, the 225 × 225 angular resolution of the STATIC d1 data product is insufficient to accurately resolve the narrow solar wind distribution. In addition, ions scattering off STATIC's upper deflector support posts lead to a few percent of the true distribution appearing in other angular sectors and energies. Together with true variability in the solar wind direction and speed, these factors contribute to the apparent wide antisunward H-ENA distribution created from the neutralization of H⁺ in the upstream undisturbed solar wind. The instantaneous antisunward H-ENA distribution would appear much narrower in space angle and energy to an ENA instrument capable of resolving the distribution. The uncertainties in the CEX cross sections are comparatively small, below ±10% and ±15% for energies over and under 500 eV, respectively (Lindsay & Stebbings 2005).

Ignoring electron‒ion recombination as an ENA production term (e.g., O⁺ + e⁻ → O⁰) could hypothetically imply that ENA fluxes emanating from the lowest altitudes are somewhat underestimated. Considering, however, that recombination cross sections are small and to first order depend on electron energy as $\propto 1/{E}_{{e}^{-}}$ (Nahar 1999), recombination can only be a significant production term in the cold and dense lower ionosphere (∼130 km altitude), which is optically thick to ENAs. In the hot and relatively tenuous top-side ionosphere, from where the observable ionospheric ENAs are sourced, we take recombination rates to be negligible.

The largest source of systematic uncertainty is likely inherited from the coronal model by Valeille et al. (2010) and Modolo et al. (2016), mainly considering how unconstrained the exospheric structure was in 2010. In particular, the parameterized model is spherically symmetric, while in reality the exosphere has a three-dimensional structure (Chaffin et al. 2015) and also varies with season (Rahmati et al. 2018). In effect, ENA fluxes emanating from the near-planet nightside are likely overestimated due to the associated higher production rates, at least relative to the dayside.

Inaccuracies in the neutral densities at high altitudes would influence the total column density and thus the fluence of antisunward H-ENAs formed from undisturbed solar wind protons. Here, the estimated fluence of this population, ∼3 × 10⁵ cm⁻² s^−1, increases to ∼7 × 10⁵ cm⁻² s⁻¹ close to the planet, indicating that about 1% of the solar wind protons are converted to H-ENAs, an estimate consistent with those found by Kallio et al. (1997) and Halekas et al. (2017).

Overall, the H-ENA fluxes on the order of up to 3 × 10⁵ cm⁻² s⁻¹ sr⁻¹ found here to emanate from the subsolar magnetosheath are largely consistent with Mars Express NPD observations (Futaana et al. 2006; Wang et al. 2019), particularly considering the unconstrained differences in upstream conditions and the instantaneous nature of the observations, as opposed to this time-average distribution. From a vantage point in the planet's shadow, the H-ENA fluxes from the magnetosheath flanks surrounding the planet are generally in the range of 10⁴–10⁵ cm⁻² s⁻¹ sr⁻¹, which is also largely consistent with the MEX observations (Galli et al. 2008).

Synthetic ENA distributions are useful in a number of aspects. As described earlier, the ENAs carry information about both the plasma and neutral distributions along the LOS; however, actual ENA measurements are required to collect this information. Designing an investigation to measure ENAs strongly benefits from some understanding of the environment to be measured in order to optimize the orbit, operations, and instrumental properties, such as the energy range, sensitivity, and dynamic range. For example, we can clearly see in Figures 4 and 5 that the bulk of O-ENAs anywhere around the planet have energies <100 eV. Relatively weak fluxes of 0.1–1 keV O-ENAs appear at vantage points over the electric field pole (∼90° in Figure 5), likely originating as pick-up O⁺ undergoing the initial phase of acceleration. The most energetic (1–10 keV) O-ENA fluxes are only found at about 180°–240° in the MSE XZ-plane (Figure 5), originating near the planet in the −Z_MSE hemisphere. For both of these energetic populations, the typical flux is 10³–10⁴ cm⁻² s⁻¹ sr⁻¹; as such, it is perhaps not particularly surprising that the NPD on Mars Express with its lower-energy limit of 100 eV and 10⁴ cm⁻² s⁻¹ sr⁻¹ sensitivity threshold did not conclusively detect any O-ENAs (Galli et al. 2008).

With past or future measured ENA distributions in hand, comparisons with synthetic ENA distributions provide constraints on either the plasma or neutral environment. If the average ion phase-space distribution during actual measurements can be constrained, then exospheric models could be fitted such that the synthetic ENA distributions match the measured distributions. In this respect, we should consider that the ENA distributions derived here represent the estimated time-averaged near-Mars ENA environment, which is inherently dependent on the prevailing solar wind and ionospheric conditions, as well as the exosphere models used. While the average ion distribution (Figure 2) is produced using nearly all currently available and calibrated MAVEN STATIC data, this time period (2014 November–2020 October) was dominated by a deep solar minimum and any accurate comparisons with past or future ENA measurements may require constraints on specific upstream or seasonal conditions. ENA observations during other parts of the solar cycle will reflect the cycle's combined effect on the Martian plasma and neutral environments. Solar maxima are associated with increased solar extreme ultraviolet (EUV) ionizing radiation, which increases ionization rates in the upper atmosphere, in turn populating the hot oxygen corona through the dissociative recombination of ${{\rm{O}}}_{2}^{+}$ (Lillis et al. 2017). The increased ion production due to higher EUV levels also supplies more ions for outflow, though decreases the characteristic energy of the tailward outflowing ions (Ramstad et al. 2017). As such, we would expect low-energy O-ENA fluxes to increase drastically under high-EUV solar maximum and otherwise nominal solar wind conditions. The effect on high-energy O-ENAs is uncertain due to the competing effects of reduced solar wind–ionosphere coupling and increased CEX rates in the denser exosphere. Nevertheless, solar wind parameters are also far more variable during solar maxima and likely to yield at least periods of both stronger and higher-energy H-ENA and O-ENA fluxes as fast solar wind streams and interplanetary shocks interact with the planet.

These synthetic ENA distributions are also unable to reproduce dynamics and features that are not driven by the solar wind electric field, such as ENAs created by ion flows organized around the crustal fields (Wang et al. 2014) or other geographic features, which are smeared over time in the MSE reference frame. Additionally, there are other potentially significant ENA sources that do not produce ENAs through ion‒neutral CEX, such as backscattering (Kallio & Barabash 2000) and sputtering (Shematovich & Kalinicheva 2020), which would not be represented in these synthetic ENA distributions.