On‐the‐Fly Calculation of Absorbed and Equivalent Atmospheric Radiation Dose in A Water Phantom with the Atmospheric Radiation Interaction Simulator (AtRIS)

Computation of the ionization and radiation dose in arbitrary (exo‐) planetary atmospheres due to energetic particles is recently becoming more important due to several reasons that are either correlated with the detection of trace gases for life on exoplanets or with computing dose rates at arbitrary altitudes in the Earth atmosphere. We previously presented Atmospheric Radiation Interaction Simulator, a new Geant4‐based code tailored specifically to enable parametric studies of radiation propagation through exoplanetary atmospheres (Banjac et al., 2019 https://doi.org/10.1029/2018JA026042). Therein, the calculation of ion‐electron pair production rates, which are a mandatory input for chemical and atmospheric modeling, has been presented and validated against Earth measurements and also other, similar, but solar‐system‐specific Geant4‐based codes (PLANETOCOSMICS). In addition to providing input for atmospheric modeling of exoplanets, with AtRIS we aim to directly characterize the habitability by calculating the absorbed dose. In this technical validation study, after showing a detailed analysis of the secondary particles contributing to the atmospheric radiation, we describe a feature of the code which makes direct parametric studies of the interrelation of incident radiation and the resulting absorbed dose throughout the atmosphere possible. In a validation case study configured using an atmospheric model obtained with NRLMSISE‐00 and a primary proton and helium GCR flux calculated using a recent improvement of the force‐field approach, we have compared simulation results with measurements obtained with the Flight Radiation Environment Detector (FRED). We show that Atmospheric Radiation Interaction Simulator (AtRIS) can reproduce the measured dose rate dependence on altitude.


Introduction
When energetic particles pass through matter, they interact in several ways. While the most dominant interaction mechanism is of electromagnetic nature, especially at higher energies, nuclear interactions become more significant. The effects range from harmless excitation and elastic scattering to DNA-wrecking ionization and nuclear spallation. In addition, radiation influences indirectly the conditions of a biosphere through a range of chemical and climate-related processes. From a biological point of view, these can be both harmful and beneficial. While the harmfulness of radiation to biological systems is known, little is known about the ongoing discussions related to the role of radiation as a mutation mediator and its evolutionary consequences. A good review of the status quo has been provided by Ari and Melot (2014).
Atmospheric Radiation Interaction Simulator (AtRIS) was developed to support the ongoing interdisciplinary action to investigate the habitability of exoplanets; however, it is ideally suited to enable parametric studies of the absorbed and dose equivalent at different positions in the Earth atmosphere. Many aspects of the code have been previously described and validated . To avoid duplication, we limit the following exposition to new results. The reader should consider the referred publication for additional information. In AtRIS the particle transport is based on a Monte-Carlo approach utilizing GEANT4 (Agostinelli et al., 2003;Allison et al., 2006Allison et al., , 2016. Internally, AtRIS calculates an Atmospheric Response Matrix (ARM) of the planet, which describes how different primaries distribute their energy throughout the planetary model. For the on-the-fly calculation of dose rates, we will employ an approach involving precalculated look-up tables describing the response of a specific phantom (see below). While providing several reference physics lists, GEANT4 also gives considerable flexibility in defining the physics of the simulation.
In order to achieve its objectives, that is, to calculate the atmospheric ion-electron pair production rate as well as the absorbed and equivalent radiation dose, the primaries and the resulting secondaries are tracked throughout the planetary model. The simulation code takes care to retain information about all energy losses that these particles suffer. Ultimately, one gets a response matrix of the atmosphere, which describes how different primary particles distribute their energy in the form of ionization throughout the atmosphere (see Banjac et al., 2019, Figure 2 therein). Still, in order to quantitatively compare the radiological conditions of different atmospheres, we need to investigate all the individual contributions to the net ionization.
One of the outputs of AtRIS is to provide detailed information about secondary particles crossing the interfaces between sensitive volumes (which are typically concentric shells). In principle, this information is sufficient to completely describe the radiation field, at least at that particular set of discrete altitudes. To calculate dosimetric quantities, one would have to determine the energy and angular distributions of all secondary particles and use this as input for a separate simulation in which a phantom is irradiated. This approach, however, has many disadvantages such as, for example, (i) the AtRIS simulation (mode 0, see Banjac et al., 2019, Appendix A) that provides the necessary information about secondary particles could result in several terabytes of data, (ii) the post processing would include many error-prone steps, and (iii) for each atmospheric shell interface, a separate simulation would need to be performed. Several studies based on existing particle transport codes similar to AtRIS have successfully calculated the radiation conditions on Mars and Venus (Dartnell et al., 2015;Guo et al., 2018;Matthiä et al., 2016). These are often either not publicly available, or are only applicable to a specific planet, or calculate the dose in air instead of the dose in a meaningful phantom. Finally, a change of input parameters (e.g., the primary particle spectrum) would typically make a new simulation necessary. With AtRIS, however, once we have the Absorbed Dose Atmospheric Response Matrix (aARM), this is not necessary, as will be shown below in section 4, where our model is described. Before that, in sections 2 and 3 we will investigate the atmospheric secondary particles and describe how the phantom look-up tables are calculated. Finally, in section 5 we will evaluate the performance of AtRIS by comparing simulation results and measurements.

Choosing A Reference Phantom for Astrobiology
The main goal of AtRIS is the calculation of the ion-electron production rate, which is typically expressed in ions·cm −3 ·s −1 , in order to enable atmospheric chemistry modeling. Ion-electron pair production rate is calculated by dividing the ionization rate (expressed in eV·cm −3 ·s −1 ) with the so called w-value, which represents the average energy required to ionize an atmospheric molecule/atom, with 35 eV being typically used for Earth's atmosphere (Rees, 1963). From the ionization rate we can quickly calculate the absorbed dose rate (expressed in J·kg −1 , i.e., in Grays) in air by dividing with the ambient air density. Hereby, air has a general meaning and does not necessarily correspond to a composition typically found in Earth's atmosphere. Examples of use of the described procedure can be found throughout the literature (see, e.g., Dartnell et al., 2015). Still, knowing the absorbed or even dose equivalent in air is not sufficient when assessing the radiation impact on biological systems. This becomes obvious immediately when the dependence of various cross sections on target material composition and density is considered.
A proper approach is to calculate the absorbed dose and dose equivalent rates that a given, standardized phantom, whose composition better reflects biological systems, would suffer. For this purpose, the International Commission on Radiation Units (ICRU) has defined the ICRU sphere: "A sphere of 30 cm diameter made of tissue equivalent material with a density of 1 g/cm 3 and a mass composition of 76.2% oxygen, 11.1% carbon, 10.1% hydrogen, and 2.6% nitrogen. Used as a reference phantom in defining dose equivalent quantities" (Quantities & Units, 1980). The composition of the ICRU phantom has been chosen to reflect life as we know it. Still, it is doubtful that the use of Tissue Equivalent Material is meaningful when it comes to investigating the radiation conditions of an exoplanetary atmosphere. Furthermore, in recent times more human-like phantoms have become popular. However, it is also doubtful whether such complexity is necessary to study exoplanets. Therefore, we propose a more general phantom with equal dimensions composed of liquid water. Thus, we distance our investigations from the radiation protection side in order to investigate the atmospheric radiation environment of exoplanets from a radiological point of view. It pays to mention

Secondary Particles and Their Primary Parent Particles
Galactic Cosmic Rays (GCRs) and Stellar Cosmic Rays are the driving force behind the atmospheric ionization, causing particle cascades, which are also known as a cosmic ray showers. A Cosmic Ray (CR) will typically indirectly produce millions or billions of secondary particles at various altitudes. Together with the primaries, these secondary particles constitute the radiation field, which is directly causing the ionization of atmospheric matter. Thus, in a first step we need to investigate the particle flux in order to determine what particles contribute to Earth's atmospheric radiation. For that, we have performed a simulation with an atmosphere obtained using the NRLMSISE-00 model (Picone et al., 2002). The primary spectrum was composed of protons uniformly sampled from the energy range between 100 MeV and 100 GeV. This case study is not related to a particular time period or a location but rather qualitatively describes the relation between primary particles and the resulting secondary particles. Thus, a uniform, that is, flat primary proton spectrum covering the most significant GCR proton energy range from 100 MeV to 100 GeV is used. Figure 1 displays a unique double 2-D histogram in which, 1. Each column corresponds to a certain particle type, as indicated. 2. Each row corresponds to an energy bin of the primary flux. The energy bin edges are indicated as ticks on the left y-axis. 3. Each orange row has been normalized to 1 and shows how primary protons from the corresponding energy bin produce different types of secondaries throughout the atmosphere. Thus, we see that protons between 100 and 125 MeV produce mostly -rays, neutrons, protons, electrons and positrons. In addition, some muon anti neutrinos, deuteron and -particles, as well as relatively few triton, helium-3, oxygen and nitrogen ions and are produced, Thus if one wants to see what kind of particles protons between 100 MeV and 125 MeV produce in the atmosphere, one needs to consider only the bottom orange row. 4. Each blue column shows the normalized polar angle distribution of the corresponding particle type throughout the atmosphere (for all primary energies), without compensating for the angular sensitivity of the detectors (spherical shells). The polar angle is measured relative to the azimuth and angles between 90 • and 180 • correspond to particles moving toward lower altitudes, that is, downward. Thus, if one wants to see the angular distribution of all -rays, one needs to consider the first blue column from the left. 5. Histogram fields enclosed in red rectangles have a low count number. A dashed edge indicates a count number below 100, and a full edge corresponds to a count number below 30. 6. The color bars use a logarithmic scale to indicate the relative distributions. 7. The top panel shows in purple both the absolute (left y-axis) and the relative (right y-axis) production of different kinds of secondaries, across the whole energy range of the primary protons sampled from a uniform distribution. 8. In the top panel we also show with green points how the relative production of different kinds of secondaries is changed when the primary protons are sampled from a realistic GCR spectrum corresponding to a solar modulation potential = 600 MV (for more details see section 5.1 and Figure 6). The calculation of these relative abundances has been derived using data shown in Figure 2, and it does not include scarse particles due to insufficient statistics for the convolution procedure (see section 4). Accross the whole atmosphere, most abundant are -rays (69.8%), followed by neutrons (17.4%), protons (8.7%, including primary protons), electrons (2.7%), and positrons (1.1%).
Several conclusions about the secondary particle flux can be drawn from the data shown in Figure 1. For example, we see that a primary proton needs to have a kinetic energy higher than 3 GeV in order to generate a kaon and also that kaons are exclusively directed downwards. As can be seen the main components of cosmic ray showers are gammas, neutrons, electrons, and protons. A significant amount of different kinds of neutrinos is produced; however, these do not contribute to the radiation dose. Also noticeable in the case of pions, kaons, and muons is the fact that in order for a given secondary to be produced, the primary energy of the primary proton must exceed the energy threshold. Furthermore, these particles that result from head-on collisions have a dominantly downward momentum. The angular distributions are nearly isotropic or have a dominant downward direction with a considerable spread. Still, no particle species has a focused angular distribution. Finally, hydrogen and helium ions are relatively abundant, but the higher Z ions are observable only in very low numbers corresponding to a relative abundance of less then 1 part in a million.
Particles that had a relative abundance of less then 1 in 10 8 were not shown in Figure 1, since they cannot significantly affect the calculated dose rate. Practically, the dose rate calculation can be even further limited to the particle species indicated in the top panel of Figure 1 with an abundance higher than one particle in 10 6 (this threshold is indicated with a red strip), less the neutrinos, which are weakly interacting with matter and do not contribute to the dose rate. By considering other related model uncertainties and the expected resolution of transit spectroscopy, we conclude that this threshold value is sufficient even for particles that The projection to the x-axis( i.e., energy axis), shown in red, has as dimension flux·energy. While they-axis projection is flux· sheet thickness and thus corresponds to the particle number density as a function of altitude.
have a higher biological effectiveness. This leaves out the particles (rays, in the case of ) that are significant for the calculation of dose rate: , e ± , ± , ± , n,n, p ± , d, t, He 3 and particles.
Additional insight about the factors determining the atmospheric absorbed dose can be gained by analyzing the altitude dependence and the energy distribution (spectrum) of the relevant types of secondary particles. However, in order for such an investigation to be meaningful, we have switched from a uniform to a more realistic primary particle spectrum consisting of primary GCR proton spectrum calculated with the Force-Field approximation (Usoskin et al., 2011) with a modulation potential = 600 MV and a cutoff rigidity R C = 0.45 GV (Shea & Smart, 2001) and extending to 1 TeV. Since this is a qualitative analysis, these parameters have been chosen at random, and it is not necessary to consider primary GCR helium nuclei. An example proton GCR spectrum calculated with a slightly different modulation parameter ( = 623 MV) can be seen in Figure 6. For each type of secondary particle, a 2-D histogram, with the kinetic energy on the x axis and the altitude on the y axis, is shown in Figure 2. The particle species is indicated with magenta in the top left corner of each panel. Therein, we can better identify the particle species that contribute the most to the altitude-dependent absorbed dose rate. As can be seen, the surface flux is mostly composed of electrons, muons, neutrons, and gammas. We also observe that d, t, He 3 and particles do not penetrate down to the surface. Figure 2 also includes projections of the 2-D histograms to the x-axis and to the y-axis. The x-axis projection corresponds to flux energy across the whole atmosphere (indicated with red and with a scale on the right hand side), while the y-axis projection corresponds to an altitude profile of the particle flux (indicated with black and with a scale on top of the figure).

Ionization Efficiency and the Precalculated Phantom Look-up Table
We proceed to quantify how efficient particles are at ionizing our phantom. For each particle species, we have simulated billions of particles with energies between 1 eV and 10 TeV. The energy range has been split in logarithmically equidistant bins using 50 bins per decade. We then proceeded to calculate the ratio between deposited energy E d and the particle's kinetic energy E i . This ratio describes the relative ionization efficiency, where T designates the particle species (type). Results are shown in the two top panels of Figure 3. The left panel contains stable particles. A relative ionization efficiency of 1, indicated by the green dotted line, reflects the case that the primary particle and all its secondaries have been stopped within the phantom. As can be seen, rays are the most penetrating ones. The critical stopping energy, below which most particles Note. In some cases ( + , + ), the resulting secondaries are again unstable. Energy distribution in a three-body and two-body decay has been properly taken into account. The most dominant process is indicated, except for p-p and n-n annihilation. Energy carried away by neutrinos is considered to be lost.
are stopped, is smallest for gammas (around 50 keV). Electrons are the second most penetrating particles, while heavier particles, like alpha particles, are least penetrating due to the high stopping power.
The right panel shows the same quantity for unstable particles and antiparticles. It is obvious that at lower energies the relative ionization efficiency exceeds 1 that the deposited energy exceeds the kinetic energy of the particle. This is due to the fact that through the decay or the annihilation process, additional energy is released inside the phantom. A further perspective on this phenomena can be obtained by perceiving the absolute ionization efficiency, that is, the average energy deposited as ionization, E d , which is shown in the bottom panels of Figure 3. In the case of stable particles (bottom left), the line with slope 1, indicated by the green dotted line, corresponds to particles, which are completely absorbed. The stable particles approach this line from higher towards lower energies. In the case of unstable and antiparticles (right hand side), the reference line is crossed, and the corresponding absolute ionization values show a plateau-like leveling on the left side of the reference line. Thereby, the height of the plateau is determined by the average energy deposited in the phantom following a decay or annihilation process.
Still, if we consider the involved decay processes, we see that this average energy deposited following a decay or annihilation is somewhat less than the energy being released through the process itself (see Table 1). When a positron annihilates with an electron, two 0.511-MeV gammas are created, which have a relatively high probability to pass through the phantom. Therefore, only ≈ 18.8% of the theoretical maximum available energy, which is calculated by solving an energy balance equation, gets deposited on average within the phantom. The annihilation of antiprotons and antineutrons is much more complicated and typically results in various mesons, which consequently produce gammas, electrons, positrons, and neutrinos. Energy carried away by neutrinos is completely lost. Also, depending on the decay type (two-or three-body), the energy of a neutrino can be as high as 50% of the available net energy. Also, + and K + , whose main decay product is a high energy +, are extremely inefficient in depositing energy following a decay inside the phantom. This has to do with the fact that the resulting muon interacts with electrons, which are much lighter. The antimuon moves without experiencing much directional changes and thus exits the phantom following an almost straight line. The muon however, interacts with the positive protons in the nucleus. Since protons are heavier, the muon is more likely to change its direction resulting in an increase in average path length and a corresponding decrease of the effective range within the phantom. While protons and antineutrons have the greatest absolute contribution, their relative ionization effectiveness is lowest, since the resulting secondaries are highly energetic and are more likely to leave the phantom.

10.1029/2019JA026622
By exploiting the spherical symmetry of the ICRU phantom making it possible to define the phantom in the same manner that a planet is specified by using the Planet Specification Format (PSF) file format, which is a custom format developed specifically to enable the quick and flexible specification of (exo) planets within AtRIS (see section 2.2 and Figure 1 in Banjac et al.,2019, and the following section, for more details). We calculated with AtRIS (mode -2, see Banjac et al., 2019, Appendix A) the relative ionization efficiency of each significant particle. This information is stored in form of an external asci look-up table. Knowing the ionization efficiency enables a quicker determination of the absorbed dose rate. Let us consider for example a 1-GeV proton incident upon the phantom. Instead of simulating all the interactions that such a proton and all of its secondaries would suffer in the phantom, we simply look up the corresponding ionization efficiencies ( Figure 3) and find that, on average, for a 1-GeV proton ≈ 6.395% of its kinetic energy would get deposited in form of ionization within the phantom. An additional significant advantage of this approach is the increased statistical significance. If we were to simulate the interaction of any particular particle with the phantom, we would waste time without gaining statistically significant results. In our approach, the number of particles per energy bin (100 million) was sufficient to account even for very rare processes, as is confirmed by the smoothness of the graphs that have been shown in Figure 3. Since the phantom look-up tables have not been hard-coded, but are rather processed at the start of a simulation from an external text file, AtRIS can be configured to work with different phantoms. We provide several phantoms for the user to choose from and also instructions on how new phantoms can be generated. Details on how to access the precalculated look-up tables describing the phantom response is given in the acknowledgments. Note that in this paragraph we wanted to point out the intricacy of the interaction of unstable particles and antiparticles with the phantom material. The reader should keep in mind that the relative contribution of such particles at energies where the decay and annihilation processes are probable is still not significant. This has to do with the relative scarcity of such particles and with the scale of the involved energies. For example, while a 10-keV neutron will deposit approx 128 keV to the phantom through decay and the consequent interaction of the resulting secondaries with the phantom material, a 7-MeV neutron will deposit 7 MeV on average, which corresponds to a factor of ≈ 55.

AtRIS: Model Description
In literature, we can find particle transport codes implementing several approaches to model the dose rates in atmosphere: for example, PHITS (Sato et al., 2013) and HZETRN/OLTARIS (Slaba et al., 2010). A typical approach to this is to count the secondary particles, recreate their spectra, and use this information to set up a separate simulation around the phantom, or alternatively, to use precalculated flux-to-dose functions.
The approach that we have adopted with AtRIS can be described as a self-consistent automation of this procedure.
To facilitate the quick specification of planet models, we have developed the PSF, a custom interface providing the necessary flexibility at a fraction of the complexity of a Geometry Description Markup Language (GDML) file. A provided script can convert the PSF file to a GDML file, which is Geant4 compatible. This approach reduces the complexity of the internal GDML interface while providing all the flexibility that might be necessary for investigations involving (exo-) planets. In addition to a macrofile configuring the simulation, the user is required to provide a .bins file specifying the energy binning that shall be used. Instead of sampling primary particles from a given flux, AtRIS simulations are performed with a logarithmically uniform energy spectrum. The main output of AtRIS is the so-called ARMs that describe the mapping between primary particle energy, altitude and quantities like ionization, absorbed dose rate, and dose rate equivalent. Using provided post processing scripts, a convolution with a realistic primary flux can be performed afterwards. Such an approach was necessary in order to make AtRIS suitable for parameter studies. Once the ARMs have been calculated, there is no need to repeat the calculation if the input spectrum changes. In the context of astrobiology, the benefit of this approach is obvious if we consider how little is known about the radiation environment of rocky exoplanets around M-dwarfs and K-stars . For a more detailed description of AtRIS the reader is referred to Banjac et al. (2019).
To describe the implementation and how exactly the absorbed and dose equivalent ARMs are calculated, let us assume that a primary proton has been generated with a starting (vertex) energy E 0 . Let j be the index, which connects the E 0 to the user-provided energy binning, that is, E 0 ∈ [E j , E j+1 ] with j = 0, 1, … . Furthermore, we index the interfaces between sensitive volumes with i, whereby i = 0 corresponds to the volume between the core and the sensitive crust, i = 1 corresponds to the interface between the crust and the

Journal of Geophysical Research: Space Physics
10.1029/2019JA026622 first sensitive air volume and so on. The index i is called the Sensitive Detector ID (SDID). As AtRIS simulates the propagation of this proton through the atmosphere, at each volume interface, and for all significant secondary particles, the following steps are performed: 1. Identify particle type T and its kinetic energy E i . 2. Retrieve the ionization efficiency  R (T, E i ) of the penetrating particle from the internal look-up table.
3. Calculate the average absorbed dose rate that the phantom would suffer if the given particle with its kinetic energy would pass through the phantom: 4. In an internal 2-D matrix, representing the aARM, add the so calculatedD to the existing value of the entry (i, j), whereby i corresponds to the SDID, while j describes the energy bin of the primary proton. 5. For particles with a radiation weighting factor W R larger than 1, determine the factor. For neutrons, the factor is calculated using the kinetic energy E i (see, e.g., ICRP 60,1991 ICRP 103, 2007). 6. Calculate the average dose rate equivalent that the phantom would suffer if the given particle with its kinetic energy would pass through the phantom according to (see, e.g., ICRP 60, 1991): 7. In an another internal 2-D matrix, representing the Dose Equivalent Atmospheric Response Matrix (eARM), add the so calculatedH to the existing value of the entry (i, j), whereby i corresponds to the SDID, while j describes the vertex energy bin.
The ionization efficiency is a precalculated step function, which describes the relative kinetic energy that the particle would on average deposit within the phantom. Since the relative ionization efficiency has been calculated discreetly with 650 bins spanning the energy range between 1 eV and 10 TeV, the second and third steps represent an interpolation procedure that is meant to minimize the discretization error. Please note that significant secondary particles are those that have previously been identified as being a significant factor for the determination of atmospheric radiation (see Figure 1).
The resulting ARMs have the dimension of energy instead of the needed Grays per particle (or in the case of dose equivalent, Sieverts per particle). AtRIS performs several normalization steps in order for each entry of the aARM and eARM to describe the effect that a single primary particle of that particular starting kinetic energy has at a particular interface designated with the SDID. These are 1. Divide each column with index j by the number of primary particles generated in the j-th energy bin. 2. Multiply each row with a factor: where A ICRU = 4 · (15cm) 2 is the sensitive area of the phantom, while A SDID = 2 · 4 R 2 SDID is the surface of the atmospheric interface corresponding to that particular SDID, with the added factor 2, due to the fact that particles can penetrate the boundary from both sides. 3. Finally, we divide all matrix entries with 14.1372 kg (the mass of our r = 30 cm water phantom that has been introduced in section 1.1).
The normalization described in the second step accounts for the smaller surface of the phantom compared to the available surface of the interface between two atmospheric sheets. The same approach can also be followed with a nonspherical phantom by replacing the corresponding surface term of the phantom in Equation (4). Whether the secondary particle flux is isotropic or not, is of no consequence, as long as the atmospheric boundary interface and the phantom have similar orientations.
After the completion of the described normalization steps, the only task remaining is to scale the results so that they correspond to a real flux and to add the contributions of the different vertex energies together. The procedure is illustrated in Figure 4 using a primary GCR flux composed only of particles calculated with the force-field approach (Usoskin et al., 2011) while using a modulation parameter of 629 MV and a cutoff rigidity of 0.45 GV. Please note that since this is a qualitative study, we do not need to include photons. particles are used here as an example. Let us consider a simulation performed with K energy bins with an that describes how − particles contribute as a function of primary energy to the absorbed dose across the whole altitude range from ground to 100 km. By multiplying [ M A ] with the number of primary −particles per energy bin [ n ] that is displayed in the upper left panel, we find that primary − particles with energies between ≈ 1 and ≈ 10 GeV are the main contributors to atmospheric absorbed dose rate (lower right panel). By projecting to the y-axis, we get the altitude absorbed dose rate height profile (upper right panel).
atmosphere composed of N sensitive shells. Let the real number of particles in the j − th energy bin be n j , where j ∈ [0, … , K − 1]. The vector [ n ] is shown in the upper left panel with a red line. It has been obtained by calculating a third degree spline representation of the differential flux j, and then integrating the spline while using the bin edges as integration limits. We designate the aARM, of dimension K×N, with M A (shown in the bottom left panel). We can now formulate the described procedure mathematically: where D i is an N-dimensional vector describing for the whole range of altitudes the absorbed dose, as shown with the orange line in the upper right panel. Alternatively, we can multiply each column of with the corresponding number of primary particles n j to obtain the Matrix [ D A ] that is shown in the lower right panel. Note that the upper right panel showing the altitude profile of the absorbed dose rate represents the final step, that is, the result of the calculation. It is obtained either by following the procedure described in Equation (5) or by summing the matrix D A accros its columns, that is, if we project the matrix on to the y-axis, we obtain the dose rate altitude profile D i (Equation 5). Although high energy particles produce high dose rates per particle (lower left panel), after scaling with a realistic primary flux (upper lef panel), we see that the primary contributor to atmospheric dose rate for the case of particles comes from the energy range of ≈ 1 to ≈ 10 GeV (lower right panel). This shows that particles are not relevant when it comes to calculating the absorbed dose at altitudes below 5 km.
Having the data available in a matrix form essentially describes how primary particles of different energies distribute their energy and generate secondaries throughout the whole atmosphere (see Figures 2 and 4). If we combine this with the calculated phantom response (Figure 3), we can obtain useful perspectives in to the cause of the radiation hazard throughout the atmosphere. For example in Figure 5, we see how different types of secondary particles relatively contribute to the net absorbed dose altitude profile, which is shown in the bottom right panel of Figure 4. While this is not a new result, it demonstrates the flexibility of AtRIS. At altitudes above 35 km, primary alpha particles are the dominant cause of absorbed dose. Below 35-km protons begin to take precedence, and at 10 km we see a significant electron and positron contribution. The contribution of neutrons is highest at an altitude of ≈ 7 km, while at surface muons are the dominant contributors to absorbed dose. We can clearly observe that antiprotons and antineutrons are not significantly contributing to the atmospheric absorbed dose at any altitude.

Comparison Against Stratospheric Balloon Measurements
On 28 September 2011, at 16:00 UTC, the Flight Radiation Environment Detector (FRED), a particle telescope developed to measure dose rates due to neutral and charged particles at high altitudes, was flown on the stratospheric balloon BEXUS 13. The experiment was performed at the Esrange Space Center in Kiruna, Sweden (67 • 53'38"N, 21 • 06'25"E). Of the total 4-hr flight time, more then 2 hr were spent at a floating altitude of ≈ 25 km, where the corresponding cutoff rigidities varied between 0.45 and 0.6 GV. The experiment also included a GPS sensor, thermometer, and a barometer, so that it was possible to measure the relation between dose rate and pressure as well altitude. Extensive description of the instrument and the experiment and data postprocessing was published by Möller (2013). FRED, which is closely related to DOSTEL (Beaujean et al., 1999), is built as a stack of four rectangular segmented semiconductor detectors with a field of view of ≈ 120 • , which are for dosimetry operated as single counters. Assuming an ideal instrument measuring each particle (its energy deposition and species), using the radiation weighting factors W R (see ICRP 60, 1991), which depend on the Linear Energy Transfer (LET), the dose equivalent could be easily calculated with where D T is the absorbed dose in grays (Gy) by radiation type T. However, since with FRED it is not possible to measure the LET of each particle (and thus determine W R ). Instead, in order to calculate the dose equivalent H, Using coincidence measurements and thanks to the detector segmentation, which has been introduced to limit the uncertainty in the path length of the incident particles, FRED is able to measure the LET spectrum in silicon. From this, Möller (2013) first calculated the LET spectrum in water and used the ICRP 60, 1991 quality-factors function to approximate the dose equivalent (see section 6.2.3 and Figures  6.29-6.32 therein). The absorbed dose measured by Möller (2013) is shown in the top panel of Figure 8 using orange error bars, while the bottom panel contains the calculated dose equivalent (green error bars). In the following, we will describe our simulation setup and compare simulation results to the measured absorbed dose rate and the calculated dose equivalent of Möller (2013).

Primary Spectra
The primary particle flux, here composed of GCR protons and helium nuclei, was modeled based on the Force-Field approximation (see, e.g., Caballero-Lopez & Moraal, 2004;Moraal, 2013). Based on the Local Figure 6. Differential spectra calculated for the Flight Radiation Environment Detector validation case study using the Force-Field approximation, as described in the text. The proton differential spectrum J is shown with orange, while the helium nuclei spectrum with green. To compare with literature it might be necessary to convert the spectrum from ·cm −2 to ·m −2 by multiplying with a factor of 10 4 . Interstellar Spectrum LIS model by Herbst et al. (2017) and by assuming the modified solar modulation parameter approach by Gieseler and Heber (2017), a solar modulation parameter of = 623 MV for the time of the flight was derived. The simulations were performed using logarithmically uniform energy bins covering the energy range between 100 MeV (corresponding to a cutoff rigidity of 0.45 GV at Kiruna) and 1 TeV, with 20 bins per decade. The calculated differential spectra are shown in Figure 6.

Atmosphere Definition
Based on the NRLMSISE-00 atmospheric model by Picone et al. (2002) we configured our atmosphere as a series of 130 concentric spherical shells. Location and time used as input for the NRLMSISE-00 model were chosen to match the measurement flight. For the altitude range between 0 and 30 km we have used 0.5-km-thick shells in order to have finer altitude steps over the altitude range corresponding to the balloon flight. Between 30 and 100 km, 1-km-thick shells were used. This custom atmospheric model was obtained by combining the output of two NRLMSISE-00 model runs. The model was set to produce an atmosphere composed out of molecular nitrogen and oxygen, and atomic helium, argon, hydrogen, nitrogen, and anomalous oxygen.

Compensating for Phantom/Detector Differences
While the 30-cm water phantom developed for AtRIS is a good match for the calculation of the radiation hazard in the context of astrobiology, it cannot be directly validated since no analogous measurements exist. Typical dosimetric experiments are performed using detector stacks most often consisting out of semiconductor discs or slabs, as is the case with FRED, with a typical thickness that is less than 1 mm; that is. measurements, which are available for validation, have been performed using instruments that differ in size, shape, and composition from our phantom.
In order to facilitate the validation of the procedures implemented in AtRIS, we have resolved to the development of a custom phantom mimicking FRED. The spherical shell phantom was replaced by a geometry corresponding to FRED's detector stack consisting of four square semiconductor detectors (side 22 mm, thickness 300 m) enclosed within a 0.8-mm-thick magnesium casing. The calculation of the phantom response, which is shown in Figure 7, was identical to the procedure that has been already used and described in section 1.1. Compared to the (ICRU) phantom, the FRED phantom (Figure 7) responds differently to radiation in several aspects: 1. The presence of the 0.8-mm-thick magnesium shielding results either in partial absorption of the particle's kinetic energy or in complete stopping. Here, with partial absorption we understand the case where some energy of the primary particle is carried over to the sensitive detector through the shielding either by the primary particle itself or by some secondary particle or particles, which were created in the interaction of the primary and the magnesium shielding. In other words, the deposited energy in the phantom is not necessarily a consequence of the primary particle interacting with phantom material, but it could well be due to some secondary managing to tunnel through. Also note that the shielding is particularly effective against heavy particles like He 3 nuclei alpha particles. 2. Particles are more likely to penetrate the FRED phantom. For example., whereas the (ICRU) phantom could completely absorb ≈ 700 MeV -particle, the FRED phantom suffers on average only ≈ 5.6 MeV in the form of ionization from such a particle. This causes the (ICRU) sphere to suffer significant absorbed dose from energetic particles, which would otherwise penetrate the thin semiconductor detectors of FRED without depositing much energy in form of ionization. 3. The decaying and annihilating particles deposit significantly less energy. For example, a moun that would have deposited on average ≈10.95 MeV within the ICRU phantom deposits only ≈0.45 MeV within the shielded FRED phantom. 4. Our ICRU phantom is completely unshielded and thus lacks the protection provided by FRED's magnesium shielding against lower-energy radiation. On the other hand, the absolute contribution to absorbed  Figure 3). The shielding stops low energy particles and absorbs a portion of the energy of particles that manage to penetrate to the sensitive detector matter. The much thinner geometry of Flight Radiation Environment Detector phantom manifests itself in lower absolute ionization at higher energies across all particle species. In particular, particles that previously deposited significant amounts of energy through decay and annihilation processes are now much less significant.
dose of such radiation in the case of the (ICRU) sphere is low (see Figure 3). Even more so, if we consider the large mass of the sphere. Figure 8 shows the calculated absorbed (top panel) and dose equivalent rate (bottom panel) using the new FRED-specific phantom (purple). Results obtained with the (ICRU) phantom are shown in green. The use of the custom FRED phantom delivers good agreement between simulation and measurement, in particular when comparing simulated and measured absorbed dose. The simulated dose equivalent using the FRED phantom is somewhat lower than the results published by Möller (2013). However, we need to keep in mind that this quantity was not directly measured but rather indirectly calculated.

Analysis
The lower-then-expected dose equivalent and the small discrepancy at high altitudes between simulated absorbed dose and several measured data points could be explained with an underestimated primary helium flux and the fact that we did not include primary Z>2 GCR nuclei. The force-field approach that we are using to generate the primary spectra suggests a helium-to-proton ratio of 0.05 in LIS (Usoskin et al., 2011), which is derived from recent measurements of the ratio at energies above 10 GeV/nuc −1 that yielded 0.05±0.002 as a result. It is uncertain whether extrapolating this factor to the whole energy range remains a valid procedure. Furthermore, due to limitations of the force-field model, which is used to generate the primary spectra and is limited to protons and helium nuclei, we were not able to directly investigate our Z > 2 hypothesis. Yet, in order to get some insight in to these effects, we have also calculated the absorbed dose and dose equivalent due to a primary helium flux that has been increased by 50%. While the obtained results, which are shown in Figure 8 with red dashed and dotted lines, in some aspects better fit the data (e.g., dose equivalent), they raise other concerns. For example, the agreement between simulation and measurement in the case of the In the case of absorbed dose, the agreement between simulation and measurements is outstanding. Except for the uppermost data point at ≈25 km, the simulation is within the measurement uncertainties. The agreement in dose equivalent (bottom panel) is somewhat worse. The results obtained using the (ICRU) phantom, shown in green, are typically higher by a factor of 2. The red plots correspond to the hypothetical case with an increase primary helium flux, as is described in the text.
absorbed dose at ≈ 15km worsens. This, however, might have to do with the fact that the stopping power of the atmosphere for Z>2 nuclei is higher due to their higher charge. Thus Z > 2 GCR nuclei would be more likely to deposit their energy in the higher regions of the atmosphere.
While some of the differences in phantom response (Figure 7) were rather striking, the net effect is not as dramatic. The dose rates seen by the (ICRU) phantom are at altitudes above 6 km higher by a factor of ≈2. However, part of this is due to the fact that the (ICRU) phantom is composed of water, whereas the FRED phantom is silicon. Typically, on earthbound dosimetric experiments, a factor of 1.2 is used to convert absorbed dose rates from silicon to water (see Beaujean, 2005, and references therein). The remaining 66% difference in seen absorbed dose rate is due to differences in phantom size and shape. Due to limitations of the force-field model, which is used to generate the primary spectra and is limited to protons and helium nuclei, we were not able to directly investigate our hypothesis. Furthermore, our unshielded (ICRU) also absorbs low energy particles, which are in the case of FRED stopped by the magnesium shielding.
Note that this validation case study would not have been possible if AtRIS did not provide the user the possibility to provide an external file describing the phantom. This particular feature could be of benefit to instrument builders since they can calculate the exact response of their instrument, provide it as a file to AtRIS, and estimate the absorbed and dose equivalent rate. Information on how to access the FRED look-up tables is provided in the acknowledgments.

Limitations of AtRIS
In its current version, AtRIS has two limitations that need to be mentioned. First, AtRIS cannot be used to model the absorbed dose and dose equivalent under complicated shielding scenarios. As a part of our validation effort, we have tried to reproduce the absorbed dose measured on commercial aircraft traveling across a wide range of geomagnetic cutoff rigidities (from India to Frankfurt) using a DOSTEL-like dosimeter. In several experiments, the dosimeter was placed in a case behind the last passenger seat. Such placement had for us an unfortunate consequence: It made for a highly nonuniform shielding environment, which made it impossible for us to develop a consistent phantom for AtRIS. Depending on the thickness, and the vertical cutoff rigidity, and the thickness of the shielding, we could observe either an increase or a decrease in the absorbed dose, compared to the unshielded case. Still, even so we obtained results that were within 30% of the measurements.
This validation problem is a manifestation of the fact that the radiation field is strongly modulated by the presence of the aircraft shielding. With AtRIS we are characterizing the field, while typical experiments placed dosimeters into scenarios corresponding to radiation protection conditions, and have thus included the effects of various kinds of shielding. This has made our search for adequate measurements, against which to validate our model, extremely difficult. In fact, the only other validation opportunity that we found are the measurements of the Mars Science Laboratory/Radiation Assessment Detector on the Martian surface. However, there are several considerable obstacles, which are preventing us to pursue this within the scope of this publication.
An additional limitation of AtRIS is the adopted approach for calculating the dose equivalent using a fixed radiation weighting factor of 20 for particles heavier than protons and 2 for protons ICRP 60, 1991. The so assigned factor is adequate only for low-energy nuclear fragments. However, at high altitudes the relativistic Z > 1 GCR nuclei are not as scarce, and thus their biological impact gets greatly overestimated. In other words, AtRIS can reliably calculate the dose equivalent only at altitudes corresponding to an atmospheric depths of 20 g/cm −2 and more. The remove this limitation, in the future we plan to implement a phantom LET response in to AtRIS and thus make the direct calculation of the radiation weighting factors possible.

Discussion, Summary, and Outlook
In this paper, we have shown how AtRIS, a new Geant4-based atmospheric particle transport code, that has previously been described and validated , can be used to quickly and reliably calculate the absorbed dose rate and the dose rate equivalent, while primarly targeting the calculation of the atmospheric ion-pair production rate for the purpose of atmospheric climate and chemistry modeling. In this approach, a response of a standardized phantom is precalculated and used to compute the absorbed and dose equivalent on-the-fly without prolonging the simulation by more then 5% and without requiring any additional steps. Furthermore, the user can, by replacing the file describing the response of the phantom, set up custom phantoms. Thus, AtRIS uniquely combines several output capabilities into a single package: the calculation of atmospheric ion-electron pair production rates, the calculation of secondary particle spectra, and the direct calculation of absorbed and dose equivalent rate. Furthermore, thanks to the general nature of the output of AtRIS, the capability to perform parametric investigations is inherent. Thanks to this feature we are able to quickly recalculate, without repeating the simulation, the dose rates for different vertical cutoff rigidities, and different solar modulation parameter values. This opens the door for gaining further insight into the effects that the primary spectrum shape, a planetary magnetic field or a dense astrosphere might have on the habitability of an exoplanet.
Our validation case study supports the conclusion that AtRIS is able to calculate the absorbed dose rate and dose rate equivalent. When using the default (ICRU) water phantom with a radius r = 15 cm, relative to existing measurements, AtRIS greatly overestimates both the absorbed dose and the dose equivalent. We have shown convincing evidence that this effect is due to the considerably larger sensitive mass and different material of the water spherical phantom. We were able to implement a custom phantom approximating the instruments. The observed disagreements are rather consequences of our choice of reference phantom than systematic errors. Furthermore, while codes that concentrate solely on the calculation of the atmospheric radiation exist, AtRIS represents the first successful implementation of an atmospheric particle transport code, which is, in addition to electron-ion production rate and secondary particle spectra, able to determine dosimetric quantities relevant to life without extensive postprocessing. Previous solutions were able to calculate the absorbed dose rate in ambient air or followed an approach based on using precalculated flux-to-dose functions, which were to be applied in a separate step (see, e.g., Guo et al., 2018).
The techniques developed and presented hold the potential to studying the habitability and radiation environment of not only exoplanets but also solar system planets and planetary moons for which the authors provide as a convenience to the user interfaces to existing atmospheric models (NRLMSISE-00; Picone et al., 2002;MCD,;Forget et al., 1999;Millour et al., 2015;HAMMONIA;Schmidt et al., 2006) and a model of the Venusian atmosphere (Keating et al., 1985;Seiff et al., 1985). For the specification of exoplanets, the authors have developed a flexible PSF. Finally, the community stands to benefit from the openness of AtRIS, which has been published under a GNU GPLv2 licence and for which an open wiki-based documentation is being made available. Access instructions are provided in the Acknowledgments.