© The Authors 2026

1 Introduction

Atmospheric escape has played a fundamental role in shaping the current state and long-term evolution of Mars (Lammer et al. 2013; Jakosky 2021). Geological and isotopic evidence indicates that early Mars likely possessed a much denser atmosphere, capable of sustaining surface liquid water, which could have supported a more habitable environment (e.g., Wordsworth 2016; Jakosky et al. 2017; Ramirez & Craddock 2018; Kite 2019). Over billions of years, however, the gradual loss of volatiles transformed Mars into the cold and arid planet we observe today. A variety of mechanisms have been proposed to account for this loss, including thermal escape, photochemical escape, ion pickup, bulk ion outflow, and atmospheric sputtering (Jakosky et al. 2018, and references therein), with the latter recognized as a particularly effective process. In this mechanism, incident high-energy particles collide with neutral atmospheric species, initiating cascades of secondary collisions, some of which transfer sufficient kinetic energy to allow particles to escape the planet’s gravitational potential Johnson (1994). Unlike other processes that predominantly affect light species, sputtering can eject heavier molecules, particularly CO₂ and Ar, thereby significantly altering the atmospheric composition and enhancing atmospheric loss (e.g., Wang et al. 2014, 2015; Jakosky et al. 2017; Leblanc et al. 2018; Curry et al. 2025).

Despite its recognized importance, the quantitative contribution of sputtering to atmospheric escape on Mars remains uncertain. Early modeling studies suggested that sputtering could account for the removal of a substantial fraction of Mars’ ancient atmosphere (e.g., Luhmann et al. 1992; Johnson et al. 2000). More recent investigations, however, have reported widely varying escape rates, reflecting the strong sensitivity of sputtering escape to factors such as solar wind energy spectra, atmospheric density profiles, and adopted collision cross sections (e.g., Wang et al. 2014; Leblanc et al. 2018). A central unresolved issue is the identification of the dominant drivers of sputtering escape. While solar wind protons are generally considered the primary source of momentum transfer (e.g., Leblanc et al. 2018; Shematovich 2021), magnetosheath-heavy ions and energetic neutral atoms (ENAs) generated through charge exchange may also play critical roles under certain solar wind conditions (e.g., Fang et al. 2013; Ma et al. 2025).

It is important to distinguish these sputtering sources. Many previous studies have focused on sputtering driven by pickup O⁺ ions, which can transfer momentum more efficiently to heavy atmospheric species because of their larger masses and, in many cases, higher characteristic energies (e.g., Luhmann et al. 1992; Johnson et al. 2000; Leblanc et al. 2018). Such particles are therefore expected to be more effective in ejecting heavy species such as CO₂ and Ar. By contrast, sputtering driven by precipitating solar wind protons and hydrogen ENAs represents a lighter but highly variable source that is directly linked to the solar wind interaction with Mars. The situation is further complicated by the highly variable interaction between the solar wind and the Martian magnetosphere: Mars’ induced magnetosphere is dynamic and strongly modulated by upstream solar wind conditions, giving rise to pronounced temporal and spatial variability in particle precipitation and sputtering efficiency (e.g., Dubinin et al. 2006; Brain et al. 2010; Wang et al. 2018). These uncertainties make it difficult to reconcile model predictions and to place robust constraints on the long-term contribution of sputtering to Martian atmospheric loss.

The simultaneous operation of China’s Tianwen-1 mission (Zou et al. 2021) and NASA’s Mars Atmosphere and Volatile Evolution (MAVEN) mission (Jakosky et al. 2015) provides an unprecedented opportunity to reduce these uncertainties and advance our understanding of sputtering escape on Mars. MAVEN has been a cornerstone of Mars aeronomy, delivering comprehensive measurements of the solar wind environment, magnetospheric configuration, and upper atmospheric structure, while Tianwen-1 complements this dataset with independent in situ plasma and ENA observations. The ability to analyze coincident measurements from two spacecraft makes it possible to disentangle the relative roles of ENAs and precipitating solar wind protons in driving sputtering escape and evaluate how their contributions vary under different conditions.

This study builds on the Monte Carlo sputtering framework and Tianwen-1 ENA context developed by Ma et al. (2025) but addresses a different and more specific objective. Here, we combined MAVEN-derived precipitating solar wind proton fluxes with Tianwen-1 hydrogen ENA observations and employed a test particle Monte Carlo model to quantify the sputtering escape driven by these two light incident particle populations. This dual-mission approach enables a quantitative assessment of proton and hydrogen ENA-induced sputtering at Mars, providing an observationally constrained estimate of their contribution to atmospheric escape rather than a complete description of all sputtering sources. These results help clarify the role of lightparticle-driven sputtering in Martian atmospheric evolution and its implications for the long-term habitability of the red planet.

Fig. 1 Precipitating solar wind protons and hydrogen ENAs at Mars during December 2021, shown for each MAVEN and Tianwen-1 orbit. (a) Energy spectrum of precipitating solar wind protons. (b) Hydrogen ENA velocity. (c) Hydrogen ENA flux.

2 Data and methods

2.1 Precipitating solar wind proton observations

Precipitating solar wind protons can interact with the Martian atmosphere and drive sputtering escape. To quantify this process, information on the precipitating proton flux is required. Ideally, this flux should be evaluated near the Martian exobase, where the transition from collisionless to collisional behavior occurs and sputtering becomes effective. This region is often taken to lie near 200 km (e.g., Fu et al. 2020). However, ion-neutral collisions can still occur above this altitude. Chaufray et al. (2007) reported that the charge-exchange probability can reach nearly 50% above 200 km. Consequently, if measurements are restricted to the exobase, the total particle energy ultimately deposited into the atmosphere may be underestimated, because a fraction of the incident energy has already been lost through neutralization and thermalization. For these reasons, many studies adopt an upper altitude limit of ~350 km. For example, Leblanc et al. (2015) assumed that ions with downward velocity vectors lying within a 75° cone around the radial direction would ultimately impact the atmosphere.

In this study, we instead sampled a higher-altitude range of 450–550 km. To ensure that the selected particles can still reach the atmosphere, we adopted a more stringent angular criterion and restricted the allowable angle to 60°. Two instruments – the Solar Wind Ion Analyzer (SWIA; Halekas et al. 2015) and the SupraThermal and Thermal Ion Composition instrument (STATIC; McFadden et al. 2015) – on board MAVEN are capable of measuring the precipitating proton flux. For SWIA, the 3D measurements comprise 64 angular bins. We classified an anode as observing precipitating ions when, within 450–550 km, the angle between the center of the anode field of view (FOV) and the local zenith direction is smaller than 60°. We then selected all anodes satisfying this criterion throughout the altitude interval and summed their counts to obtain the precipitating flux. To ensure adequate spatial coverage, we retained only the intervals in which the FOV coverage, defined as the sum of the solid angles of the selected SWIA bins divided by the solid angle of a cone with a 60° apex angle, 2π(1 – cos60°) ≈ 3.14 sr, exceeds 65%. For each time and energy bin, we summed the counts from all selected anodes to compute the precipitating flux. Due to its limited FOV, STATIC was not used for flux calculations. Figure 1a presents the energy spectrum of the precipitating protons as a function of date, clearly showing that the proton flux exhibits pronounced temporal variability.

2.2 Hydrogen ENA observations

Tianwen-1, China’s first Mars exploration mission, was launched on July 23, 2020, and began its orbital science observations in November 2021. Hydrogen ENA observations used in this study were obtained from the Mars Ion and Neutral Particle Analyzer (MINPA; Kong et al. 2020) on board the Tianwen-1 orbiter. For the MINPA ENA detection, the ambient neutral atoms were ionized upon impact at an Al₂O₃ conversion surface and subsequently analyzed by an electrostatic analyzer (ESA) and time-of-flight (TOF) system. This design enables ENA measurements over an energy range from 57 eV to 3.5 keV with an energy resolution of approximately 68%, and MINPA is capable of resolving energetic neutral hydrogen, helium, and oxygen.

The instrument’s FOV is 360° × 9.7°. Advanced by the orbiter attitude and the MINPA mounting position, this FOV can cover the Sun-pointing direction. In this study, we focused on the hydrogen ENA measurements. By leveraging the advanced Extreme Ultraviolet (EUV) absorption technique (Wang et al. 2023), MINPA can detect ENAs directly from the Sun-pointing direction, enabling direct observations of penetrating hydrogen ENAs (Ma et al. 2025). For our analysis, we used the dataset from December 2021, when Tianwen-1 spent a substantial fraction of its orbit upstream of the Martian bow shock, providing favorable conditions for upstream solar wind hydrogen ENA observations. Due to the low count rates, the original 16 s data of hydrogen ENA fluxes and velocities were averaged over 10 minutes.

The hydrogen ENA count rates measured by MINPA were converted into physical particle fluxes using a dedicated data processing procedure that accounts for the instrument response and calibration characteristics. The performance of the MINPA hydrogen ENA detection unit, including the energy response, angular scattering at the conversion surface, and geometric factors, was characterized through extensive ground calibrations conducted at the MEsskammer fu¨r FlugzeitInStrumente und Time-Of-Flight (MEFISTO) facility at the University of Bern. Based on these calibration results and in-flight measurements, the conversion from count rates to hydrogen ENA fluxes was performed using a parameterized spectral fitting approach. The observed hydrogen ENA count rate spectra in individual viewing sectors were fitted to derive characteristic energies and count-rate amplitudes. The fitted energies were corrected for energy losses at the conversion surface, from which the bulk velocity of hydrogen ENAs was obtained, while the amplitudes were converted into physical fluxes using the corresponding geometric factors. Figures 1b and 1c show the corresponding hydrogen ENA velocities and fluxes during this period, revealing pronounced variability in both parameters.

2.3 Monte Carlo model

To obtain the sputtering escape rates of various atmospheric species, information on the sputtering yield is required, which is defined as the average number of escaping atmospheric molecules produced per incident high energetic particle (Johnson 1994). Given the stochastic nature of particle collisions and the involvement of multiple recoils, sputtering yields are most reliably evaluated with Monte Carlo simulations. This approach has been widely applied in studies of sputter-induced escape on solar system bodies (e.g., Shematovich et al. 2003; Leblanc et al. 2018; Snowden & Higgins 2021; Huang et al. 2024; Curry et al. 2025) and more recently extended to exoplanets (Gu et al. 2023a). In this study, we employed the test particle Monte Carlo model previously developed for the same planet (Ma et al. 2025), which is briefly summarized here for clarity.

The background atmosphere, consisting of CO₂, CO, N₂, and O, was adopted from Gu et al. (2025) over the altitude range 100–500 km and was assumed to be spherically symmetric. Below 100 km, the mean free path is sufficiently short that particles rapidly lose energy and thermalize to the local temperature within a distance not exceeding 10 m, consistent with local thermalization. At 500 km, the collision probability drops to ~1%, therefore, the atmosphere above this altitude is not expected to exert an appreciable influence on the derived sputtering yields. The Martian background atmosphere is discretized into 80 vertical layers, each with a depth of 5 km.

For a proton or hydrogen ENA entering the Martian upper atmosphere with a specified incident energy and incidence angle, we followed its trajectory step by step under Martian gravity until a collision occurred. Collisions with ambient atmospheric particles were implemented using three independent random numbers uniformly distributed between 0 and 1. The first (R₁) was used to sample the altitude of the next collision, z′, according to $${\int }_{z_0}^{z^{\prime }}\sum _i^{}{\sigma }_i{n}_i\mathrm{d}z=-\cos \beta \ln (1-R_1),$$(1)

where z is altitude, z₀ is the initial location, β is the angle between the particle’s velocity vector and the local vertical (downward) direction, and σ_i and n_i are the total collision cross section and number density of species i (CO₂, CO, N₂, and O), respectively.

Once a collision occurs, a second random number (R₂) was generated to select the collision target among CO₂, CO, N₂, and O. To this end, we defined the selection weight for species i as w_i = σ_iN_i∕Σ(σ_iN_i), where σ_i is the total collision cross section and N_i is the vertical column density of species i between z₀ and z′. The incident particle was taken to collide with CO₂ if R₂ < w_CO2, with O if w_CO2 ≤ R₂ < w_CO2 + w_O, and with the remaining species by continuing the same cumulative-sum procedure.

After the target species was specified, the scattering angle of the incident particle in the laboratory frame, θ_LAB, was determined using a third random number (R₃): $${\theta }_{LAB}=[R_3({\theta }_{max}^{{\alpha }_E}-{\theta }_{min}^{{\alpha }_E})+{\theta }_{min}^{{\alpha }_E}{]}^{1/{\alpha }_E}.$$(2)

Here, θ_min is the minimum scattering angle (taken to be 0.01 °), and θ_max is the maximum scattering angle, given by $$\begin{array}{rl}{\theta }_{max}& =\pi ({\rho }_m<1),\\ & =\arcsin (\frac{1}{{\rho }_m})({\rho }_m\ge 1),\end{array}$$(3)

where $${\rho }_m=\frac{m_a}{m_b\sqrt{1-\frac{m_a+m_b}{m_b}\frac{|Q|}{\epsilon }}}.$$(4)

In Eq. (4), m_a and m_b are the masses of the incident and target particles, respectively, ε is the kinetic energy of the incident particle, and Q is the kinetic energy converted into internal energy in an inelastic collision. In the present model, ionization, excitation, and dissociation were included only as inelastic energy-loss channels. That is, they reduced the kinetic energy available to the post-collision particles, but the chemical products of these processes were not explicitly followed. For example, when CO₂ dissociation occurred, the subsequent CO, O, or C fragments were not tracked as separate particles in the Monte Carlo simulation. The model therefore follows the trajectories and energies of the incident particle and recoiling atmospheric target particles, while treating internal excitation, ionization, and dissociation as energy losses. We computed Q following the formalism of Firsov (1959): $$Q(eV)=\frac{4.3\times {10}^{-8}v(Z_A+Z_B{)}^{5/3}}{[1+0.3b(Z_A+Z_B{)}^{1/3}{]}^5},$$(5)

where b is the impact parameter in units of Bohr radius, v is the relative velocity in units of centimeters per second, and Z_A and Z_B are nuclear charge numbers of the incident and target species, respectively.

The scattering angle in the center-of-mass frame, θ_CM, is related to θ_LAB through $$\begin{array}{rl}{\theta }_{CM}=\arccos [& -{\rho }_m+{\rho }_m{\cos }^2({\theta }_{LAB})\\ & +\cos ({\theta }_{LAB})\sqrt{1-{\rho }_m^2{\sin }^2({\theta }_{LAB})}].\end{array}$$(6)

The post-collision velocities of the two collision partners were determined from conservation of momentum and energy.

The above procedure was repeated until one of the following termination criteria was met: (1) the incident proton reaches the lower boundary, or its kinetic energy falls below the local escape energy through a cascade of collisions, in which case it is no longer tracked; (2) the incident proton reaches the upper boundary, where it is assumed to escape the atmosphere if its kinetic energy exceeds the local escape energy, or elastically reflected downward otherwise. To ensure statistical robustness, each simulation included 100 000 incident protons or hydrogen ENAs. A total of 171 simulations were performed, spanning 11 energy levels for hydrogen ENAs and eight for protons, with nine incidence angles per energy level uniformly sampled from 0° to 80°.

We next describe the cross sections adopted to characterize elastic scattering and charge exchange between energetic particles and ambient neutral species. Elastic scattering cross sections for H colliding with CO₂, CO, N₂, and O were adopted from Lewkow & Kharchenko (2014). Elastic scattering cross sections for H⁺-O and H⁺-N₂ were taken from Noël & Prölss (1993), and for H⁺-CO and H⁺-CO₂, we assumed the same elastic scattering cross sections as for H⁺-N₂. For other collision pairs, elastic scattering cross sections below 10 eV were adopted from Gu et al. (2020) and are were to energies above 10 eV by scaling their relative magnitudes to match the energy dependence reported in Lewkow & Kharchenko (2014). Charge-exchange cross sections for H⁺ with O and N₂ were adopted from Lindsay & Stebbings (2005). Similarly, we assumed that the charge-exchange cross sections for H⁺-CO and H⁺-CO₂ are identical to those for H⁺-N₂. The scattering angle distribution was parameterized by α_E in Eq. (2). This parameter was obtained by fitting laboratoryframe differential scattering cross sections using the reduced coordinates ρ and τ (Noël & Prölss 1993): $$\log (\rho )={\alpha }_0+{\alpha }_E\log (\tau ),$$(7) $$\rho (\tau )={\theta }_{LAB}\sin ({\theta }_{LAB})\frac{\mathrm{d}\sigma }{\mathrm{d}\mathrm{\Omega }}({\theta }_{LAB},\epsilon ),$$(8) $$\tau =\epsilon {\theta }_{LAB},$$(9)

where ε is the kinetic energy of the incident particle, θ_LAB is the scattering angle in the laboratory frame, σ is the total scattering cross section, and dσ/dΩ is the differential scattering cross section per unit solid angle Ω. Specifically, at an incident energy of 1 keV, α_E for elastic scattering, α_E is generally in the range −1.2 to −1.0, corresponding to strongly forward scattering.

3 Results

3.1 Sputtering yields

Figure 2 displays the sputtering yields of CO₂ and O induced by precipitating solar wind protons and hydrogen ENAs as functions of incidence angle and incident velocity (or energy). Here, an incidence angle of 0° corresponds to precipitating particles enter the atmosphere perpendicularly. The figure reveals several interesting features.

First, our calculations yield sputtering yields of 2.8×10⁻⁵−2×10⁻² for CO₂ and 8×10⁻⁶−3.7×10⁻³ for O induced by hydrogen ENAs, and 3×10⁻⁶−5×10⁻² for CO₂ and 1×10⁻⁶−6.5×10⁻³ for O induced by precipitating solar wind protons. Overall, the sputtering yields induced by precipitating solar wind protons exhibit a broader range than those induced by hydrogen ENAs. This behavior arises from the wider incident energy distribution of solar wind protons relative to hydrogen ENAs. The latter mainly originate from the charge exchange between solar wind protons and hydrogen atoms in the Martian exosphere, and thus their velocities are approximately equal to the solar wind speed, with a typical thermal broadening of about 10 eV. In contrast, precipitating solar wind protons are decelerated while passing through the Martian magnetosheath, resulting in a wider range of energies, 10¹–10⁴ eV, as illustrated in Fig. 1. Therefore, in the subsequent escape rate calculations, we adopted a representative single energy for hydrogen ENAs, whereas for precipitating solar wind protons, the flux and sputtering yields were integrated over their entire energy range.

Second, both CO₂ and O sputtering yields increase with increasing incidence angle and incident velocity (or energy). For example, the CO₂ sputtering yield induced by hydrogen ENAs rises from 2.8×10⁻⁵ to 2.3×10⁻⁴ as the incident velocity increases from 280 to 800 km s⁻¹ at a fixed incidence angle of 0° Similarly, at an incident velocity of 540 km s⁻¹, the CO₂ sputtering yield increases from 1.3×10⁻⁴ to 1.3×10⁻² as the incidence angle increases from 0° to 80°. This behavior reflects that particles with higher kinetic energies generate more atmospheric recoils during energy degradation (e.g., Johnson et al. 2000; Gu et al. 2019, 2023a). Moreover, oblique incidence angles confine hydrogen ENAs and recoil trajectories to higher altitudes, where the lower atmospheric density enhances the escape probability (Johnson et al. 2000; Huang et al. 2024, 2025).

Third, the sputtering yields of CO₂ are generally larger than those of O, with this contrast becoming more pronounced at higher incident velocities (or energies). Specifically, at a fixed incidence angle of 0°, the CO₂ and O sputtering yields induced by hydrogen ENAs are 2.8×10⁻⁵ and 8×10⁻⁶, respectively, for an incident velocity of 280 km s⁻¹, while the corresponding values are 2.3×10⁻⁴ and 4.7×10⁻⁵ at 800 km s⁻¹. This behavior can be explained by the altitude-dependent atmospheric composition. At oblique incidence angles and low velocities, hydrogen ENAs and recoil trajectories are confined to higher altitudes where the O/CO₂ mixing ratio is relatively high. Conversely, at higher speeds, they penetrate deeper into the atmosphere, where the ratio decreases rapidly above the Martian homopause.

For illustration, Fig. 3a shows the production yields of all ambient particles for hydrogen ENAs with an incident angle of 50° and an incident velocity of 400 km s⁻¹, while Fig. 3b presents the production yields of above-threshold O and CO₂ recoils, defined as particles whose nascent kinetic energies exceed the local escape energies of 1.8 and 4.8 eV, respectively. Figure 3c further shows the altitude distribution of the ultimately escaping O and CO₂ particles, i.e., particles that retain sufficient kinetic energy after subsequent collisional transport and finally escape through the upper boundary. Figure 3a demonstrates that the production yield decreases rapidly with decreasing altitude as the background gas density increases approximately exponentially, indicating that the maximum penetration depth is below 100 km.

The comparison between Figs. 3b and 3c highlights the important distinction between above-threshold recoil production and actual escape. Although a substantial fraction of abovethreshold O and CO₂ recoils is produced below the exobase, these particles do not necessarily escape. Particles generated in dense regions may undergo additional collisions, lose kinetic energy, change direction, and eventually become thermalized before reaching the exobase. As a result, the production yield of ultimately escaping particles shown in Fig. 3c is much smaller than the above-threshold production yield shown in Fig. 3b. The altitude dependence in Fig. 3c indicates that the final escaping population is preferentially associated with higher altitudes, where collisional attenuation is weaker. This result demonstrates that a nascent energy exceeding the local escape energy is only a necessary condition for escape, not a sufficient one. The final escape probability is controlled by both the initial recoil energy and the subsequent collisional history. Thus, the escape rates reported below were not obtained by directly integrating the above-threshold production in Fig. 3b, but by following particles through the full Monte Carlo transport process until they either escaped, thermalized, or were lost through downward transport.

Furthermore, the calculated total production yield is on the order of unity, comparable to the results reported in Lewkow & Kharchenko (2014). Figures 3b and 3c indicate that collisional filtering is strongest at lower altitudes, whereas above roughly 150 km a large fraction of above-threshold recoils can survive subsequent collisions and escape in this particular simulation. Thus, the escape probability increases with altitude as the probability of further collisions decreases, rather than being controlled only by production near the nominal exobase.

Figure 4 shows the nascent energy distributions of O and CO₂ recoils in the 1–10 eV range produced at altitudes of 100– 500 km. The distributions are dominated by low-energy recoils, but a weak high-energy tail extends beyond the local CO₂ escape energy. Although this tail represents only a small fraction of the recoil population, it is the component most relevant to the production of above-threshold CO₂ recoils and to the residual escaping fraction.

Fig. 2 Distributions of the O and CO2 sputtering yields induced by precipitating lar wind protons (upper row) and hydrogen ENAs (bottom row) as functions of incidence angle and incident velocity (or energy).

Fig. 3 Production and escape yield profiles for hydrogen ENAs with an incident angle of 50° and an incident velocity of 400 km s−1. (a) Total recoil production yield of all ambient particles per incident H ENA per unit altitude. (b) Production yield of above-threshold O and CO2 recoils, defined as particles whose nascent kinetic energies exceed the local escape energy. (c) Production altitude distribution of ultimately escaping O and CO2 particles after subsequent collisional transport. The difference between panels b and c illustrates the strong collisional filtering of hot recoils produced in dense atmospheric regions.

Fig. 4 Energy distributions of nascent O and CO2 recoils in the 1–10 eV range. The vertical dashed lines indicate the corresponding escape energies.

3.2 Comparison with Lewkow and Kharchenko (2014)

It is useful to compare our results with those of Lewkow & Kharchenko (2014), who performed Monte Carlo simulations of atmospheric sputtering on Mars induced by precipitating hydrogen ENAs. They concluded that heavy species such as CO₂ and Ar have no realistic probability of escape, whereas our calculations reveal a nonzero CO₂ sputtering yield under the precipitating proton and hydrogen ENA conditions considered here. This apparent difference is important because it concerns whether light incident particles can contribute to the escape of heavy neutral species from Mars.

A key point in this comparison is the distinction between above-threshold recoil production and actual escape. In a collision cascade, some CO₂ recoils may be produced with nascent kinetic energies exceeding the local escape energy. However, this condition alone does not guarantee escape. If these particles are produced deep in the atmosphere, they may undergo subsequent collisions, lose kinetic energy, change direction, and become thermalized before reaching the exobase. Therefore, the production of above-threshold CO₂ recoils should not be interpreted as the CO₂ sputtering yield.

Our Monte Carlo calculation follows each energetic incident particle and recoil through its subsequent collisional history until it reaches the lower boundary, becomes energetically unable to escape, or crosses the upper boundary with kinetic energy exceeding the local escape energy. The CO₂ sputtering yield reported here therefore represents only the subset of recoils that survive this collisional transport process. This distinction is illustrated directly in Fig. 3. Figure 3b shows the production of above-threshold O and CO₂ recoils, whereas Fig. 3c shows the much smaller subset that ultimately escapes. The large reduction from Figs. 3b to 3c demonstrates that most hot recoils, especially those produced at low altitudes, are removed by subsequent collisions and thermalization before they can escape. This behavior is consistent with the physical argument of Lewkow & Kharchenko (2014) that heavy recoils generated deep in the atmosphere have a very low escape probability.

The main difference between the two calculations appears to arise before collisional transport, namely in the nascent CO₂ recoil energy distribution. In Lewkow & Kharchenko (2014), the production of nascent CO₂ molecules with energies exceeding the escape energy is essentially negligible. By contrast, our Monte Carlo treatment produces a small but finite high-energy tail. Figure 4 supports this interpretation by showing that the nascent CO₂ recoil distribution is dominated by low-energy particles, while a weak high-energy tail extends beyond the local CO₂ escape energy. Although this tail represents only a small fraction of the total nascent CO₂ population, it is the component relevant to the above-threshold recoil production and to the final escaping fraction after collisional attenuation.

This difference is likely related to the adopted angular distribution of the differential scattering cross section. Rare large-angle scattering events can transfer substantially more kinetic energy to CO₂ than the mean collision, and the probability assigned to this angular tail can therefore strongly affect the production of above-threshold CO₂ recoils. Figure 5 compares the differential cross sections for 1 keV H–CO₂ collisions in the center-of-mass frame. The distribution from Lewkow & Kharchenko (2014) was calculated using their analytical formula, whereas the distribution used in the present study was estimated from 1 000 000 randomly sampled scattering events using Eq. (2), with a characteristic value of α_E = −0.85. Although the two distributions are broadly similar over most scattering angles, the sampled distribution gives a slightly larger probability at large scattering angles. This small difference is important for heavy-species escape because the escaping CO₂ population is controlled by rare, high-energy recoil events rather than by the mean energy transfer.

Several model differences may also contribute to the quantitative discrepancy. The present study uses precipitating solar wind proton fluxes derived from MAVEN/SWIA and hydrogen ENA fluxes derived from Tianwen-1/MINPA, whereas Lewkow & Kharchenko (2014) adopted prescribed incident hydrogen ENA conditions. Differences in the background atmosphere, incident energy distribution, angular distribution, and collision treatment may also affect the altitude distribution and energy spectrum of nascent recoils. However, because CO₂ escape is controlled by the small high-energy tail of the recoil distribution, the treatment of energy transfer during individual collisions, especially the angular distribution of the differential cross section, is likely a particularly important factor.

Finally, we emphasize that our results should be viewed as complementary to, rather than in direct contradiction with, Lewkow & Kharchenko (2014). Both studies indicate that hydrogen ENA/proton-induced escape of heavy species is strongly limited by collisional thermalization. The main difference lies in the rare high-energy tail of the nascent CO₂ recoil distribution: in our calculation, a small fraction of CO₂ recoils can acquire energies above the local escape threshold and, after subsequent collisional filtering, survive to escape under the observationally constrained precipitating particle conditions examined here. Accordingly, the CO₂ sputtering yield reported in this work should be regarded as the final escaping component of this rare high-energy recoil population after thermalization and collisional attenuation, rather than as the total production of energetic or above-threshold CO₂ recoils.

Fig. 5 Comparison of differential cross sections for 1 keV H–CO2 collisions in center-of-mass frame. The values indicated by the solid blue line were derived from the analytical formula of Lewkow & Kharchenko (2014), while the values indicated by the dashed orange line were estimated from 1 000 000 randomly sampled scattering events using Eq. (2), with a characteristic value of αE = −0.85.

3.3 Sputtering-induced escape fluxes

By combining the calculated sputtering yields of CO₂ and O with observations of precipitating solar wind protons and hydrogen ENAs, we can estimate the corresponding escape rates driven by atmospheric sputtering. For a given species i, its hemispheric escape rate L_ij induced by energetic particles j is given by $$L_{ij}=2\pi {\int }_{0^{\circ }}^{{90}^{\circ }}(R_{Mars}+z_{upper}{)}^2{Y}_{ij}(\theta ){\mathrm{\Phi }}_j\sin \theta d\theta ,$$(10)

where R_Mars = 3400 km is the solid body radius of Mars, z_upper = 500 km is upper boundary, θ is the incidence angle, Y_ij is the sputtering yields of species i induced by j, and Φ j is the incident flux of j.

Figure 6 presents the CO₂ and O escape rates induced by precipitating solar wind protons (panel a) and hydrogen ENAs (panel b) during December 2021. According to our calculations, the median CO₂ and O escape rates are 1×10²² s⁻¹ and 1×10²¹ s⁻¹ for precipitating solar wind protons, and 1.4×10²² s⁻¹ and 8× 10²¹ s⁻¹ for hydrogen ENAs. The figure clearly shows that for hydrogen ENAs, the sputter-induced escape rate of CO₂ generally exceeds that of O, whereas for precipitating solar wind protons, the escape rates of the two species are comparable. However, during intervals of strong precipitating proton flux at low energies, the O escape rate becomes comparable in magnitude to that of CO₂. Specifically, on December 3, when the low-energy precipitating proton flux is significantly enhanced, the daily averaged CO₂ and O escape rates driven by precipitating proton sputtering are 3×10²² s⁻¹ and 7×10²¹ s⁻¹, respectively. This trend is consistent with the fact that precipitating solar wind protons possess higher average incident energies than hydrogen ENAs.

It is also necessary to compare the relative importance of these two sputtering sources. This comparison is shown in Figs. 6c and 6d, which present the daily averaged CO₂ and O escape rates induced by the two sources described above. The figures show that precipitating solar wind protons generally produce stronger CO₂ escape, whereas hydrogen ENAs drive stronger O escape. However, hydrogen ENAs can produce comparable or even greater escape rates during periods of abnormally enhanced ENA flux, as indicated by the spike-like features in Fig. 1c on December 9, 16, and 17.

The sputter-induced escape rates derived here were compared with those reported in previous studies. Specifically, using a precipitating hydrogen ENA flux of 9.2×10⁷ cm⁻² s⁻¹, which is one order of magnitude larger than the incident flux considered in this study, Lewkow & Kharchenko (2014) obtained an O escape rate of ~8×10²² s⁻¹. More recently, Shematovich (2017) estimated an O escape rate of 8×10²³ s⁻¹ to 4.5×10²⁴ s⁻¹ for comparable hydrogen ENA precipitation conditions. Using a modeled spectrum of precipitating solar-wind protons (Shematovich 2017), which is larger than the SWIA-based spectrum adopted here, that study also reported an O escape rate of 1.1×10²³ s⁻¹ driven by precipitating solar-wind protons. Given that we adopted a substantially smaller hydrogen ENA and proton flux, these published estimates are broadly consistent with our results for O loss. Interestingly, previous studies have largely focused on sputtering-driven O loss and did not explicitly quantify CO₂ escape, whereas our calculations indicate that CO₂ escape is generally stronger than O under the same precipitation conditions.

4 Discussion and concluding remarks

This study quantifies atmospheric sputtering on Mars driven by hydrogen ENAs and precipitating solar wind protons, using Monte Carlo simulations constrained by Tianwen-1 and MAVEN observations. Both particle populations contribute to the nonthermal escape of CO₂ and O, though their relative importance depends strongly on the incident energy spectrum and geometry. The calculated sputtering yields of CO₂ and O increase with both incidence angle and incident velocity (or energy), consistent with enhanced recoil production and greater escape probabilities in the rarefied upper atmosphere. Precipitating solar wind protons, characterized by a broader and higher-energy spectrum, generally produce stronger CO₂ and O escape, whereas hydrogen ENAs can generate comparable or even higher escape rates during episodes of abnormally enhanced ENA fluxes. This highlights the dynamic coupling between solar wind variability and atmospheric erosion efficiency.

The sputtering source examined here should be distinguished from sputtering driven by pickup O⁺ ions, which has been the focus of many previous studies. Because of their larger mass and, in many cases, higher characteristic energies, pickup O⁺ ions can transfer momentum to atmospheric species more efficiently and are therefore expected to be more effective at ejecting heavy neutrals such as CO₂ and Ar. By contrast, the present study focused specifically on precipitating solar wind protons and hydrogen ENAs. These particles are lighter and less efficient per incident particle, but their precipitation is directly linked to the variable solar wind interaction with Mars and can therefore provide a complementary, condition-dependent sputtering source. Accordingly, the escape rates reported here should not be interpreted as a complete estimate of all sputtering escape at Mars, but rather as an observationally constrained estimate of the proton/hydrogen ENA-induced component.

Our results also show that heavy-species escape by pro-ton/hydrogen ENA-induced sputtering is strongly regulated by collisional thermalization. Although CO₂ recoils with nascent energies exceeding the local escape energy can be produced, only a small fraction ultimately survive subsequent collisions and escape. This distinction is important when comparing with Lewkow & Kharchenko (2014), who emphasized that heavy species such as CO₂ and Ar have very low escape probabilities under hydrogen ENA precipitation. Our calculations are consistent with the general physical picture that heavy recoils produced deep in the atmosphere are efficiently thermalized. The nonzero CO₂ escape rates obtained here represent the residual escaping fraction after strong collisional filtering, rather than the total production rate of energetic or above-threshold CO₂ recoils.

Our calculations suggest equivalent total O and C escape rates of 10²²–10²³ s⁻¹, driven by hydrogen ENA and precipitating solar wind proton sputtering. For comparison, photochemical escape, primarily through dissociative recombination of O₂⁺, is the dominant present-day source of atomic oxygen loss on Mars, yielding escape rates on the order of 10²⁵–10²⁶ s⁻¹ (e.g., Fox & Hac´ 2009, 2014, 2018; Cravens et al. 2017; Lillis et al. 2017). Photochemical escape of carbon, produced mainly by photolysis of CO₂ and CO, contributes a total C escape rate of approximately 10²³–10²⁴ s⁻¹ (e.g., Fox & Bakalian 2001; Cui et al. 2019; Lo et al. 2021; Gu et al. 2023b; Ramstad et al. 2015). Ion pickup and ion outflow processes, which transport charged species along open magnetic field lines, provide additional losses of 10²³– 10²⁵ s⁻¹ for O⁺ or O₂⁺ and 10²²–10²⁴ s⁻¹ for CO₂⁺ (e.g., Barabash et al. 2007; Nilsson et al. 2010; Maes et al. 2021; Li et al. 2023; Song et al. 2025). Although the sputtering escape rates for O inferred in this study and previous work (e.g., Lewkow & Kharchenko 2014; Leblanc et al. 2018; Ma et al. 2025) are generally smaller, typically on the order of 10²²–10²³ s⁻¹ under moderate solar activity, this process may at times dominate the escape of C on Mars by driving CO₂ loss. Moreover, sputtering remains uniquely important because it can eject heavy neutral species, particularly Ar (e.g., Curry et al. 2025), that are largely inaccessible to photochemical or ion processes. Over geological timescales, this mechanism may have contributed a substantial fraction of the ancient Martian atmospheric loss.

Overall, our results indicate that precipitating solar wind protons and hydrogen ENAs make comparable contributions to nonthermal escape on Mars. Coordinated Tianwen-1 and MAVEN observations provide an effective means to constrain these mechanisms simultaneously. The present work therefore complements previous studies of pickup-ion sputtering by isolating the contribution from light incident particles, namely proton and hydrogen ENAs, under observationally constrained conditions. Future missions equipped with higher-sensitivity ENA imagers and ion analyzers will be critical for quantifying the variability, efficiency, and cumulative impact of sputtering on Mars’ atmospheric evolution and long-term habitability.

Fig. 6 CO2 and O escape rates driven by precipitating solar wind protons and hydrogen ENAs. (a) CO2 and O escape rates driven by precipitating solar wind protons for each MAVEN orbit. (b) CO2 and O escape rates driven by hydrogen ENAs for each Tianwen-1 orbit. (c) Comparison of daily averaged CO2 escape rates driven by precipitating solar wind protons and hydrogen ENAs, with the dashed line indicating equal escape rates. (d) Same as panel (c) but for O escape.

Acknowledgements

The authors acknowledge supports from the National Natural Science Foundation of China through grants 42241112, 42522803 and 42475134.

References

All Figures

	Fig. 1 Precipitating solar wind protons and hydrogen ENAs at Mars during December 2021, shown for each MAVEN and Tianwen-1 orbit. (a) Energy spectrum of precipitating solar wind protons. (b) Hydrogen ENA velocity. (c) Hydrogen ENA flux.
In the text

	Fig. 2 Distributions of the O and CO2 sputtering yields induced by precipitating lar wind protons (upper row) and hydrogen ENAs (bottom row) as functions of incidence angle and incident velocity (or energy).
In the text

Fig. 3 Production and escape yield profiles for hydrogen ENAs with an incident angle of 50° and an incident velocity of 400 km s−1. (a) Total recoil production yield of all ambient particles per incident H ENA per unit altitude. (b) Production yield of above-threshold O and CO2 recoils, defined as particles whose nascent kinetic energies exceed the local escape energy. (c) Production altitude distribution of ultimately escaping O and CO2 particles after subsequent collisional transport. The difference between panels b and c illustrates the strong collisional filtering of hot recoils produced in dense atmospheric regions.

In the text

	Fig. 4 Energy distributions of nascent O and CO2 recoils in the 1–10 eV range. The vertical dashed lines indicate the corresponding escape energies.
In the text

	Fig. 5 Comparison of differential cross sections for 1 keV H–CO2 collisions in center-of-mass frame. The values indicated by the solid blue line were derived from the analytical formula of Lewkow & Kharchenko (2014), while the values indicated by the dashed orange line were estimated from 1 000 000 randomly sampled scattering events using Eq. (2), with a characteristic value of αE = −0.85.
In the text

Fig. 6 CO2 and O escape rates driven by precipitating solar wind protons and hydrogen ENAs. (a) CO2 and O escape rates driven by precipitating solar wind protons for each MAVEN orbit. (b) CO2 and O escape rates driven by hydrogen ENAs for each Tianwen-1 orbit. (c) Comparison of daily averaged CO2 escape rates driven by precipitating solar wind protons and hydrogen ENAs, with the dashed line indicating equal escape rates. (d) Same as panel (c) but for O escape.

In the text