Dose calculations at high altitudes and in deep space with GEANT4 using BIC and JQMD models for nucleus–nucleus reactions

Radiation exposure of aircrew is more and more recognized as an occupational hazard. The ionizing environment at standard commercial aircraft flight altitudes consists mainly of secondary particles, of which the neutrons give a major contribution to the dose equivalent. Accurate estimations of neutron spectra in the atmosphere are therefore essential for correct calculations of aircrew doses. Energetic solar particle events (SPE) could also lead to significantly increased dose rates, especially at routes close to the North Pole, e.g. for flights between Europe and USA. It is also well known that the radiation environment encountered by personnel aboard low Earth orbit (LEO) spacecraft or aboard a spacecraft traveling outside the Earth's protective magnetosphere is much harsher compared with that within the atmosphere since the personnel are exposed to radiation from both galactic cosmic rays (GCR) and SPE. The relative contribution to the dose from GCR when traveling outside the Earth's magnetosphere, e.g. to the Moon or Mars, is even greater, and reliable and accurate particle and heavy ion transport codes are essential to calculate the radiation risks for both aircrew and personnel on spacecraft. We have therefore performed calculations of neutron distributions in the atmosphere, total dose equivalents, and quality factors at different depths in a water sphere in an imaginary spacecraft during solar minimum in a geosynchronous orbit. The calculations were performed with the GEANT4 Monte Carlo (MC) code using both the binary cascade (BIC) model, which is part of the standard GEANT4 package, and the JQMD model, which is used in the particle and heavy ion transport code PHITS GEANT4.


Introduction
It is well known that the radiation environment encountered by personnel aboard low Earth orbit (LEO) spacecraft, such as the NASA Space Shuttle and the International Space Station (ISS), or personnel aboard a spacecraft traveling outside the Earth's protective magnetosphere on interplanetary missions is different from that on the ground. These personnel are exposed to radiation from both galactic cosmic rays (GCR) and energetic solar particle events (SPE). In the case of an ISS-type orbit (51.6 • inclination, ∼450 km altitude), it has been estimated that GCR contribute about 50% [1] to the total dose equivalent rate received by astronauts and cosmonauts. The remaining part of the dose equivalent rate originates mainly from protons trapped by the Earth's magnetic field. The relative contribution to the dose from GCR when traveling outside the Earth's magnetosphere, e.g. on Moon or Mars missions, is even greater. The interactions of neutrons are of also of great importance, especially in LEO, and estimates of the neutron contribution to an astronaut's total dose equivalent inside a spacecraft at LEO range from 30 to 60% [1]- [3]. To investigate the radiation-shielding qualities of proposed vehicle models and to simulate the radiation environment, as well as to estimate the biological effects of cosmic radiation, e.g. the organ and tissue equivalent doses in the astronaut's or cosmonaut's organs, so-called particle and heavy ion transport codes are required. These codes have to be validated carefully to make sure that they fulfill preset accuracy criteria, e.g. to be able to predict particle fluences and energy distributions within a certain accuracy. The code chosen and supported by ESA for such studies is the Monte Carlo software package GEANT4 [4], which has been used in Europe and the US for studies of space radiation effects. It has, e.g., recently been used for studies of the Columbus/ISS radiation environment [5]- [10] within the ESA-sponsored project 'Dose Estimation by Simulation of the ISS Radiation Environment' (DESIRE) [11]. A longterm goal of the DESIRE activity has been to develop a European software package that can be used to predict the radiation risks for human space flights, both in LEO as well as to the Moon and Mars. However, the DESIRE project identified problems with existing GEANT4 physics models required for GCR simulation [6]. The problems were concerned with both GCR ions 3 heavier than carbon for all energies as well as all GCR ions with energies above 10 GeV n −1 . This resulted in misleading simulation results where the dose equivalent rates were only about 30% of the experimental values [8].
Since the issue of ICRP Publication 60 [12], in which the exposure of aircrew is recognized as an occupational hazard, radiation protection of crew members has also been widely discussed. One possibility is to calculate particle spectra with the PLANETOCOSMICS Monte Carlo package [14] developed by Desorgher. This package provides a layered geometry of the Earth's atmosphere based on GEANT4 [4], and interactions of the particles within the atmosphere can be simulated with the GEANT4 standard models for electromagnetic interactions, the socalled quark gluon string precompound (QGSP) binary cascade (BIC) HP model for hadronic interactions and the BIC model for inelastic interactions of heavy ions with matter. However, this model is supposed to be applied to light ion reactions below 10 GeV n −1 and is not recommended for ions heavier than carbon. For heavier projectile-target combinations and for energies above 10 GeV n −1 , no validated model exists in GEANT4 yet.
To be able to use the advantages of PLANETOCOSMICS, GEANT4 and the detailed geometry models in Columbus and ISS included in the DESIRE project, etc, and at the same time improve the accuracy in the estimations of ambient dose and dose equivalents, it must be possible to call hadronic and heavy ions reaction models from GEANT4, e.g. JAM [16] and Jaeri quantum molecular dynamics (JQMD) [17], or other transport codes, e.g. the particle and heavy ion transport code system (PHITS) [18] optimized for heavy ion reactions. Koi et al [19] have therefore developed interfaces between GEANT4, and JAM, JQMD and PHITS. Another possibility is to use the BIC model also for ions heavier than carbon and for higher energies, but the results from such a calculation must be benchmarked and validated.
In this paper, it is shown that an interface between GEANT4 and PHITS has only a very minor influence on the accuracy when calculating neutron and proton fluxes at aviation and high altitudes using PLANETOCOSMICS and GEANT4, which is consistent with the fact that only a minor part of the neutrons and the protons in the atmosphere are secondary particles from heavy ion reactions. However, as the altitude increases to LEO and in deep space, where the dose contribution from the heavy charged particles from GCR is important, the need for accurate heavy ion models increases. This is also important when estimating the risk for radiationinduced failures in electronic devices, such as single-event upsets (SEU), multiple-bit upsets (MBU), etc at LEO and in deep space. By interfacing GEANT4 to JQMD or extending the use of BIC beyond its recommendations, better prediction of the radiation risks both for humans and electronic devices can be achieved.

Background
USA and Russia have successfully performed many manned space flights for more than four decades, and when China sent the first Taikonaut to space in 2003, it became the third country in the world to send humans to space. As a part of Europe's strategy for space exploration, endorsed by the European Union Council of Research, the European Space Agency (ESA) set up the Aurora Programme in 2001. The aim of this program is to explore the solar system and the universe, stimulate the development of new technology, and inspire young people of Europe to take a greater interest in science and technology. A second objective is to search for life beyond the Earth. Future missions under the program are supposed to carry sophisticated exobiology payloads to investigate the possibility of life forms existing on other worlds within the solar system. The space environment is, however, complex with both spatial and temporal variations. The composition and energy spectra of the radiation from GCR and SPE are quite specific and must be considered separately. The GCR consist mainly of protons and ions with energies up to several hundreds of GeV, with their peaks ranging from several hundreds of MeV up to around 1 GeV. Even if ions heavier than α particles only comprise around 1% of the GCR's flux, their contribution to the dose is significant, since dose is related to deposited energy, which is proportional to the square of the ion's charge Z . In deep space, more than 50% [1] of the total equivalent dose rate received by astronauts and cosmonauts is induced by GCRs. The highest contribution to the equivalent dose is caused by heavy charged particles and it has been shown that more than 60% of the equivalent dose contribution of the GCR component evaluated at the skin, assuming exposure during solar minimum behind 5 g cm −2 aluminium shielding [20], comes from primary and secondary particles with Z > 2. It is therefore essential to use a particle and heavy ion transport simulation tool that is able to calculate the dose contribution from the heavy charged particles with a satisfying accuracy when estimating the radiation risks for human space flights, both in LEO as well as on interplanetary missions.
The radiation environment on the ground is radically different from that in deep space owing to the protective geomagnetic field and the shielding provided by the atmosphere. The radiation environment in the atmosphere is a result of the interaction of the GCR and SPE with the particles in the atmosphere, and therefore depends on the altitude. The ionizing environment at standard commercial aircraft flight altitudes consists mainly of secondary particles, including protons, neutrons, photons, electrons, muons and pions. Secondary neutrons have the major contribution to the dose equivalent rates and cause an enhanced radiation level compared with that on Earth. Muon fluxes do not change very much with altitude in comparison with other particles because of their strong penetrability in the atmosphere. The ambient radiation dose increases with altitude by approximately 15% for each increase of around 2000 feet (∼600 m), depending on the latitude, so the dose rates at a conventional flight altitude are approximately 100 times higher than that at sea level [21]. As a result of this discussion following ICRP Publication 60 [12], many countries have issued regulations or recommendations to set annual dose limitations for aircrews. To certify that these limitations are followed, the use of an aviation-dose calculation code is indispensable since the doses depend on altitude, vertical cutoff rigidity and solar modulation potential along the flight routes in a complicated manner, and it is impractical to measure the doses for all flight conditions. Several calculation codes, e.g. EPCARD [24], CARI-6 [25] and PCAIRE [26], have already been developed in order to estimate the aircrew dose. These codes can calculate route-doses-the dose en route between two airports-by specifying the flight conditions such as departure and destination airports. The accuracy of the calculation has been well verified by experimental data. However, the codes lack information on dose contribution from SPE and the dose at very high altitudes, where contributions from alpha particles and heavier ions might play a role. It is therefore worthwhile to develop a new route-dose calculation code that explicitly determines the atmospheric cosmicray spectra on flight routes, and combines these data with user-specified fluence-to-dose conversion coefficients. Establishment of a reliable model for calculating the particle spectra under any global conditions is the key issue in the development. One example of such a code is the PHITS-based analytical radiation model in the atmosphere (PARMA) [27].  [28,29], even though CORSIKA and COSMOS are mainly aimed at simulations of air showers.

PLANETOCOSMICS
PLANETOCOSMICS is a Monte Carlo package developed by Desorgher. This package allows computation of the hadronic and electromagnetic interactions of cosmic rays with the Earth, Mars and Mercury environment based on the GEANT4 [4] Monte Carlo code. For each planet it is possible to take into account the presence of the planetary magnetic field, atmosphere and soil. For each planet, different magnetic fields and atmospheric models are available. The main applications of the code are in: 1. the computation of flux of particles resulting from the interaction of cosmic rays with a planet environment at user-defined altitudes, and/or atmospheric depths; 2. the computation of the propagation of charged particles in the planet magnetosphere; 3. the computation of cut-off rigidity (mainly for the Earth) at a given position on a planet and for different directions of incidence; and 4. the visualization of magnetic field lines, and the trajectories of primary and secondary particles in the planet environment.

GEANT4
GEANT4 is a well-known three-dimensional particle and heavy ion Monte Carlo transport code, supported by a large collaboration [4], which uses object-oriented C + + as the programming language. To allow for maximum flexibility, GEANT4 lets the user make a selection of models of particle interactions relevant for the simulation scenario. This physics configuration is done by choosing among users' applications called 'physics lists'. Processes can be assigned to particles, and models to processes, until all relevant interactions for all relevant particles in a simulation have been defined. This approach allows for e.g. multiple models of an interaction to exist in GEANT4 and to be easily replaceable in user code depending on the desired application. The geometry in which particles are propagating in GEANT4 is constructed by placing constructive solid geometry (CSG) primitives, such as boxes or cylinders, with associated material information and geometry transformations inside each other. The result is a geometry tree originating from a 'world'-volume. The geometry primitives can be combined with Boolean operations. In addition to CSG primitives, geometries constructed from their boundary surfaces can also be used. Logical movement of a particle in the simulation from one volume to another requires evaluation if the particle trajectory intersects another volumea computationally expensive operation if done naively. In general, various optimization techniques are used. GEANT4 uses a voxelization algorithm which allows it to perform well, even with a flat (i.e. non-hierarchical) geometry tree. GEANT4 has successfully been used for a number of applications, ranging from simulations of large complex high-energy experiments to microdosimetry.

Hadronic physics models in GEANT4
The hadronic models in GEANT4 follow three different basic approaches: data-driven, parameterization-driven and theory-driven modeling; in the overall framework they offer both complementary and alternative options. The main data-driven models in GEANT4 deal with neutron and proton-induced isotope production, and with the detailed transport of neutrons 6 at low energies. The codes for neutron interactions are generic sampling codes, based on the ENDF/B-VI data format [30], and evaluated neutron data libraries such as ENDF/B-VI [31], JENDL3.2 [32], JENDL-HE [33,34] and FENDL2.2 [35]. The approach is limited by the available data for neutron kinetic energies up to 20 MeV, with extensions to 150 MeV for some isotopes. The data-driven isotope production models that run in parasitic mode to the transport codes are based on the MENDL [36] data libraries for proton-and neutron-induced production. They complement the transport evaluations in the sense that reaction cross sections and final state information from the transport codes define the interaction rate and particle fluxes, and the isotope production model is used only to predict activation. The low-and high-energy parameterized models (LEP and HEP) in GEANT4 are based on the GEANT3 GHEISHA package [37] and include fission, capture, inelastic and elastic scattering processes. The energy regimes of the models extend up to 25 GeV for the low-energy version, and from 10-20 GeV up to 10 TeV for the high-energy version. The assumptions of the parameterized model do not hold in the sub-GeV energy regime [38]. In GEANT4 there is also a leading particle biased model, which is a partial rewriting of the MARS code system [39] called G4Mars5GeV.
Theory-based modeling is the basic approach in many hadronic models. It includes a set of different theoretical approaches to describe hadronic interactions, depending on the addressed energy range and computing performance needs. One important theory-based model is the intra-nuclear cascade (INC) model. If the wavelength of a projectile incident on a nucleus is shorter than the average distance between the constituent nucleons, and also shorter than the mean-free path inside the nucleus, the particle-nuclei reactions can be viewed as a sequence of particle-nucleon reactions which results in an INC of particles inside the nucleus. It is generally considered that this condition is fulfilled inside the nucleus in the energy region 0.2-1 GeV [40]. Since any secondary particles created in the cascade would also have to obey these restrictions, the INC model has its limitation. However, the majority of the reactions are peripheral reactions where the density is lower and therefore the mean free paths longer. One discriminant between INC and quantum molecular dynamics (QMD) models is that in QMD the nucleons move in a self-consistent, continuously evolving nuclear field, but in INC models they have a frozen motion and position in the nuclear mean field. Additionally, Pauli-blocking of secondary particles can occur, limiting the creation of low-energy particles. In practice, INC models are used for particle energies down to about 100 MeV. Two INC models are available in GEANT4: the Bertini cascade (BERT) model [41] and the BIC model [42] 5 . In the BERT model, the nucleus is described with three concentric spherical regions, approximating the nuclear matter density distribution and Fermi energy. In the BIC model, each nucleon has its own Fermi momentum and at the beginning of a reaction an impact parameter selection is made, and an initial interaction point is determined. If the projectile crosses the nucleus without interacting, a new selection is made. Interactions are then calculated between cascade particles and nucleons, or the nuclear medium, according to free-space cross sections modified with respect to Fermi motion. The nuclear potential is taken into account. Reactions generating secondaries violating the Pauli principle (i.e. if the secondary particle has a momentum below the local Fermi momentum) are suppressed. The cascade progresses until various limitations on the remaining energy in the nucleus are reached. Modeling of the remaining nucleus after the cascade is done according to the sequence of pre-equilibrium calculations, 7 evaporation/break-up and final nucleus de-excitation. The BERT model contains its own implementations of these models, whereas the BIC model uses the default GEANT4 implementations described below. At the end of the cascade, several cascade particles are in general still propagating through the nucleus. These are predominantly of a lower energy than what is suitable according to INC assumptions. The nucleus is therefore handed over to the precompound model [37], which provides an exciton-based statistical model of the emission of nucleons and fragments. At the end of this stage, the remaining excitation energy is shared by a large number of nucleons. Nucleons and light fragments can 'evaporate' from the nucleus at this stage, or the nucleus might undergo fission. The remaining nucleus is then de-excited by emission of gamma-rays. Several models for calculating heavier nucleus-nucleus interactions in GEANT4 are under testing and benchmarking, but no validated model has yet been included. A description of these models is beyond the scope of this paper.
For the ion physics, GEANT4 also includes the abrasion-ablation model used in the NUCFRG2 model [43], G4WilsonAbrasionModel and G4WilsonAblationModel. This abrasion model is a simplified macroscopic model for nucleus-nucleus interactions based largely on geometric arguments, and a nuclear ablation step to provide a better approximation for the final nuclear fragment from an abrasion interaction.

PHITS
The PHITS [18] is a three-dimensional Monte Carlo code which can simulate not only the trajectory of the charged particles in a magnetic field but also the collisions and the ionization process at the same time, in the same way as GEANT4. PHITS can calculate fluxes, doses, energy-deposition distributions and many other observables, and is used in a wide field of applications, e.g. semiconductor soft errors, BNCT, shielding around accelerators, radiotherapy, and radioprotection of personnel on space missions. Both the combinatorial geometry (CG), which is used widely by the HETC, NMTC and LAHET, and the general geometry (GG), which is used by the MCNP code system, are available using the same description for all geometries as in the PHITS geometry input. PHITS is currently maintained by a cooperation between RIST (Research Organization for Information Science and Technology, Japan), JAEA (Japan Atomic Energy Agency, Japan), KEK (High-Energy Accelerator Research Organization, Japan) and Chalmers University of Technology (Sweden). Validation of calculations with heavy ion measurements performed at GSI (Gesellschaft für Schwerionenforschung) in Germany, HIMAC (Heavy Ion Medical Accelerator in Chiba) in Japan, and NSRL (NASA Space Radiation Laboratory) at BNL (Brookhaven National Laboratory) in the USA are ongoing and results have been presented in e.g. [44]- [55].

Hadronic interaction models in PHITS
In PHITS, neutrons can be transported from thermal energies up to 200 GeV, and the same method as in the MCNP4C code [56] is employed for neutrons with energies between 1 and 20 MeV based on Evaluated Nuclear Data such as the ENDF/B-VI [57], JENDL-3.3 [58] and LA150 libraries [59]. Above 20 MeV, the simulation model JAM (jet AA microscopic transport model) developed by Nara et al [60] is used. JAM is a hadronic cascade model, which explicitly treats all established hadronic states including resonances with explicit spin and isospin as well as their anti-particles. All hadron-hadron cross sections are parameterized based on the resonance and string model by fitting the available experimental data. Below the energy in the center-of-mass (c.m.) system E < 4 GeV, the inelastic hadron-hadron collisions are described by the resonance formations and their decays, and at higher energies, string formation and their fragmentation into hadrons are assumed. For protons and other hadrons, JAM is used above 1 MeV, but for charged particles below 1 MeV only the ionization process is considered until the particles are stopped. The total reaction cross section, or the lifetime of the particle for decay, is an essential quantity in the determination of the mean free path of the transport particle. In PHITS either the NASA systematic, developed by Tripathi et al [63]- [65], or the Shen formula [66] can be used for the total nucleus-nucleus reaction cross section. When calculating the transport of a nucleus in materials using PHITS, only the ionization process is taken into account below 10 MeV n −1 , but above 10 MeV n −1 nucleus-nucleus collisions up to 100 GeV n −1 are described by the simulation model JQMD developed by Niita et al [69]. In the QMD model, the nucleus is described as a self-binding system of nucleons, which are interacting with each other through the effective interactions in the framework of molecular dynamics. One can estimate the yields of emitted light particles, fragments and excited residual nuclei resulting from the heavy ion collision. The QMD simulation describes the dynamical stage of the reactions. At the end of the dynamical stage, excited nuclei are created and must be forced to decay in a statistical way to get the final observed state. In PHITS, the GEM model [70] is employed for light particle evaporation and fission process of the excited residual nucleus.

Interfacing PHITS and its physics models to GEANT4
The integration of the PHITS, JQMD and JAM interfaces into the GEANT4 hadronic framework is shown schematically in figure 1. The G4HadronicProccess class is the base class of hadronic interactions and it has two important functionalities. One is calculating the cross section of a reaction and proposing the mean free path to the GEANT4 kernel, and the other is generating final state particles from the reaction. Hence, the cross section of a reaction and 9 its final state generator are completely divided in the framework. Each concrete class of the G4HadronicProccess has cross section data class(es) and final state generator class(es). These final state generators are called models in the framework. G4HadronicInteraction in the figure is the base class of these models. PHITS2G4InelasticModel is a concrete class of G4HadronicInteraction, in which the PHITS library is interfaced. JQMD2G4InelasticModel and JAM2G4InelasticModel are the classes which contain interfaces to JQMD and JAM, respectively. G4HadronicDiscreteProcess in the figure is a derived class of G4HadronicProccess and becomes the base class of the process for usual point-like hadronic interactions handling models and cross sections. Native cross section classes of GEANT4 based on Shen [66], Kox [61], Sihver [62] and Tripathi [63]- [65] formulae are used for calculating mean free paths in geometries. It should be mentioned here that final state generators also need cross sections while calculating a reaction, for example, sampling impact parameter and so on. For this kind of usage of cross sections in JQMD, JAM and PHITS, their own cross sections are used. To establish the re-productivity of Monte Carlo events, the same random number generator is used both in GEANT4 and in PHITS, JQMD and JAM.

Neutron fluences in the atmosphere
A comparison of calculated and measured particle fluxes and effective neutron dose rates in the atmosphere has recently been published [13]. The calculation in this paper was performed with the PLANETOCOSMICS [14] package based on the GEANT4 [2] and two local interstellar models of galactic cosmic rays. To better understand how much the results differ when using the native BIC implemented in GEANT4 with the JQMD/JAM interface a comparison of the secondary particle fluences calculated with these models is presented below. The calculations of the interactions of primary ions with molecules in the air allow an estimation of the influence of heavy and high-energy ions of cosmic origin on the radiation environment in the atmosphere. The BIC model implemented in GEANT4 is only recommended to be used for inelastic interactions of ions with Z < 6 and energies below 10 GeV n −1 , since it does not consider participant-participant collisions. The JQMD model should be valid for all heavy ions up to around 100 GeV n −1 .
To overcome the energy restriction for alpha particles which are the second most important primary source for secondary particle fluences, two free protons and two free neutrons (2p2n) with a fourth of the energy have been used to replace one alpha particle over the whole energy range up to 100 GeV n −1 in order to evaluate the importance of alpha particles with an energy above 10 GeV n −1 . Besides this replacement, the same physics models were used for the 2p2n approximation and the BIC model, which also means that secondary ions produced in collisions are transported using the BIC model. Additionally, the interface by Koi et al [19] as described in section 5, was used in GEANT4, which adds the alternative to choose the JQMD and JAM models for heavy ion physics implemented in the PHITS code.
Altogether, three different scenarios were tested against each other where the differences only apply to the physics of ion projectiles with A > 1: • native GEANT4 BIC model for all ions with A > 1 and E < 10 GeV n −1 ; • approximation of an alpha particle by two free protons and two free neutrons and the BIC model for ions heavier than alpha; • JQMD/JAM interface with the JQMD model used for all ions with A > 1. The CREME96 model [67] was chosen to describe the primary particle fluences with an isotropic angular distribution for ions from hydrogen up to iron and energies from 1 MeV n −1 to 100 GeV n −1 . Ions heavier than carbon have been considered only in the case where the JQMD model was used and the energies are restricted to E < 10 GeV n −1 for the binary cascade model. The geometry consists of a layered atmosphere and was created with the PLANETOCOSMICS toolkit applying the NRLMSISE00 [68] model to set the density, temperature and composition of the atmosphere. The source of the primary particles was set to an altitude of 180 km above ground with an isotropic angular distribution. Figures 2-4 show the differential fluence rate of secondary neutrons at different altitudes for the three different approaches to treat the transport of ions in the atmosphere and it can be seen clearly that in all cases the BIC model falls short of both the JQMD model and the case where the 2p2n approach was used. The differential fluence rate dJ/dE, which is shown in figure 2-4, is the number of particles per area, time and kinetic energy.
For high altitudes, the agreement between the JQMD model and the 2p2n method is good only for energies below a few tens of MeV, whereas for higher energies above a few hundreds of MeV the primary neutrons which have not yet suffered a significant energy loss result in a very high fluence rate when free nucleons are used. These free neutrons are obviously unphysical and result from the approximation of an alpha by four free nucleons. The result changes for mid-altitudes, where it can be assumed that all primaries have undergone at least one major interaction with molecules from the air and the majority of primary alphas are fragmented. The agreement between the 2p2n approximation and the JQMD model at these altitudes is very good. Obviously, heavier ions than carbon do not play an important role and the 2p2n approximation can be used to calculate secondary neutron fluences; calculations for other secondaries show similar results.
Further down in the atmosphere the situation is different. In figure 4, at an altitude of 10 km, good agreement between JQMD and the 2p2n approximation is observed for energies above 10 MeV, whereas for lower energies a difference can be seen. This behavior is a result of the lower neutron fluences mainly originating from primary protons and not just from heavier or more energetic primary ions, which can be recognized in figure 4 from the contribution of primary protons to secondary neutrons. This effect may be explained by the differences in the transport of secondary ions produced by collisions of high energetic primary protons with molecules in the air. These secondary ions are then transported by the ion physics model used in the corresponding calculation giving rise to the differences in the neutron fluences. This is evidence for the fact that the selection of an ion physics model can have an impact on secondary particle fluences if the transport of primary protons is calculated.
The difference in the secondary neutron fluence is not very significant for radiological issues though, as the biological effectiveness decreases with lower energies (see conversion coefficients in [73]) and is at least one order of magnitude lower for neutrons with E < 1 MeV compared with a few hundreds of MeV.
One can conclude that primary ions heavier than alpha particles do not play an important role in the calculation of secondary particle fluences and the radiation exposure at aviation altitudes in the atmosphere, but contributions of alpha particles with energies above 10 GeV n −1 cannot be neglected. The approximation to use four free nucleons (two protons and two neutrons) to substitute an alpha particle can be used for altitudes below ∼20 km but fails for higher parts of the atmosphere.

Irradiation from GCR in space
The earth's atmosphere and the geomagnetic field provide a very good shielding against ionizing radiation. It was shown in the last section that the contribution of ions heavier than He to the radiation environment in the atmosphere is marginal due to their limited range. The situation is different at ISS altitudes or in deep space, where heavier ions contribute a considerable amount to the radiation exposure. When estimating the radiation risks in space, it is therefore crucial to use a transport code which includes models for nucleus-nucleus interaction in a sufficiently broad energy range. The code chosen by ESA for calculating the radiation environment in spacecraft is GEANT4 [4]. However, problems have been identified with the existing GEANT4 physics models required for simulating the radiation environment on the ISS [6], where heavy ion interactions from GCR with both the materials of the ISS and the human body are important. One way to solve this problem and to increase the accuracy of the calculations of particle spectra after the GCR has penetrated different shielding materials at LEO and in deep space is to interface the GEANT4 code to other existing heavy ion reaction models and heavy ion transport codes, e.g. JQMD and PHITS. To certify the accuracy of the heavy ion physics models in PHITS, an extensive benchmarking is currently ongoing [44]- [55]. A neutron spectrum measured inside the Space Shuttle by the Bonner Ball Neutron Detector (BBND) during the STS-89 flight [71] has also been compared with the corresponding data calculated [72] using PHITS. The agreement between the calculations and the measurements has been satisfactory for space-shielding applications.
To compare the abilities of the native GEANT4 binary cascade model with the JQMD/JAM interface, the irradiation of an 'imaginary' spacecraft has been performed for primary ions up to iron using input spectra from the CREME96 [67] model during solar minimum in a geosynchronous orbit. The simulations were performed by generating particles on the inside surface of a source sphere with a diameter of 14 m. In the middle of the sphere, the imaginary spacecraft was simulated as a simple closed aluminium cylinder with length 10 m and a diameter of 3 m. The thickness of the walls of the cylinder was varied between 0 and 10 mm both on the long and short sides. The primary ions were transported through the aluminium shielding of the spacecraft, and both the absorbed and the equivalent dose rate in a water sphere placed at the center of the spacecraft were calculated for different shieldings. The sphere was divided into shells with a thickness of 1 cm, and the dose was determined for each layer. The different layers of the water sphere in the spacecraft provide additional shielding, and by looking at each shell individually, the dose can be estimated for increasing shielding thickness.
The results were obtained using two different approaches for the transport of ions with Z > 1: • BIC model for Z = 1-26; • JQMD/JAM interface with the JQMD model used for all ions with Z = 1-26. Moreover, the contributions from the different ions were investigated to check if the recommended restriction in the BIC model to Z > 6 delivers acceptable results. The results for the dose equivalent calculated at different depths for a 1 mm aluminium hull using GEANT4 with the JQMD model on the one hand, and the binary cascade model on the other, are compatible with each other (figure 5). The dose from primary iron ions contributes up to 25% of the total dose and all ions with Z > 6 add up to a total contribution of more than 60%. This shows the importance of including the contribution from ions heavier than carbon even if the binary cascade model is not optimized for such heavy ion interactions. From figure 5, it is also obvious that the contribution from heavier ions to the total dose equivalent decreases with increasing shielding. This is due to the breakup of heavier ions into fragments. The fraction of the dose equivalent from primary iron ions even exceeds the dose deposition caused by protons close to the surface of the sphere. The calculated dose equivalent is between ∼65µSv h −1 (∼1.56 mSv d −1 ) at the surface of the sphere and ∼30 µSv h −1 (∼0.72 mSv d −1 ) at the center of the sphere, i.e. after an additional shielding of about 25 g cm −2 of water. The quality factor is calculated according to ICRP60 [12] and is a function of the linear energy transfer (LET) and mirrors the biological effectiveness of a particle (for details see [12]). The quality factor averaged over the shells is Q ≈ 2.4 using the JQMD and the BIC model including ions not heavier than iron and Q ≈ 1.5 if the BIC model is used in its suggested limits. The decrease of the quality factor from Q ≈ 4 at the surface to Q ≈ 2 at the center of the sphere that is shown in figure 6 reflects the attenuation of heavy ions with a high quality factor in water. Of course this effect cannot be observed if the ions under consideration are restricted to Z 6, where the quality factor is constant over all depths (Q ≈ 15). The decrease in the dose equivalent is directly related to the change of the radiation field and the resulting change in the quality factor and does not mirror a change in the absorbed dose, which does not change significantly.   In figure 7, we show the dose equivalent calculated at different depths for a 10 mm aluminium hull. The BIC model predicts absorption rates that are higher than those predicted by JQMD for depths larger than about 10 g cm −2 , which reflects the differences in the details of the two fragmentation models. As the shielding thickness increases, higher-generation particle production becomes more important; the discrepancy observed here might be connected with the failure of the BIC model in the description of reaction for ions heavier than carbon or more energetic than 10 GeV n −1 .

Conclusions
We have compared the calculations of neutron fluence distributions in the atmosphere with the PLANETOCOSMICS Monte Carlo package based on the GEANT4 Monte Carlo code, using both the BIC model and JQMD, included in the PHITS code, for the inelastic interactions of heavy ions with matter. We have also calculated the total dose equivalents and the quality factors at different depths in a water sphere in an imaginary spacecraft during solar minimum in a geosynchronous orbit. The BIC model implemented in GEANT4 is only recommended to be used for inelastic interactions of ions with Z < 6 and energies below 10 GeV n −1 , since it does not consider participant-participant collisions. In this paper, it has been shown that the use of four free nucleons (two protons and two neutrons) to substitute an alpha particle gives good accuracy when calculating the neutron fluence distributions with the BIC model at altitudes below ∼20 km, but fails for higher parts of the atmosphere where there are more highly energetic alpha particles. It has also been shown that if the recommended limitations of the BIC model are followed, there will be a large difference between the results when using the BIC model and the JQMD model. The magnitude of this difference depends on the shielding material and the thickness due to the restriction of E < 10 GeV n −1 for the BIC model. Magnetic shielding could also influence the results when lower energetic ions are reflected and energies above 10 GeV n −1 become more important (as we see in the differences for high altitudes in the atmosphere). So even at LEO and at lower latitudes, i.e. higher cutoff, there might be a difference between the results from BIC and JQMD. This should be further investigated. However, if both the BIC and JQMD models are used for nuclei with 1 Z 26, with energies E < 100 GeV n −1 for JQMD and E < 10 GeV n −1 for BIC, the difference seems to be minor. The main conclusion is therefore that there is no major difference in the results when calculating the neutron fluence distributions in the atmosphere and dose equivalents in space with GEANT4 using the BIC model if it is used for nuclei 1 Z 26, with energies E < 10 GeV n −1 , compared with using the JQMD model for the same nuclei and energies below 100 GeV n −1 .