Production rates of cosmogenic nuclides in extraterrestrial material using GEANT4 software

We present a model for the calculation of the production rates of cosmogenic nuclides in extraterrestrial material. The model is based on the Monte Carlo simulation software Geant4. Using this software an application simulating the irradiation of the spherical body with predefined chemical composition by galactic cosmic-ray protons was developed. The fluxes of secondary neutrons and protons generated within the body of the sample are calculated. These are further used for the calculation of the production rates for cosmogenic nuclides. The plausibility of the presented model was tested in case of the well-studied meteoritic sample Knyahinya for which the production rates of cosmogenic nuclides were previously measured and also calculated using MCNP simulation software. The production rates for 3He, 10Be, 21Ne, 22Ne, 26Al, 36Cl, 38Ar, 39Ar, 41Ca and 53Mn have been calculated in Knyahinya meteorite and for 10Be, 26Al and 36Cl in Apollo 17 sample 73,002.


Introduction
Cosmic rays represent subatomic particles (protons 89%, helium nuclei 10% and heavier nuclei ≈ 1%) reaching energies from a vast interval reaching up to tens of EeV (1 EeV = 10 18 eV) region.Different sources of cosmic rays are distinguished.In the region of lowest energies up to several GeV the Sun is the primary source of cosmic-ray particles.Afterwards there are galactic and extragalactic cosmic rays.
These particles represent the natural environment to which bodies of the Solar system are exposed.Interaction of these high-energy particles with the matter leads to the generation of secondary particles-protons, neutrons, pions, α-particles etc. [1].Generated cascades of these particles further lead to the spallation reaction within the body and to the generation of the cosmogenic radionuclides.Investigation of the depth dependent production rates of the cosmogenic nuclide within the irradiated body can provide information about its characteristics or enables us to calculate cosmic-ray exposition (CRE) time [2][3][4].Also, the shielding of the matter can be studied.In case of meteorites, it is possible to estimate the fall time on the Earth of the sample or to estimate the dimension of the parent body.Nowadays, there exist useful simulation software, with which it is possible to simulate the interaction of cosmic-ray particles with the body of predefined chemical composition such as Geant4 [5,6] or MCNP [7,8].Recently, there existed several already published works such as [1,[9][10][11][12][13] which are focusing on the simulation of the bombardment of different meteoritic samples or other minor planets of the Solar system by cosmic-ray particles.These samples differ in the dimension, density, or chemical composition.These are the parameters which can influence the final production rates of cosmogenic nuclides.
Our aim was to perform calculations using Monte Carlo simulation toolkit.
Geant4 which is nowadays vastly used in different fields of particle physics or astrophysics.We present a model based on Geant4 software which simulates a sphere of predefined chemical composition bombarded by cosmic-ray protons and which further calculates the production rates of selected cosmogenic nuclides.The applicability of the presented model has been tested in case of the calculation of the cosmogenic nuclides in the sample of meteorite Knyahinya and further extended by calculating the rates inside Moon using Apollo 17 samples.In order to calculate the production rates, we have used the theory described in Method section.In Geometry section, we focus on the description of the geometry implemented in Geant4 environment for our analysis.The analysis has been divided in 3 separate steps which are discussed in Analysis section.The results of our simulations for Knyahinya meteorite and Apollo 17 sample together with discussion form the part of Results.Last section is the Conclusion of our work.

Method
The energy spectrum of cosmic rays is steeply decreasing.Origin of these particles differ according to the energy region.Particles with the lowest energies (up to the GeV region) originate in the Sun.Their trajectories and flux are significantly influenced by solar magnetic field and are also dependent on the solar cycle.Their energies are too low and therefore they cannot initialize nuclear reaction in the irradiated body.Some of most energetic solar particles can start nuclear reaction leading to the production of cosmogenic nuclides, however only in outermost layers of irradiated body.These layers are removed during the burning of meteoroid in Earth atmosphere or are influenced by gardening effect, in the case of planetary bodies and therefore they are not subject of our interest.Particles with energies from GeV region up to 10 15-16 eV, where the so called "knee" is observed are of galactic origin.These represent the major interest of our work as these are particles significantly contributing to the production rates of different cosmogenic nuclides caused by nuclear-spallation reactions.
In the vicinity of the Earth the energy spectrum of galactic cosmic-ray protons can be approximated by the refined formula adopted from [14] In Eq. 1, E represents proton kinetic energy in MeV, E 0 is the rest energy of a proton (938 MeV), θ is the solar modulation parameter in MeV.We have adopted θ = 550 MeV corresponding to the calm phase of the solar cycle.C = 1.244 × 10 6 cm −2 s −1 MeV −1 , K = 780•exp (− 2.5 × 10 −4 × E) and γ = 2.65.The presented galactic cosmic-ray energy spectrum was implemented in our simulations.
Galactic protons bombard bodies presented in the Solar System leading to the nuclear reactions within these bodies.They generate secondary particles such as protons, neutrons and α-particles.These have sufficient energy in order to ignite reactions through different mechanisms and leading to the production of isotopes of nuclides in different depths.Different types of cosmogenic isotopes can be considered.Firstly, they are stable isotopes.These are nuclides which do . not undergo further decay and considering that nuclei cannot escape from the sample, their concentration is increasing with time.Afterwards, there are short-lived radionuclides which have the lifetime on the order of hours, days or several years.These can be used for the calculation of the time from the fall of the sample on the Earth ś surface.With sufficient exposition time in the Solar system, the concentration of the generated nuclei reaches a value which represents the balance between their production rate and decay.As the sample falls on Earth, it is not anymore exposed to the cosmic-ray bombardment and the concentration of the nuclides begins to decrease.The last group of cosmogenic nuclides are longlived nuclides with the lifetime on the order of 100 000 of years and more.As the sample is exposed to cosmic rays in the extraterrestrial environment for long-enough time and the lifetimes of these nuclides are also long-enough, these nuclides reach saturation where it is possible to neglect the decay of the nuclides.There are several interesting isotopes which are generally measured and studied in the samples of the meteorites.The model we are presenting has been tested for long-lived nuclides 10 Be, 26 Al, 36 Cl, 41 Ca and 53 Mn.Their half-lifes can be seen in Table 1.We have calculated also production rates in case of 39 Ar which is radioisotope with significantly shorter half-life on the order of 268 years.The resulting production of the specified cosmogenic nuclide is dependent on two main parameters.Firstly, it is the flux of the secondary particles (protons, neutrons and α-particles) generated due to the interaction of the primary cosmic particles.In our work we are focusing only on primary protons generating secondary particles.Primary α-particles as a contributor to the production rate of the desired nuclide are neglected in our simulations.The fluxes can be simulated using different Monte Carlo based simulation software.After calculation of the fluxes for different depths, it is necessary to consider the excitation functions for the desired cosmogenic nuclide.These represent the major source of the inaccuracies in the calculations of the Ca have been adopted from http:// www.nucle ide.org/ DDEP_ WG/ DDEPd ata.htm library. 10Be half-life was adopted from [15]. 53Mn halflife was adopted from [16] Isotope Half-life (years) final production rates as these values must be empirically measured or calculated using various models.However, the measurements for different channels are usually carried out for a limited range of energies and for a limited number of values which means the interpolation, or the extrapolation of the values are needed.In some cases, the measurements for the desired channel are missing.The resulting production rate P j can be calculated as [9] In Eq. 2, index j represents a type of cosmogenic nuclide within the spherical sample of radius R. N A is Avogadro constant, A i the atomic mass of target element i, c i is the concentration of the target element i in g/g and k is the index for the reaction particle type (primary proton, secondary proton and secondary neutron).σ i,j,k (E) represents the excitation function (or cross section) for the production of nuclide j from target element i by the reactions induced by particles of type k.J k is the differential flux density of particles of type k.P j can be calculated directly in S.I. units.In the case of stable isotopes ccSTP g −1 My −1 (ccSTP is cubic-centimeters-at-standard-temperature-and-pressure) is used in order to quantitatively represent the number of atoms generated over a millions of years.In the case of radioactive nuclides more frequently dpm.kg −1 (decays per minute per kilogram) is used.

Geometry
The aim of the presented work has been to use Monte Carlo simulation software Geant4.Using this software environment, a spherical body is created of predefined chemical composition which is bombarded by the flux of galactic cosmic-ray protons corresponding to the flux described by Eq. 1.The inner concentric shells are created within the body in which the generation of secondary protons and neutrons is tracked.The flux of secondary particles is further used for the calculation of the production rates of different cosmogenic nuclides.Between concentric shells, we have assumed constant step of decreasing diameter which is parameter directly set by the user.The whole spherical body represents in this analysis the minor body of the Solar System of homogenic density corresponding with its chemical composition to the chemical composition of the known sample.

Geant4
For our simulations we have been using Geant4 (GEometry ANd Tracking) Monte Carlo simulation toolkit.Geant4 is free for download and available at the webpage site https:// (2) geant4.web.cern.ch/.This software simulates the propagation of high-energy particles through matter [5,6,17], Geant4 is very versatile and multifunctional software which enables to the user to develop his own geometrical object and the environment exposed to the flux of high-energy particles.It is written in C + + and enables to create the resulting files in ROOT [18].
Using Geant4, we have created the studied geometry as described earlier.We further prepared the scripts where user can directly specify chemical composition and dimensions of the object.Geant4 offers a choice of the so called Physicslists.Physicslist is one of the three mandatory classes provided by the user to the G4RunManager.It specifies the particles that can be used in the simulation and the list of the interaction processes of these particles.In the presented analysis, we are focusing mainly on the use of Physicslist QGSP_BIC_AllHP (it use Quark gluon string model for energies > ~ 20 GeV and Binary Cascade Model for energies < ~ 10 GeV).The reason for choosing this Physiscslist, is the tracking of neutron energies with high precision which it offers.In our simulation, there is a significant number of secondary neutrons which are one of our aims of interest.For more information about Physicslists see the official webpage of Geant4 software https:// geant4-userd oc.web.cern.ch/ Users Guides/ Physi csLis tGuide/ html/ index.html.In this analysis we decided to use only one Physicslist as our intention is to compare the validity of the presented methodology in comparison with the previously published works.In the presented model, we have used Geant4 only for the calculation of the fluxes of secondary particles.As a next step, the dataset of cross section, adopted from [19], have been implemented for the calculation of production rates for cosmogenic nuclides.

Analysis
We have divided the analysis into 3 separate steps.In the first step the irradiation of the investigated body by the flux of cosmic-ray protons is simulated.We have assumed the integral cosmic-ray proton flux to be 4.8 protons cm −2 s −1 as it was assumed in the work [20].The energy of the injected protons follows the energy spectrum according to the Eq. 1.
The next step is the calculation of the number of secondary protons and neutrons for each shell.The numbers are further stored in the separate files containing histograms for each shell.The energy and the absolute number of particles for each energy bin are stored.Files are created separately for protons and neutrons.
As a part of the last step, the information about the excitation functions for different nuclides are included.For this purpose, separate files have been created for different interaction channels containing the energy in MeV and experimentally measured cross sections expressed in milibarns.
From Eq. 2, there are 3 parameters which define the occurring spallation reaction: target nucleus, secondary particle (neutron or proton in our case) and outcoming cosmogenic radionuclide.Tables for each interaction have been created.
In the final step we numerically calculate the production rates using Eq. 2 for specific nuclide.The integral in Eq. 2 was approximated as a discrete calculation using the trapezoidal rule for numerical integration.
The output of these calculations is the depth profile for specific nuclide-the production rate of studied nuclide for each depth.In order to compare our calculations with experimental data, we have performed normalization to the flux of galactic cosmic-ray protons, real chemical composition of studied object and to ccSTP g −1 My −1 in case of stable isotopes and to dpm.kg −1 in case of cosmogenic radionuclides.

Excitation functions
As we have already mentioned, the major inaccuracy in the calculations of the production rates of cosmogenic nuclides in general represents the sparse knowledge of excitation functions of the spallation reactions needed for the calculation of the production rates.In order to define the value of cross section for each energy bin of the secondary particles, it is necessary to perform linear interpolation of excitation function values.We have used the simplest linear interpolation method during every calculation.
The used cross section [19] contains the information about the measured values of cross section for all 3 parameters (target nucleus, incident secondary particle and outcoming nuclide).In work [10], the authors used the same dataset of cross sections for the calculation of production rates of different cosmogenic nuclides for Knyahinya meteorite.For this purpose, we have assumed the results of this work as reference results to compare.
The values of excitation functions we have been using are limited for energies up to 10 GeV.However, the energies of secondary particles are reaching also higher energies as we have adopted the maximum energies for primary protons to be up to 20 GeV following the distribution according to the Eq. 2. For this purpose, the extrapolation of the values was needed.The simplest assumption has been made that the excitation function keeps constant trend after reaching the last experimentally measured value.

Knyahinya meteorite
We have performed the simulation of the bombardment of Knyahinya meteorite by galactic cosmic-ray protons.This meteorite has been selected as it was well studied by previous works, mainly [10] so it serves as a valuable crosscheck of the validity of the presented model.Its chemical composition is known and can be seen in Table 2.
We have performed simulations using the presented model assuming the integral galactic cosmic-ray flux to be 4.8 protons cm −2 s −1 .The spherical body with diameter 90 cm and homogeneous chemical composition according to the Table 2 with bulk density 3.35 g.cm −3 was simulated.This meteorite is LL5 chondrite.These input parameters copy the parameters input in the calculations of [10].The step for concentric shells was set to be 2.5 cm which represents in total 18 shells over the whole body.This step was chosen to represent sufficient resolution to study the depth profile for different cosmogenic nuclides.
15.1 millions of primary galactic protons were injected in the simulation with the distribution in energies according to the Eq. 2 with energies from interval 0 MeV up to 20 GeV.We have assumed that this represents a statistically large enough sample.In Fig. 1, we present the energy spectrum in absolute values of 15.1 millions of injected primary galactic protons.
The calculated fluxes for both protons and neutrons for different depths can be seen in Fig. 2. As a comparison, the fluxes from [21] near the surface and halfway to the center of the sample are also introduced.The last bin was omitted as this represents the very center of the meteorite.
Using the analysis process described in Analysis section, the production rates of five noble-gas isotopes, namely 3 He, 21 Ne, 22 Ne, 38 Ar and 39 Ar.Besides 39 Ar, production rates for other five radioactive cosmogenic nuclides, namely 10 Be,  and 4).In Fig. 3, calculated production rates for 7 nuclides are plotted compared with the measurements from [22].In Fig. 4, we introduce predictions for three nuclides, 36 Cl, 39 Ar and 41 Ca, for which measurements are not introduced.

Apollo 17 sample 73,002
We have performed calculations of the production rates of three nuclides, namely 10 Be, 26 Al, 36 Cl, inside the upper layer of lunar surface (see Fig. 5).The chemical composition was chosen to be similar to the drill 73,002 of Apollo 17, see Table 3.We have simulated actual dimension of the Moon with diameter 3474 km and density 1.8 g.cm −3 .Incress 2.5 cm was assumed between different layers.However, we visualize the final results in g.cm −2 in order to be capable to perform the comparison with results from other drills.15 millions of primary galactic protons have been simulated.Integral flux 4.8 protons cm −2 s −1 was assumed.

Discussion
Fluxes of secondary neutrons and protons have been calculated for different depths of meteorite, see Fig. 2. Comparing the calculated fluxes for both protons and neutrons with the calculations based on MCNP simulation code presented in [21], it is possible to see that the energy spectra are in satisfying accordance with the published predictions.Furthermore, it is possible to see that the relative flux of protons is decreasing with increasing depth in the meteoritic sample.This is also in agreement with the predictions from [21].The relative flux of neutrons is increasing with increasing depth what has been also expected.It is possible to assume that our model correctly calculates the fluxes of secondary neutrons and protons for different layers.
In the next step, the productions rates for ten different cosmogenic nuclides for Knyahinya meteorite have been calculated, Figs. 3 and 4. The presented model gives results which are in good agreement with the published production rates calculated in work [10].Comparing results with the measurements from [22], it is possible to see that our predictions are supported by measured data for the majority of studied cosmogenic nuclides for which there were accessible data.
There is observed slight discrepancy between the calculated and measured values of the production rates in case of 26 Al and 38 Ar.We assume that this can be caused by the set of used cross sections during the calculations.It would be beneficial to use different set of cross sections or refined cross sections and observe the influence on the final production rates.
The presented model was further applied for the calculation of the production rates of 10 Be, 26 Al and 36 Cl in Moon.For this purpose, the actual dimension of the Moon was simulated with the chemical composition corresponding to the Apollo 17 lunar samples.It is possible to directly see the Fig. 1 Energy spectrum of injected primary galactic cosmic-ray protons Fig. 2 Fluxes of neutrons (left) and protons (right) in different depths of Knyahinya meteorite.Fluxes near the surface and halfway to the center of the sample from [21] are introduced Fig. 3 Depth profiles of calculated production rates for Knyahinya for different cosmogenic nuclides using the presented model.Data are adopted from [10].Measurements from [22] are plotted.First row: exponential suppression of the production rates with increasing depth.In case of 36 Cl the contribution from neutrons is dominant over the whole depth profile.On the other side for 10 Be, 26 Al the protons represent the dominant contribution near the very surface of Moon and are quickly suppressed.
In total, the presented model proposes the calculations of the production rates of the studied sample which are in good agreement with previously published works based on different Monte Carlo simulation software and are also in good agreement with the measured data.This means it could be used as a valid model for the future calculations for different samples either of meteoritic material or other types of extraterrestrial material.

Conclusion
In the proposed work, we present the newly developed model for the calculations of the production rates of the isotopes of cosmogenic nuclides based on the Monte Carlo simulation software Geant4.The model has been tested for meteoritic sample Knyahinya calculating the production rates for 10 cosmogenic nuclides.For 7 of these nuclides, there were accessible data to compare with our simulations.We further used our model for the calculation of the production rates of 3 cosmogenic nuclides inside Moon based on the chemical composition of Apollo 17 samples.We conclude that the

Fig. 4
Fig.4 Depth profiles of calculated production rates for Knyahinya for different cosmogenic nuclides using the presented model.Data are adopted from[10].First row: 36 Cl (left) and 39 Ar (right), second row: 41 Ca

Table 2
Percentual abundance of different elements in the Knyahinya meteorite sample [10] Al, 36 Cl, 41 Ca and 53 Mn, have been calculated (see Figs.3 26

Table 3
Percentual abundance of different elements in the Apollo 17 sample 73,002