Implementation and validation of the GEANT4/AtRIS code to model the radiation environment at Mars

A new GEANT4 particle transport model -- the Atmospheric Radiation Interaction Simulator (AtRIS, Banjac et al. 2018a. J. Geophys. Res.) -- has been recently developed in order to model the interaction of radiation with planets. The upcoming instrumentational advancements in the exoplanetary science, in particular transit spectroscopy capabilities of missions like JWST and E-ELT, have motivated the development of a particle transport code with a focus on providing the necessary flexibility in planet specification (atmosphere and soil geometry and composition, tidal locking, oceans, clouds, etc.) for the modeling of radiation environment for exoplanets. Since there are no factors limiting the applicability of AtRIS to Mars and Venus, AtRIS' unique flexibility opens possibilities for new studies. Following the successful validation against Earth measurements Banjac et al. 2018, J. Geophys. Res., this work applies AtRIS with a specific implementation of the Martian atmospheric and regolith structure to model the radiation environment at Mars. We benchmark these first modeling results based on different GEANT4 physics lists with the energetic particle spectra recently measured by the Radiation Assessment Detector (RAD) on the surface of Mars. The good agreement between AtRIS and the actual measurement provides one of the first and sound validations of AtRIS and the preferred physics list which could be recommended for predicting the radiation field of other conceivable (exo)planets with an atmospheric environment similar to Mars.


Introduction
There are mainly two types of energetic particles in the heliosphere that may impose risks for deep space and planetary missions: Galactic Cosmic Rays (GCRs) and Solar Energetic Particles (SEPs). GCRs are energetic charged particles, comprised of 2% electrons and 98% atomic nuclei with the later contributed by 87% protons, 12% helium, and about 1% heavier nuclei (Z ! 3) (Simpson, 1983). They have energies from less than 1 MeV nuc À1 up to hundreds of TeVs with a power-law distribution at energies above~GeV nuc À1 (Allkofer, 1975). SEPs on the other hand are energetic particles (mainly protons and electrons) emitted from the Sun and accelerated by solar flares and/or Coronal Mass Ejection (CME) associated shocks. SEPs generally have energies from a few keV up to several hundreds MeV (occasionally even reaching 1 or 2 GeV or further above) and can reach significantly higher fluxes at these energies compared to background GCRs.
It is likely that the first human-visited planet will be our neighbour planet Mars. It has a thin atmosphere with its surface pressure less than 1% of that at Earth's surface making it much easier for high energy particles to reach the Martian surface. Therefore, the assessment of the Martian radiation environment is necessary and fundamental for (a) mitigating radiation risks for nearfuture robotic and crewed missions and (b) better understanding the impact of energetic particles on the preservation of organic biosignatures on Mars. Primary GCRs and SEPs passing through the Martian atmosphere may undergo inelastic interactions with the ambient atomic nuclei losing their energies and also creating secondary particles via spallation and fragmentation processes. These secondary particles may further interact with the atmosphere as they propagate downwards and even with the Martian regolith, finally resulting in very complex spectra including both primaries and secondaries at the surface of Mars (e.g., Saganti et al., 2002).
There are various particle transport codes such as High charge (Z) and Energy TRaNsport (HZETRN, Slaba et al., 2016;Wilson et al., 2016), Particle and Heavy Ion Transport code System (PHITS, Sato et al., 2013) and GEometry And Tracking (GEANT4, Agostinelli et al., 2003;Allison et al., 2016) which can be employed for studying the particle spectra and radiation dose at Mars. Various studies have also combined these particle transport codes with different GCR and/or SEP spectra for estimating the radiation exposure on the surface of Mars (e.g., Simonsen et al., 1990;Simonsen & Nealy, 1993;Saganti et al., 2004;Keating et al., 2005;De Angelis et al., 2006;Ehresmann et al., 2011;McKenna-Lawlor et al., 2012;Guo et al., 2017a). Articles collected in a special issue at Life Sciences in Space Research  have included most recent studies to model the radiation environment at Mars which is also compared with in situ measurement at the surface of Mars (Ehresmann et al., 2017;Guo et al., 2017b) by the Radiation Assessment Detector (RAD, Hassler et al., 2012) on board the Curiosity rover belonging to the Mars Science Laboratory (MSL, Grotzinger et al., 2012).
In particular, PLANETOCOSMICS (http://cosray.unibe. ch/~laurent/planetocosmics/) is a toolkit based on GEANT4 with a specific application purpose to simulate particle transport in planetary magnetic fields as well as interactions when passing through planetary environments and creating secondary particles (Desorgher, 2005). Modeling the radiation environment on the surface of Mars using PLANETOCOSMICS has been carried out by previous researchers (e.g., Dartnell et al., 2007;Ehresmann et al., 2011;Gronoff et al., 2015;Matthiä et al., 2016;Guo et al., 2018c) and has been validated by Matthiä et al. (2016) when compared to energetic charged and neutral particle spectra on the surface of Mars measured by MSL/RAD.
Since the landing of MSL at the Gale Crater of Mars on August 6, 2012, RAD has been providing the first in situ detection of the radiation environment at the surface of Mars . These novel measurements provide evaluations of the Martian radiation level at an atmospheric depth around 22 g cm À2 . It has been found that the measured GCR radiation dose rate is anti-correlated with the surface pressure (proportional to the atmospheric column depth) which changes both daily and seasonally Guo et al., 2015Guo et al., , 2017a up to ±25%. Meanwhile the varying heliospheric conditions modulate the GCR fluxes and thus RAD measured radiation is anti-correlated with solar activity both in the long term evolution  and in the short term due to e.g., interplanetary CMEs and their associated shocks passing Mars (Witasse et al., 2017;von Forstner et al., 2018;Winslow et al., 2018;Guo et al., 2018b). RAD also measures the flux spectra of energetic charged particles, such as protons, deuterons, tritons, 3 He and 4 He, carbon, nitrogen, oxygen and iron ions which penetrate downwards and stop within the detector set with energies up to about 100 MeV nuc À1 (Ehresmann et al., , 2017. With an inversion method exploiting the response matrix of the neutral particle detection efficiency in the scintillators, RAD can also measure spectra of neutral particles on Mars from~10 to 900 MeV (Kohler et al., 2014;Guo et al., 2017b).
The last solar cycle (no. 24) has been rather quiet and RAD has only detected a few SEP events at the Martian surface and most of them are rather insignificant apart from the September 10, 2017 event Zeitlin et al., 2018;Guo et al., 2018a). This was the first Ground Level Enhancement (GLE) resulting from a very intense SEP detected at the surface of two different planets: Earth and Mars. In general, SEP events are sporadic, often rather impulsive and could be extremely hazardous especially when the observer has a direct magnetic connection to the particle acceleration region and injection site at the Sun. In the case of Mars, there is no global magnetosphere that could shield the atmosphere from energetic particles. In addition, its thin atmosphere can only stop charged particles below~150 MeV nuc À1 (Guo et al., 2018c). Therefore it is important to reliably model SEP events and their induced radiation in order to give immediate and precise alerts for future human missions to Mars. To do so, it is essential to well understand the interactions of atmospheric molecules with incoming particles.
Recently, a new GEANT4 particle transport modelthe Atmospheric Radiation Interaction Simulator (AtRIS, Sect. 2.2)has been developed by Banjac et al. (2018) in order to model the interaction of radiation with planets. The upcoming instrumentational advancements in the exo-planetary science, in particular transit spectroscopy capabilities of missions like JWST and E-ELT, have motivated the development of this particle transport code with a focus on providing the necessary flexibility in planet specification (atmosphere and soil geometry and composition, tidal locking, oceans, clouds, etc.) for the modeling of exoplanets. The application of AtRIS to model the radiation environment on Mars has been realized and validated for the first time in this study. The Mars Climate Database (MCD, http://www-mars.lmd. jussieu.fr) offers the possibility to access Martian atmospheric properties, such as temperature, density and composition, for different altitudes, seasons and even the time of the day on Mars. MCD has been developed using different Martian atmospheric circulation models which are further compared and modified by the observation results from past and current Mars missions (Lewis et al., 1999). Therefore, it provides a Martian atmospheric environment which can be implemented into the planetary particle transport toolkit for the purpose of simulating high energetic particles interacting with the Martian atmospheric atoms.
In our atmospheric setup for the AtRIS particle transport model, we use the composition, density and temperature profiles from MCD at Gale Crater on Mars where MSL's rover Curiosity landed on August 6, 2012 (coordinate: 4.5°S, 137.4°E). The elemental composition of the Martian atmosphere consists of C, O, N, Ar, and H with more than 95% of the molecules being CO 2 . The atmospheric condition is set to be "clim aveEUV" for climatology scenario with average solar EUV radiation. The Martian solar longitude is set to be 200.5°when the surface pressure is close to the annual average pressure at Gale Crater measured by MSL~840 Pa (e.g., Guo et al., 2015). Figure 1 shows the atmospheric structure, i.e., density and pressure versus the atmospheric altitude used for the AtRIS simulations.

Model description: AtRIS
GEANT4 is a Monte Carlo approach widely used for simulating the interactions of particles as they traverse matter (Agostinelli et al., 2003).
The Atmospheric Radiation Interaction Simulator, AtRIS, is a GEANT4 based particle transport code developed to simulate the propagation of energetic particles through planetary atmosphere and regolith. AtRIS allows rather flexible geometry and composition definitions of the planet. This flexibility and ease of use has been achieved through a custom interface called the planet specification format, that has been optimized for the specification of planets and their atmospheres. For Earth, an interface based on the so-called NRLMSISE-00 (Picone et al., 2002) model is provided within AtRIS. Similarly, for Mars, the MCD interface has been implemented (see Sect. 2.1).
AtRIS can calculate ion and electron pair production rates, secondary particle distributions (as a function of energy, directionality, planet altitude and so on), as well as absorbed and equivalent dose rates for a 30 cm diameter ICRU sphere phantom composed of water. The tracking of charged particles through magnetic fields is not implemented in AtRIS. As Mars lacks a global magnetic field, it serves as a good validation object for particle propagation modeling of AtRIS. Control over hadronic and electromagnetic processes is provided via the standard GEANT4 messenger and compounded physics list naming scheme (Allison et al., 2016;Geant4_Collaboration, 2017). Different physics lists tested and compared in this study will be discussed in detail in Section 3.1.
The main features of AtRIS are (i) the planet specification format as explained above, (ii) the Atmospheric Response Matrices (ARMs), quantifying the relation between primary energy and altitude-dependent ionization as well as equivalent/absorbed dose rate, and finally, (iii) the spectrum folding procedure used to calculate net quantities like the electron-ion pair production rate, by implementing a convolution of a measured spectrum and the ARM. A more detailed description is given in Banjac et al. (2018), where the results of AtRIS were compared with and validated against different kinds of Earth measurements (i.e., ion pair production and secondary particle fluxes, absorbed dose rate and dose equivalent).
In the current AtRIS setup for the Martian environment as shown in Table 1, we use a sphere with a radius of 3390 km (approximately the average radius of Mars) representing Mars. The soil composition is approximated as 50%Si, 40%O, and 10%Fe (mass fractions) and the crust (soil) sheet is 100 m thick. The maximum height of the atmosphere from the MCD model is 100 km which is divided into 500 m thick layers. The accumulated column depth at the surface is about 22 g cm À2 corresponding to a surface pressure of 830 Pa (Fig. 1). Table 1 provides the main features of the atmospheric and regolith structures used for the AtRIS simulations.

The matrix realization of AtRIS
Recently, Guo et al. (2018c) have developed a generalized approach based on the GEANT4/PLANETOCOSMICS transport code and the MCD Martian environment setups to quickly model the Martian surface radiation level of any given incoming proton/helium ion spectra. Such an approach called "Planetomatrix" can be mathematically described by ARMs and visually illustrated as 2-D histograms (Figs. 1 and 2 of Guo et al., 2018c).
Similar to this Planetomatrix approach, we construct the "response function" of the Martian atmosphere based on simulated particle spectra from AtRIS. We can describe the statistical transformation of the atomic and nuclear interaction process for a particle spectrum (of particle type i as a function of energy E 0 ) above the Martian atmosphere resulting in a particle spectrum (of particle type j as a function of energy E) on the Martian surface in a matrix M ij (E 0 , E). As particles could also interact with the Martian atmosphere and regolith and produce albedo particles contributing to the upward fluxes, we have also considered the generation of such upward particles in our simulations. Since the energy spectra of upward-and downward-traveling particles are dissimilar, we have separately constructed the upward and downward directed matrices for each primarysecondary case. Once such a matrix is constructed, it can be folded with any incoming SEP or GCR spectra within the energy range of the simulated particles for calculating the surface spectra without rerunning the particle transport simulations. The upward and downward secondary particle spectra generated by different types of primary particles can be combined. The following steps are used for constructing the matrix and folding it with a given incoming spectra to obtain the surface spectra: 1. Simulate primary particles (type i) with energy in the range of 1 and 10 5 MeV through the Mars AtRIS model described in Section 2.2. A flat spectrum over 50 energy bins distributed uniformly in logarithmic scale is used. There are N i (E 0 ) particles simulated in each energy bin. 2. Based on the simulation results create a matrix M ij (E 0 , E j ) for each secondary particle type j with certain directions, e.g., downward-directed protons on the surface of Mars generated by primary protons. Each column of the matrix is a histogram of secondary particles H i (E j ) representing the number of particles at energy E j created by N i (E 0 ) primary particles with energy E 0 . 3. Divide each histogram H i (E j ) in M ij (E 0 , E j ) by the number of simulated particles located in each incoming energy bin N i (E 0 ) to generate the normalized histogram due to a single primary particle with energy E 0 . The normalized histograms constitute the normalized matrix which is . For a given GCR/SEP spectrum f i which is often an isotropic flux in deep space (e.g., with units of particles/ sec/m 2 /sr/MeV), first interpolate this spectrum using the incoming energies with the unit of (particles/sec/m 2 /sr). 5. Calculate the total primary particle count rate (particles/ sec) taking into account of the integrated geometric factor of particles arriving at the planet by multiplying the above Here A in is the area of which primary particles were fed into the simulation. When using a sphere covering the top of the planetary atmosphere as the source sphere, A in is 4pR 2 with R (meter) being the radius of this sphere (i.e., R top in Table 1). X in is the integration of the cosine of zenith angle of the source particles and is p when they were injected isotropically inward from the source sphere. 6. Fold such scaled incoming particle count rate (energy dependent) with the normalized matrix to obtain the histogram of secondaries (particles/sec) corresponding to each column of the matrix (at certain incoming primary energy). 7. Scale the secondary particle histogram at each column into a differential spectra per geometric factor by dividing it with X out A out W(E j ). Here W(E j ) is the energy bin width (MeV) of the secondary particle histogram. A out is the surface area of the planet where the output secondaries are counted and it is 4pR 2 surf with R surf (see Table 1) in centimeter. And the X out is the integration of the cosine of zenith angle of secondaries crossing the surface plane.
In this study, we consider three cases of secondary particle directions from AtRIS simulations: (a) all secondaries propagating downwards with the flux-weighted solid angle X out of p, (b) all secondaries propagating upwards with X out = p and (c) secondary downward directed particles with a zenith angle within 36°, i.e., Q out = 1.1 sr. This corresponds to the inner cone angle used for counting downward propagating charged particles stopping inside the RAD detector (Hassler et al., 2012).
Finally the surface secondary spectra F ij (E j ) induced by primary particle spectrum F i (E 0 ) obtained through step 1-7 has unit of (particles/sec/cm 2 /sr/MeV). The above procedure can be also mathematically summarized as: where M ij E 0 ; E j À Á is the atmospheric response matrix which takes into account the scaling of the energy bin widths of the output histograms and the geometric factors used in the simulation. The surface spectrum of particle type j resulting from primary particle type i is the multiplication of the matrix with the primary particle spectrum and this operation is essentially the sum product over the second axis (E 0 ) of M ij and F i , i.e., The simulation is often set up with the areas A in and A out equal to each other where the planet is approximated as a rectangular box which has the same size on top of the atmosphere and on the regolith surface. In the current AtRIS/MCD setup, A in is only slightly bigger than A out as R top is the radius of the planet including the atmospheric layer and it is slightly larger than R surf . Given three cases of solid angle of secondaries detected on the surface, i.e., downward, upward and within 36°z enith angle (RAD inner cone), a matrix for each case can be constructed for certain secondary particles induced by given primary particle type.
Although the construction of each matrix is time-consuming, the multiplication of different input spectra F i (E 0 ) with such a matrix to generate different surface spectra F j (E j ) is very much simplified. The spectra of the surface secondary particle are the combined spectra resulting from different types of primary particles which arrive at the top the Martian atmosphere: here F dn j ðE j Þ, F up j ðE j Þ and F RAD j ðE j Þ represent the energy spectra of secondary particle type j, resulting from various primary particle types, averaged in the downward, upward and within 36°zenith angle respectively. As a first step of the model verification, we have simulated protons and 4 He ions as primary particles because they constitute the majority of GCR particles. We considered the secondary charged particles (type j) whose spectra are also measured by RAD on the surface of Mars, such as protons, deuterons, tritons, 4 He and 3 He ions as output particles.

Validation of AtRIS
GEANT4 (version 10.4.p02 used here) offers a wide variety of models for handling physical processes within different energy ranges. It is not yet entirely clear what is the most accurate and efficient physics list describing high energy (from hundreds of MeV to tens of GeV) particles integrating with the Martian atmosphere . One of the goals of MSL/RAD is to help validate the appropriate physics list which could precisely model the high energetic cosmic ray interaction with the Mars atmosphere (Hassler et al., 2012).
In this study, we employ four different physics lists (Table 1) FTFP_BERT_HP in model C is recommended by Geant4 collaboration "for cosmic ray applications where good treatment of very high energy particles is required" (Geant4_Collaboration, 2017). It contains all standard electromagnetic (EM) processes by default. It uses Bertini-style cascade for hadrons <5 GeV and the FTF (Andersson et al., 1987;Nilsson-Almqvist & Stenlund, 1987) model for simulating the interaction of mesons, nucleons and hyperons in the 3 GeV-100 TeV energy range. Finally and in model D, the Liege Intra-nuclear Cascade model INCL++ has been recently extended in GEANT4 to handle reactions between 3 and 15 GeV incident energies (Mancusi et al., 2014) and extensive benchmarks have shown the INCL++ model has a very good predictive power for the particles related to neutron production in spallation reactions (Leray et al., 2011). A map of physics models in different energy ranges used for the INCLXX physics list can be found here http://irfu.cea.fr/ dphn/Spallation/physlist.html.
We have validated the accuracy of the AtRIS approach when applied with different physics lists (named model A, B, C and D) to the Martian environment, especially for the generation of charged particle spectra by primary protons and 4 He ions. We first compare the ARMs of each primary-secondary pair when using different physics lists as shown in Section 3.1. Then we compare such modeled surface secondary spectra based on different physics lists with the MSL/RAD measurements as explained in Section 3.2.

Matrices constructed via AtRIS
Using the matrix approach described in Section 2.3, we compare the ARMs derived from above four different models: A) QGSP_BIC_HP, B) QGSP_BERT_HP, C) FTFP_BERT_HP and D) FTFP_INCLXX_HP. To ease the comparison between different simulation results, normalized matrices M N ij E 0 ; E j À Á (step 3 in Sect. 2.3) have been plotted and discussed. The interaction of primary protons with the Martian atmosphere has been modeled and ARMs have been constructed for five different secondary types including Hydrogen isotopes (proton, deuteron and triton) and Helium isotopes ( 3 He and 4 He). Each secondary type is differentiated as downward-and upward-directed with the respective matrix. A total of ten matrices have been generated for primary proton interactions. Similarly ten other matrices represent the primary 4 He ion interaction with the Martian atmosphere. Although there are other secondaries such as electrons, muons, positrons which also contribute to the radiation environment on Mars, we focus on the five hydrogen and helium isotopes as their spectra have been well measured on the surface by MSL/RAD ) and a direct model-observation comparison is possible. Figure 2 shows the matrices which describe the surfaced downward secondary protons which result from the primary protons under an atmospheric configuration as described in Section 2.1. The general statistics of the matrices constructed from four different physics lists are rather similar as shown in the top row of the figure. We have also accounted for the Poisson uncertainty in the Monte Carlo simulations in each bin of the matrix and scaled it with the statistics in the corresponding bin to obtain the normalized uncertainties as plotted in the second row of the figure. As shown, the normalized uncertainties are relatively small (below 20%) in most bins apart from the low energy output bins where the bin widths are smaller (energies are binned into logarithmic scale) with lower statistics.

Primary proton matrices
To better quantify the differences between models, we have calculated the (relative) image difference between matrices of different models as shown in the last two rows of Figure 2, which we define as Image Differencing Matrices (IDM). It is visible that there are some systematic differences between different models and boundaries of patterns are formed at energies where physics processes are switching in the model.

For instance, A-B IDM reflects the difference between BIC and
BERT models for incoming particles below~10 GeV and the consistency between two physics lists at higher energies where QGSP is used for both models. Alternatively, B-C IDM shows the difference between QGSP and FTFP at high energies and the agreement of BERT at lower energies. In comparison, C-D shows a good agreement between FTFP_BERT_HP and FTFP JNCLXX_HP apart from low output energy ranges where BERT seems to have a higher efficiency. This is because when simulating p-p interactions at energies below 20 MeV, FTFPJNCLXX_HP uses the models from the HP library. INCL is turned on only when one of the involved particles has the atomic number larger than one.
As shown in the bottom row, the mean relative differences are about 20.9% between A and B, 23.2% between B and C and 17.7% between C and D which are not much larger in comparison to the mean relative uncertainties of each model (second row). Therefore the relative IDM indicates a general good agreement between different physics lists for protons generating secondary protons on the surface of Mars. Besides, a 2D Gaussian filter has also been applied to the IDM of models C and D, shown in the last column of the last two rows, in order to smooth out features resulting from fine bins with low statistics. Figures 3 and 4 show the matrices of primary protons generating secondary deuterium and tritium hydrogen isotopes detected as downward particles on the surface of Mars.  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins.
Compared to the generation of secondary protons, the p-d and p-t reaction probability is much lower as can be seen in the matrices of the top panels which have the same color bar scale as the proton-proton matrices. The deuterium and tritium matrices lack the high energy diagonal component, while most secondaries are due to spallation reactions of high energy protons (above a few GeV) with the atmosphere nucleus. Model D is clearly more efficient in generating deuterium and tritium particles via such reactions as is readily seen in the model D matrices. Model A in the range with the binary cascade model is least efficient in such productions. Indicated by the IDM between different models, in the energy range below~1 0 GeV, deuteron and triton generation efficiency of different models follow: BIC < BERT < INCL; at higher energies, FTFP has a slightly higher efficiency than QGSP. Figure 5 shows the matrices of primary protons resulting into surface downward 4 He ions. These reactions are even more rare than the above proton-deuteron and proton-triton cases. Despite of the low statistics of such reactions, the IDM between A and B indicates a slightly higher efficiency of BERT and the IDM between B and C reflects a higher production rate in FTFP compared to QGSP. And it is also noticeable from the IDM between C and D that INCL predicts a higher probability of such spallation interactions in the relevant energy range.  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins. Finally, Figure 6 shows surface downward 3 He ions induced by primary protons. The features of matrices in this case are very similar to those of deuterium, tritium and 4 He ions in Figures 3-5 indicating similar cascading processes in generating these different types of secondaries. However 3 He ions have a slightly lower statistics than 4 He ions which may suggest some of the 3 He ions are from fragmented 4 He ions first generated in the atmosphere.
Primary protons and their secondaries generated in the atmosphere may reach the Martian surface and generate upward particles in the regolith. As an example, we only show the matrices constructed for the case of upward protons as shown in Figure 7. Similar to the downward proton case in Figure 2 the IDM of A-B and B-C show higher efficiency of BERT in the middle energy range and FTFP in the high energy part. The IDM of C-D shows features that are more similar to those of downward protons rather than those of deuterium, tritium and  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins. This is also suggested by the similar atmospheric cutoff energies of primary protons when comparing the downward and upward secondary proton matrices.
The atmospheric cutoff energy is about 160 and 200 MeV in Figures 2 and 7 respectively which approximately corresponds to the minimum energy a proton needs to traverse through 22 g cm À2 of atmosphere.

Primary 4 He ion matrices
Similar to the ARMs of primary protons, the primary 4 He ion interaction with the Martian atmosphere and the generation of secondaries have been mapped to matrices with some examples shown in Figures 8-11. Figure 8 displays the matrices for 4 He generated surface protons. The general shape from four physics lists looks rather similar and also very comparable to the proton-proton matrices in Figure 2, indicating that a substantial amount of 4 He ions fragment into protons in the atmosphere and follow the same energy loss processes as in the proton-proton case. It is noticeable that model D predicts a slightly higher efficiency in the primary energy range from a few GeV to~20 GeV corresponding to the INCL model in the physics list of model D.  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM. m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins. Figure 9 shows the 4 He-4 He matrices from four physics lists which have similar structures. It is visible that there are three components in these matrices. The first is represented by the top-right diagonal component showing some 4 He particles with energies J 2-3 GeV arriving at the Martian surface as 4 He ions with similar energies, i.e., these relativistic particles propagate through the 22 g cm À2 of atmosphere without any interactions. The second component is shown as a vertical branch connecting to the lower-left part of the diagonal structure. These are primary ions with energies between~650 MeV and~3-4 GeV which traverse through the atmosphere with substantial ionization energy loss. The third component is the blob structure at the lower-right part of the matrices (similar to the p-d, p-t matrices in Figs. 3 and 4) which represent high energy primary particles interacting with the atmosphere and generating secondary 4 He ions via spallation process. There is a gap beneath the diagonal component while there is an enhanced production of 4 He produced protons at this range (Fig. 8). This suggests that 4 He particles easily fragment into protons as they traverse through the atmosphere.
As shown by the IDM in the lower panel, all four models result in very small differences especially for the first component. The difference between different models is mostly due to statistics as the IDM shows rather similar features compared  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM. m represents the average value in each matrix (while bins with zero statistics have been excluded) and a is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins.
to the uncertainty matrices and the mean values of IDM are comparable to the mean matrix uncertainties. Model D seems to predict a larger effect of the third component, especially in the energy range where the INCL model is applied. Figure 10 represents the 4 He-deuteron case and also shows the three-component structure in four models. Here the first and second components correspond to 4 He ions fragmenting into 2 H particles as they propagate through the atmosphere. It is also shown by the IDM of A-B that Bertini cascade is more efficient and Binary cascade. The third component predicted by model D is much more enhanced compared to model A and B and the transition between the second and the third components in model D is much smoother. The average difference between model D and C is significant and as large as 100%. Similar features of the matrices are also shown in Figure 11 for the 4 He-3 He case with model D predicating an enhancement of 3 He production especially at the range between the second and third components. This might be due to 4 He being more efficient in fragmenting into 3 He particles.

Comparison of modeled spectra with RAD proton measurements
In order to generate surface secondary particle spectra, we fold the above matrices with the primary proton and helium GCR spectra calculated from a standard GCR modelthe  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins. Badhwar-O'Neill 2010(BON10, O'Neill, 2010. The BON10 spectra have been computed using the spherically symmetric Fokker-Planck equation, with the Local Interstellar Spectrum (LIS) at the boundary of the heliosphere (~100 AU) as a boundary condition. The model uses an input parameter of the solar modulation U, in units of volt, as derived from the International Sunspot Number. The modulation parameter U is a depiction of solar activity and its modulation magnitude on the GCR particle fluxes. The modeled spectra at Earth have been fitted and adjusted to measurements obtained by a multitude of instruments, including all available observations of GCR particle energy spectra from 1955 to 2010. Note that we use the near-Earth GCR spectra to approximate those at Mars with a distance of~1.5 AU to the Sun. The radial gradient of GCRs in the inner heliosphere has been estimated to be about 3% per AU for protons with energies around 1 GeV (e.g., Gieseler & Heber, 2016). This translates into about 1.5% of increase of GCR flux from Earth and Mars. However, this correction is much smaller compared to other systematic and statistical uncertainties in the models employed in the current work and therefore we omit the radial gradient of the heliospheric modulation of GCRs.
The solar modulation parameter used here as input for the BON10 GCR model is equal to the average value, 550 mV,  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins during the periods of the selected RAD measurements (Ehresmann et al., 2017). Primary protons and 4 He ions from the GCR model have been folded with the matrices to generate surface secondary particles. The resulting secondary spectra of protons, 4 He and 3 He ions, deuterons, tritons in the energy range of 10-100 MeV nuc À1 in the RAD view cone have been compared with MSL/RAD surface measurements as shown in Figure 12. Although RAD is not always positioned perpendicular to the surface as the rover body is inclined occasionally, Wimmer-Schweingruber et al. (2015) have found from the measurement that the variability of downward particle fluxes with the zenith angle is small. Therefore a separate set of matrices with surface particles constrained within the RAD view cone (<36°of downward zenith angle) has been constructed. Besides, we have also propagated the Poisson uncertainties from the simulations through the matrix multiplication for obtaining the uncertainties of the modeled surface spectra (Eq. (3)).
The proton spectra modeled by four different physics lists all agree very well with each other and also the RAD measurement. The modeled fluxes are slightly lower than the observed spectra which is reasonable as GCR particles heavier than 4 He have not been considered in the model. Despite the large  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins. uncertainties due to low statistics, it is clearly shown that the fluxes for both deuterium and tritium particles from model A, B and C are considerably lower than the measured spectra while those predicted by model D match well with the data. As explained in Section 3.1 and shown in Figures 3, 4, and 10, model D has a higher efficiency in generating secondary particles via spallation process (or the third component) which are located at the energy range RAD measures.
For Helium isotopes, all four models show a good agreement of 4 He ions, although with a slightly higher flux, compared to the RAD data. For 3 He case, model D again has a better prediction in comparison to the RAD measurement in the energy range of 10-100 MeV/nuc. The higher production rate of 3 He in model D is clearly shown in Figures 6 and 11 as explained in Section 3.1.
We have also compared the modeled spectra in a larger range beyond the energy range of the RAD data as shown in Figure 13(a) and (d) for H and He isotopes respectively. For protons and 4 He ions, all models agree with each other at different energy ranges. However for deuterons, tritons and 3 He ions, model D shows significant enhancement of the flux compared to other models for particles up to~500 MeV/nuc. At higher energies the discrepancy between different models is much smaller as also shown in the deuterium, tritium and 3 He matrices.  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins.
The RAD measurement locates within the energy range where predictions from different physics lists differ the most and the current comparison strongly supports the FTFP JNCLXX_HP model for GCR particles transported through the Martian environment.
The angular dependence of the secondary particles on the surface of Mars is also very important for future human exploration of Mars and understanding radiation effects on Mars. The downward spectra within the RAD view cone shown in Figure 13(a) and those averaged over the whole downward-directed solid angle shown in panel (b) have very similar levels of flux. This is consistent with the results found by Wimmer-Schweingruber et al. (2015) from the RAD measurement that within 15°of the Rover tilt angles, the radiation field is mostly isotropic.
However the upward spectra are rather dissimilar from the downward spectra and all upward spectra lack the high energy component ( J GeV) compared to downward spectra. Particle fluxes of protons and 4 He ions are particularly lower than in the downward direction where more primary GCR particles contribute to the surface downward radiation. As the upward flux, mostly produced in the Martian regolith, is much lower than the downward flux, this suggests the effectiveness of shielding using the Martian regolith for potential human habitat on Mars.  (Table 1). Second row: Normalized statistical uncertainty of the matrices. Third row: Image Differencing Matrices (IDM) of model A and B (1st column), B and C (2nd column), C and D (3rd column) which is also visualized with a 2D Gaussian filter applied (4th column). Fourth row: the normalized IDM.
m represents the average value in each matrix (while bins with zero statistics have been excluded) and r is the standard deviation. All matrices are shown in the energy range from 1 MeV to 100 GeV (in logarithmic scale) for both input bins and output bins.

Summary and conclusion
We have applied the novel tool AtRIS (Banjac et al., 2018) with a specific implementation of the Martian atmospheric and regolith structure to model the radiation environment at Mars.
We have validated the accuracy of the AtRIS model when applied with four different physics lists to the Martian environment, especially for the generation of charged particle spectra by primary protons and 4 He ions. We first visualized the atmospheric response matrix of each primary-secondary chain via a  (Table 1). The spectra are averaged within 36°of downward zenith angle corresponding to the RAD view cone. (a) Input: GCR H + 4 He, output: downward H isotopes, (b) input: GCR H + 4 He, output: downward He isotopes.
2-D histogram where the energy-dependent efficiency of primary particles generating secondaries is nicely shown. We compared these matrices obtained via four different physics lists of GEANT4 by calculating and visualizing their IDMs. The IDM results suggest that in general the Bertini cascade model is more efficient than the Binary cascade model while INCL model is most productive especially for deuterium, tritium and 3 He particles. For two primary particle types of protons and 4 He ions, four different physics lists agree with each other reasonably well. For other hydrogen and helium isotopes, model D FTFP_INCLXX_HP generates significantly more secondaries due to spallation processes of primary particles in the energy range of~GeV and 20 GeV where the Liege Intra-nuclear Cascade model is applied, as shown in Figures 3, 4, 6, 10 and 11.
To benchmark the AtRIS model based on different GEANT4 physics lists when applied to the Martian environment, we compared the AtRIS predicted spectra with the energetic particle spectra recently measured by MSL/RAD on the surface of Mars (Ehresmann et al., 2017). We folded the above ARMs with BON GCR spectra of protons and 4 He particles to obtain the secondary particle spectra of Hydrogen and Helium isotopes: protium, deuterium, tritium, 3 He and 4 He particles. This process has summed up the same secondary type generated by both proton and 4 He primaries. But we have ignored the contribution by other heavier GCR particles. This may have caused the slightly smaller flux of proton in all four models in comparison with the MSL/RAD measurement. The comparison shows that all models have an agreed prediction of the proton and helium flux in comparison to the RAD measurement in the energy range of 10-100 MeV/nuc. However, model D has a higher prediction of deuterium, tritium and 3 He flux and a better match with the RAD data. This is due to the same reason we have already observed in the ARMs.
In general, the good agreement between AtRIS and the actual measurement extends the validation of AtRIS to non-Earth systems. For the specific case of a very thin atmosphere, we have analyzed how different physics lists may impact the generation of secondary particles through the Martian atmosphere. While we have seen that FTFPJNCLXX_HP provides a better agreement with the MSL/RAD data for the specific cases that were examined, a further study is needed to validate this for heavier primary ions and other secondary particles as well as for planets with thicker atmospheres and different atmospheric compositions (Earth-like and Venus-like). For instance, secondary muons are a major surface radiation component at the largest moon of Saturn-Titan. Heavier ions and neutrons have more enhanced biological effectiveness and should also be carefully examined and included when using AtRIS to predict the equivalent radiation dose.
In order to make AtRIS easily accessible by the community, we are preparing an online documentation including examples for AtRIS in the form of a wiki page. In 2019, AtRIS will be published under a GNU GPL licence.