1 Introduction

Mars is enveloped by an atomic hydrogen exosphere detectable at Lyman-$$\alpha $$ wavelengths (e.g., Bhattacharyya et al., 2015, 2020; Chaffin et al., 2015, 2018; Chaufray et al., 2008; Chaufray, Gonzalez-Galindo, et al., 2021; Clarke et al., 2017; Mayyasi et al., 2017). In the absence of a global dipole magnetic field, newly ionized hydrogen atoms within the exosphere can be accelerated by the solar wind convective electric field, a process known as ion pickup, enabling their eventual escape to space (e.g., Barabash et al., 1991; Bertucci et al., 2013; Connerney et al., 2015; Dubinin et al., 2006; Rahmati et al., 2017; Romanelli et al., 2016; Q. Xu et al., 2020; Yamauchi et al., 2015). Aside from the pickup processes occurring in the solar wind and magnetosheath, acceleration via the $$\mathbf{J}\times \mathbf{B}$$ force in the magnetotail provides another major pathway for planetary ion escape (Dubinin et al., 2011; Lundin, 2011; Y. J. Ma et al., 2019). These ion escape mechanisms, coupled with neutral loss processes such as Jeans escape (Gu et al., 2022), may have collectively contributed to the long-term climate evolution of Mars (e.g., Baker, 2001; Jakosky & Phillips, 2001; Villanueva et al., 2015). Currently, the major mechanism for atmospheric hydrogen loss on Mars is likely to be thermal (Jeans) escape of neutral atoms (Jakosky, 2021; Jakosky et al., 2018), with the Mars Atmosphere and Volatile EvolutioN (MAVEN; Bougher et al., 2015; Jakosky et al., 2015; Lillis et al., 2015) mission quantifying the neutral hydrogen escape rate and its temporal variability (e.g., Halekas, 2017; Mayyasi et al., 2023; Rahmati et al., 2018). Nevertheless, few observational studies focused on the $${\mathrm{H}}^{+}$$ escape rate (He et al., 2024) due to the difficulty of distinguishing planetary protons from solar wind protons in spacecraft data.

Efforts from model calculations have established a basic understanding of the $${\mathrm{H}}^{+}$$ escape rates on Mars. Quantitatively, Jarvinen et al. (2016), Sun et al. (2023) and Chaufray et al. (2024) simulated the $${\mathrm{H}}^{+}$$ loss driven by solar wind interactions and found that the $${\mathrm{H}}^{+}$$ loss rate is 1–2 orders of magnitude lower than the neutral escape rate, implying a secondary role under current conditions. However, under the impact of the 2017 September space weather events, Romanelli et al. (2018) modeled a tenfold increase in the global escape rate of planetary $${\mathrm{H}}^{+}$$, matching the contemporaneous thermal escape rate of neutral hydrogen (Lee et al., 2018; Mayyasi et al., 2018). Notably, while studies by Romanelli et al. (2018) and Chaufray et al. (2024) incorporated the charge exchange reactions between planetary neutrals and solar wind protons, which is a key ionization mechanism for the extended hydrogen corona (Rahmati et al., 2018), this process was omitted in the works of Sun et al. (2023) and Jarvinen et al. (2016).

The relative importance of different heavy ion escape channels has been extensively investigated through simulations (Curry et al., 2013; Curry, Luhmann, Ma, Dong, et al., 2015; Curry, Luhmann, Ma, Liemohn, et al., 2015; Jakosky et al., 2018; Jarvinen & Kallio, 2014; Sakai et al., 2023; Sakakura et al., 2022; Sakata et al., 2024). However, the dominant escape mechanisms for planetary $${\mathrm{H}}^{+}$$ remain poorly constrained. To address this gap, the objectives of this study are to assess the relative importance of the two primary $${\mathrm{H}}^{+}$$ escape channels (pickup ion escape and ion outflow through the magnetotail), determine the primary source of the escaping $${\mathrm{H}}^{+}$$, and identify the key factors controlling $${\mathrm{H}}^{+}$$ escape rates.

2 Model Descriptions

To avoid artificial effects from the outer boundary, we employ a large computational domain spanning $${-}\mathrm{40}\ {\mathrm{R}}_{\mathrm{M}}\le \mathrm{X},\mathrm{Y},\mathrm{Z}\le \mathrm{40}\ {\mathrm{R}}_{\mathrm{M}}$$ in the Mars Solar Orbital (MSO) coordinate system, in which the X-axis is directed from Mars toward the Sun, the Z-axis aligns with Mars' rotational north pole, and the Y-axis completes the right-handed coordinate system. A nonuniform spherical grid is adopted, with a resolution of $$3{}^{\circ}$$ in both longitude and latitude. The inner boundary of the domain is set at an altitude of $$\mathrm{100}\ \mathrm{k}\mathrm{m}$$ above the surface of Mars, with a radial resolution of $$\mathrm{9.5}\ \mathrm{k}\mathrm{m}$$ near this boundary to resolve the variations of ion densities with scale heights (Y. Ma et al., 2004). Above $$\mathrm{600}\ \mathrm{k}\mathrm{m}$$, the grid spacing increases radially outward to balance resolution and computational efficiency. For all simulated cases in Table 1, the upstream solar wind temperature and the interplanetary magnetic field (IMF) are held constant at $$\mathrm{3.5}\times \mathrm{1}{\mathrm{0}}^{\mathrm{5}}\ \mathrm{K}$$ and a $$\mathrm{3}\ \mathrm{n}\mathrm{T}$$ Parker spiral configuration, respectively. The crustal field is fixed on the dayside, with the subsolar point positioned at 180$${}^{\circ}$$W, 0$${}^{\circ}$$N, except in Case 3. For each fluid, the ion escape rate is estimated by integrating the ion flux across a spherical surface at $$\mathrm{12}\ {\mathrm{R}}_{\mathrm{M}}$$ (Y.-J. Ma & Nagy, 2007; Y. Ma et al., 2014; Y. J. Ma et al., 2015, 2017; Y. Ma, Russell, et al., 2018; M. Wang et al., 2024).

3 Results

Figure 1 displays the simulated distributions of magnetic field $$(\mathbf{B})$$, solar wind convective electric field ($${\mathbf{E}}_{\mathbf{S}\mathbf{W}}=-{\mathbf{U}}_{\mathbf{S}\mathbf{W}}\times \mathbf{B}$$, where $${\mathbf{U}}_{\mathbf{SW}}$$ represents the solar wind proton velocity) and ion densities for Simulation Case 1. The transition from solar wind proton-dominated flows to the regime governed by planetary ions is marked by the ion composition boundary (ICB; Halekas et al., 2018; J. Wang et al., 2020; S. Xu et al., 2023), which is a critical layer that governs ionospheric outflow and solar wind penetration, as demarcated by the black line in Panel (c). The $${\mathrm{E}}_{\mathrm{S}\mathrm{W}}$$ is $${\sim} \mathrm{1}\ \mathrm{m}\mathrm{V}/\mathrm{m}$$ upstream of the bow shock (Panel b), and intensifies severalfold within the magnetosheath, except in the subsolar region. Upon crossing the ICB, the $${\mathrm{E}}_{\mathrm{S}\mathrm{W}}$$ drops abruptly due to the deflection of solar wind protons $$\left({\mathrm{H}}_{\mathrm{S}\mathrm{W}}^{+}\right)$$, as seen in Panels (b) and (c). As a fundamental plasma interface, the dynamics of ICB location regulate the efficiency of pickup ion escape, as the ion motions above the ICB location are governed by the $${\mathbf{E}}_{\mathbf{SW}}$$ as shown in Panels (d-f). For $${\mathrm{O}}_{2}^{+}$$ and $${\mathrm{O}}^{+}$$, the asymmetric density distributions along the Z-axis arise from the large curvature radii ($${\sim} \rho $$, $$\rho $$ is Larmor radius) of their $${\mathbf{E}}_{\mathbf{SW}}$$-driven cycloidal motion: $${\mathrm{n}}_{\mathrm{i}}{\mathrm{q}}_{\mathrm{i}}\left({\mathbf{E}}_{\mathbf{S}\mathbf{W}}+{\mathbf{U}}_{\mathbf{i}}\times \mathbf{B}\right)$$, where $${\mathrm{n}}_{\mathrm{i}}$$, $${\mathrm{U}}_{\mathrm{i}}$$ and $${\mathrm{q}}_{\mathrm{i}}$$ are density, velocity and charge for the ions, respectively. In contrast, $${\mathrm{H}}^{+}$$ exhibit tailward escape trajectories aligned with solar wind flow (Panel d), as their pickup gyroradii are smaller than the characteristic size of the Martian induced magnetosphere. The red streamtraces reveal that $${\mathrm{O}}_{2}^{+}$$ are efficiently accelerated by the strong $${\mathrm{E}}_{\mathrm{S}\mathrm{W}}$$ at the ICB, driving the formation of an ion plume (Y. Dong, Fang, et al., 2015; Qiao et al., 2024) and confirming their ionospheric origin, while the pickup $${\mathrm{H}}^{+}$$ originate outside the ICB, because $${\mathrm{H}}^{+}$$ accelerated at the ICB can recirculate back into the ionosphere due to their small gyroradii. For pickup $${\mathrm{O}}^{+}$$, the origins include both ionospheric ions and the ionized exospheric particles (Panel f). Although the pickup ions outside the ICB may be partially deflected into the magnetotail by local electromagnetic fields (Y. Dong et al., 2023; S. Xu et al., 2021), the bulk transport within the magnetotail remain governed by dense, slow-moving ionospheric outflows.

The escape fluxes for Case 1 are presented in Figure 2, where the ICB demarcates two distinct escape regions: magnetosheath and solar wind dominated by pickup ions and the magnetotail dominated by ionospheric outflow. The escaping $${\mathrm{O}}_{2}^{+}$$ and $${\mathrm{O}}^{+}$$ are predominantly concentrated in the magnetotail, consistent with observational studies (e.g., Barabash et al., 2007; Y. Dong et al., 2017; Inui et al., 2019; Nilsson et al., 2023; Ramstad et al., 2015; Schnepf et al., 2024). In contrast, $${\mathrm{H}}^{+}$$ fluxes exhibit an asymmetric spatial distribution, driven by the $${\mathbf{E}}_{\mathbf{S}\mathbf{W}}\times \mathbf{B}$$ drift in $${-}$$Y direction, and are predominantly concentrated within the magnetosheath and solar wind regions, indicating that pickup acceleration is the dominant escape channel. The $${\mathrm{H}}^{+}$$ escape rates through the two regions separated by the ICB are calculated by integrating ion fluxes at the cross-sectional plane of $$\mathrm{X}=-\mathrm{3}\ {\mathrm{R}}_{\mathrm{M}}$$. The low escape rate of $${\mathrm{H}}^{+}$$ through the magnetotail $$\left(\sim \mathrm{0.5}\times \mathrm{1}{\mathrm{0}}^{\mathrm{24}}\ {\mathrm{s}}^{-\mathrm{1}}\right)$$, which contribute merely $${\sim} \mathrm{5}\mathrm{\%}$$ of the total escape rate in Table 1 Case 1, further confirm that ionospheric outflow plays a minor role in hydrogen escape under current solar wind conditions.

To identify the dominant source of escaping $${\mathrm{H}}^{+}$$, we analyze a scenario where the solar wind charge exchange with hydrogen exosphere $$\left({\mathrm{H}}_{\mathrm{S}\mathrm{W}}^{+}-\mathrm{H}\right)$$ is deactivated (Case 2). The omitted $${\mathrm{H}}^{+}$$ production rate by $${\mathrm{H}}_{\mathrm{S}\mathrm{W}}^{+}-\mathrm{H}$$ charge exchange (Equation 1), mapped in Figure 3b, is calculated based on the solar wind flux (Figure 3a) and the neutral density (Figure S1 in Supporting Information S1). Within a spherical volume of $$10\ {\mathrm{R}}_{\mathrm{M}}$$, the $${\mathrm{H}}_{\mathrm{S}\mathrm{W}}^{+}-\mathrm{H}$$ reaction yields a total $${\mathrm{H}}^{+}$$ production rate of $$\mathrm{6}\times \mathrm{1}{\mathrm{0}}^{\mathrm{24}}\ {\mathrm{s}}^{-\mathrm{1}}$$. As expected, the simulated escape rate drops by $$\mathrm{6}\times \mathrm{1}{\mathrm{0}}^{\mathrm{24}}\ {\mathrm{s}}^{-\mathrm{1}}$$, demonstrating a high efficiency for charge exchange-generated $${\mathrm{H}}^{+}$$ to escape, with negligible recycling into the ionosphere.

The results of Table 1 Case $$3-5$$, where $${\mathrm{H}}^{+}$$ escape rates vary under nightside crustal field, southern summer dust storm, and extreme solar wind conditions, demonstrate that $${\mathrm{H}}^{+}$$ loss is influenced by the interplay between the exospheric hydrogen density and the factors controlling exospheric ionization. For Case 3, the slight increase of escape rate is mainly due to the decline of dayside ICB location when the crustal field locates on the nightside (Halekas et al., 2018), because more ions can be produced and picked up at low altitude. Case 4 examines the effect of elevated exospheric hydrogen density (Figure S1 in Supporting Information S1) resulting from dust storms at southern summer (Qin, 2021), while Case 5 simulates the interplanetary coronal mass ejection (ICME) condition with intensified solar wind fluxes (Y. Ma, Fang, et al., 2018; Romanelli et al., 2018). In Figures 4b and 4c, the $${\mathrm{H}}^{+}$$ escape fluxes increase significantly compared to Figure 2a, highlighting the sensitivity of $${\mathrm{H}}^{+}$$ escape to both internal and external drivers.

4 Discussions and Conclusions

In this study, the multifluid model of Najib et al. (2011) is used to quantify the contribution of the pickup escape channel to the total ion loss rate. The simulation results reveal that, unlike heavy ions such as $${\mathrm{O}}_{2}^{+}$$ and $${\mathrm{O}}^{+}$$, planetary $${\mathrm{H}}^{+}$$ does not predominantly escape through the magnetotail. Instead, over $$\mathrm{90}\mathrm{\%}$$ of the $${\mathrm{H}}^{+}$$ loss occurs within the magnetosheath and solar wind, with the pickup process identified as the dominant escape channel. Remarkably, the pickup $${\mathrm{H}}^{+}$$ is mainly sourced from solar wind charge exchange reactions with planetary hydrogen in the exosphere, in contrast to heavy ions (e.g., $${\mathrm{O}}_{2}^{+}$$, $${\mathrm{O}}^{+}$$), which forms a plume originating from the ICB, or ionopause (S. Xu et al., 2023), due to the finite Larmor radius effect.

Romanelli et al. (2018) simulated a significant increase in $${\mathrm{H}}^{+}$$ escape rates during an ICME event, while Zhao et al. (2023) observed a signature suggestive of enhanced $${\mathrm{H}}^{+}$$ loss during the passages of solar wind stream interaction regions (SIRs). These findings point to enhanced solar wind flux as the key driver during space weather events, as it directly amplifies the $${\mathrm{H}}^{+}$$ production rate by $${\mathrm{H}}_{\mathrm{S}\mathrm{W}}^{+}-\mathrm{H}$$ charge exchange reaction. Our simulations confirm this relationship: under the intense solar wind condition (Table 1 Case 5), the global $${\mathrm{H}}^{+}$$ escape rate surges. In addition, seasonal variations in hydrogen exospheric density (Chaffin et al., 2018; Chaufray et al., 2015; Clarke et al., 2017; Halekas, 2017; Rahmati et al., 2018), primarily driven by dust storms via atmospheric heating and water vapor transport (Chaffin et al., 2021; Heavens et al., 2018; Stone et al., 2020), can significantly modulate $${\mathrm{H}}^{+}$$ escape rate. Although our MHD calculations constrain the dust storm-driven increase in $${\mathrm{H}}^{+}$$ escape rates to less than an order of magnitude, the adopted exobase H density $$\left(\mathrm{2}\times \mathrm{1}{\mathrm{0}}^{\mathrm{5}}\ \mathrm{c}{\mathrm{m}}^{-\mathrm{3}}\right)$$ may underestimate dust storm effects. MAVEN observations reported perihelion-averaged exobase densities of $$(\mathrm{4}-\mathrm{9})\times \mathrm{1}{\mathrm{0}}^{\mathrm{5}}\ \mathrm{c}{\mathrm{m}}^{-\mathrm{3}}$$ (Mayyasi et al., 2023), suggesting potential underestimation of both exospheric hydrogen density and consequent $${\mathrm{H}}^{+}$$ escape enhancement during extreme dust events. We have also examined the influences of crustal magnetic field locations and solar activity on $${\mathrm{H}}^{+}$$ escape rates. Under nightside crustal field conditions, magnetotail $${\mathrm{H}}^{+}$$ fluxes increase slightly due to reduced dayside magnetic shielding, while pickup $${\mathrm{H}}^{+}$$ escape rates rise as the ICB lowers, enhancing ion exposure to solar wind acceleration.

In summary, our results highlight the critical role of pickup escape channel when estimating the global $${\mathrm{H}}^{+}$$ loss rate. The modeled $${\mathrm{H}}^{+}$$ escape rate under present-day conditions is approximately one order of magnitude lower than the Jeans escape rate of neutral hydrogen (Chaffin et al., 2018), while the influence of solar wind conditions on the $${\mathrm{H}}^{+}$$ escape rate implies a substantially higher hydrogen loss rate in the early history of Mars. Future refinements will focus on utilizing observational data, such as the energetic neutral atoms (ENA) detected by MINPA instrument (Kong et al., 2020) onboard Tianwen-1 spacecraft (Zou et al., 2021), to better characterize the pickup $${\mathrm{H}}^{+}$$ escape rate at Mars.