Comparison of Monte Carlo Tools for Displacement Calculations for Protons and Heavy Ions Irradiation on Iron and EUROFER97

Investigation of neutron-induced damage is a field of research relevant for present and future designs of nuclear reactors. Neutron exposure is known to induce degradation in the mechanical properties of reactor components [1] and an evaluation of the amount of this phenomenon is crucial to determine operational safety limits, lifespan, and life extension of nuclear facilities. While radiation damage is normally studied through neutron irradiation in nuclear or material test reactors, an alternative approach offering several advantages [2] consists of performing ion irradiations at accelerator facilities [3]. Benefits of ion irradiation include the possibility to reduce the activation of the specimens and shorten the irradiation time at a fraction of the cost of an equivalent neutron irradiation one. On the other hand, due to the shallow penetration of charged particles into matter, the investigation of the induced ion damage requires micro-mechanical and surface-sensitive testing methods. For example, the ion range in the case of this research spans from one to a few tens of micrometers: this shallow layer is the region of interest where mechanical properties are expected to change.

The designing phase of an ion irradiation experiment is usually based on displacement calculations carried out with Monte Carlo codes. Since the codes capable of performing the transport of energetic ions and the related damage profiles are numerous, few comparative studies aimed at testing results consistency have been carried out by different authors [4–6]. In such studies, though, displacement calculation performance has been mainly tested for materials and ion energies common in high-energy physics. Simulations of nuclear reactor environments require irradiations with incident particles of considerably lower energy; hence, the agreement between the different codes may not be valid anymore. In this research, we present a comparative study, focused on displacement calculations, with the widely popular general-purpose Monte Carlo codes cern-fluka 4-2.1 [7,8], mcnp 6.1 [9], phits 2.88 [10], and srim 2013 [11]. These codes have been selected due to their broad adoption for ion transport calculations in matter, while also being easy to obtain or readily available to the paper's authors.

The concept of displacements per atom or dpa was developed from the theories characterizing the generation of primary point defects in materials [12]. A displacement intuitively happens when an atom is moved from its lattice position to an interstitial one, creating a Frenkel pair or primary point defect. However, due to diffusion and recombination phenomena, displacement calculations should not be considered directly representative of damage production but rather as a radiation exposure measure [13]. To be representative of the actual induced damage, and not just the instantaneous creation of defects, damage profiles should be combined with models considering the thermal diffusion of vacancies and interstitials, the presence of defect sinks, and all kinds of mechanisms that influence the defects' concentration, which cannot be done by the selected codes. The advantage of the dpa quantity over the simple energy deposition (dose) is that, by definition, it calculates out the energy dissipated through electronic stopping mechanisms, which are known not to cause damage in most metallic materials. Moreover, employing dpa offers an advantage over fluence units, since comparison of the dpa estimates enables the correlation of irradiations performed on different materials, incident beam energies, and different kinds of particles, such as electrons, ions, and neutrons, as long as the main damage production mechanism is due to atomic recoils [14].

Displacement calculations play a key role in estimating the effects related to the exposure of components in nuclear environments, especially so for fusion reactors, where dpa levels are foreseen to be ten to hundred times higher than in conventional fission reactors [19]. This aspect prompted the need to develop alloys capable of retaining mechanical properties in the extreme fusion reactor environment while trying to reduce the activation and radiation-induced swelling effect [20,21]. The target material choice for the case study was the EUROFER97 steel, due to its relevance in irradiation studies as a structural alloy for fusion reactors [22]. The beam energy of past studies (5 MeV) regarding the irradiation of EUROFER97 with heavy Fe ion beams has been adopted as a starting point for this work [23]. For the case of proton irradiation instead, the beam energy has been selected as close to the highest value which also avoids target activation above the free handling limit [24,25]. Beam energies for both types of incident particles are summarized in the first line of Table 2. Displacement calculations have also been performed with pure iron G379 as the target material, to provide the reference results needed to investigate density and composition effects. In the simulations, the original composition of EUROFER97 has been reasonably simplified by considering only alloy elements present above 1% by weight. Relevant material characteristics of EUROFER97 and G379 are displayed in Table 1, including the ASTM-recommended values for the displacement energies [18].

Two sets of simulations have been carried out for each target material with the four codes. The first set, performed at a single beam energy, simulates the damage profiles obtained with Fe ions and protons beams; the second set estimates the (average) total number of displacements induced by an incident particle in the target materials, for seven increasing beam energies. The setup parameters shared among the four codes can be found in Table 2 for the damage profile simulations and in Table 3 for the total induced displacements versus energy simulations.

For the cern-fluka, mcnp, and phits codes, the incident beam has been shaped with a circular cross section diameter (d) of 5 mm, together with a flat energy and current profile along the beam radius. The detector region, responsible for measuring the dpa quantity, is positioned within a larger cylindrical target with a d = 6 mm and a thickness (t) of 1 mm. This setup enables the recoil cascades to scatter within the detector's boundaries. For the damage profile runs, the detector is a cylinder with a diameter of 5 mm, and a depth slightly larger than the particles' range in the material, ensuring that the displacement profile remains fully enclosed within its geometry (see Table 2).

For simulations comparing total displacements against energy, where spatial details of displacement generation are unnecessary and only a collective value of displacements matters, we have employed an infinite geometry approximation. This approximation assumes detector and target dimensions that significantly exceed the maximum range of incident particles (R_p) in the target material, ensuring complete energy deposition across the geometry for all beam energies. For cern-fluka, mcnp, and phits codes, displacement calculation output files are in units of dpa/incident ion. To convert them to the (average) total number of displacements/incident ions, it is essential to multiply the output results by the number of atoms within the detector volume. However, the srim code provides this value by default in the TDATA output file.

In all simulation runs, 10⁵ primary particles (N_pp) were utilized, resulting in a relative error of the results (E_rr) approximately equal to 0.003, with $Err=1/N$ being the typical relative error in a Monte Carlo simulation based on the number of trials N.

Tables 1–3 present the parameters utilized in running the simulations, which are shared across all four codes. The following subchapters will outline software-specific setup parameters and provide recommendations on their implementation.

In srim runs, the protocol outlined by Stoller et al. [14] was adhered to. This involves utilizing the “Ion Distribution and Quick Damage Calculation” mode instead of the “Full Cascade” mode (FC), setting surface and lattice binding energy to 0, and applying the NRT-dpa formula to the collisional energy distribution obtained. This approach yields the number of displacements, which can be converted—if necessary—to dpa values by normalizing with the appropriate number of atoms in the irradiated volume. Stoller recommends this approach for ensuring consistency with the NRT model's theory and for mitigating overestimations observed when using the FC calculation mode [14].

mcnp code is a software especially devoted to the simulation of neutron transport; to use this code for charge particle transport, some adjustments were needed to obtain sound results:

1. In mcnp, proton calculations have been possible through the adoption of displacement cross sections developed by the Institute of Neutron Physics and Reactor Technology at the Karlsruhe Institute of Technology (KIT) [26]. Unlike the other codes, which are provided with their own implementation of the displacement cross sections, in mcnp, it is necessary to outsource them. Displacement cross sections are a crucial element in displacement calculations since their convolution with the particle fluence yields the number of displacement estimate. The KIT displacement cross section was supposed to adhere to the format of the mcnp dosimetry libraries but their compilation was found to be incorrect; hence, they could not be directly used as tally modifiers (FM card). A workaround for this problem was found by introducing the cross section data directly in the input file as an energy-dependent point-wise response function (DE and DF cards) of an F4 tally. Since the KIT displacement cross sections were available for protons only, Fe ion displacement calculation simulations have not been performed with mcnp.

2. The mcnp simulation results, which the code provides in terms of dpa for each element of the compound separately, have been combined to obtain the total number of displacements in the irradiated material following the procedure also known as Bragg's rule [27], as shown here below

   $N=10−24∑jcjNj$

1. For the proton beam simulations, the proton energy cutoff has been set equal to 1 keV for consistency with all the other software and with the calculated cross sections that are not measured below this threshold.

1. In the cern-fluka code, the DEFAULTS card recalls specific configurations optimized for different uses. This card was set to DAMAGE since is the developers' recommendation for dpa calculations [8]. A choice that enables the transport of charged particles down to the lowest admissible value (1 keV/n for hadrons and light nuclei). The DAMAGE configuration also enables simulation of the energy straggling phenomenon.

2. To obtain dpa values coherent with the NRT model, the DPA-NRES scoring mode was selected.

3. The fluka code expects a single value of the displacement energy (Ed) for the target material, so for the EUROFER97 alloy, a suitable average one was provided by weighting the compound elements' displacement energies with the respective atomic fractions in the alloy, as reported here

   $Ed=∑jcjEd,j$

4. c_j is the atomic fraction and E_d,j is the displacement energy of the jth component of the alloy. By this formulation, one obtains a value of 40.17 eV, close to the one of pure iron, owing to the limitedness of the Tungsten amount, the only alloy component with a different value of E_d (see Table 1). A different formulation for the average displacement energy, derived from experimental testing, can be found in Ref. [28]. For this study's case, adopting either one of the formulations yields almost equal values.

1. In the [P a r a m e t e r s] card, the energy straggling for charged particles was enabled by setting nedisp = 1. By setting this value to 1, the straggling model implemented according to the Landau and Vavilov theory is chosen.

2. In the [P a r a m e t e r s] card of the input file, it is possible to lower the energy transport threshold. For proton runs, it is set to a value (1 keV) which is consistent with the transport model validity and the other software here used. For Fe ion irradiations, the emin(19) parameter, corresponding to the nucleus (transport) cutoff energy, was set to 0.001 MeV/n. Leaving it to the default value would cutoff the incoming ion right at the start, producing no dpa whatsoever.

3. The [T-D P A] tally, which performs the displacement calculation according to the NRT model, has been chosen.

Damage profiles display the distribution of the primary point defects induced in the target material and will be plotted in terms of dpa/incident ion as a function of the particle penetration depth. To correctly read the graphs it is necessary to point out that the dpa levels are obtained from the normalization of the number of displacements in the segment of the discretized domain by the number of atoms included in the corresponding node. For all the simulations a discretization of 100 nodes has been adopted, meaning that the number of atoms in each segment is equal to 1/100 of the total number of atoms in the whole irradiation volume.

Figures 1(a) and 1(b) show the distribution of the induced damage due to 5 MeV protons for the EUROFER97 and the G379 target materials respectively. The codes seem to be in accordance when foreseeing the shape of the profile, with all the trends characterized by an initial low damage region (due to the protons being still very fast and thus interacting mainly by electronic collisions) and a sharp displacement peak located at the end of the particle track (nuclear collisions are more probable at lower energies). Points of agreement between the codes can be observed in the form of a coherent prediction of the depth of the peaks, suggesting a similar implementation of the electronic stopping power, which is the main mechanism that defines the range of the protons. However, inconsistency in computing the peaks, and in general, the trend-line heights, is noticed, with values differing up to around 100% between the lowest mcnp and highest cern-fluka simulated dpa levels. This could be due to differences in the values of the scattering cross sections adopted by the codes which affect the fraction of the incident particle energy that goes into the generation of displacements [29].

When comparing Fig. 1(a) with Fig. 1(b) some effects due to the different material properties can be noticed. Displacement peaks are located at slightly higher depth in EUROFER97 with respect to G379, owing to the difference in the density values of the two target materials. Since the stopping power is a linear function of the atomic density, this reflects in a shorter range of the beam particles in the denser target material. Accordingly, assuming that for the two targets, the protons deposit very closely the same amount of energy into recoils, the density difference also translates into taller trend lines for the G379 target. An additional factor contributing (though in a very limited way) to the lower dpa levels for the EUROFER97 target is the higher value of E_d, owing to its Tungsten component, which reduces the number of displacements produced by a given collisional energy.

Analogous considerations can be drawn for the case of 5-MeV Fe ions, with Fig. 2(a) corresponding to the profile induced in the EUROFER97 target and Fig. 2(b) in the G379 one. Results for the mcnp code are not present due to the lack of cross sections to perform the displacement calculation. In Fig. 2, the profiles show a smoother trend with respect to the ones obtained in Fig. 1 for proton irradiations. This is related to the more pronounced fluctuations in energy (Landau–Vavilov straggling) that characterize the slowing down of heavier ions in matter. Still, while peak depths and profile shapes show marked differences, dpa levels can be seen going to zero at a depth of around 1.8 μm for all three codes, meaning that the codes output coherent predictions of the Fe ions range in the targets.

With respect to proton irradiations of Fig. 1, in Fig. 2, the effects due to target density difference (different peak height and particle range values) are less noticeable. The electronic stopping power (which is proportional to the target density) has in fact a lower stake in slowing down energetic heavy ions with respect to lighter protons.

The sharp profile associated with the fluka code instead is to be related to the code's proprietary implementation of the energy straggling phenomenon. To highlight this, in Fig. 2(a), the profile from fluka is compared to two profiles from phits code, corresponding to two different runs where the energy straggling package has been enabled and disabled. In phits, the package is implemented according to the Landau–Vavilov theory, and its activation is recommended by the developers [30]. On the other hand, fluka relies on an original approach to the straggling problem and has been run with the DAMAGE configuration. The DAMAGE mode is the developer's recommendation for displacement calculations, which among several things also enables the simulation of the ionization fluctuations (i.e., energy straggling) [31]. It can be noticed the that fluka-DAMAGE profile is quite similar to the one obtainable with phits if one were not to activate the Landau–Vavilov energy straggling package, suggesting that the implementation of the straggling phenomenon is optimized for much higher energies [32].

In Fig. 3, it is possible to observe the influence of mass and energy of the beam's ions on the total displacement production. For each code, seven lumped (i.e., with no spatial distribution) displacement calculations have been performed with beam energies ranging from 100 keV to 100 MeV. These calculations have been carried out with both proton and Fe ion beams and results are all shown in Fig. 3. Since increasing the beam energy will intuitively produce a larger number of displacements, monotonically increasing trends were expected to show up. This happens with all the codes, except for Fe ions with phits, which show a decrease at about approximately 5 MeV bombarding energy.

In the case of proton irradiations toward higher energies, the dpa growth is less steep for the srim simulations with respect to the other codes and could be attributed to the lack of modeling of nuclear reactions (non-elastic reactions) which might contribute to the increase of displacement production through the release of extra energy as kinetic energy of secondary products (Q-factor). On the other hand, in the case of Fe ions, no nuclear reactions are produced, as they are not energetically viable.

The Fe ion trends show a decrease in steepness when moving toward higher beam energies, almost reaching a saturation level. This can be explained by the fact that the ratio between nonionizing energy loss to the total energy loss of a particle impinging on a target is strongly energy dependent. The nonionizing energy loss fraction is defined by Lindhard's partition function and corresponds to the energy available for damage production (also termed damage energy). Given the damage efficiency ratio T_dam/E_PKA, where T_dam is the damage energy and E_PKA is the energy of the primary knock-on atom, according to Lindhard's theory, the ratio decreases as E_PKA grows [29]. For example, in the case of 500-keV Fe ions on EUROFER97, about 50% of ion's kinetic energy is converted into T_dam, with this fraction reducing to 1% for a 100-MeV ion.

Overall, the Fe ions induce several more displacements per incident particle with respect to protons, due to the more favorable mass ratio between the target atom and the incident ion [29], (ratio close to 1:1 for Fe ions on Eurofer97 or G379, while ratio approximately 60:1 for protons on Eurofer97/G379).

A second plot for G379 target material has not been displayed, since, as demonstrated in the next section, an almost identical figure would be obtained.

A final set of plots in Figs. 4 and 5 is meant to investigate differences in the generation of displacements due to the composition of the G379 and EUROFER97. The total number of displacements induced by a proton and an Fe ion has been calculated for seven beam energies, from 100 keV up to 100 MeV.

The results are depicted in such a manner that the (x,y) coordinates represent the displacements generated at a fixed energy, with EUROFER97 simulation results plotted against those of G379. This visualization allows us to observe whether certain materials produce a higher number of displacements (and thus dpa) under the same irradiation conditions compared to others. Figures 4 and 5 illustrate that simulation points for all codes align closely with the 45-deg line. Due to the similarity in composition between the two target materials, only minor deviations from the 1:1 ratio of induced displacements are expected. Indeed, some slight deviations from the line have been observed, albeit too small to be discernible on a logarithmic graph. Such deviations have always kept equal to about a 2% value of relative difference, and just in one case went up to 5%. In other words, simplifying the composition of an iron-based alloy (such as EUROFER97) by replacing it with a pure iron, will only cause a minor approximation in terms of estimated dpa, as long as the alloy and iron share similar average atomic mass number and average displacement energy. Moreover, the density parameter, which might differ considerably between different iron-based alloys, will not affect the generation of displacements (and so the total dpa estimate) but only their spatial distribution. This approximation could be useful when studying the radiation damage of an iron-based alloy of unknown composition or if one wants to speed up the setup of a simulation by avoiding specifying all the alloy components of the target materials.

The study here presented analyzed the simulation results without discussing the root causes of disagreements between the codes: to perform this quite challenging task, one would need the access to the implementations of the physical models in the source codes, which are not released to the public. Moreover, the absolute accuracy of the results can only be checked with the feedback from empirical studies, and the majority of those found in the literature focus on a highly different energy range or on irradiations in which the temperature plays a role in the damage and displacement profile, for these reasons results cannot be extrapolated to our own case [33,34].

Some lessons on the code performances can still be learned:

• Although srim has traditionally served as the benchmark code for such studies, it possesses both limitations and advantages. Notably, srim lacks the capability to simulate nuclear reactions, leading to an underestimation of the dpa at energies pertinent to these reactions. As the contribution of nuclear reactions to the dpa escalates with increasing energy, so does the extent of underestimation. However, srim incorporates tuning parameters, such as surface and lattice binding energy, which are not present in other codes. These parameters enable more precise estimates of displacements.

• mcnp relies on external cross sections for displacement calculations involving ions, restricting its applicability to the data provided by other studies. Consequently, in this work, it was unfeasible to conduct Fe ion displacement calculations using mcnp due to the unavailability of the required displacement cross sections.

• The energy straggling implementation in fluka for heavy ion simulations appears to be less refined compared to other codes (see Fig. 2(a)), particularly at the energies selected for this study. It's important to consider this factor if a more precise distribution of displacements is desired.

• The phits code exhibits a tendency to underestimate displacement generation when subjected to impinging Fe ions with energies exceeding 5 MeV. The authors were able to pinpoint a more precise switching value, which stands at approximately 2.75 MeV. This underestimation becomes more pronounced with higher impinging energies.

• Each of the four codes provides a distinct estimate of the dpa for identical irradiation conditions. For instance, if one were to plan an experimental study targeting 10 dpa with a 1-MeV proton, on an EUROFER97 sample, the resulting irradiation times would vary significantly, as demonstrated in Table 4. Without conducting a refined empirical study, it remains uncertain which result is the most accurate. However, it is still feasible to consider the differences using scaling coefficients.

A comparative study between the fluka, mcnp, phits, and srim Monte Carlo codes, aimed at evaluating the respective capabilities of performing displacement calculations, has been presented. The study focuses on proton and Fe ion irradiations of EUROFER97 and G379 specimens, two popular irradiation techniques that try to replicate the material damage induced in reactor environments. The simulation setup has been described, together with the additional specific software's recommendations needed to get coherent results among the codes.

Analyzing the software's output, several differences have been highlighted (e.g., displacement peaks' heights, sharpness and locations, particle straggling influence, particle penetration depth). Moreover, the codes exhibit discrepancies in calculating displacement estimates, with variations sometimes exceeding a factor of two or more. It has been demonstrated that different irradiation times can be obtained from various codes when aiming for a fixed required damage level. Without an empirical study encompassing the same range of particle energies, there is no evidence indicating which code is the most reliable. Such a study would be highly beneficial for the scientific community and the organic progression of this research. The four codes coherently showed that the chemical composition of the two different iron-based alloys just marginally affects the generation of displacements, implying that in a dpa calculation replacing the alloy composition with just its main constituent (i.e., iron) the error one commits is limited, and in any case smaller than the discrepancies among different codes.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.