Modeling the Lunar Radiation Environment: A Comparison Among FLUKA, Geant4, HETC‐HEDS, MCNP6, and PHITS

Radiation transport codes have been an increasingly important tool for studying the space radiation environment, which includes high‐energy and high‐nuclear‐charge particles. The unique advantage of transport models lies in covering a wider range of particles, energies, and angles than would be attainable in a laboratory or measurable by an instrument. However, since there are several transport codes developed by different teams that have contributed heavily to the literature, differences are expected between such codes. In this work, we use five such radiation transport codes (FLUctuating KAscade, GEometry ANd Tracking, High‐Energy Transport Code‐Human Exploration and Development in Space, Monte Carlo N‐Particle, and Particle and Heavy Ion Transport code System) to study the radiation environment near the Moon, specifically the lunar “albedo” radiation, which is the radiation emitted by the lunar surface through interactions with incident galactic cosmic rays and solar energetic particles. The primary goal of this paper is to provide a general characterization of the lunar albedo radiation and to find the areas where the modeled transport codes agree and disagree by using almost identical input parameters. The results of this work show overall good agreement between the codes. However, some areas of discrepancies exist that are reported herein. Thus, this paper equips the space weather and radiation biology communities with a comparison between popular radiation transport models applied to lunar albedo radiation. The overall agreement and, in some cases, discrepancies between these transport codes provide fundamental insight necessary for assessing code reliability and identifying where further study and improvements are needed to advance our understanding of lunar albedo radiation.

Various missions to the Moon have provided measurements of albedo radiation, mainly gamma rays, neutrons, and protons. Radiation transport codes have also been a valuable resource for studying lunar albedo particles, especially with their unique advantage of providing spectra of a wider range of particles, energies, and angles than would be attainable in a laboratory or measurable by an instrument. Different studies using Monte Carlo (MC) transport codes have focused on selected albedo species and energy ranges, with the codes including MCNPX Prettyman et al., 2006), Geant4 (Looper et al., 2013;Yamashita et al., 2008), FLUKA (De Angelis et al., 2007), PHITS (Ota et al., 2011), and MCNP6 (Zaman, 2020;Zaman et al., 2021). This paper provides a generalized characterization of lunar albedo species by comparing results from various MC transport codes created using almost identical input parameters. Transport code comparisons help in improving radiation modeling by uncovering the areas where the codes might agree or disagree, which provides a future reference for any work studying the lunar albedo radiation environment or using the modeled radiation codes for space radiation protection purposes. In addition, the energy spectra are important for calculating the biological risk from the radiation near the lunar surface.
In a recent paper, , which was a follow-up study to Matthiä et al. (2016), examined the Martian radiation environment by providing comparisons between measurements from the Radiation Assessment Detector of the Mars Science Laboratory and calculations using five different radiation transport codes: HETC-HEDS (de Wet & Townsend, 2017), PHITS , Geant4 (Matthiä & Berger, 2017), MCNP6 (Ratliff et al., 2017), and HZETRN (Slaba & Stoffle, 2017). Results presented in  indicated a factor of two to three agreement among codes. In this work, we compare the spectra of albedo particles generated using FLUKA, Geant4, HETC-HEDS, MCNP6, and PHITS for the Moon, which, unlike Mars, does not have an atmosphere. Except for the geometry and physics models, all the input parameters were identical between all the codes in an effort to produce the most comparable results. These identical parameters include: the incident energy and angular distributions of the GCR source term, the elemental composition and density of the material representing the lunar regolith, and the output tallies.
Section 2 briefly discusses each MC transport code utilized in this work and provides a description of the input parameters used in each code. Section 3 is divided into two parts: the energy spectra of all the modeled species, and the angular distributions of selected particles. The conclusions and future work are presented in Section 4. Last, Appendix A provides a summary of the transport and nuclear physics models used in each code and Appendix B presents additional plots of the results obtained in this work.

Transport Models
Each MC code used in this work is described below. Detailed diagrams of the transport and physics models used in each code are available in Appendix A, apart from HETC-HEDS which is described in this section. It is important to note that the diagrams in Appendix A inevitably provide only approximations and do not necessarily cover the all the processes involved.

FLUKA
FLUKA (FLUctuating KAscade) MC transport code, developed by the European Organization for Nuclear Research and the Italian Institute for Nuclear Physics, is "a general-purpose Monte Carlo code for the interaction and transport of hadrons, heavy ions, and electromagnetic particles from few keV (or thermal energies for neutrons) to Cosmic Ray energies in arbitrary materials" (Battistoni et al., 2015). It endeavors to perform its calculations using complete physical models, duplicating interactions at the microscopic level. The FLUKA calculations are calculated using FLUKA version 4-2.2 and the Flair graphical interface (https://fluka.cern; Ahdida et al., 2022;Andersen et al., 2004;Battistoni et al., 2015;Bohlen et al., 2014;Fedynitch, 2015;Roesler et al., 2001;Vlachoudis, 2009).

Geant4
Geant4 (GEometry ANd Tracking) is an open-source MC radiation transport toolkit, which enables the user to model the transport of arbitrary energetic particles through a 3D material geometry. It tracks primary particles, and any secondary particles that are generated, individually until they stop, interact, decay, or escape the geometry, giving the user access to all properties for logging purposes at any point along their trajectories. The version used in this work is Version 10.6 of the code with patch level 2 (geant4.10.06. p02) on Linux clusters, using the reference physics list Shielding_EMZ. This physics list is based on a Bertini style cascade below several GeV and the Fritiof model above, with improved models for neutron interactions and for collisions between two heavy nuclei, both of which are important for the present work (Allison et al., 2016).

HETC-HEDS
HETC-HEDS (High-Energy Transport Code-Human Exploration and Development in Space) is a three-dimensional MC particle transport code modified by Townsend et al. (2005) to include high-nuclear-charge and high-energy transport. It is the only code in this work that was developed primarily for space radiation studies. "Each particle in the cascade (generated by HETC-HEDS) is followed until it eventually disappears by escaping from the geometric boundaries of the system, undergoes nuclear collision or absorption, comes to rest due to energy losses from ionization and excitation of atomic electrons, or, in the case of pions and muons, decays. Neutrons produced below a given cut-off, usually 20 MeV, and photons produced in the cascade or from de-excitation of gammas are not transported. Instead, the information is stored for transport by (other) codes" . In this work, the lower energy limit for neutrons was set to 1 MeV. For heavy particle fragmentation, HETC-HEDS uses a modified version of the NUCFRG2 model, which was developed at NASA Langley (Wilson et al., 1994). The modification made for use in HETC-HEDS was to replace the light ion yields with a stochastic model that accounted for particle energies and directions of travel after the collision . For Coulomb collisions, the code utilizes the continuous slowing down approximation (Turner, 2004).

MCNP6
The Monte Carlo N-Particle (MCNP) code, developed by Los Alamos National Laboratory, is "a general-purpose, continuous-energy, generalized-geometry, time-dependent, Monte Carlo radiation-transport code designed to track many particle types over broad ranges of energies" (Werner, 2017). In this work, we use MCNP version 6.2, the current version of MCNP6, which is the merger of the previous codes: MCNPX and MCNP5. The recommended default physics model options are chosen for all MCNP simulations. These include the Cascade-Exciton Model (CEM03.03) (Mashnik & Sierk, 2012), ISABEL (Yariv & Fraenkel, 1979) for light ions with energies <940 MeV/n, and the Los Alamos version of the Quark Gluon String Model (LAQGSM03.03) (Mashnik, 2011) for heavy ions, light ions with energies >940 MeV/n, pions with energies >2,500 MeV/n, and protons and neutrons with energies above 3,500 MeV/n (Werner, 2017;Werner et al., 2018).

PHITS
Particle and Heavy Ion Transport code System (PHITS) is a MC code developed by Japan Atomic Energy Agency, Research Organization for Information Science and Technology, High Energy Accelerator Research Organization, and several other institutions. The code calculates the reaction and transport of numerous particles (Sato et al., 2018). In this work, we used PHITS version 3.22 with the Intra-Nuclear Cascade of Liége (INCL4.6) and JAM for the nuclear reaction models and the GEM model for the evaporation model (Boudard et al., 2013;Nara et al., 1999). The continuous slowing down approximation is applied using the ATIMA code (Scheidenberger & Geissel, 1998). Electrons, positrons, and photons were transported using the EGS5 algorithm (Hirayama et al., 2005). For neutron transport below 20 MeV, event generator mode combined with evaluated nuclear data was utilized (Nitta et al., 2008).

Computational Methods
The source term used in this work is identical among all the codes and is generated using the Badhwar-O'Neill 2014 Galactic Cosmic Ray Flux Model for a solar modulation parameter of 500 MV (O'Neill et al., 2015), approximating conditions near a solar minimum. During such periods with low solar activity, the radiation environment near the lunar surface is dominated by GCRs. Here we use the most abundant isotope for each incident GCR ion from Z = 1 to Z = 28. Figure 1 illustrates the energy spectra of the source term used in this work, which had an isotropic directionality in each MC code input.
Due to the variations in the computational capabilities for each code, the number of source particle histories ran differed between the codes. FLUKA ran 6.3 million histories, PHITS and HETC-HEDS ran 3 million histories, and MCNP6 ran 1 million histories, all without the use of any special biasing techniques in an analog transport. However, the Geant4 runs were performed on a Linux cluster at the Aerospace Corporation, with about a thousand cores at a time running for several weeks. Consequently, the Geant4 calculations benefitted from the ability to simulate many more incident particles than any of the other codes, which is reflected in the smoothness of the Geant4 curves compared to others in the figures that follow. For each cosmic-ray primary ion species from Z = 1 to 28, isotropic primary particles in the energy range E from 0.1 MeV/n to 50 GeV/n were sampled with a differential spectrum proportional to 1/E below 10 MeV/n (equal number of particles per logarithmic energy interval) and proportional to 1/E 2 above 10 MeV/n (number of particles per logarithmic energy interval proportional to 1/E). This distribution of number of simulated primary particles versus energy was determined ad hoc because individual primary particles' shower histories take much longer to run at higher primary energies, and the total runtime for all particles had to be kept within reasonable bounds. No additional biasing on secondary-particle species, energy, direction, etc., was performed. The albedo results were then tabulated for 903 logarithmically spaced primary-particle energy bins, and reweighted by the GCR spectrum for each species. About 6E10 protons were simulated from 0.1 to 10 MeV/n and about 1E10 from 10 to 5E4 MeV/n; 1/10 these numbers of helium ions, 1/100 these numbers of carbon and oxygen ions, and 1/1,000 these numbers of other species were simulated. The geometry used in FLUKA, Geant4, MCNP6, and PHITS consisted of a sphere with a radius of 1,737 km, representing the radius of the Moon. Due to geometrical constraints for very large objects in HETC-HEDS, a cylinder deep and wide enough to include almost all the reactions (10-m thick and 10-km wide) was used instead.
Surface roughness was not included under the assumption that it would be smoothed out in observations made at satellite altitudes. The elemental composition of the material representing the lunar regolith is taken from Zaman et al. (2020) for "the average elemental composition of 10 different samples taken from the six Apollo missions, the three Soviet sample return missions, and ferroan anorthosite (FAN), a material widely used in simulating the lunar regolith." Table 1 presents the exact elemental composition used by all MC codes for this work.
The current (in units of # per GCR source particle per MeV of albedo-particle energy) of the following particles is tallied as a function of their energy as they leave the modeled lunar surface: neutrons, gamma rays, electrons, positrons, negative and positive muons, negative and positive pions, protons, deuterons, tritium ions, helium-3 ions, and alphas. For HETC-HEDS, the energy spectra of only the albedo protons and neutrons were tallied in the present work, with other particles, such as pions and light ions, to be investigated in future work. Finally, the angle of the albedo neutrons and protons as . The division between Z = 1, Z = 2, and Z = 3-28 is shown for demonstration purposes, but individual abundances and spectra for each ion species were used in the source term for all codes.   they leave the lunar surface is tallied with respect to the surface normal. The next section presents the resulting energy spectra and angular distributions for the different MC codes, with some of the plots that are discussed being relegated to Appendix B to keep the main body of the paper focused on the most significant results.

Neutrons
The energy spectra of albedo neutrons emitted off the modeled lunar regolith are illustrated in Figure 2. The neutron results for all five codes are in good agreement over the entire energy range considered. Initially, MCNP6 showed a relative underestimate of neutrons at energies below ∼1 MeV, being lower than other codes by around one order of magnitude at 0.1 MeV and two orders of magnitude at ∼0.01 MeV. Further investigation showed that this deviation was caused by (a) the default disabling of light-ion recoil physics and ions from neutron capture in MCNP6, and (b) the usage of nuclear models, instead of tabulated data, for all neutron interactions regardless of energy. This was resolved by controlling the proton and neutron physics cards and specifying ENDF/B-VII.1 evaluated data to handle low-energy neutrons. The new results show a much better agreement between MCNP6 and other codes as presented in Figure 2.

Gamma Rays
For gamma ray spectra, the results are presented in Figure 3. The general trends from all codes agree reasonably well, with the highest differences in the modeled continua reaching a factor of ∼3 at energies above 1 GeV. One key point in the MCNP6 results is the effect of neutral pions. By default, neutral pions are turned off in MCNP6, and the absence of neutral pion decay leads to the underestimate of high-energy gamma rays as noted by Ratliff et al. (2017) and . The agreement of MCNP6 results with those of other codes in this work was only reached after manually turning on neutral pions in the MCNP6 run.
A closer look at Figure 3 for energies below 10 MeV shows some variations between the codes. This is available in Figure B1 in Appendix B. Gamma rays generated by the radiative capture reactions and induced by low-energy neutrons, such as Fe photopeak at ∼7.6 MeV (Yamashita et al., 2013), appeared in the results only after removing any cutoff energy for neutrons. However, in PHITS, such photopeaks were only visible after turning on the event generator mode. FLUKA did not produce sharp photopeaks like other codes, although similar binning was used by all codes. This is caused by the group-wise treatment of some gamma ray production reactions in FLUKA. Last, although the codes agree on some photopeaks like the annihilation peak at 0.511 MeV and O peak at in the y-axis stands for galactic cosmic ray source particle). Note that the lower energy limit of neutrons in High-Energy Transport Code-Human Exploration and Development in Space is 1 MeV.
6.13 MeV (Reedy, 1978), the locations of other photopeaks vary among the codes and require a deeper look into their differences.

Electrons and Positrons
The energy spectra of albedo electrons and positrons are presented in Figures 4 and 5 respectively. For albedo electrons, the codes are in good agreement apart from Geant4 at lower energies, which deviates significantly at energies below ∼6 MeV and reports an electron current that is higher than other codes by a factor of ∼6 at 0.1 MeV. The curves for Geant4 are also smoother than those given by other codes for these particles, and the same is true for other rarer albedo species discussed below. This is due not to any differences between the codes themselves, but rather to the fact that we ran Geant4 on Linux clusters with thousands of cores as mentioned.
A slight increase reaching a factor of ∼1.5 in also observed in the current obtained from MCNP6 between 5 and 100 MeV relative to other codes. Initial investigation of the Geant4 deviation at lower energies shows that this might be due to electrons created by Compton scattering. In general, electrons are created through several processes including particle decay which dominates at higher energies, pair production, ionization, and photoelectric effect. Examining the electrons created through Compton effects in the different MC codes is the subject of future work. For albedo positrons, all the codes are within good agreement, as shown in Figure 5. Since decay processes dominate at higher energies, both albedo electrons and positrons exhibit a similar behavior at such  energy regions. This is illustrated in Figure B2 in Appendix B, which combines electron and positron spectra for each code. Every code has an energy range below which the current of albedo positrons starts to drop due to the dominance of ionization processes at lower energies.
Such symmetry was also seen in the modeled electron and positron spectra at the Martian surface. Compared to , the results of MCNP6 in this work are in better agreement with other codes. This, again following the work of Ratliff et al. (2017), is due to turning on neutral pions in MCNP6, which are turned off by default. Such agreement further confirms that, as stated by the authors in , not accounting for neutral pion decay caused the drop in MCNP6 electron and positron spectra in their work. In addition, differences in the level of agreement between the codes in this work and between the codes in  are expected for two main reasons: the teams responsible for running different codes in  were given more freedom in input parameters; and the presence of an atmosphere at Mars creates a modified source of incident radiation near the Martian surface due to interactions between the incident particles and the Martian atmosphere. Thus, the results in this work are expected to show better agreement due to the simpler conditions of the comparison.

Muons
Unlike the surfaces of Earth and Mars, the absence of an atmosphere at the Moon makes the decay of pions created in the lunar regolith the main if not the sole source of muons. This is reflected in the energy spectra of negative and positive muons in Figures B3 and B4 in Appendix B, where their magnitude and poor statistics relative to other albedo particles imply the relative low production of such particles in our simulations. Although the statistics are not good enough to arrive at a conclusion on how the energy spectra would be shaped (with MCNP6 e.g., running 1 million source particles), two unique features can be seen in the spectra of positive muons: a step at ∼4 MeV, and a step at ∼150 MeV. One possible reason, given the absence of lunar muon measurements, is that such features are artifacts of the nuclear models. A sophisticated study of lunar muons is outside the scope of this work, especially with the requirement of large computational resources for all the codes to improve the statistics of rare albedo species. However, further investigation of the step at 4.12 MeV showed that the likely reason is the decay of a positive pion at rest to a positive muon of 4.11981 MeV and an oppositely directed muon neutrino of 29.7919 MeV. The asymmetry in the negative muon spectra and the absence of the 4.12 MeV peak is probably due to the negative pion interaction with the target material and consequent absorption before it decays. As for the step at 150 MeV, it is not clear why such feature exists, and it requires further investigation into the underlying processes and nuclear models involving lunar muons.

Pions
Negative and positive albedo pion spectra are illustrated in Figures B5-B7 in Appendix B. Recall that any albedo pion is created in the lunar regolith. For negative pions, FLUKA reports a current higher by factor of ∼3 at higher energies compared to other codes, while MCNP6 starts to deviate slightly from Geant4 and PHITS, which are in good agreement throughout, at energies below ∼100 MeV. A similar behavior can be seen in the positive pion spectra, where FLUKA is also an outlier by a similar factor. Here, however, MCNP6 seems to be in better agreement with the other two codes. Figure B7 shows how symmetric the production of negative and positive pions is in each code. A symmetry can be seen in the spectra generated by FLUKA, and a slight increase is observed in the spectra of positive pions from Geant4, which has the best statistics compared to other codes as previously mentioned. For MCNP6 and PHITS, the statistics are affected negatively by the low production of albedo pions in regolith relative to other particles.

Protons
The study of albedo protons has been of great importance recently after the launch of the Cosmic Ray Telescope for the Effects of Radiation (CRaTER) instrument aboard the Lunar Reconnaissance Orbiter (Spence et al., 2010). As shown in Figure 6, all the codes are within good agreement when it comes to albedo protons, with the maximum difference being a factor of ∼1.5. The only exception is the current reported by MCNP6 at energies above ∼3,500 MeV, where its results are higher than those from other codes by a factor of ∼4. This is most likely due to the shift in nuclear models beyond 3,500 MeV in MCNP6 from CEM03.03 to LAQGSM03.03.

Light Ions
The comparisons of light ion spectra are presented in Figures B8-B11 in Appendix B. In general, the codes seem to agree reasonably well in their general trends. However, the convergence statistics in the MCNP6 and PHITS results are not sufficient, especially with He-3 and alpha ions, due to the relative low production of such particles.

Angular Distributions
The angular distributions of albedo neutrons and protons are presented in Figures B12 and B13 in Appendix B. Note that the tallied emission angle is off the surface normal, which means that zero degrees in Figures B12  and B13 corresponds to the particles leaving the lunar surface heading straight upward. For the angular behavior of albedo neutrons, all codes agree within a factor of ∼2, with the current decreasing at an increasing rate as the emission angle off surface normal increases. This is expected since first, all neutrons are created in the lunar regolith, and second, the angular behavior is in a direct relationship with the distance required to reach the lunar surface. From the same point in the lunar regolith, a particle escaping the lunar regolith in the straight up direction traverses a minimal distance to reach the surface, while a particle escaping at a shallow angle traverses a longer distance through the medium to reach the surface.
Investigating the angular behavior of albedo protons also shows an agreement by a factor of ∼2. Unlike albedo neutrons, a slight increase in the current of albedo protons up to angles between 70 and 80 degrees off surface normal can be noticed. This is due to the presence and dominance of protons in the GCR source term, specifically GCR protons striking the lunar surface at shallower incident angles. Zaman (2020) showed that GCR protons striking the lunar surface straight down produce an albedo proton angular distribution with a shape similar to that of albedo neutrons, where the current decreases at an increasing rate with increasing emission angle, unlike the GCR protons striking the lunar surface at shallower angles. In general, all the codes in this work show similar behaviors for albedo neutrons and protons, with the visible effect of the incident GCR protons on the angular distributions of albedo protons.
Albedo protons and neutrons were selected to show the asymmetric behavior of such particles. Most of the other species show a trend similar to that of albedo neutrons. Other energetic albedo charged particles, which are created closer to the surface through fragmentation and spallation processes, show trends similar to that of albedo protons. However, due to the dominance of the lower energy particles for such albedo species, the general trend incorporating all energies show a behavior similar to the one seen in the angular distribution of albedo neutrons. The angular distributions of HETC-HEDS were not included due to differences in the geometry model used.

Conclusions and Future Work
In this work, we compared the simulated lunar albedo radiation environment between five MC radiation transport codes: FLUKA, Geant4, HETC-HEDS, MCNP6, and PHITS. We attempted to make the input parameters as close to identical as possible to investigate regions of agreement and disagreement between the codes. In general, the codes were within good agreement for energy spectra of gamma rays, neutrons, electrons, positrons, pions, and protons, and light ions. However, some areas of disagreements were noticed, which include: • photopeaks in the gamma ray spectra for all the codes • a relative overestimate of low-energy electrons, probably those created by Compton scattering, in Geant4 • overall disagreement in muon energy spectra • a relative overestimate of pions in the high-energy range in FLUKA • a relative overestimate of albedo protons with energies above 3,500 MeV in MCNP6 (where the physics models change by default) In addition, the codes show two steps in the albedo muon spectra that require further investigation. When it comes to the angular behavior of albedo protons and neutrons, all the codes agree on the general trends of such particles.
Future work will include direct comparisons of the modeled codes with measurements taken from the CRaTER instrument, further investigation of each albedo species studied in this work, especially in areas where the codes disagree, and the inclusion of deterministic radiation transport codes (e.g., HZETRN) to the comparisons. Whereas this work studied the albedo radiation emitted from a hydrogen-free regolith as a generalized material for the lunar surface, future work will also include a comparative MC study of the albedo spectra emitted off a hydrogenated regolith, due to the importance of hydrogen for current and future missions to the Moon.
The prospects of sending humans back to the Moon highlight the importance of understanding the lunar radiation environment. Remarkably, recent observations by CRaTER (e.g., Schwadron et al., 2018, Zaman et al., 2021 suggest that the albedo radiation encodes information about the composition and hydrogenation of the regolith. This highlights the need to apply models to improve our understanding of the albedo radiation. Further, the albedo radiation presents a hazard to future human explorers. The model comparisons presented here equip the space weather and radiation biology communities with a better understanding of the results of popular radiation transport models applied to the lunar albedo radiation problem. The overall agreement and, in some cases, discrepancies between these transport codes provide new insight necessary for assessing code reliability and identifying where improvements are needed to advance our understanding of lunar albedo radiation and the space radiation environment in general.

Appendix A: Summary of Transport and Physics Models
This appendix provides an approximation of the transport and physics models used in the codes. Since each code uses extensive and complex algorithms covering varying and overlapping ranges of energies, reactions, and species, creating such simplified diagrams provides an incomplete picture of all the processes involved. Thus, this appendix is to be used only for a general comparison among the processes involved based on the available resources. For a complete picture of the models involved, please contact the respective teams that developed each code ( Figures A1-A4). Figure A1. A graphical approximation of the transport and physics models used in GEometry ANd Tracking (Wright, 2012).              Figure B12. The angular distributions of the albedo neutrons (above 10 keV) emitted off the modeled lunar regolith. Figure B13. The angular distributions of the albedo protons (above 10 keV) emitted off the modeled lunar regolith.

Data Availability Statement
The energy and angular data set used for FLUctuating KAscade, GEometry ANd Tracking, High-Energy Transport Code-Human Exploration and Development in Space, Monte Carlo N-Particle and Particle and Heavy Ion Transport code System spectra and distributions in the study are archived at Cosmic Ray Telescope for the Effects of Radiation's official website via https://crater-web.sr.unh.edu/Zaman2022/index.html.