# Effects of Airless Bodies' Regolith Structures and of the Solar Wind's Properties on the Backscattered Energetic Neutral Atoms Flux

Sébastien Verkercke, Jean-Yves Chaufray, François Leblanc, E. M. Bringa, D.

Tramontina, L. Morrissey, A. Woodson

► To cite this version:

Sébastien Verkercke, Jean-Yves Chaufray, François Leblanc, E. M. Bringa, D. Tramontina, et al. "Effects of Airless Bodies' Regolith Structures and of the Solar Wind's Properties on the Backscattered Energetic Neutral Atoms Flux." The Planetary Science Journal, 2023, 4, pp.197. 10.3847/PSJ/acf6bd. insu-04267570

HAL Id: insu-04267570

https://insu.hal.science/insu-04267570

Submitted on 1 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Effects of Airless Bodies' Regolith Structures and of the Solar Wind's Properties on the Backscattered Energetic Neutral Atoms Flux

S. Verkercke¹, J.-Y. Chaufray¹, F. Leblanc², E. M. Bringa³, D. Tramontina³, L. Morrissey⁴, and A. Woodson⁵

¹LATMOS/CNRS, Université Versailles Saint Quentin, Guyancourt, France; sebastien.verkercke@latmos.ipsl.fr

²LATMOS/CNRS, Sorbonne Université, Paris, France

³CONICET and Facultad de Ingeniería—Universidad de Mendoza, Mendoza, 5500, Argentina

⁴Memorial University, Canada

⁵University of Virginia, Charlottesville, VA, USA

Received 2023 May 31; accepted 2023 August 29; published 2023 October 31

## Abstract

The surfaces of airless planetary bodies, such as the Moon or Mercury, are covered with regoliths, which interact with the solar wind. The solar protons can either be absorbed by the surface or neutralized and reflected as hydrogen energetic neutral atoms (ENAs). The ENA flux is thought to depend mostly on the structure of the upper regolith layer. By using a model combining a Monte Carlo approach to describe a solar proton's journey through the lunar surface with molecular dynamics to characterize its interactions with the regolith's grains, we highlight the surface roughness as a key parameter that influences the backscattered H ENA flux. By considering spherical silica grains, the lunar regolith's structure is described using the open-source code Large-scale Atomic/Molecular Massively Parallel Simulator (or LAMMPS), which allows a realistic description of grain-on-grain contacts. The roughness of the modeled regolith, characterized by the roughness ratio, is shown to be dictated by the surface energy and the grain-size distribution. This work shows that a rougher surface favors deeper penetration of the protons inside the regolith, which increases the number of collisions and thus decreases their reflected fraction. The angular distribution of the backscattered H ENAs is influenced by both the surface roughness and the solar zenith angle. We show that the angular distribution of the backscattered ENAs is anisotropic and is influenced by the regolith's structure, which is consistent with Chandrayaan-1 measurements. This work aims for a better understanding of the interactions ongoing at this interface and intends to look into the possibility of deducing information on the surface structure solely from ENA flux measurements. Highlighting the key structural parameters influencing the ENA backscattering will also help the development of models of surface-bounded exospheres.

Unified Astronomy Thesaurus concepts: Planetary surfaces (2113); Regolith (2294)

## 1. Introduction

Airless planetary bodies' surfaces, such as the Moon's or Mercury's, are composed of complex regolith structures that interact directly with the impinging solar wind (McComas et al. 2009). In the case of the Moon, this incoming flux of solar protons has been observed to be partially neutralized and backscattered as hydrogen energetic neutral atoms (ENAs). The first observation of these ENAs was reported in 2009 by McComas et al. (2009), who used Interstellar Boundary Explorer (IBEX) measurements to derive a reflection coefficient of ~0.1. Since then, other studies based on IBEX measurements have reported other values of this reflection coefficient, with Rodríguez Moreno et al. (2012) suggesting a value of 0.09 ± 0.05 or Funsten et al. (2013) reporting reflection coefficients ranging from 0.08 to 0.18 depending on whether the solar wind is fast or slow, respectively. Indeed, protons with larger energies can penetrate deeper into the intro-granular structure when colliding with a regolith grain, increasing their probability of absorption (Leblanc et al. 2023). Other missions, like Chandrayaan-1 or Chang'E 4, have also performed measurement of the backscattered H ENAs. Chandrayaan-1 observations suggest an ENA reflection coefficient ranging

between 0.1 and 0.2, without any correlation to the solar-wind velocity (Futaana et al. 2012; Vorburger et al. 2013). More recently, Zhang et al. (2020) reported a reflection coefficient of 0.32 based on the measurements of Chang'E 4. Such a large range of reflection coefficients underlines the lack of understanding of the lunar regolith's structure and its influence on the particles impacting it. Additionally, the multiple missions monitoring the H ENA reflection of the Moon each have very different observation conditions: IBEX is on an highly elliptic Earth orbit, Chandrayaan-1 on a polar lunar orbit, and the rover Chang'E 4 is monitoring directly at the surface (Futaana et al. 2012; Allegrini et al. 2013; Zhang et al. 2020). Further characterization of the angular distribution of the reflected ENAs is therefore necessary to accurately reconstruct the global reflection coefficient (Schaufelberger et al. 2011; Vorburger et al. 2013). The ENAs' angular distribution function has already been derived from Chandrayaan-1 measurements by Schaufelberger et al. (2011), in which they found that four distinct features emerged with increasing solar zenith angle (SZA): less azimuthal uniformity, a bigger ratio of sunward versus antisunward flux, a decrease in scattering angle with respect to the local, macroscopic surface normal (i.e., shallower scattering), and a decrease in amplitude. However, the origins of these features remain unknown. To the best of our knowledge, no modeling of this angular distribution has ever been done, which stresses the necessity of understanding how the surface's structure influences the ENA reflection through modeling.

<figure><img src="image_3.png" /></figure>

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Some models have already been developed in order to describe the solar wind-lunar regolith interaction. Hodges (2011) used the Lunar Exosphere Simulator toolkit to describe the solar ions' interaction with the regolith. In this model, the absorption coefficient of H is computed as being energy dependent, but the energy loss of particles was derived by fitting to match the measurements from Wieser et al. (2009). More recently, Szabo et al. (2022) modeled the interaction of solar protons with different regoliths characterized by different porosities using the ion-solid-interaction software SDTrimSP-3D. They concluded that with increasing porosity of the lunar regolith there was a reduction of the reflected fraction of H ENAs and that the structure of the regolith, especially its uppermost layer, has to be considered in order to understand the ions' impacts on airless bodies. While each of these studies focuses primarily on retrieving a global reflection coefficient for a fixed regolith structure, this current study varies the regolith structure to ascertain its influence on the angular distribution of backscattered H ENAs. Additionally, we also propose the roughness of the surface as one of the major structural parameters determining the ENA backscattering process rather than the porosity.

The open-source code Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS; Thompson et al. 2022) was used to describe the structure of the lunar regolith, which allows a simulation of thousands of micron-sized grains that were allowed to settle under the lunar gravity to simulate a realistic 3D packing of the regolith (Sarantos & Tsavachidis 2020; Szabo et al. 2022). This regolith model, described in Section 2.1, was then exposed to a simulated solar wind. This latter simulation, described in Section 2.2, uses Monte Carlo modeling to reconstruct the path of the protons in the regolith, while using molecular dynamics (MD) to describe the interactions between protons and regolith grains (Leblanc et al. 2023). This allows us to simulate the backscattered ENA flux of the lunar surface with different structural properties and under different solar-wind conditions. In doing so, we highlight the surface roughness as one of the main structural parameters which directly influence the interactions between the regolith and its exosphere (Section 3). We discuss these results in Section 4.

## 2. Models

In order to realistically simulate the interactions of a solar-wind proton with lunar regolith, we used two models: one describing the regolith structure based on a granular contact model, and the other describing the trajectory and the interactions of the protons within the regolith. Both of these models rely on the LAMMPS MD simulator, which, by solving Newton's equations for each particle and grain, allows an individual description of each proton-grain or grain-grain interaction.

### 2.1. Regolith Structure Model

In order to model realistically a regolith structure, a contact-stress model is needed to account for all the grain-on-grain interactions. To do so, a realistic description of grain-on-grain contacts was considered using a Johnson-Kendall-Roberts (JKR) contact model. This latter model incorporates the effects of adhesive contact forces between grains close together into the classic Hertz contact model (Johnson et al. 1971). The JKR model is an extension of the Hertz contact model, which describes the normal and shear forces based on theoretical

analysis of the deformation of smooth, elastic spheres in frictional contact (Mindlin & Deresiewicz 1953; Elata & Berryman 1996). In the case of the JKR model, an additional attractive force is introduced to the classic Hertz model, as it was observed experimentally. When considering spherical grains, this attractive force acts on the circular contact area, which can be described by a radius, a, and can be expressed as

$$F_{\text{adh}} = 2\pi a^2 \gamma, \quad (1)$$

with γ being the surface energy per unit area in joules per square meter. Using the LAMMPS MD simulator, and based on this contact model, we generated between 1200 and 800,000 spherical silica grains with radii ranging between 10 and 100 µm in a 1 mm × 1 mm simulation box that settled under a lunar gravity. SiO₂ being the most abundant component of the lunar regolith (McKay et al. 1991), the lunar surface was considered as solely composed of SiO₂. The grains were generated with a random spatial distribution with a small randomized initial velocity, which was incorporated to avoid the formation of ordered structures. The floor of the simulation box along the z-axis acted as a trap for the grains, as any grain falling onto it sticks there for the rest of the simulation. Enough grains were used so that the bottom layer of the simulated regolith would not affect the rest of the structure (Sarantos & Tsavachidis 2020). We used periodic boundaries on the x and y sides of our simulation box so that the structure would not be affected by the boundary's conditions (Figure 1). The number of simulated grains depends on the grain size used for the simulation, but the total mass generated was constant through all the simulations so that all generated structures could be compared. Grain-size distributions for the lunar regolith are known from the samples returned from the Apollo missions. The size of these grains ranges between a few microns to a few hundred microns, with a mean size between 40 and 80 µm, depending on the sample from which the distribution was computed (McKay et al. 1974; Tsuchiyama et al. 2022). Even though these distributions are known, it is computationally advantageous to build a given simulated regolith using grains of a single size. Thus, most of our simulations employ a regolith composed solely of 10, 20, 30, 50, or 100 µm grains. The novelty to our approach comes from using a realistic 3D packing with a realistic JKR contact model in order to simulate the lunar regolith and the grain-on-grain interactions.

### 2.2. Proton-Grain Interaction Model

As we are here interested in the influence of the structure of the upper regolith's layer and of the solar wind's characteristics' effects on the backscattered ENA flux, the journey of a solar proton through the regolith and its interactions with the grains need to be fully modeled. In light of this, the fate of a solar proton in the lunar surface can be described by combining two models, as in Leblanc et al. (2023): a calculation of the particle's trajectory in a layer of grains, and a calculation of the interaction between a particle and a grain (using MD).

To model the solar-wind interaction with the surface, we followed 100,000 test particles launched toward the regolith. The incident angle of the solar protons was defined by the chosen SZA, where SZA = 0° is taken to be the direction normal to the floor of the simulation box. This solar wind was
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Graphic representation of slices of regolith models created with LAMMPS while using a JKR contact model with $\gamma = 0.025 \text{ J m}^{-2}$ (a) and $\gamma = 1.5 \text{ J m}^{-2}$ (b). Both are formed of 50 µm grains. The brighter grains are closer to the simulation box sides, while the darker ones are occluded by surrounding grains. The simulation box is represented in black.</figcaption></figure>

simulated using the solar-wind conditions measured by Cluster between 00:00:00 and 10:00:00 UT on 2009 February 20 and between 14:00:00 and 22:00:00 UT on 2009 February 21. At this time, the solar-wind velocity, $V_{sw}$, ranged between 320 and 380 km s$^{-1}$, the solar-wind density, $n_{sw}$, ranged between 9 and $10 \text{ cm}^{-3}$, and the measured temperature, $T_{sw}$, varied between $0.6 \times 10^5$ and $2.8 \times 10^5$ K. These solar-wind conditions were chosen to be similar to the ones observed during Chandrayaan-1/CENA observations as reported by Futaana et al. (2012) and by Vorburger et al. (2013). At each time step, we evaluated if the particle encountered a grain, in which case various parameters derived from MD were used to determine the energy loss and direction of the particle if re-emitted (Leblanc et al. 2023). As the lowest energy range of CENA or IBEX is 10 eV (Wieser et al. 2009; Futaana et al. 2012), we removed from the simulation all test particles re-emitted with energies below 10 eV. Over the energy range considered here (10 eV to a few keV), the main source of energy loss comes from the particle-grain collisions, while the acceleration due to gravity of the re-emitted protons is negligible compared to their velocities, and that due to electromagnetic forces is negligible on the dayside where the surface potential is usually below a few tens of V (Halekas et al. 2011). Periodic boundaries are considered on the sides of the regolith so that a particle leaving on one side of the simulation box re-enters it through the opposite side. The particles passing through the bottom layer of the regolith are considered lost and are removed from the simulation. The particles which are not lost are followed until they intersect the top face of the simulation box. At this point, the given particle's weight, velocity and charge are saved. The time step of the simulation of ~0.3 ps was optimized to speed computation while accurately describing all the proton-grain

collisions as well as the position of each impact at the grain surface.

For each collision with a grain, the model computes the position of the impact at the grain surface and the incident impact angle with respect to the normal. It also determines the charge of the test particle, where we assumed a probability of 91% of the protons to be neutralized after each interaction with a grain (Hodges 2011). Using the prediction from MD, finally the likelihood that the particle is absorbed by the grain and, if not absorbed, the re-emission direction and energy of the particle after collision with the grain are evaluated. This description of the proton-grain encounters accounts for the interaction of the proton with the molecules composing the regolith grain. In this description, the interactions between the atoms are described by the reactive force field (or ReaxFF; Van Duin et al. 2001; Yu et al. 2016), which accounts for the screened Coulomb interaction and various bonding scenarios. Previous studies have validated this approach in the analysis of H interactions with $SiO_2$ (Morrissey et al. 2022). The incident particles are approximated as neutral hydrogen atoms since these particles move at high speeds, which makes the small lunar electromagnetic force negligible. These atoms are also considered to undergo many collisions with atoms in the first few nanometers of the sample. Since the energy of the incident protons is low, electronic stopping is neglected. Random impact points are selected on the sample for each bombardment, leading to different local environments around the impact and various outcomes, including the possibility of reflection. Some 500 events are simulated for each (energy, angle) combination. This ensures reflection coefficients with an error of less than 10% (Leblanc et al. 2023). By doing so, we inferred functions describing the absorption probability, the energy loss per collision, and the resulting velocity after each
collision for every incident angle and energy. These functions and the sample characteristics are described in more detail in Leblanc et al. (2023).

### 3. Backscattered Hydrogen Energetic Neutral Atoms

#### 3.1. Regolith Structures and Their Effects on the Backscattered H Fraction

The lunar regolith's structural properties are believed to influence the backscattered ENAs from the lunar surface (Wieser et al. 2009; Allegrini et al. 2013; Szabo et al. 2022). Most of the parameters associated with the composition of the grains, here SiO₂, which contribute to the JKR contact model (e.g., density, Young's modulus, Poisson ratio) are well constrained by different experimental measurements (Pabst & Gregorová 2013). However, the surface energy of silica grains is uncertain, as the actual measured values range over an interval of two orders of magnitude between 0.025 and 2.2 J m⁻² (Brunauer et al. 1956; Chokshi et al. 1993; Kimura et al. 2015). Since the surface energy of the grains influences directly the grains' interactions, the structure of the regolith should be different depending on the surface energy value that is used in the JKR model. In this work, we tested three different values for the silica grains' surface energy, γ: 0.025, 0.25, and 1.5 J m⁻², respectively. For each of these surface energies, we tested different grain sizes (10, 20, 30, 50, and 100 µm) in order to explore the influence of the grain size on the regolith's structure. We also modeled regolith structures using the grain-size distribution derived from the Apollo 17 sample 71241 described by McKay et al. (1974). In doing so, we obtained different regolith structures depending on the grain size and on the surface energy used. The global structure of a regolith is mainly described by the arrangement of its grains. This can be characterized by the porosity, defined as

$$P = 1 - V_{\text{reg}}/V_{\text{w}}, \qquad (2)$$

with $V_{\text{reg}}$ the volume occupied by the regolith and $V_{\text{w}}$ the total volume (grains + voids; Yan 2006; Szabo et al. 2022). This definition describes a mean porosity even if the porosity of the regolith is not uniform with depth, and that there is a general trend of increased porosity right at the surface. We noticed that at equal surface energy, a regolith composed of larger grains exhibits a smaller porosity than a regolith composed of smaller grains. The two main forces governing the grains' contact are the gravity and the adhesive force (Equation (1)), which depend on the cube and the square of the grains' radius, respectively (Johnson et al. 1971). Larger grains thus favor a tighter packing as they are more readily pulled downward by gravity into gaps and—depending on the Poisson ratio—compacted under the weight of those on top. Meanwhile, regoliths generated using a smaller grain surface energy exhibit a more compact packing (Figure 1(a)) than those with a larger value of γ because a large surface energy (Figure 1(b)) favors adherence between grains, resulting in a highly porous, fairy-castle-like structure for the regolith (Szabo et al. 2022). This mechanical adhesive force can be interpreted in the case of a contact between two spherical elastic grains, as a result of the equilibrium that obtains between the stored elastic energy and the lost surface energy (Johnson et al. 1971).

In order to identify the structure's influence on the reflection coefficient, we exposed our regolith models to the simulated solar wind described in Section 2.2. Simulations were performed at intervals of 10° in SZA from 0° to 80°. This allowed us to compute the reflection coefficient as a function of SZA, as was done in Szabo et al. (2022). Figure 2 presents the reflection coefficients of three different regoliths as a function of the SZA and compares these with the model of Szabo et al. (2022) and Chandrayaan-1 observations reported by Vorburger et al. (2013). One can see that, as presented by Szabo et al. (2022), the modeled reflection coefficient seems independent of the SZA until 60°, after which it starts to increase slightly with the SZA. As the SZA increases, some protons' incident angles with respect to the local normal of the impacted grain become blocked by the regolith structure, which causes the distribution in local incident angle to slightly shift to larger values. As shown by the flat surface curve (black solid line in Figure 1) representing the reflection probability of a proton after a collision with a grain as predicted by MD, larger incident angles favor reflections (Leblanc et al. 2023). Overall, this characteristic is consistent with the results of Szabo et al. (2022). To compare our results with Szabo et al. (2022), we characterize our regoliths by their mean porosity. For regolith with a porosity of 0.83 (orange), our results are consistent until 50° SZA, after which our results exhibit smaller dependency of the reflection coefficient with increasing SZA. Most likely, this effect is solely caused by the fact that SDTrimSP predicts a larger dependency of the reflection coefficient with increasing SZA than MD does, as shown in Figure 2. This difference might come from the use of a binary collision approximation for each proton's impact in SDTrimSP, whereas many-body collisions are accounted for in MD while including complex low-energy attractive interactions in the regolith. However, for regolith with a porosity of 0.55 (blue), our model predicts a lower reflection coefficient than Szabo et al. (2022) over all SZAs. This would suggest that, while ENA backscattering depends on the regolith porosity, this structural parameter has limited control of the backscattering of ENAs. Moreover, one can observe that two regoliths characterized by different porosities (red and orange solid lines with $P=0.71$ and $P=0.83$, respectively) can exhibit a similar reflection coefficient at all SZAs. This further suggests that the mean porosity is limited in its ability to describe the ENA backscattering.

Since the surface layer influences the reflection of ENAs (Schaufelberger et al. 2011; Szabo et al. 2022), we propose to use the roughness ratio (R) of the regolith as another key structural parameter governing this interaction. The roughness ratio, also known as the Wenzel roughness ratio, is computed as the textured surface area of the regolith divided by the projected surface area of the regolith (equivalent to the surface of the top of the simulation box in our case; Wenzel 1936; Hu & Sun 2016; Bai et al. 2021):

$$R = \text{Area}_{\text{Textured}}/\text{Area}_{\text{Projected}}. \qquad (3)$$

Rougher surfaces thus exhibit larger roughness ratio, and a flat surface exhibits a roughness ratio of 1. To compute the textured surface area of the regolith, we used a ray-tracing routine to map the regolith surface. The rays originate from a grid set above the regolith and propagate normally toward its surface, allowing one to evaluate the height at which each ray encountered a regolith grain. The height is computed with respect to the bottom of the simulation box, and the grid is set
<figure><img src="image_5.png" /></figure>

**Figure 2.** Reflection coefficient as a function of SZA. The blue, red, and orange solid curves correspond, respectively, to regoliths with porosities of 0.55 (Figure 1(a)), 0.71, and 0.85 (Figure 1(b)). The blue and red regoliths are composed of 50 and 10 µm grains, respectively, with a surface energy $\gamma = 0.025 \text{ J m}^{-2}$, while the orange regolith is composed of 50 µm grains with a surface energy of $1.5 \text{ J m}^{-2}$. The blue and orange dashed curves are the reflection coefficients retrieved by Szabo et al. (2022) for regoliths of porosities 0.55 and 0.83, respectively. The numbers in brackets correspond to the numbers in Figure 3. The black solid and dashed lines represent, respectively, the reflection coefficient dependency on SZA for a flat surface as predicted by MD (Leblanc et al. 2023) and SDTrimSP (Szabo et al. 2022). The observed reflection coefficients derived from CENA observations by Vorburger et al. (2013) are also plotted in gray. The error bars on the solid curves are based on the standard deviation of the reflection coefficient of protons launched from different random initial conditions.

<figure><img src="image_6.png" /></figure>

**Figure 3.** Reflection coefficient as a function of the regoliths' roughness ratio at (left) SZA 0° and (right) SZA 70°. The numbers in brackets correspond to the regoliths of Figure 2. The stars serve as indicators in Section 3.2. The error bars are based on the standard deviation of the reflection coefficient of protons launched from different random initial conditions. The different colors stand for the different types of grains composing the regoliths. The “hexagonal lattice” regoliths were created based on a hexagonal lattice to achieve compact packing.

fine enough so that the regolith's computed surface area converges. To ensure convergence, we always used more than 40,000 rays per mm². This method allows us to better characterize the surface of the regolith that is exposed directly to the solar protons.

Since our model shows that for any regoliths the average number of collisions by an incident solar-wind proton is lower than 3, it seems reasonable to assume that the roughness ratio, which characterizes the texture of the surface reached by the incident protons, is one of the most relevant structural parameters. Figure 3 presents the reflection coefficient as a function of the roughness ratio for all the modeled regoliths at

SZAs of 0° and 70°. We chose these SZAs to compare two extreme cases that could still be compared to observational data, as Chandrayaan-1 produced no relevant data beyond SZA 70° (Vorburger et al. 2013). The reflection coefficients presented by Vorburger et al. (2013) are represented by the gray shaded area for each SZA. One can observe that the reflection coefficient follows an exponential-like decrease with increasing roughness ratio. A larger roughness ratio implies that the regolith's surface hit by the solar protons on first impact extends deep into the regolith's structure. This favors deeper penetration of the protons in the regoliths, increasing the average number of proton-grain collisions and thus their absorption (Szabo et al. 2022). Our
results thus show that the absorption is greater for regoliths exhibiting larger roughness ratios. The three regoliths presented in Figure 2 (1–3), which are characterized by porosities of 0.55, 0.71, and 0.85, exhibit an average number of collisions by an incident solar-wind proton of 2.29, 2.62, and 2.6, respectively. The regoliths modeled using the grain-size distribution published by McKay et al. (1974, black dots) are well contained in the observation range.

Utilizing the roughness ratio also highlights that regoliths with similar roughness ratios can also exhibit slightly different reflection coefficients, with an increasing dispersion in reflection coefficients with increasing roughness ratio. For a regolith model to be realistic, its reflection coefficient should be in the range of Chandrayaan-1's measurements at any SZA. Figure 3 shows that the regoliths respecting this condition have roughness ratios between 5.4 and 19.2. Even if the roughness ratio tends to become larger with the porosity of a regolith, this constraint gives hardly any information with regard to the lunar regolith's porosity, since our range of regolith porosities spans from 0.5 to 0.86. However, Figure 3 shows that roughness ratios between 10 and 20 can exhibit similar reflection coefficients, which fails to achieve a more accurate correlation between the structure and the ENA reflections than when using the porosity. Nevertheless, using the roughness ratio as a regolith structural parameter is relevant in ENA backscattering analysis. More generally, in studies of surface-bounded exospheres, one could choose between porosity or roughness ratio depending on the specific nature of the study, as porosity is more relevant in the description of gaseous or thermal diffusion in the regolith (Yan 2006), whereas roughness ratio could be more relevant when studying sputtering (Brooks & Ruzic 1990).

## 3.2. Regolith Structures' Effects on the Energetic Neutral Atom Angular Distribution

Extensive observations of the H backscattering from the lunar surface by Chandrayaan-1 have been used to describe their angular distribution function (Schaufelberger et al. 2011). From these measurements, four distinct features emerged with increasing SZA: less azimuthal uniformity, a greater ratio of sunward versus antisunward flux, shallower scattering, and a decrease in amplitude. In order to verify that our model captures these features, we divided the hemisphere above the regolith into solid angles defined by intervals of 10° in both elevation and azimuth angles. The elevation angle, $\phi$, is defined as the angle between the backscattered ENA velocity vector and the normal to the surface of the regolith (i.e., the normal to top face of the simulation box). The azimuth, $\theta$, of each particle was also monitored and is computed as the angle between the backscattered ENA velocity projected on the surface and the projected Sun direction. By choosing two regoliths (marked by stars in Figure 2) exhibiting roughness ratios representative of the interval of realistic values established previously, and by tracking the solid angles in which each reflected particle is located, we simulated the ENA flux. In order to study the distribution, we define the sunward and antisunward directions as the particles having an azimuth $|\theta| < 90°$ or $> 90°$, respectively, and as introduced by Schaufelberger et al. (2011). Our model predicts that the elevation angular distribution of ENAs backscattered from a flat SiO$_2$ sample (based on MD) is a Gaussian distribution centered around 45° with an FWHM of 18°, for incident angles smaller or equal to 60°. For incident angles greater than 60°, the reflected angle follows a Rayleigh

distribution centered on the incident angle and with a FWHM following a relation like $\sigma(i) = -0.0163 i + 1.48$, with i the incident angle with respect to the surface normal (Leblanc et al. 2023). The model predicts a uniform azimuthal angular distribution for particles reflected on a flat SiO$_2$ sample (see the Appendix).

Additionally, for the regolith models validated previously by observations, we compute the average roughness, which is defined as (Myers 1962; Thomas 1981)

$$S_a(\text{SZA}) = \iint_A |Z(x, y)| \, dx dy, \text{ with } Z(x, y) = z(x, y) - \hat{z}, \qquad (4)$$

where x and y are the plane coordinates, $z(x, y)$ is the height of the surface at (x, y) with respect to the bottom of the simulation box, and $\bar{z}$ is the mean height of the surface. This value varies with the angle under which it is computed as, from the perspective of the incoming protons, the surface structure induces shadows (Sarantos & Tsavachidis 2020). The average roughness was computed with a ray-tracing routine using rays propagating at a certain SZA. This value tends to be maximal at SZA = 0° and minimal at SZA = 90°. The average roughness is reported in grain size units. If the regolith is composed of a grain-size distribution, the size of the largest grains is taken as the reference. We justify this choice of unit as, from the point of view of an approaching proton, only the geometrical arrangement of the grains matter, not their size.

Figure 4 presents the sunward–antisunward fluxes ratio at (a) SZA 20° and (b) SZA 80° as a function of $S_a$ for the different regoliths studied above. Figures 4(c)–(f) present the full angular distribution maps (elevation and azimuth) for the two extreme regolith structures modeled, marked with stars in Figure 3, and at two different SZAs. These regoliths are composed of 50 and 10 µm grains, respectively, and were modeled using $\gamma = 0.025$ and $1.5 \text{ J m}^{-2}$. By comparing the values of the sunward–antisunward fluxes ratio presented in Figures 4(a) and (b), one can clearly see that with increasing SZA, this ratio increases for any regolith structure, which is consistent with the results reported by Schaufelberger et al. (2011). The angular distribution is more uniform in azimuth at SZA = 20° than SZA = 80° for the two selected regoliths (Figures 4(c)–(f)). The peak of the angular distribution moves to larger elevation angles with increasing SZA, which causes the distribution to be shallower, as observed by Schaufelberger et al. (2011). Moreover, regoliths which exhibit a larger average roughness also present a larger sunward–antisunward fluxes ratio, as shown in Figures 4(a)–(b). This difference is increased at larger SZA, as clearly shown by Figures 4(c)–(f). This suggests that the increase of sunward–antisunward fluxes ratio with SZA observed by Chandrayaan-1 could be partially caused by the roughness of the regolith. The regoliths modeled using the grain-size distribution published by McKay et al. (1974, black dots) exhibit a sunward–antisunward fluxes ratio well within the range of the sunward–antisunward fluxes ratio from regoliths modeled using a single grain size. Additionally, using the largest grains of the distribution to describe the average roughness shows a similar trend of sunward–antisun-ward fluxes ratio with the average roughness than for the other regoliths, showing that the regolith structure mainly depends on these larger grains.
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Sunward vs. antisunward fluxes ratio as a function of the average roughness evaluated for multiple regolith structures at SZA 20° (a) and 80° (b). Normalized angular distribution maps for two extreme regolith structures at SZA 20° (c) and (d), and at SZA 80° (e) and (f). The bottom half of each map represents the sunward direction hemisphere and the upper half the antisunward hemisphere. The error bars are based on the standard deviation of the ratios of protons launched from different random initial conditions.</figcaption></figure>

However, our model always overestimates this ratio when compared to the observations of Chandrayaan-1 (Schaufelberger et al. 2011). At SZA 20° and 80°, respectively, Schaufelberger

et al. (2011) derived a sunward–antisunward fluxes ratio of 1.33 ± 0.06 and 1.38 ± 0.08, whereas we predict respective average values of 1.61 ± 0.08 and 2.72 ± 0.18. On the one hand,
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Left: field of view of the CENA instrument mapped on the lunar surface. The red and blue colored bands show the seven channels of the field of view (FOV), which is within the sphere of the lunar surface at the equator, and the spacecraft trajectory is shown in green. Right: projection of the FOV of each channel on the lunar surface at the subsolar point as reported by Bhardwaj et al. (2005).</figcaption></figure>

this discrepancy may be due to a geometrical factor neglected in the analysis of the CENA observations. Indeed, CENA's measurements, performed at lower elevation angles, imply that the lunar surface intercepted by CENA's field of view got larger (Bhardwaj et al. 2005), by up to ~12° of longitude taken at lunar equator for the channels at the extremities for observations performed at 100 km altitude (Figure 5). For these observations through the lateral channels (CH-0, -1, -5, -6), an interval of SZAs should be considered to derive accurately the surface reflection coefficient rather than a single SZA value. CENA measurements of the ENA backscattered flux at large elevation angles from the antisunward direction are increased by a contribution from the surfaces exposed to a lower SZA, whereas the observations made from the sunward direction are decreased by a contribution from surfaces exposed to a larger SZA. As the backscattered flux decreases as cos(SZA), these contributions would lead to the measurement of larger antisunward flux and smaller sunward flux, leading to a smaller sunward versus antisunward fluxes ratio at large SZA.

The distribution function derived by Schaufelberger et al. (2011) is based on a mean angular distribution taken over 15° intervals of SZA and does not take the latter geometrical effect into account. As shown by our model and by CENA observations (Schaufelberger et al. 2011), as the SZA increases the elevation angular distribution gets shallower, meaning that it is at these angles, where the observations are most affected by the geometry, that most of the ENA flux is emitted. Additionally, the number of observation at low elevation angles becomes smaller with increasing SZA (Schaufelberger et al. 2011; Vorburger et al. 2013), leading to larger uncertainty over the observed distribution function. At SZA 80°, CENA observations miss most of the $|\theta| < 30°$ and $150° < \theta < 210°$ intervals (Schaufelberger et al. 2011). Removing the flux passing through these regions in our model decreases the modeled sunward-antisunward fluxes ratio at SZA 80° by ~24%. These missing observations could also lead to an underestimation of the measured sunward-antisunward fluxes ratio at large SZA. On the other hand, the protons in our simulations were all approaching the regolith surface at the

same SZA angle. This does not include any velocity distribution due to a thermal component, which would tend to decrease the sunward-antisunward fluxes ratio by only a few percent. The sputtering contribution is also neglected in our model, which would mainly contribute to the antisunward flux as sputtering yield is preferentially forward oriented (Oen & Robinson 1978; Bay et al. 1980). However, it would also be negligible as the sputtering flux is small compared to the reflected ENA flux (Leblanc et al. 2023), and the sputtered particles have energies around 10 eV, which are barely detectable by CENA (Bhardwaj et al. 2005; Futaana et al. 2012). Moreover, our model only considers the influence of the micrometric structure of the regolith on the ENA backscattering, but large-scale topography and shadows should be included to compare our model with the observations.

From Figures 4(c)–(f), one can also immediately see that the angular distribution of the backscattered ENAs of the lunar surface depends on the SZA at which the protons are arriving. At SZA = 20°, for both regolith structures, the ENA angular distributions are mostly centered around the normal to the surface. As the SZA increases, the angular distribution peak moves to larger elevation angles, as observed by Schaufelberger et al. (2011). The decrease in amplitude of the backscattered ENA flux with increasing SZA cannot be seen in Figure 4 as the normalized fluxes are presented. However, since the backscattered flux depends on the solar-wind flux, which itself depends on the cosine of the SZA, the backscattered ENA flux is obviously smaller at large SZA. In our case, for any regolith structure the total flux reflected at SZA = 80° corresponds to only ~23.2% of the total reflected flux at SZA = 20°, meaning that the amplitude of the angular distribution does decrease with greater SZA, as observed by Schaufelberger et al. (2011).

## 4. Discussion and Conclusions

Atmosphereless bodies' surfaces, such as Mercury's or the Moon's, interact directly with the solar wind. Several missions have reported measurements of solar-wind protons being
backscattered from the lunar surface (Futaana et al. 2012; Funsten et al. 2013). Different reflection coefficients have been reported and attributed to structural properties of the regolith without much evidence of this influence or by attributing it solely to the porosity of the regolith (Wieser et al. 2009; Allegrini et al. 2013; Szabo et al. 2022). Using a Monte Carlo method combined with LAMMPS Molecular Dynamics Simulator for both the description of the solar proton-grain interactions (Leblanc et al. 2023) and of the grain-on-grain contacts, we tried to build a global physical description of the regolith's structure and of its interaction with the solar wind. This work allowed us to highlight the roughness ratio as one of the main structural parameters controlling the backscattering of the ENAs of the lunar surface. This parameter could also be important in the surface-bounded exosphere dynamics. The roughness ratio depends on different grain parameters such as size and surface energy (in the JKR contact model). We report that, at a similar surface energy and while considering regolith composed of same-sized grains, the grain size can influence the structure of the regolith and thus the ENA reflection. This means that different surface roughness can arise from the use of a granular contact model in the regolith structure description. Using a realistic grain-size distribution based on an Apollo 17 sample results in a regolith roughness ratio similar to the ones obtained when using same-sized grains. This suggests that these structures are shaped by the larger grains, whereas smaller grains only slightly fill the voids of these larger structures.

We report results which are consistent with previous models, such as the one presented by Szabo et al. (2022). However, we suggest that the roughness ratio also exhibits a relation with the ENA reflection coefficient, similarly to the porosity. A larger roughness ratio was shown to decrease globally the ENA backscattered fraction of the lunar regolith. The irregularities of the surface increase the average number of proton-grain collisions, thus reducing the global reflection coefficient. The range of roughness ratio which exhibits ENA reflection coefficients matching the reflection coefficient measured by CENA on Chandrayaan-1 spans from 5.4 to 19.2. The regolith structures based on the grain-size distribution published by McKay et al. (1974) present reflection coefficients which are consistent with the observations (Futaana et al. 2012; Vorburger et al. 2013). The lunar regolith structure and its roughness appear to be key parameters in the ENA backscattering. Further analysis of these influences could consider either the porosity or the surface roughness of these structures, depending on the specific needs of surface-bounded exosphere studies.

Through this work, we also analyzed the angular distribution of the backscattered H ENA and found a good agreement with the observations reported by Schaufelberger et al. (2011), notably regarding the anisotropy of the angular distribution and its dependency with SZA. This angular distribution is shown to greatly depend on the micrometric structure of the regolith and on the SZA. Regoliths presenting a larger average roughness exhibit a larger sunward-antisunward fluxes ratio. For any regoliths, we also observe that the angular distribution loses its azimuthal uniformity with increasing SZA and that the elevation angle distribution gets shallower. This global behavior of the angular distribution with SZA is consistent with observations from CENA (Schaufelberger et al. 2011). At SZAs of 20° and 80°, our model predicts mean sunward-antisunward fluxes ratios of 1.61 ± 0.08 and 2.72 ± 0.18, respectively, which is not consistent with the ratio derived from the observations by Schaufelberger et al. (2011). The lack of

agreement of our modeled ratio with the values of this ratio predicted by the measured angular distribution is not yet well understood, but could result from a measurement bias due to the geometry of the observations and/or to the lack of observations at large elevation angles for large SZA. Our model neglects the thermal contribution to the incident protons' velocities, the sputtering flux, and the large-scale topography effects, which could also partially explain these discrepancies. Finally, our regolith structures were considered stable at all times. However, on the dayside of the Moon, regolith cavities can accumulate positive charges and generate much higher electrostatic potential than 10 V (Nakazono & Miyake 2023). These larger potentials for the regolith's microcavities could lead to an instability in the topmost grains. This effect might remove fairy-castle-like structures, which would reduce the roughness of the regolith, and thus the sunward-antisunward fluxes ratio, if accounted for in our model. Comparing these modeled and observed ratios accurately could allow one to better constrain the lunar surface's structure. This work underlines the important role of the lunar regolith's structure, especially of the topmost regolith layer that constitutes the surface roughness, on the ENA backscattering. It also shows the influence of the roughness of the regolith on the ENA angular distribution, which needs to be accounted for when comparing observations performed under different spatial conditions. This study could be extended to Mercury by accounting for the influence of the planet's magnetic environment.

## Acknowledgments

S.V., F.L., and J.-Y.C. acknowledge the support of ANR, France of the TEMPETE project (grant No. ANR-17-CE31-0016) and the support of CNES, France for the BepiColombo mission. This research was supported by the International Space Science Institute (ISSI) in Bern, through ISSI International Team project #559 "Improving the Description of Exosphere Surface Interface." The authors also acknowledge the support of the IPSL CICLAD data center for providing access to their computing resources and data. D.T. and E.M.B. thank the support from PICTO-UM-2019-00048. Data may be obtained upon request from S.V. (email: sebastien.verkercke@latmos.ipsl.fr).

## Author Statement

S.V.: conceptualization, methodology and software of the LAMMPS regolith part. F.L and J.-Y.C: conceptualization, methodology and software of the MC part. D.T. and E.B.: conceptualization and software of the MD part. L.M. and A.W.: supervision and validation.

## Declaration of Competing Interests

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## Data Availability

Data will be made available on request.

## Appendix

Figure A1 shows the angular distribution maps for reflected H ENA of a flat surface as predicted by MD for SZA 0°(right) and 80°(left).
<figure><img src="image_9.png" /><figcaption>**Figure A1.** Angular distribution maps for reflected H ENA of a flat surface as predicted by MD for SZA 0° (right) and 80° (left).</figcaption></figure>

## ORCID iDs

S. Verkercke https://orcid.org/0000-0002-1966-6553

F. Leblanc https://orcid.org/0000-0002-5548-3519

D. Tramontina https://orcid.org/0000-0001-5356-6719

L. Morrissey https://orcid.org/0000-0001-7860-9957

## References

Allegrini, F., Dayeh, M. A., Desai, M. I., et al. 2013, P&SS, 85, 232

Bai, P., Leping, Z., & Xiaoze, D. 2021, ICHMT, 120, 105036

Bay, H. L., Bohdansky, J., Hofer, W. O., et al. 1980, ApPhy, 21, 327

Bhardwaj, A., Barabash, S., Futaana, Y., et al. 2005, JESS, 114, 749

Brooks, J. N., & Ruzic, D. N. 1990, JNuM, 176, 278

Brunauer, S., Kantro, D. L., & Weise, C. H. 1956, CaJCh, 34, 1483

Chokshi, A., Tielens, A. G. G. M., & Hollenbach, D. 1993, ApJ, 407, 806

Elata, D., & Berryman, J. G. 1996, Mech. Mater., 24, 229

Funsten, H. O., Allegrini, F., Bochsler, P. A., et al. 2013, JGRE, 118, 292

Futaana, Y., Barabash, S., Wieser, M., et al. 2012, JGRE, 117, E05005

Halekas, J. S., Saito, Y., Delory, G. T., et al. 2011, P&SS, 59, 1681

Hodges, R. R. 2011, GeoRL, 38, L06201

Hu, H., & Sun, Y. 2016, IJHMT, 101, 878

Johnson, K. L., Kendall, K., & Roberts, A. 1971, RSPSA, 324, 301

Kimura, H., Wada, K., Senshu, H., et al. 2015, ApJ, 812, 67

Leblanc, F., Deborde, R., Tramontina, D., et al. 2023, P&SS, 229, 105660

McComas, D. J., Allegrini, F., Bochsler, P., et al. 2009, GeoRL, 36, L12104

McKay, D. S., Fruland, R. M., & Heiken, G. H. 1974, LPSC, 5, 887

McKay, D. S., Heiken, G., Basu, A., et al. 1991, in Lunar Sourcebook, ed. G. Heiken et al. (Cambridge: Cambridge Univ. Press), 285

Mindlin, R. D., & Deresiewicz, H. 1953, JAM, 20, 327

Morrissey, L. S., Pratt, D., Farrell, W. M., et al. 2022, Icar, 379, 114979

Myers, N. O. 1962, Wear, 5, 182

Nakazono, J., & Miyake, Y. 2023, JGRE, 128, e2022JE007589

Oen, O. S., & Robinson, M. T. 1978, JNuM, 76, 370

Pabst, W., & Gregorová, E. 2013, Ceramics-Silikáty, 57, 167

Rodríguez Moreno, D. F., Saul, L., Wurz, P., et al. 2012, P&SS, 60, 297

Sarantos, M., & Tsavachidis, S. 2020, GeoRL, 47, e2020GL088930

Schaufelberger, A., Wurz, P., Barabash, S., et al. 2011, GeoRL, 38, L22202

Szabo, P. S., Poppe, A. R., Biber, H., et al. 2022, GeoRL, 49, e2022GL101232

Thomas, T. R. 1981, Precis. Eng., 3, 97

Thompson, A. P., Aktulga, H. M., Berger, R., et al. 2022, CoPhC, 271, 108171

Tsuchiyama, A., Sakurama, T., Nakano, T., et al. 2022, EP&SS, 74, 172

Van Duin, A. C., Dasgupta, S., Lorant, F., et al. 2001, JPCA, 105, 9396

Vorburger, A., Wurz, P., Barabash, S., et al. 2013, JGRA, 118, 3937

Wenzel, R. N. 1936, Ind. Eng. Chem., 28, 988

Wieser, M., Barabash, S., Futaana, Y., et al. 2009, P&SS, 57, 2132

Yan, N. 2006, Doctoral Dissertation, Université Pierre-et-Marie-Curie

Yu, Y., Wang, B., Wang, M., et al. 2016, JNCS, 443, 148

Zhang, A., Wieser, M., Wang, C., et al. 2020, P&SS, 189, 104970