Magnetic Evidence for an Extended Hydrogen Exosphere at Mercury

Remote observations by the Mariner 10 and MErcury Surface, Space ENvironment, GEophysics, and Ranging (MESSENGER) spacecraft have shown the existence of hydrogen in the exosphere of Mercury. However, to date the hydrogen number densities could only be estimated indirectly from exospheric models, based on the remotely observed Lyman‐α radiances for atomic H, and the detection threshold of the Mariner 10 occultation experiment for molecular H2. Here, we show the first on‐site determined altitude‐density profile of atomic H, derived from in situ magnetic field observations by MESSENGER. The results reveal an extended H exosphere with densities that are ∼1–2 orders of magnitude larger than previously predicted. Using an exospheric model that reproduces the H altitude‐density profile allows us to constrain the so far unknown H2 density at the surface, which is ∼2–3 orders of magnitude smaller than previously assumed. These findings demonstrate the importance of (a) in situ measurements supporting remote observations of Mercury's exosphere that will be realized in the near future by the BepiColombo mission and (b) that dissociation processes play a crucial role in Mercury's exosphere.

To date, however, the neutral H exosphere has been measured only on basis of remote optical detections. In this study, we determine for the first time the local H density profile from in situ (magnetic field) observations. We survey the magnetic field data of MESSENGER in the solar wind for so-called ion cyclotron waves (ICWs), which were specifically generated by freshly ionized H atoms. As the neutral H atoms from Mercury's exosphere become photoionized they start to gyrate around the background magnetic field (IMF) and get picked up by the solar wind. Since the velocity of the new born planetary protons (couple of km/s) is very different from the solar wind velocity (hundreds of km/s), the solar wind plasma becomes unstable to different plasma waves via resonant and nonresonant instabilities (Gary, 1991). As has already been shown at Mars and Venus, and since recently also at Mercury, most prominently ICWs are excited by this instability (Delva et al., 2008;Mazelle et al., 2004;Russell et al., 2006;Schmid et al., 2021). Figure 1 shows a schematic illustration of the ICW generation mechanism. ICWs are transverse electromagnetic waves near the proton cyclotron frequency. They propagate nearly parallel to the background magnetic field and are either left-or right-hand elliptically polarized: Theory suggests that left-hand polarized waves are produced by a perpendicular pick-up geometry (background magnetic field perpendicular to the plasma flow). Right-hand polarized waves are produced by a parallel pick-up geometry (background magnetic field parallel to the plasma flow) (Wu & Davidson, 1972;Wu et al., 1973). In the solar wind, mainly parallel pick-up takes place due to the small angle between the interplanetary magnetic field and the solar wind streaming direction, which is typically ∼30° at Mercury's orbit (James et al., 2017;Schmid et al., 2021). The right-hand polarized waves propagate in sunward direction with a phase speed on the order of the Alfvén velocity ( A = √ 0 , with B the magnetic field strength, μ 0 the permeability of free space, and ρ the (1) The atomic hydrogen (H, gray dots) get ionized by photons from the sun (purple line); (2) the newborn ions (H + , blue) start to gyrate around the interplanetary magnetic field (IMF, green lines) and get picked up by the solar wind flow; (3) In the solar wind frame of reference, the freshly picked up ions form a secondary distribution in velocity space that is highly unstable to the cyclotron wave instability and ion cyclotron waves (ICWs, orange lines) are excited.
3 of 14 mass density of the charged particles in the plasma). That speed is slower than the solar wind velocity. Hence, the waves are carried in antisunward direction over the spacecraft, thereby reversing the right-hand polarization by the anomalous Doppler shift (Mazelle & Neubauer, 1993). Consequently, in the spacecraft frame a left-hand polarization is observed. In that frame, the waves are shown to be always observed at the local ion gyrofrequency, since the newborn exospheric ions have a negligible velocity relative to the spacecraft. This immediately excludes confusion with ICWs generated at the bow shock by back-streaming solar wind protons because those waves will be observed at the spacecraft with frequencies very different from the local proton gyrofrequency due to their large velocity relative to the spacecraft (Delva et al., 2008(Delva et al., , 2011.
In Section 3, we identify ICWs that are specifically generated by the pick-up of freshly ionized planetary hydrogen. Thereto, we search for time intervals with large transverse magnetic field fluctuations, which are left-hand polarized in the spacecraft frame, and close to the local proton gyrofrequency. From the observed wave power of the identified ICWs, we are able to derive the local atomic H density, necessary to produce the observed waves and obtain the first in situ determined altitude density profile of hydrogen around Mercury in Section 4. Based on the determined H density profile, we also introduce an exospheric model that includes hydrogen photochemistry that can explain the origin of the discovered extended H exosphere and which allows us to constrain the so far unknown H density at the surface. In Section 5, we discuss the findings and make a quantitative statement on the reliability of the estimated hydrogen densities.

Ion Cyclotron Wave Identification Criteria
Our starting point is the recently published pick-up ICW event list that consists of 5455 events in the space environment around Mercury (Schmid et al., 2021). To identify the pick-up generated ICWs, they used the 20 Hz magnetic field observations by MESSENGER (Anderson et al., 2007;Solomon et al., 2007) between March 2011 and April 2015 and applied the following steps to a ∼100 s long sliding interval: 1. Within each interval the magnetic field data are transformed into a mean-field-aligned (MFA) coordinate system, where the parallel component, ̂‖ = 0∕| 0| , is given by the average magnetic field, B 0 = [B x,0 , B y,0 , B z,0 ], and the perpendicular components in this coordinate system are chosen to be ̂⟂ 2 =̂‖ × [0, 0, 1]| and ̂⟂ 1 =̂⟂2 ×̂‖ . 2. Each interval (2048 data points) is split into seven subintervals of ∼30 s (512 data points) with 50% overlap. The magnetic field data of each subinterval are Fourier transformed and the power spectral density matrix is evaluated. 3. The diagonal elements of the matrix give the in-phase power densities, parallel (P ‖ ) and perpendicular ( ⟂ = 1 2 ⋅ (P ⊥1 + P ⊥2 )) to the mean magnetic field (B 0 ). The off-diagonal elements of the matrix yield the out-of-phase cross powers, that is, the field rotation sense around the mean field. The complex off-diagonal elements of the spectral matrix are used to determine the ellipticity and the handedness of the observed wave (Arthur et al., 1976;Fowler et al., 1967;Means, 1972;Samson & Olson, 1980). Negative/positive signs refer to left/right-handed polarization of the wave in the spacecraft frame. 4. To evaluate the coherency between the input signals in a particular frequency range and to obtain how stable the components are in phase, the degree of polarization (DOP) of each subinterval is determined. Hundred percentage indicates a pure state wave and values less than 70% indicate noise (Samson & Olson, 1980).
The arithmetic means of the obtained power densities and ellipticities of the seven subintervals are calculated. A crucial condition for ICWs generated by local ion pick-up is that the observed wave frequency in the spacecraft frame is the same as in the plasma frame (no Doppler shift) and thus close to the local proton gyrofrequency (Delva et al., 2008). To provide a reliable identification, we calculate the proton gyrofrequency c,H + = qB 0 / (2πm) and error range Δ c,H + = qσ B /(2πm), with proton mass m, charge q, and the average and standard deviation of the magnetic field magnitude B 0 and σ B , for each ∼100 s time interval and apply the following criteria in the frequency range ΔF 4 of 14 • The power density per component is integrated in the frequency range ΔF to account for power maxima just below the calculated gyrofrequency. The ratio between the integrated perpendicular E ⊥ and parallel fluctuations E ‖ is evaluated and needs to be larger than 5: E ⊥ /E ‖ > 5. • Within ΔF the ellipticity ϵ should be smaller than −0.5, to ensure a left-handed polarization of the observed wave. • The DOP of all subintervals is required to be larger than 0.7 within ΔF, to maintain large coherency of the observed wave and that the signal-to-noise ratio is high. • The maximum of the perpendicular fluctuating field P ⊥ is within the limits of ΔF, to ensure that the observed wave is dominated by the ion cyclotron mode.  Figure 2a shows the magnetic field observation in MFA coordinates. The two perpendicular components (red and blue) are coherent and their fluctuations dominate over the parallel magnetic field variations (green). This can be also seen in Figure 2b where the perpendicular component of the power spectral density (red) prevails over the parallel component (green), indicating that the observed wave is rather transverse than compressional around the proton cyclotron frequency c,H + = 0.43 Hz (marked as solid black line). The area between the two dashed lines illustrates the integration frequency range ΔF, used to evaluate the power densities, ellipticity, and DOP. Figure 2c shows the hodogram in each plane of the MFA coordinate system for the time interval from Figure 2a. The observed wave is almost circularly left-hand polarized with an estimated ellipticity of ∼−0.70.
For this study, only those time intervals when MESSENGER was located in the solar wind are preselected. Utilizing an extended boundary data set (Philpott et al., 2020;Winslow et al., 2013), 3969 (of 5455) time intervals were allocated in the solar wind.

Injection Velocity Estimate
Due to the limitation of the plasma measurements of MESSENGER, we use a solar wind propagation model (Tao et al., 2005) provided by the Automated Multi-Dataset Analysis database that was successfully applied for Mercury (Schmid et al., 2021) to get an approximate estimate of the plasma density n SW and velocity V SW of the solar wind during the ICW observation period. As the injection velocity (V inj ), we use the aberrated solar velocity (V SW ), modified by the orbital motion of Mercury (V Mercury ) as provided by the Navigation and Ancillary Information Facility (NAIF; Acton, 1996) The injection velocity vector is also used to determine the solar zenith angle (SZA), which is used to select events dayside of the terminator (SZA < 90°), to maintain that the observed waves are freshly generated from local ion pick-up and to omit waves that might already diminish (Delva et al., 2009). From the 3969 preselected time intervals, 2247 pertain to observations dayside of the terminator. These 2247 ICWs finally constitute the data set for this study.

Results
As mentioned in Section 2, the ICWs should be generated locally through initial ionization of the neutral atomic H. To test this, we transform the observation locations into electromagnetic coordinates to examine possible asymmetries with respect to the convection electric field. ICWs propagate with speeds that are on the order of (or lower than) the local Alfvén velocity (V A ), which should be much lower than the injection speed (V inj ) (see e.g., Delva et al., 2009). Thus, we also evaluate the V A /V inj ratios during the ICWs detection. Since the initial velocity of the neutral H before ionization is negligible in the planetary frame, we assume that the injection velocity of the newborn exospheric ions directly corresponds to the aberrated solar wind velocity in the solar wind frame of reference.   Figure 3a that the observed ICWs are generated locally: (a) ICWs occur at large positive X MBE , far from the planet, suggesting that the ICWs should have propagated against the solar wind flow with a velocity faster than the solar wind speed, which is very unlikely. (b) ICWs are evenly distributed between ±Z MBE and since there is no known mechanism to move ions across the magnetic field against the electric field into the negative motional electric field region, the ICWs needed to be generated locally (Delva et al., 2008). Figure 3b depicts the normalized occurrence rate of the estimated V A /V inj ratios. The histogram confirms the second assumption: The Alfvén velocity is significantly smaller than the injection velocity of the newborn exospheric ions into the background (solar wind) plasma. From the results in Figure 3, we conclude that the underlying assumptions for reliable density estimations are fulfilled.
From the observed wave power of the ICWs, we derive the required pick-up proton densities in the same way as previous studies did at Venus (Delva et al., 2009): The total free energy, E free , which is required to excite cyclotron waves from a pick-up ion ring-beam distribution is approximately given by (Huddleston & Johnstone, 1992): Here, m i and H + are the mass and density of the pick-up ions, V A is the local Alfvén velocity, calculated from the modeled solar wind density and the in situ magnetic field measurements, V inj is the plasma injection velocity, and α(V inj , B 0 ) is the pitch angle between the plasma injection velocity and background magnetic field. Inverting Equation 1 for H + yields the pick-up ion density, under the assumption that the entire free energy of the ring-beam distribution is transferred to the wave and corresponds to the observed ICW energy, E free = ∫ ΔF P ⊥ df (Delva et al., 2009); H + can thus be understood as a lower limit for lower energy transfer rates. Hybrid simulations have shown that the ICWs grow rapidly until a quasi-steady level is reached after 60 − 100 ion gyrations (Cowee et al., 2012). It should be mentioned that simulations assume specific seed particle distributions, which develop in time in an isolated system. Under quasi-stationary conditions, the freshly produced (through photoionization) and lost (to the solar wind) ions are in equilibrium; however, the characteristic time (2π c,H + ⋅ t) of 60-100 ion gyrations until the ICWs fully develop can be directly applied to an open system as is the case here. Based on the (conservative) assumption that the full energy transfer from the ions to the waves takes 100 gyro periods, the pick-up ion density in this time should be balanced by the ion production rate, which can be estimated by multiplying the neutral H density n H with the photoionization rate ν: with H + being the estimated pick-up ion density from Equation 1 and H + the gyrofrequency. The photoionization rate varies significantly with the solar activity and the radial distance of Mercury to the Sun. Therefore, we modified the ionization rate according to the Flare Irradiance Spectral Model (FISM-P; Chamberlin et al., 2008) of Mercury, which reflects the solar activity, and the radial distance of Mercury to the Sun during the ICW observations: First, we normalized the FISM-P irradiance between [0,1], where 0 corresponds to 0.028 Wm −2 s −1 nm −1 and 1 to 0.1 Wm −2 s −1 nm −1 (which in turn corresponds to the minimum and maximum spectral irradiance index at 121.5 nm during solar cycle 24). Based on the normalized FISM-P irradiance index, we interpolate the corresponding photoionization rate at Earth's orbit between minimum 7.26 × 10 −8 s −1 (quiet Sun = 0) and maximum 1.72 × 10 −7 s −1 (active Sun = 1) (Huebner & Mukherjee, 2015). In the last step, we rescale this ionization rate from Earth's orbit (1 AU) to the square of the distance of Mercury during the ICW observation (0.31-0.47 AU). Figure 4 shows the radial dependence of the estimated H number density. The radial distance is given from the planet's center in Mercury radii (R M = 2,440 km) (top axis) and from the surface (bottom axis). The blue boxes represent the median number densities of H and the upper/lower quartiles. The black error bars are the maximum and minimum densities within the 0.5 R M bins. The obtained medians decrease from ∼100 cm −3 to ∼10 cm −3 between 2 R M and 7.5 R M . The relationship is approximately logarithmic between number density and radial distance. Although the number densities from individual ICW events can vary considerably within the bins (see error bars), the differences between the upper and lower quartiles are small. Interestingly, the obtained H density profile follows the trend of exospheric Monte Carlo models, which have successfully been applied in previous studies (Pfleger et al., 2015;Wurz & Lammer, 2003;Wurz et al., 2019) but are 1-2 orders of magnitude larger and lie in-between the atomic (H) and molecular (H 2 ) profiles of these models (see e. g., Figure 3 in Wurz et al. (2019)). This might indicate that dissociation processes, which have not been 7 of 14 included in the models, may play an important role in Mercury's exospheric composition. The dashed green line in Figure 4 shows the result of the H 2 profile that is modeled with a kinetic Monte Carlo model (see Appendix A) for the commonly assumed surface number density of 1.4 × 10 7 cm −3 , given by the upper limits derived from the detection threshold of the Mariner 10 occultation experiment (Hunten et al., 1988; R. M. Killen & Ip, 1999). From this H 2 profile, we model the corresponding ionization, dissociation, and recombination products that result in H + 2 , H + and H atoms (details of the model are given in Appendix A). The solid green line shows the resulting H density profile originating from the dissociated H 2 molecules, which are most likely the major source for H atoms at altitudes that are >1.5 R M . Although the modeled density profile reproduces the trend of the altitude H density profile obtained from the in situ ICW observations, the modeled densities are still higher, indicating that the H 2 surface density should be lower than ∼1.4 × 10 7 cm −3 . The best results that reproduce the inferred H densities from the ICW observations are obtained for a H 2 surface density of ∼8 × 10 4 cm −3 (dashed red line). The red lines correspond to the H 2 and H profiles based on this lower H 2 surface density which is ∼2-3 orders of magnitude lower than the previously assumed surface density of ≤1.4 × 10 7 cm −3 .

Conclusions and Discussion
So far, only remote measurements of H Lyman-α emissions have been used to evaluate the H exosphere at Mercury. In this study, we present for the first time the number density profile of H in Mercury's exosphere, based on in situ magnetic field measurements of ICWs by the MESSENGER spacecraft. For this study, we assume that the observed ICW wave energy exactly corresponds to the free energy in the ring-beam distribution of the pick-up ions (Huddleston & Johnstone, 1992). Simulations, however, show that the energy transfer from the ions to the waves might be 3-4 times smaller because the energy is distributed between wave growth and ion heating (Cowee et al., 2007). Consequently, the number densities obtained in this study might be underestimated. The variability of the energy transfer efficiency might also result in the broad distribution of estimated number densities at similar radial distances from Mercury (see black error bars in Figure 4). However, these simulations only consider a rather perpendicular pick-up geometry with only small parallel relative drift velocities. Due to the small cone angle between the interplanetary magnetic field and solar wind velocity, the ICWs used in this study are rather generated under quasi-parallel pick-up configurations (Schmid et al., 2021). Theoretical studies on the growth rates for parallel and perpendicular pick-up geometries suggest that perpendicular picked up ions produce ICWs with low growth rates and that parallel picked up ions generate ICWs with large growth rates (Wu & Davidson, 1972;Wu et al., 1973). Therefore, the observed ICWs are most likely fully developed in less than the assumed 100 ion gyro periods. Together with the assumption that all of the free energy from the pick-up ions is transferred to the ICW, the derived H density may thus be understood as a lower limit but is still larger than predicted by exospheric models, which assume a thermal H atom population with subsolar H atom surface densities between 23 cm −3 and 230 cm −3 (Hunten et al., 1988;R. M. Killen & Ip, 1999;Wurz & Lammer, 2003). It should be noted that it is likely that the H profile is not strictly isotropic in latitude and local time and may have a variability along the year. To achieve the extension of our study, a comparison with a subset of measurement at, for example, only the subsolar region would be more correct, but the statistical data amount is not sufficient to reach this goal at the present stage. We expect that additional measurements by the two spacecraft BepiColombo mission further increase the number of records in the data set upon its arrival in 2025 (Benkhoff et al., 2010). With the Mercury Magnetospheric Orbiter (Mio), one of the spacecraft, we may obtain a higher rate of statistic around the planet due to the short period of Mio (as compared to MESSENGER), which will allow us to concentrate on specific local time effects in the Hermean exosphere.
In order to make a quantitative statement on the reliability of the estimated H densities in this study, we transform them to Lyman-α radiances and compare them to the previously detected radiances of MESSENGER and Mariner 10. The radiance R is given in Rayleigh by 4π R = g ⋅ N/10 6 , where N is the integrated column density along the LOS in atoms cm −2 and g the photon scattering coefficient (referred to as g-value). The g-value for H at Mercury is 5.3 × 10 −3 photons atoms −1 s −1 (Hunten et al., 1988). When the exosphere of Mercury is viewed externally along the LOS of the spacecraft that is tangent to a spherical shell of radius r (i.e., radial vector form the center of the planet), the column density N at this position can be expressed in terms of the local scale height h and the local density of scattering H atoms n H with = H ⋅ √ 2 ⋅ ℎ ⋅ (Chamberlain & Hunten, 1987, Chapter 6). h can be directly obtained from Figure 4 by e-folding of the H density n H at position r. Based on this simple approach, we are able to transform the estimated H densities n H to their Lyman-α radiances at various altitudes and thus obtain the expected airglow around Mercury. The blue line in Figure 5 shows the Lyman-α radiances obtained from the median (blue dots), upper, and lower quartile (blue error bar) of the observed H densities in Figure 4. The gray and black dots indicate the Lyman-α observations of MESSENGER and Mariner 10 during their flybys (Ishak, 2019;Vervack et al., 2011). Although the H densities derived from the ICW observations are more than one order of magnitude larger than estimated by previous models, they are still in good agreement with the upper limits of the Lyman-α radiances measured by MESSENGER and Mariner 10, suggesting that the derived densities are in the correct order of magnitude and that dissociation processes play an important role in Mercury's exosphere.
Several exospheric models have been developed to explain the complex, and interrelated, source and loss processes of Mercury's exosphere (Jones et al., 2020;R. Killen et al., 2007;Mura et al., 2007;Wurz & Lammer, 2003;Wurz et al., 2010). Previous studies suggested that H 2 molecules originate from H + ions from the solar wind or the exosphere (through ionization) that are backscattered to surface (e.g., Anderson et al., 2011). Part of these ions diffuse at the surface, where they recombine and degas as H 2 O molecules (Hunten et al., 1988;Potter, 1995;Tucker et al., 2019). Another part will produce H2O from reactions with OH groups within the H-saturated regolith. The grains are saturated with solar wind H, thus the H 2 recombination might happen in the grains as well, and the H 2 will diffuse out, since it is chemically less bound to the mineral than the H atom (Jones et al., 2020). Thus, H 2 formation competes with the production of OH and H 2 O in the regolith. In experiments, it is found that on the Moon only ∼2% of the implanted H + is released as H 2 (Crandall et al., 2019). However, on the Moon the OH is only observed at higher latitudes, where the temperature is low enough. On Mercury's dayside the temperature 9 of 14 is too high that this process becomes important but close to the terminator the temperature is cooler where this might happen. Although the solar wind proton flux is higher at Mercury's orbit compared to the Moon at 1 AU, the planet's magnetosphere protects large areas from solar wind precipitation, which may indicate that the H 2 formation in the surface is also a not very efficient process (Hunten et al., 1988). Therefore, these studies suggested that photolysis of exospheric OH and H 2 O molecules, stemming from the bombardment of micrometeoroids (Broadfoot et al., 1976;Hunten et al., 1988;R. M. Killen & Ip, 1999;Killen et al., 1997), may be a much more efficient source of H 2 near the surface than chemical reactions in the surface (Jones et al., 2020). Note that these micrometeoroids are smaller and more common in comparison to the large meteoroids, which might yield transient enhancements of the exosphere at high altitudes (Jasinski et al., 2020;Mangano et al., 2007;Pokorný et al., 2018).
One can estimate the H 2 O density at Mercury's surface by using the estimated flux of 1 × 10 8 cm −2 s −1 of H 2 O from the vaporization of micrometeoroids (Killen et al., 1997). If we scale the H 2 O photodissociation time for average solar activity at 1 AU (Huebner & Mukherjee, 2015) to Mercury's average orbital distance of 0.38 AU, we obtain ∼10 4 s. The average micrometeoroid-related column density of H 2 O is then 1 × 10 12 cm −2 . By using an average temperature of 4000 K for the ejecta gas/water vapor (Wurz & Lammer, 2003), one obtains a scale height for micrometeoroid-related H 2 O vapor of ∼500 km. This would yield a surface number density of 2 × 10 4 cm −3 . This value is lower than the upper limit of possible H 2 O surface density of 1.5 × 10 7 cm −3 , estimated by earlier studies (Hunten et al., 1988). However, there will also be a thermal H 2 O population on Mercury, that is produced by surface reactions and evaporation from ice deposits on the nightside or from the planet's interior (Deutsch et al., 2019;Hunten et al., 1988;Killen et al., 1997Killen et al., , 1997Moses et al., 1999). H 2 O molecules will be dissociated in H and OH and ∼13% will yield H 2 and O atoms (Gombosi et al., 1986;Hunten et al., 1988). Photochemical reactions most likely enhance the lifetime of H 2 O molecules near Mercury's surface by ∼8 times (Hunten et al., 1988).  (Ishak, 2019;Vervack et al., 2011). 10 of 14 Moreover, a fraction of OH will be adsorbed at the surface where it also reacts with H so that H 2 O can be recycled near the surface too. Previous studies showed that the number density of thermally released H 2 O molecules quickly drops off to negligible values above 1.3 R M (Wurz & Lammer, 2003). The hotter micrometeoroid-related H 2 O population reaches 2 R M with a number density of a few cm −3 . Because of the decreasing availability of OH and O molecules at these distances, H atoms that will be produced via dissociation of H 2 molecules will not be efficiently removed by photochemical processes with these molecules. Although, H atoms (that are produced from H 2 O molecules that originate from vaporized micrometeoroids) may contribute to the atomic H number densities, which we have inferred from the ICW observations upstream of the bow shock in the solar wind, we expect that the main source of our derived H number densities between 2-8 R M is the dissociation of H 2 molecules.
The simulation output that best reproduces the observationally derived H density at distances above >1.5 R M yields a H 2 surface density of ∼8 × 10 4 cm −3 . Such an H 2 surface density yields a modeled atomic H density of 100-10 cm −3 between 2-8 R M (solid red line in Figure 4) with an escape rate of dissociated H atoms of ∼6 × 10 25 s −1 . From our analysis, we can therefore constrain the so far unknown and overestimated H 2 surface number density to ∼8 × 10 4 cm −3 . This value allows us to study in the future the details of the solar wind implantation into Mercury's regolith that leads to H, H 2 , OH, and H 2 O production and exospheric release as well as H 2 O photochemistry in the exosphere. It will give us the opportunity to investigate and separate the H 2 O sources and sinks on the innermost planet of the solar system. We expect that future measurements by the BepiColombo mission, in particular by the SERENA package (Orsini et al., 2021a(Orsini et al., , 2021b onboard the Mercury Planetary Orbiter spacecraft will help refine our knowledge about Mercury's exosphere. While the SERENA-STROFIO instrument will measure the H 2 and OH components and probably also the H component, the SERENA-PHEBUS observations of Lyman-α and other spectral lines can give the column densities of H, H 2 , and OH. Finally, the same analysis presented in this study can be repeated with the magnetic field measurements of MPO-MAG (Glassmeier et al., 2010;Heyner et al., 2021) and Mio-MGF (Baumjohann et al., 2019(Baumjohann et al., , 2020, which can be compared to the observations of SERNA-PICAM and Mio-MSA (Saito et al., 2010) that will respectively detect the low energy H + components nearby and far away from Mercury. SCHMID ET AL. The H, H + , and H + 2 profiles were obtained after integrating the system of ordinary equations for ionization, dissociation, and recombination, which can be written as follows: Although, atom recombination has only a small influence in the results, in accordance to previous works (Erkaev et al., 2016(Erkaev et al., , 2017Tian et al., 2008;Yelle, 2004) the process is included and treated as a three body reaction H + H + M → H 2 + M. Equation A2 were solved in normalized quantities and the distance r from the planetary center is normalized to R p . L is normalized to ( , where T 0 = √ kB 0∕ H , N 0 is the H 2 number density at the planetary surface. The total density ρ is normalized to (N 0 H 2 ). The normalized rates with the normalization factor in brackets can be written as follows: ions, α 1 , α 2 , and α 3 are the coefficients of recombination (H + + e − → H), , and (H + H → H 2 ), respectively, and ν i1 , ν i2 , and ν d are rates of photoionization of H, H 2 , and dissociation, respectively. Note that the ionization and dissociation rates of H and H 2 correspond to the average values over all 2247 ICW events. The H 2 photoionization and dissociation rates were calculated for each event separately according to the Flare Irradiance Spectral Model (FISM-P, (Chamberlin et al., 2008)) and radial distance from Mercury to the Sun during the event observation (see also section: Hydrogen density estimate from ICW observation), and that the recombination coefficients are taken from previous publications (Johnstone et al., 2015;Yelle, 2004).
From the studied hydrogen-bearing species, molecular hydrogen is considered to be the major constituent and thus the total mass density is assumed to be approximately equal to that of the H 2 density.
We further determine the density ρ-function by interpolating the H 2 profile from Monte Carlo simulation H 2 results. Using the determined loss rate L and density ρ, we integrate the system of Equation A2 with respect to quantities H + , H + 2 , and s, which finally yields the radial distributions of the dissociated atomic H and ionized particles.

Data Availability Statement
The magnetic field (MAG) data from the MESSENGER spacecraft are publicly available at the NASA Planetary Data System (PDS) and can be retrieved on their website (Korth, 2021): https://pds-ppi.igpp.ucla.edu/ search/view/?f=yes&id=pds://PPI/mess-mag-calibrated/data/mso. The numerical results of the model described in Appendix A and presented in Figure 4 are accessible via the Open Science Framework (OSF) repository (Schmid, 2022): https://osf.io/qtxh7/?view_only=4195a6a20d984d78b919155e37c4b736. The solar wind density and velocity data were obtained from the Automated Multi-Dataset Analysis (AMDA; Génot et al., 2021) database. All data are open-access and can be downloaded on their website via the Workspace Explorer under Solar Wind Propagation Models/Mercury/Tao Model/SW/Input OMNI: http://amda.cdpp.eu/. The orbital motion of Mercury were retrieved from the Navigation and Ancillary Information Facility (NAIF; Acton, 1996), publicly accessible on the NASA Jet Propulsion Laboratory (JPL) webpage: https://wgc.jpl.nasa.gov:8443/webgeocal-12 of 14 c/#StateVector. The solar spectral irradiance at the orbit of Mercury that was used to determine the solar activity during the event observations was obtained from the Flare Irradiance Spectral Model for Mercury (FISM-P), provided by the LASP Interactive Solar iRradiance Datacenter (Chamberlin et al., 2008), publicly accessible on their webpage: https://lasp.colorado.edu/lisird/data/fism_p_ssi_mercury/.