Calculation and benchmark of fluence-to-local skin equivalent dose coefficients for neutrons with FLUKA, MCNP, and GEANT4 Monte-Carlo codes

Dose equivalent limits for single organs are recommended by the ICRP (International Commission for the Radiological Protection publication 103). These limits do not lend themselves to be measured. They are assessed by convoluting conversion factors with particle fluences. The Fluence-to-Dose conversion factors are tabulated in the ICRP literature. They allow assessing the organ dose of interest using numerical simulations. In particular, the literature lacks the knowledge of local skin equivalent dose (LSD) coefficients for neutrons. In this article, we compute such values for neutron energies ranging from 1 meV to 15 MeV. We use FLUKA, MCNP and GEANT4 Radiation transport Monte-Carlo simulation codes to perform the calculations. A comparison between these three codes is performed. These calculated values are important for radiation protection studies and radiotherapy applications.


Introduction
Neutron radiation dose estimation is of high interest in radiation protection studies as well as radiotherapy applications. In particular, the assessment of local skin equivalent and absorbed dose (LSD) (ICRP116 2010) can be performed using fluence-to-LSD Conversion Coefficients referred as 'LSD-CC' in the rest of this article.
Such conversion coefficients are provided by the International Commission on Radiological Protection (ICRP), in the publication 116 (ICRP116 2010). However, the LSD-CC for neutrons are not available in the literature.
Recent papers can be found in the literature aiming at assessing the neutron dose. This dose has to be evaluated for instance for medical applications of proton therapy (Jarlskog et al 2008) as organ specific neutron dose, for instance. In Veinot et al (2020) the authors address the skin dosimetry specifically for operational quantities (extremity monitoring for neutrons) used to approximate protection quantities. In that work, the authors computed dose coefficients for the local skin dose in the slab, rod, and pillar phantoms, using the latest guidance from the ICRU report 95 (International Commission on Radiation Units and Measurement) (ICRU95 2020). The used phantoms differ from the ICRP recommended geometries for LSD assessment (ICRP116 2010). This difference originates from the quantities that are computed: ICRU95 provides the guidance to compute operational quantities while ICRP116 is the standard for protection quantities. The phantom skin compositions and dimensions are not the same to compute protection or operational quantities. The protection quantity phantom (ICRP116 2010) is described below in section 2. It is a skin cube with 10 × 10 × 10 cm 3 . The operational quantity slab phantom (ICRU95 2020) has 30 × 30 × 15 cm 3 dimensions.
In this paper, we compute the LSD-CC values for equivalent and absorbed skin dose using the irradiation scenario geometry presented in the publication 116, Annex G (ICRP116 2010), in order to complete ICRP data.
The LSD-CC values are of high importance for radiation protection purposes and, as an example, for the IAEA working group A1/A2 whose objective is to propose updated limits of the international regulations for the safe transport of radioactive material (Frosio et al 2019(Frosio et al , 2020. The IAEA A1/A2 working group defines the Q B quantity based on the local skin equivalent dose (IAEASSG26 2012).
In this document, the LSD-CC values for equivalent and absorbed skin dose are computed with three different radiation transport Monte-Carlo calculation codes: FLUKA, GEANT4 and MCNP. A comparison between the results of these three codes is performed to enhance the scientific confidence and the validity of the results.
In the section 2, we present the geometry model implemented in the three numerical simulations and we detail in section 3, the various assumptions taken into account. The next section 4 presents the calculation results of the three codes. A discussion is given detailing the physics processes to explain the results. Finally, in section 5, we present conclusions regarding the results and the comparisons between the three codes.

Calculation geometry and physics models
The geometry model is identical for the three calculation codes. It is setup as defined in the Annex G of ICRP116 (2010). A cube of skin of 10 × 10 × 10 cm 3 is located in vacuum (see figure 1). The used skin composition is defined in ICRP110 (2009) and shown in table 1. A circular source beam of 7 cm diameter is normally incident to the skin cube. Depending on the code used, the simulated source consists of a beam of mono-energetic neutrons (MCNP) or a beam of neutrons with energies samples from a uniform distribution over a given interval (group formalism for doses in GEANT4 and FLUKA). More details about the neutron groups are provided below in section 3. A cylindrical scoring volume is located between 50 and 100 µm depth in the skin cube, with a normal surface of 1 cm 2 . The scored quantity is the deposited energy in the scoring volume. This quantity is transformed into the absorbed dose and then to the equivalent dose using the neutron radiation weighting factors provided in ICRP92 (2003). The radiation weighting factors for neutrons are energy dependent. The multiplication of the resulting calculated doses with the source surface allows establishing the LSD-CC values for the absorbed or equivalent doses depending on the use of the neutron radiation weighting factor. The radiation weighting factors (W R ) considered in this paper are extracted from ICRP103 (2007). W R values are calculated using the energy of incident neutrons on the skin phantom. Equation (1)

Monte-Carlo simulation codes setup parameters
Neutrons are not directly ionizing particles. By their interactions with matter, secondary charged particles such as recoil protons or electrons are produced. These secondary charged particles do not deposit their energy at the location of their creation. They deposit their energy continuously when slowing down via Coulomb interactions within the material. Those secondary charged particles are the contributors to the dose in the scoring volume.

MCNP calculations
The simulated incident neutron energies correspond to the energy list defined in ICRP116 (2010) from 1 meV to 15 MeV. MCNP version 6.2 (Werner 2017, Werner et al 2018 is used. The considered nuclear libraries are derived from ENDF/B-VII (Chadwick et al 2011) for neutrons and protons and from ENDF/B-VI.8 for photons. The electron transport algorithms and nuclear libraries are derived from the ITS code version 3 (Halbleib et al 1992). However, when nuclear libraries are missing or the particle energy is outside of the range of the nuclear libraries, dedicated models are used by the code. As recommended in Reed (2007), the S(α,β) for hydrogen (thermal scattering correction) in light water is taken into account. Hence, the 'lwtr.20t' library is used.
The simulation is done in 'MODE N P E H # D T S A' to take into account all the secondary particles created in the neutron transport process (neutrons, photons, electrons, protons, heavy ions, deuterons, tritons, helions and alpha).
Most of particle transport physics parameters are set to default values (Werner 2017). For neutrons and protons, the light-ion recoil parameter is activated with a value of 1.
Concerning the particle energies transport limits, except for neutrons, the low energy cutoff for all secondary particles is set to 1 keV (default value for photon and electron, lowest value for others particles). Below this energy cutoff, the particle is no more tracked (killed) and the remaining energy is deposited locally.
The '+F6' tally is used to calculate the deposited energy (in MeV g −1 ) from all particles in the scoring volume (see figure 1). This tally has no particle designation and automatically applies to all charged particles.

FLUKA calculations
FLUKA 2011.3 (Bohlen et al 2014, Battistoni et al 2015 is used with the FLAIR version 3.1-2 graphical interface (Vlachoudis 2009). The PRECISION default setting is considered and a calculation is performed for each group of the 260 neutron groups parametrized in FLUKA, from 1 meV to 15 MeV, as FLUKA does not handle pointwise cross-sections for neutrons. The geometry described in figure 1 is used. The source is a circular beam of 7 cm diameter. The disk is uniformly emitting unidirectionally and perpendicularly to the skin cube. Heavy ions, light ions, alpha, electrons, positrons, neutrons, protons and photons are transported. The energy cutoff for both the transport and production of electromagnetic particles is lowered to 1 keV. This means that, as far the simulation codes are concerned, the electromagnetic particles deposit their energies at the location as they reach 1 keV. Thermal scattering correction for H bound in water is applied during the calculations.
FLUKA uses Kerma approximation factors also called 'first-collision dose factors' to account for the first collision charged particles doses. More precisely, secondary particles originating from primary neutrons, other than neutrons, protons (except in two cases) and photons, are not transported in FLUKA (see FLUKA's manual 3 ). The two cases in which protons are transported are the following: • Recoil protons; • Protons from (n, p) reactions.
Recoil protons originate from the elastic collision of neutrons, in particular with a hydrogenated mixture. They can induce proton reactions leading to the creation of other tertiary charged particles. The dose of these tertiaries are not included in the FLUKA Kerma approximation as they are explicitly transported. Hence, they are scored under their respective particle type dose. The same applies for the tertiary charged particles originating from secondary photons.

GEANT4 calculations
GEANT4.10.07 toolkit (Agostinelli et al 2003, Allison et al 2006 is used to compute the LSD-CC. Two patches are applied: • patch 4 for exit direction of hadronic elastic collisions perpendicular to Z-axis, • patch 5 for inelastic excitation energy set to 0 when inelastic collision energy is below 20 keV. The same energy groups of the FLUKA simulations are implemented for the primary neutron generation in order to ensure the validity of FLUKA multigroup formalism. However, GEANT4 handles the pointwise neutron cross-sections. G4NDL4.6 library is used. The QGSP_BIC model is applied with neutron 'High Precision', to use pointwise neutron cross-sections for energies below 20 MeV. Note that in GEANT4, the Doppler broadening is calculated on the fly by sampling the velocity of the neutron-collided nucleus for each interaction. This process is done at a fixed temperature in the data files for FLUKA and MCNP. G4ThermalScattering and thermal scattering correction for H bound with water is implemented as suggested in Ribeiro and Souza-Santos (2017). The authors show that the High Precision physics list with the thermal scattering correction from chemically bound atoms is recommended for neutrons in the energy range 1 meV to 15 MeV. The same particle types as for FLUKA simulations are scored. No energy threshold is applied and the standard electromagnetic physics is used for secondary particles.

Rebinning of the 260-group calculations
The 260 bin results from FLUKA and GEANT4 are finally transformed considering the ICRP116 (2010) energy points assuming a uniform distribution of the dose within each bin. We consider g 1 , . . . , g 260 the 260 FLUKA groups and h 1 , . . . , h 47 (47 energy bins from 1 × 10 −09 to 15 MeV according to ICRP116) a binning consisting of the centers of ICRP116 energy points. In order to conserve the integral under the computed LSD-CC histograms, the dose in the new energy bin h i , h i+1 (denoted D [hi,hi+1] ) as a function of the computed dose with FLUKA binning g j , g j+1 (denoted D gj∈ [hi,hi+1] ) is calculated as follows (illustration with parameter description in figure 2): gl,gl+1] (1)

Results
In this section, we present the LSD-CC calculated using MCNP, FLUKA and GEANT4. A comparison is provided between the three code results and further discussion is presented focusing on the different particle interactions occurring in the skin. Figure 3 shows the results obtained with the three calculation codes. Computation times and processing hardware differ but the calculations are tuned in order to reach a final statistical uncertainty of less than 2%.

Calculations benchmark
The top left graph provides the LSD-CC values obtained with FLUKA and GEANT4 for the 260 groups and those obtained with MCNP considering the energy points of ICRP116. The three curves show quantitative and qualitative agreement. For the MCNP values, statistical uncertainties (at 1σ) are below 1% for the whole energy range. Whereas for FLUKA and GEANT4, the statistical uncertainties (at 1σ) are below 2% from 10 keV to 15 MeV and below 10% for energies below 10 keV. These statistical uncertainties are higher compared to those from MCNP. However, the rebinning step (see section 3.4) combines results from multiple energy groups. Hence, it leads to the reduction of the statistical uncertainties of the rebinned values. The results provided by the three codes are qualitatively consistent. Geant4 results seems to slightly overestimate the LSD-CC values, in particular, at low energies.
The top right curves show proton contributions to the LSD-CC values. As it can be seen, the protons represent the main contribution, to the LSD-CC, above 1 × 10 −2 MeV and below 1 × 10 −8 MeV. From 1 × 10 −3 to 1 × 10 −2 MeV, we notice a decrease of the FLUKA proton (without Kerma) contribution compared to GEANT4. This effect is caused by part of the FLUKA proton dose that is already taken into account in the 'neutron Kerma' as protons are produced by (n,p) reactions with the skin Nitrogen (seen in the Fluka_Kerma curve). The protons produced by this latter mechanism are not recoil protons. Hence, they are not transported in Fluka. They are in fact simulated by the Kerma approximation for such reactions.
The bottom left curves represent the electron contribution to the LSD-CC. Electrons are the main contributors in the energy range between 1 × 10 −8 MeV and 1 × 10 −2 MeV. We observe slight differences at low primary neutron energies (around 1 × 10 −9 MeV) and above 10 MeV between GEANT4 and the two other codes (FLUKA and MCNP). In these two regions, the LSD-CC values are slightly higher (∼20%) with GEANT4 simulations.
The ions (shown in the bottom right curves) become strong contributors to the LSD-CC at energies above 10 MeV. As expected, we observe a difference between the FLUKA and GEANT4 curves above 1 × 10 −3 MeV as part of the ion doses are integrated in the 'Fluka_Kerma' from Kerma approximation. Below 1 × 10 −3 MeV, ions are mainly produced by recoil protons whereas above 1 × 10 −3 MeV, they mainly come from neutron reactions that produce ions. That is the reason why they are included in the 'Fluka_Kerma dose' for this energy range.
Finally, we draw the attention of the reader that alpha particles, positrons and photons contributions have also been studied. They are included in the LSD-CC values (for each code result) with contributions with less than 0.1% to the LSD-CC. The slight discrepancies observed between the three codes could originate from: • The nuclear database cross-section used for the calculation, regarding neutron transport.
Indeed, MCNP calculations rely on ENDF/BVII for neutrons and ENF/BVI.8 for photons, GEANT4 results are produced with G4NDL4.6 which is mainly based on JEFF-3.3 and FLUKA nuclear data originates from different evaluation files such as ENDF/BVII, JENDL-4.0, JEFF-3.1, etc. As FLUKA is mainly based on ENDF/BVII, this could be the reason why the obtained results are closer to MCNP. • The physics models implemented for the transport of electromagnetic and charged particles.

Local skin dose conversion coefficients
The LSD-CC values, calculated in this paper, are illustrated in figure 4 with the provided energy points of ICRP116 (2010).
The LSD-CC, computed from FLUKA and GEANT4 with 260 groups, are rebinned and centered to match the ICRP116 energy point values, as explained in section 3.4. Due to the rebinning step, the statistical uncertainties are reduced to below 2% (at 1 sigma) as the new energy bins incorporate larger particle samples. The LSD-CC values are provided for the three codes in the table 2. The maximum relative difference is observed between MCNP and GEANT4 with a value of 19% at 10 keV. The minimum LSD-CC value is reached around 1 × 10 −3 MeV due to the transition between capture and inelastic reactions for increasing primary neutron energies as detailed in section 4.3.

Discussions
To better understand the obtained results and shape of the LSD-CC curves, we investigate the reactions induced by neutrons in a large volume of skin that could be considered as infinite for all practical reasons. The used skin material has an identical elemental composition described in table 1. The objective is to increase the neutron-skin reaction probability and hence to force all neutrons to interact at least once before leaving the volume, and avoid neutron leakage. Figure 5 lists the main reactions of neutrons in the skin as a function of their incident kinetic energy, produced with GEANT4 simulations. We can see that at low energies, until 100 keV, the dose is mainly driven by electromagnetic radiations from neutron capture in hydrogen and in nitrogen. The reaction with hydrogen leads to the emission of photons or electrons and deuterons as reaction products while the reactions with nitrogen produce protons. Above 100 keV, the dose is mainly driven by recoil protons originating from elastic scattering reactions (Baiocco et al 2016) (not seen in figure 5 which only lists first inelastic and capture reactions-we refer to capture reactions only the radiative capture in this paper,    Figure 5. Main first inelastic and capture reactions (obtained using GEANT4) from neutrons in the skin as a function of primary neutron kinetic energy. The above results are calculated based on an infinite skin volume.
following GEANT4 naming convention). While the neutron energy increases, neutrons have collisions with higher atomic number atoms such as nitrogen, oxygen and carbon, producing recoil particles with enough energy to induce ionizations. Above 5 MeV, we see a change of the dominating first reactions as neutrons interact firstly with carbon and oxygen, leading to the production of alpha, photons, electrons and secondary neutrons by inelastic threshold reactions. Despite collisions are mainly happening with hydrogen in skin at low energy, the main component of the deposited energy is the nitrogen below 100 eV, which is dominating (see figure 6).
We compute the inelastic scattering length and the ratio of neutrons undergoing inelastic or capture reactions in the infinite skin volume. The results are presented in figure 7.
We can see that at low primary neutron energies, up to 1 keV, capture reactions are dominating. As the energy increases, neutrons are more likely to leak from the ICRP skin volume since the inelastic scattering length to around 840 m at 100 keV. This explains the drop in the LSD-CC values around 1 keV. Above 100 keV, the inelastic reactions (excluding capture) become dominating as the inelastic scattering length is reduced.

Conclusions
We use three Monte-Carlo simulation codes (MCNP, FLUKA and GEANT4) to compute a set of 'LSD per Fluence' conversion coefficients (LSD-CC) for neutrons in the energy range of 1 meV to 15 MeV following the ICRP116 energy points. These factors are provided for the first time as they are not given in the literature. They are necessary to calculate the protection quantity of local skin dose induced by neutron exposure. A comparison between three code results is performed. The results are both qualitatively and quantitatively consistent. However, differences up to 19% are observed between GEANT4 and MCNP at 10 keV, 18% between GEANT4 and FLUKA at 500 eV and 15% between FLUKA and MCNP at 20 keV. GEANT4 results are slightly above the two others at energies below 1 × 10 −8 MeV. Explaining these discrepancies is not the scope of this paper but they could originate from the different nuclear databases used in the various codes and from the physics models (electromagnetic particle transport, Kerma approximation). The LSD-CC values are intended to be used in the framework of the IAEA working group A1/A2, whose mandate is to provide recommendation for the improvement of the regulations for the safe transport of radioactive material (class 7). These values could also be of interest for radiation protection and medical treatment applications. In the absence of experimental values, we are unable to recommend the use of one code over another. Hence, we recommend averaging the computed LSD-CC values of the three simulation codes.