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

TL;DR

This paper presents the implementation and validation of the GEANT4/AtRIS code for modeling the radiation environment at Mars, demonstrating its accuracy through comparison with recent surface measurements and validating its applicability for planetary radiation studies.

Contribution

The work adapts and validates the AtRIS particle transport model for Mars, extending its use from Earth to Mars and demonstrating its effectiveness with experimental data.

Findings

AtRIS accurately models Martian radiation environment.

Good agreement with RAD surface measurements.

Preferred physics list identified for planetary radiation prediction.

Abstract

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.

Tables

Table 1. Table 1: Martian atmospheric and regolith properties used in AtRIS and four different physics lists of GEANT4 (more details in the text) tested in this study.

	model A	model B	model C	model D
physics lists	QGSP_BIC_HP	QGSP_BERT_HP	FTFP_BERT_HP	FTFP_INCLXX_HP
atmosphere	100 km thick spherical shell: $R_{surf}$ = 3390 km, $R_{top}$=3490 km
setup	MCD atmospheric profile at Gale Crater with a total column depth of $\sim$ 22 g/cm2
regolith	100 m thick of soil spherical shell beneath the atmosphere
setup	1.7 g/cm3, 50%Si, 40%O, and 10% Fe

Equations

$${\bar{\mathrm{M}}}_{ij}(E_0,E_j)$$

$$F_{ij}(E_j)$$

$$F_j^{dn}(E_j)$$

$$F_j^{up}(E_j)$$

$$F_j^{RAD}(E_j)$$

Full text

Following the successful validation against Earth measurements (Banjac et al., 2018), 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.

1 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$\geq$ 3) (Simpson, 1983). They have energies from less than 1 MeV/nuc up to hundreds of TeVs with a power-law distribution at energies above $\sim$ GeV/nuc (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 percent 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 near-future 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 HZETRN (High charge (Z) and Energy TRaNsport, Slaba et al., 2016; Wilson et al., 2016), PHITS (Particle and Heavy Ion Transport code System, Sato et al., 2013) and GEANT4 (GEometry And Tracking, 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 and Nealy, 1993; Saganti et al., 2004; Keating et al., 2005; De Angelis et al., 2006; McKenna-Lawlor et al., 2012; Ehresmann et al., 2011; Guo et al., 2017a). Articles collected in a special issue at Life Sciences in Space Research (Hassler et al., 2017) 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/$\sim$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; Gronoff et al., 2015; Matthiä et al., 2016; Ehresmann et al., 2011; 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 (Hassler et al., 2014). These novel measurements provide evaluations of the Martian radiation level at an atmospheric depth around 22 g/cm2. 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 (Rafkin et al., 2014; Guo et al., 2015, 2017a) up to $\pm$ 25%. Meanwhile the varying heliopsheric conditions modulate the GCR fluxes and thus RAD measured radiation is anti-correlated with solar activity both in the long term evolution (Guo et al., 2015) and in the short term due to e.g., interplanetary CMEs and their associated shocks passing Mars (Witasse et al., 2017; Guo et al., 2018b; von Forstner et al., 2018; Winslow et al., 2018). RAD also measures the flux spectra of energetic charged particles, such as protons, deuterons, tritons, 3He and 4He, carbon, nitrogen, oxygen and iron ions which penetrate downwards and stop within the detector set with energies up to about 100 MeV/nuc (Ehresmann et al., 2014, 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 $\sim$ 10 to 900 MeV (Köhler 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; Ehresmann 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 $\sim$ 150 MeV/nuc (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 model – the Atmospheric Radiation Interaction Simulator (AtRIS, Section 2.2) – has been developed by Banjac et al. 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 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.

2 Model Description and Implementation

2.1 The Mars Climate Database

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 CO2. 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 $\sim$ 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.

2.2 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 Mars Climate Database (MCD) interface has been implemented (see Section 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., ), 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/cm2 corresponding to a surface pressure of 830 Pa (Figure 1). Table 1 provides the main features of the atmospheric and regolith structures used for the AtRIS simulations.

2.3 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 atmospheric response matrices and visually illustrated as 2-d histograms (Figure 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 $\rm{\bar{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 primary-secondary 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 re-running 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:

Simulate primary particles (type $i$) with energy in the range of 1 and 105 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. 2.

Based on the simulation results create a matrix $\rm{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. 3.

Divide each histogram $H_{i}(E_{j})$ in $\rm{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 $\rm{M_{ij}^{N}(E_{0},E_{j})}\coloneqq\rm{M_{ij}(E_{0},E_{j})}/N_{i}(E_{0})$. 4. 4.

For a given GCR/SEP spectrum $f_{i}$ which is often an isotropic flux in deep space (e.g., with units of [particles/sec/m2/sr/MeV]), first interpolate this spectrum using the incoming energies ($E_{0}$) of the matrix $\rm{M_{ij}(E_{0},E_{j})}$ so that $f_{i}=f_{i}(E_{0})$. Then multiply it with $W(E_{0})$ which is the incoming bin width [MeV] of the matrix such that $F_{i}(E_{0})=W(E_{0})f_{i}(E_{0})$ with the unit of [particles/sec/m2/sr]. 5. 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 $F_{i}(E_{0})$ with $\Omega_{in}A_{in}$. 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 4$\pi R^{2}$ with $R$ [meter] being the radius of this sphere (i.e., $R_{top}$ in Table 1). $\Omega_{in}$ is the integration of the cosine of zenith angle of the source particles and is $\pi$ when they were injected isotropically inward from the source sphere. 6. 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. 7.

Scale the secondary particle histogram at each column into a differential spectra per geometric factor by dividing it with $\Omega_{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 4$\pi R_{surf}^{2}$ with $R_{surf}$ (see Table 1) in centimeter. And the $\Omega_{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 $\Omega_{out}$ of $\pi$, b) all secondaries propagating upwards with $\Omega_{out}=\pi$ and c) secondary downward directed particles with a zenith angle within 36∘, i.e, $\Omega_{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/cm2/sr/MeV]. The above procedure can be also mathematically summarized as:

[TABLE]

where $\rm{\bar{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 $\rm{\bar{M}_{ij}}$ and $F_{i}$, i.e.,

[TABLE]

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∘ zenith 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:

[TABLE]

Here ${F}_{j}^{dn}(E_{j})$, ${F}_{j}^{up}(E_{j})$ and ${F}_{j}^{RAD}(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 4He 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, 4He and 3He ions as output particles.

3 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 (Matthiä et al., 2016). 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): QGSP_BIC_HP (model A), QGSP_BERT_HP (model B), FTFP_BERT_HP (model C) and also FTFP_INCLXX_HP (model D) when applying AtRIS to model the Martian radiation environment. These physics lists are named based on the combination of various models in different energy ranges (Geant4_Collaboration, 2017):

• •

QGS stands for the Quark Gluon String model for high energy particles ($>\sim$ 20 GeV).

• •

FTF represents the Fritiof model for particles $>\sim$ 5 GeV.

• •

P stands for the Precompound model used for de-excitation process for nucleon-induced reactions below 1-2 MeV.

• •

BIC uses the Binary Cascade Model for particles $<\sim$ 10 GeV.

• •

BERT uses the Bertini Cascade Model in the middle energy range of $<\sim$ 10 GeV.

• •

INCL represents the Liége Intra-nuclear Cascade model in the middle energy range.

• •

HP option switches on the high precision neutron elastic and inelastic scattering model for neutrons below 20 MeV.

• •

At energies where different physics models are overlapping, a weighted combination of various models is considered.

Specifically speaking, model A) QGSP_BIC_HP and model B) QGSP_BERT_HP use the same Quark Gluon String (QGS) model for the high energy range while different cascade models for the lower energy range ($<\sim$ 10 GeV). 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 and 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 Liège 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 4He ions. We first compare the atmospheric response matrices 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.

3.1 Matrices constructed via AtRIS

Using the matrix approach described in Section 2.3, we compare the atmospheric response matrices 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 $\rm{M_{ij}^{N}(E_{0},E_{j})}$ (step 3 in Section 2.3) have been plotted and discussed.

The interaction of primary protons with the Martian atmosphere has been modeled and atmospheric response matrices have been constructed for five different secondary types including Hydrogen isotopes (proton, deuteron and triton) and Helium isotopes (3He and 4He). 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 4He 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 (Ehresmann et al., 2014) and a direct model-observation comparison is possible.

3.1.1 Primary proton matrices

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.

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 $\sim$ 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_INCLXX_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, FTFP_INCLXX_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. 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 $\sim$ 10 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 4He 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.

Finally, Figure 6 shows surface downward 3He ions induced by primary protons. The features of matrices in this case are very similar to those of deuterium, tritium and 4He ions in Figures 3, 4, and 5 indicating similar cascading processes in generating these different types of secondaries. However 3He ions have a slightly lower statistics than 4He ions which may suggest some of the 3He ions are from fragmented 4He 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 3He particles, indicating that these are mostly downward propagating protons reaching the soil and get re-directed upwards rather than protons produced via spallation in the regolith. 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 $\sim$ 22 g/cm2 of atmosphere.

3.1.2 Primary 4He ion matrices

Similar to the atmospheric response matrices of primary protons, the primary 4He ion interaction with the Martian atmosphere and the generation of secondaries have been mapped to matrices with some examples shown in Figures 8, 9, 10 and 11.

Figure 8 displays the matrices for 4He 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 4He 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 $\sim$ 20 GeV corresponding to the INCL model in the physics list of model D.

Figure 9 shows the 4He-4He matrices from four physics lists which have similar structures. It is visible that there are 3 components in these matrices. The first is represented by the top-right diagonal component showing some 4He particles with energies $\geq\sim$ 2-3 GeV arriving at the Martian surface as 4He ions with similar energies, i.e, these relativistic particles propagate through the 22 g/cm2 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 $\sim$ 650 MeV and $\sim$ 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 Figures 3 and 4) which represent high energy primary particles interacting with the atmosphere and generating secondary 4He ions via spallation process. There is a gap beneath the diagonal component while there is an enhanced production of 4He produced protons at this range (Figure 8). This suggests that 4He 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 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 4He-deuteron case and also shows the three-component structure in four models. Here the first and second components correspond to 4He ions fragmenting into 2H 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 4He-3He case with model D predicating an enhancement of 3He production especially at the range between the second and third components. This might be due to 4He being more efficient in fragmenting into 3He particles.

3.2 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 model – the 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 ($\sim$ 100 AU) as a boundary condition. The model uses an input parameter of the solar modulation $\Phi$, in units of volt, as derived from the International Sunspot Number. The modulation parameter $\Phi$ 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 $\sim$ 1.5 AU to the Sun. The radial gradient of galactic cosmic rays in the inner heliosphere has been estimated to be about 3% per AU for protons with energies around 1 GeV (e.g., Gieseler and 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, during the periods of the selected RAD measurements (Ehresmann et al., 2017). Primary protons and 4He ions from the GCR model have been folded with the matrices to generate surface secondary particles. The resulting secondary spectra of protons, 4He and 3He ions, deuterons, tritons in the energy range of 10 to 100 MeV/nuc 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∘ degrees 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 (Equation 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 4He have not been considered in the model. Despite the large 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 4He ions, although with a slightly higher flux, compared to the RAD data. For 3He 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 3He 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 4He ions, all models agree with each other at different energy ranges. However for deuterons, tritons and 3He ions, model D shows significant enhancement of the flux compared to other models for particles up to $\sim$ 500 MeV/nuc. At higher energies the discrepancy between different models is much smaller as also shown in the deuterium, tritium and 3He matrices. 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_INCLXX_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 ($>\sim$ GeV) compared to downward spectra. Particle fluxes of protons and 4He 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.

4 Summary and Conclusion

We have applied the novel tool AtRIS (Banjac et al., ) 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 4He ions. We first visualized the atmospheric response matrix of each primary-secondary chain via a two-dimensional 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 image-differencing matrices (IDM). 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 3He 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 $\sim$ GeV and 20 GeV where the Liège 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 atmospheric response matrices with BON GCR spectra of protons and 4He particles to obtain the secondary particle spectra of Hydrogen and Helium isotopes: protium, deuterium, tritium, 3He and 4He particles. This process has summed up the same secondary type generated by both proton and 4He 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 3He flux and a better match with the RAD data. This is due to the same reason we have already observed in the atmospheric response matrices.

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 FTFP_INCLXX_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.

Acknowledgements.

Part of this work was supported by the German Deutsche Forschungs Gemeinschaft, DFG project 282759267. The work is also supported by DLR and DLR’s Space Administration grant numbers 50QM0501, 50QM1201 and 50QM1701 to the Christian Albrechts University, Kiel. We thank Bent Ehresmann for the helpful discussions and support. J.G. and R.F.W.S. acknowledge the International Space Science Institute, which enabled part of the collaborations in this paper through the ISSI International Team 353 “Radiation Interactions at Planetary Bodies”. The authors and the editor thanks Rami Vainio and an anonymous referee for their assistance in evaluating this paper.