Current atomic displacement models

The Threshold Displacement Energy (TDE) is defined as the minimum recoil energy required to create a stable point defect. The direction-averaged TDE, denoted by E_d in the present study, is widely used in analytical atomic displacement models. Using E_d, Kinchin and Pease (KP) got a formula for computing the number of atomic vacancies (denoted by ν in the present work) induced by a PKA with the kinetic energy of E_PKA as [12]

It is noteworthy that the cut-off energy E_c is not used in more recent models nor in the present work. The main reason for this can be found in Ref. [14].

Considering electronic energy loss and a more realistic atomic collision cross section, Norgett, Robinson, and Torrens (NRT) proposed a modified KP formula based on several Binary Collision Approximation (BCA) calculations [8, 15]:

where E_d is the effective energy for atomic motion, also called damage energy, first proposed by Lindhard et al. [16]. Figure 1 shows Lindhard’s partition function (i.e., P=E_D/E_PKA) with Robinson’s analytical fitting [17] for various monatomic materials.

Fig. 1Full-screen View

(Color online) Lindhard’s partition function for selected monatomic materials.

Because some displaced atoms are recombined before reaching thermal equilibrium, the Athermal Recombination-Corrected (ARC) model corrects the NRT model for E_d>2E_d/0.8-2.5E_d[13, 18, 19]. The athermal recombination of displaced atoms cannot be simulated by BCA codes. However, it has a quite limited influence for the low PKA energy, which is the case for the present work; thus, it is not considered here.

In the KP and NRT (or NRT-based) models, one can conclude that the effective variables are E_PKA/E_d and E_D/E_d, respectively. Therefore, the present study uses the energy normalized by E_d as an essential parameter to reduce the number of variables and simplify the comparison among different materials as well as the analysis. To simplify the expressions in the following discussion, let $\widehat{\nu }$ denote the number of atomic vacancies using the normalized energy as a unique parameter, i.e., ${\widehat{\nu }}_{\text{KP}}\left(E_{\text{PKA}}/E_{\text{d}}\right)={\nu }_{\text{KP}}\left(E_{\text{PKA}}\right)$ and ${\widehat{\nu }}_{\text{NRT}}\left(E_{\text{D}}/E_{\text{d}}\right)={\nu }_{\text{NRT}}\left(E_{\text{D}}\right)$.

SRIM simulations

In SRIM-like codes, there are four methods to obtain the number of atomic displacements: number of vacancies directly from BCA simulations (vacancies.txt for SRIM) and the value calculated using the NRT formula with the damage energy from the BCA simulations for both the QC and FC options. Since the FC option is more physically reasonable, the present work is based on FC.

Because the method of using damage energy is based on the NRT formula, the direct results from collision simulations should be more reliable. Conversely, Agarwal et al. [7] recently pointed out that the latter should be incorrect according to the details of collisions and recommended using the former method. Their reasoning is absolutely convincing. It is however surprising that the results obtained with the recently developed code Iradina are consistent with those of the SRIM FC [9, 20]. In addition, because the first method using damage energy is slightly different from the NRT model calculations, and the theories behind it are well understood, the present work investigates the number of atomic vacancies from SRIM-2013 FC using the vacancies.txt file, simply referred to as SRIM FC hereinafter.

In both the KP and NRT or most other models, it is a common conclusion or assumption that only one atomic vacancy is produced for PKA energy (or damage energy) larger than E_d but smaller than ~2E_d. For a high initial energy, the atomic displacements induced by a PKA or an incident self-ion are almost identical [11]. For an initial energy comparable with E_d, a PKA is very different from an incident self-ion. Therefore, PKAs, rather than externally incident ions in SRIM simulations, are used in the present work. The original position of the PKAs is set to the center of a 10×10×10 nm³ (or larger for a few high PKA energies) cube.

Because SRIM is a stochastic code, the convergence of the Monte Carlo simulations must be ensured. Figure 2 displays the number of atomic vacancies from the SRIM FC of the neutron cascade for 50 and 80 eV Fe PKAs in pure Fe. The study of the numerical convergence is performed on the grid of 2ⁿ PKAs. One can conclude that 8192 (= 2¹³) PKAs are reasonable to ensure the convergence of SRIM Monte Carlo simulations; thus, the following studies are based on 8192 PKAs simulations. Different from the assumption that ν=1 for E_d ≤ E_PKA < 2E_d, the SRIM FC gives ν>1 for E_d < E_PKA < 2E_d, which is achievable for atomistic simulations because the TDE is direction-dependent. It is much less evident in SRIM simulations owing to the amorphism of the materials. Nevertheless, this is consistent with the case of Ni studied by Weber and Zhang [10]. Agarwal et al. [7] believed that this is due to the incorrect count of some replacements as displacements.

Fig. 2Full-screen View

Number of atomic vacancies in Fe versus the number of simulated PKAs for 50 and 80 eV Fe PKAs.

Using SRIM FC, we again compare the displacement damage induced by a PKA and an externally incident ion. Figure 3 plots the number of atomic displacements in Si induced by PKAs and externally incident Si ions (coming from one side of the simulated cube) using the initial kinetic energy and corresponding damage energy as variables. It should be noted that the PKA-induced damage energies used in the present work are the PKA energies after subtracting the ionization energies stored in the IONIZ.TXT SRIM output file. SRIM FC confirms the non-negligible differences between the atomic displacements induced by a PKA and those induced by an incident self-ion when the energies are comparable with E_d. Accordingly, the neutron cascade option must be used to study atomic displacements versus PKA energy.

Fig. 3Full-screen View

Number of atomic vacancies in Si versus the PKA (black) or incident ion (red) energy (solid symbols fitted by solid lines) and the corresponding damage energy (center-dotted symbols fitted by dashed lines) with the unit of E_d. The grey plot is the linear fitting of the vacancies versus PKA energy without PKA itself (i.e., –1), it is shown for an intuitive comparison with the vacancies induced by an incident ion.

The present work is based on selected important monatomic materials because the current analytical formula is valid only for monatomic materials. Moreover, the materials are chosen to cover a wide range of atomic numbers (from Z = 6 to 74). Fe and Ni are widely used in stainless steel, Al is used in many fission reactors [21], C, Cu, and W [22] are used for fusion applications, and Si is a necessary element for semiconductors in various applications [23]. The average TDEs of the studied elements are given in Table 1. All binding energies are set to be 0 to study the analytical atomic displacement models.

Atomic displacements from SRIM FC

Figure 4 shows the number of atomic displacements for Fe PKA (i.e., using the neutron cascade option) up to 100 keV in pure Fe by SRIM FC. Both the values from the VACANCY.TXT file and those computed with the NRT formula using the damage energy from the SRM FC are illustrated. The NRT formula is also multiplied by a factor of 2 for an intuitive comparison. The results are quite similar to the case of Ni PKA in Ni shown by Weber and Zhang [10]. For damage energy above 2.5E_d, a typical discrepancy of a factor of about 2 is found between the SRIM FC and NRT calculation. Such discrepancy has been widely recognized and analyzed [5, 6, 9–11]; therefore, the present work does not emphasize this point.

Fig. 4Full-screen View

(Color online) Number of atomic displacements versus PKA energy for Fe PKA in Fe. Vacancy is the data taken from VACANCY.TXT, whereas NRT refers to the value computed with the NRT model using damage energy computed using SRIM FC.

In the range of E_d ≤ E_D < 2.5E_d, SRIM FC show that v>1 and is a strictly increasing function of damage energy, whereas the NRT model implies that v_NRT=1. The ratio of SRIM FC to NRT increases from ~1 to ~ 2 when E_D increases from ~ E_d to 2.5E_d. This region has received much less attention because it is not crucial for the displacement damage induced by reactor neutrons and ions. However, it has a large influence on that induced by light particles (e.g., electrons, positrons, photons) (see Refs. [28, 29] for example). Therefore, the atomic displacements for energy below 2.5E_d are yet to be studied.

Atomic displacements for E_D < E_d

As the results shown in Fig. 3, atomic displacements are observed when E_D < E_d but E_PKA ≥ E_d. This is a direct consequence of the definition of E_d: a PKA with E_PKA ≥ E_d is able to produce one atomic displacement (itself or a replacement) because the energy loss in inelastic collisions occurs after the displacement of PKA [29]. Therefore, for PKA energy comparable with E_d, it is questionable to use the damage energy as the effective energy for computing the number of atomic displacements.

The author has proposed a simple modification of the NRT (mNRT) model by assuming that v=1 when (E_PKA≥E_d) & (E_D < E_d)[29]. This modification has a limited influence on the quantification of the primary damage induced by neutron and ions irradiations, of which the contribution of high-energy PKA is predominant. However, its influence on the displacement damage induced by light particles can be considerable [29]. Compared with the original NRT model, the mNRT is more consistent with the SRIM simulations for PKA energy in the range of (E_PKA ≥ E_d) & (E_D < E_d). Nevertheless, they are identical for E_D > E_d; thus, the agreement with the SRIM BCA calculations is barely improved in general.

Atomic displacements between E_d and 2.5E_d

Once E_PKA is larger than 2 times the minimum TDE, denoted by E_d,min hereinafter, it is possible to produce two atomic displacements. Because the TDE is direction dependent, it is possible that ν > 1 for E_d < E_PKA < 2E_d. In fact, for some materials, E_d > 2E_d,min. For iron, Table 2 reveals that 9 of the 11 interatomic potentials used in Ref. [30] give E_d > 2E_d,min. Because the NRT model and NRT-based models use 2.5E_d as a demarcation energy to ensure continuity, we investigate the point defects for both PKA energy and damage energy in the range of ［E_d, 2.5E_d］. Because the present work and the SRIM code use only the average TDE, it is simply denoted by TDE hereinafter if without any other statement.

Figure 5 plots the numbers of atomic displacements for various monatomic materials with PKA and damage energies between E_d and 2E_d. For the sake of simplification, they are respectively denoted by ${\widehat{\nu }}_{\text{KP}}$ (left plot) and ${\widehat{\nu }}_{\text{NRT}}$ (right plot) with the variables E_PKA/E_d and E_D/E_d. It is obviously confirmed that ν > 1 for E_PKA>E_d for all monatomic materials. Excluding the material dependence already included in E_d (and E_D), $\widehat{\nu }$ is additionally dependent on the material. For the seven monatomic materials studied in the present work, Fe, Ni, and Cu follow almost the same law. ${\widehat{\nu }}_{\text{KP}}$ of C and W are quasi-identical but smaller than that of the other five. ${\widehat{\nu }}_{\text{KP}}$ of Si is between those of Al, Fe, Ni, and Cu and the ones of C and W. However, ${\widehat{\nu }}_{\text{NRT}}$ seems to be decreasing with the increasing atomic number of the target. The main reason is that the partition function is larger for heavier atoms [14, 16] (Fig. 1). It is noticeable that ${\widehat{\nu }}_{\text{NRT}}$ of C and W are quite different, whereas their ${\widehat{\nu }}_{\text{KP}}$ are quasi-identical.

Fig. 5Full-screen View

(Color online) Number of atomic vacancies versus normalized PKA (left, noted by ${\widehat{\nu }}_{\text{KP}}$ in the text) and damage (right, noted by ${\widehat{\nu }}_{\text{NRT}}$ in the text) energies from SRIM FC of neutron cascade. The straight lines are linear fittings.

The discrepancies shown in Fig. 5 can be attributed to the different materials and different TDEs. Therefore, Figure 6 compares the results of Si with two different values of E_d. It is noted that 24 eV is the average threshold energy for a bond defect or a Frenkel pair [26]. This value is comparable with E_d=21 eV obtained by Bourgoin et al. [31]. It can be concluded that the value of E_d influences both ${\widehat{\nu }}_{\text{KP}}$ and ${\widehat{\nu }}_{\text{NRT}}$, even though they are independent of E_d in typical models (cf. Sect. 2.1). In addition, comparing the results of C (Fig. 5) and Si (Figs. 5 and 6), ${\widehat{\nu }}_{\text{KP}}$ and ${\widehat{\nu }}_{\text{NRT}}$ of Si with E_d =24 eV are larger than those of C, of which E_d=25 eV. Consequently, $\widehat{\nu }>1$ for E_PKA > E_d and $\widehat{\nu }$ depends on both the material and value of E_d. Therefore, we cannot obtain a unique simple function of E_PKA /E_d or E_D/E_d to describe the number of atomic displacements from the SRIM FC in the range of (E_PKA≥E_d) & (E_D<2.5E_d). Nevertheless, as the results shown in Figs. 5 and 6, a linear fitting of the number of atomic displacements vs. PKA energy or damage energy is suitable for each monatomic material.

Fig. 6Full-screen View

(Color online) Number of atomic vacancies versus the normalized PKA (solid points) and damage (center-dotted points) energies for Si with E_d = 36 eV (black squares) and 24 eV (red circles). The straight lines are linear fittings.

Correcting the damage energy in the displacement calculation

Robinson and Oen [32] recognized that the inelastic energy loss for atoms with kinetic energy smaller than 2.5E_d does not influence the number of atomic displacements. Thus, the inelastic energy loss when an atom slows down from 2.5E_d to 0 should be added to the damage energy for computing the atomic displacements [32]. Based on this reasoning, they obtained the effective energy as [32]

It is noteworthy that the demarcation of 2.5E_d in the NRT formula is used only to ensure the continuity of the displacement function. According to the reasoning of Kinchin and Pease [12], 2E_d is a physically crucial limit. In addition, an atom does not slow down with continuous energy loss. A collision may decrease the energy of an atom from E₁>2E_d to E₂<2E_d. Assuming the equiprobable energy distribution (i.e., hard-sphere collision [11]) for an atom slowed down to E<2E_d for the first time, one can introduce a correction factor by

Because the partition function can be considered as quasi-constant for the PKA energy from 0 to 2E_d, $E_{\text{D}}\left(E_{\text{PKA}}\right)\approx E_{\text{PKA}}\times P\left(E_{\text{d}}\right)$ when 0≤E_PKA≤2E_d. Therefore, the correction factor can be approximated using

Using Lindhard’s analytical partition function for monatomic materials, ηcan be simply calculated as

For atoms from Li to U with $E_{\text{d}}$ of several tens of eV, η≈1.2. This value is in good agreement with the experimental values of Fe and Ni summarized in the Nuclear Energy Agency (NEA) report [33], the experimental results of Cu obtained by Averback et al. [34], and many MD simulations for energy around 2E_d.

It is noticeable that the correction factor proposed by Robinson and Oen [32] is numerically close to the present one because the partition function varies insignificantly between E_d and 2.5 E_d. The difference is only a factor of 1.16 in the second term of η. Therefore, η≈1.2 for both corrections. This value leads to ν_NRT(E_D) ≈E_D/2E_d for E_D>2E_d. One obtains exactly the same formula as the KP formula by replacing of the PKA energy with the damage energy.

However, it is noticeable that κ≈0.86 or 0.8 in the formula ν_NRT(E_D) ≈κE_D/2E_d is determined by fitting the BCA calculation results [15]. Therefore, if the effective energy E_D,eff=ηE_D rather than E_D is used, the fitted constant (or widely recognized as the correction to the hard-sphere collision cross section) becomes κ’=κ/η≈0.7.

Effective energy for SRIM simulations

Because Lindhard’s partition function is slightly different from that computed by SRIM, the present correction factor for damage energy is calculated with E_D (E_d) from SRIM FC and denoted by η_SRIM. The effective energy is computed as follows:

The results shown in Fig. 5 are rescaled by η_SRIM and are illustrated in Fig. 7. Excluding W, the six monatomic materials considered in the present work (from C to Cu) almost follow the same law, which uses E_D,eff/E_d as a unique variable. Therefore, this effective energy provides the possibility to correct the current formulae and describe the results of SRIM FC. Moreover, the use of η satisfies ν(E_PKA=E_d)=1, which is consistent with the definition of E_d.

Fig. 7Full-screen View

(Color online) Number of atomic displacements versus the effective energy from the SRIM FC of neutron cascade. The straight line is the linear fitting of the six cases excluding W (R² = 0.994).

The fitted line in Fig. 7 indicates that ν_SRIM≈0.7 (E_D,eff/E_d-1)+1. Assuming its applicability up to high energies, once $v\gg \text{0}\text{.3}$, ν_SRIM≈0.7 E_D,eff/E_d=0.7η_SRIM/E_D/E_d. For most cases where η_SRIM≈1.1 or 1.2, ν_SRIM≈0.8E_D/E_d=2ν_NRT. This relationship explains the typical discrepancy of a factor of about 2 between the SRIM FC and NRT formula.

However, it should be noted that the number of atomic displacements in W still differs from the others. Moreover, the use of η_SRIM cannot make the two curves of Si in Fig. 6 coincide. In fact, it is important to indicate that the damage energy from the SRIM FC depends on the value of E_d (e.g., the example on Si shown in Fig. 8), whereas Lindhard’s damage energy is independent of E_d [16]. Using the data for E_d =36 eV as a reference, we rescale the effective energy for E_d =24 eV by a factor of 1.1 to get the same damage energy function. The rescaled data are plotted in Figure 9 together with the data versus the damage energy and effective energy. Rescaling the effective energy to eliminate the bias induced by E_d results in similar atomic displacements for Si with E_d =24 eV and E_d =36 eV.

Fig. 8Full-screen View

(Color online) Damage energy versus PKA energy from SRIM-2013 FC for Si with E_d = 36 eV (black) and 24 eV (red).

Fig. 9Full-screen View

(Color online) Number of atomic displacements versus the normalized damage and effective energies for Si using E_d = 24 eV (circle) and 36 eV (square). The straight line is the linear fitting of two data sets (R² = 0.997).

From the case of Si shown above, one can find that none of E_D,eff/E_d, E_D/E_d, and E_PKA/E_d can be the unique variable for describing the number of atomic displacements from the SRIM FC. Therefore, the number of atomic displacements versus E_PKA/E_d, E_D/E_d, or E_PKA/E_d for two arbitrary monatomic materials are not necessarily the same. The only general conclusion is that linear fitting can be used to describe the number of atomic displacements versus E_D,eff, E_D, and E_PKA for energy comparable with E_d. One can also adopt a specific value of E_d for W to reduce the difference with the other materials, as shown in Fig. 7. The difference decreases but still exists using of E_d = 55 eV [35, 36] rather than 90 eV. The difference can be further decreased by decreasing E_d; however, using an unphysical value for E_d is unnecessary.

Further comments on the damage energy

Figure 8 shows that the damage energy from SRIM FC depends on the value of E_d, whereas the Lindhard’s damage energy is independent of E_d. However, these are not physically incompatible. In fact, the original equation governing the damage energy (or atomic vacancies) includes E_d as a basic parameter [16]. The currently used Lindhard’s damage energy is TDE-independent because it is obtained according to the numerical solutions after removing E_d in the original equation (i.e., their approximation (B): E_d is negligeable when compared with kinetic energies of atoms [16]).

Table 3 summarizes the damage energies and numbers of atomic displacements for Fe PKA in Fe from SRIM FC using two different values for E_d: 40 eV and 20 eV. The results show that the damage energy depends on E_d for PKA energies up to 100 keV. For a given PKA energy, the damage energy is larger for a larger E_d. This is a consequence of knocked-on atoms having a smaller kinetic energy for larger E_d. A lower kinetic energy results in lower inelastic energy losses in subsequent collisions. Therefore, once E_d is changed, the corresponding number of atomic displacements cannot be directly predicted as the inverse proportion to E_d. Taking the Fe PKA in Fe shown in Table 3 as an example, the number of atomic displacements is reduced by a factor greater than 2 (the last column in Table 3) if E_d is doubled (20 eV → 40 eV), whereas the NRT model predicts a reduction of a factor of 2 (or smaller than 2 if the slight increase in damage energy is considered). This confirms the conclusion given in Sect. 4.2: E_D/E_d or E_D,eff/E_d cannot be the unique variable for computing the number of atomic displacements. A variation of E_d by a factor of x does not imply a variation of a factor of 1/x for the number of atomic displacements.