Low-energy atomic displacement model of SRIM simulations

Sheng-Li Chen

Abstract Radiation-induced atomic displacement damage is a pressing issue for materials.The present work investigates the number of atomic displacements using the Primary Knock-on Atom(PKA)energy EPKA and threshold displacement energy Ed as two major parameters via lowenergy SRIM Binary Collision Approximation(BCA)full cascade simulations.It is found that the number of atomic displacements cannot be uniquely determined by EPKA/Ed or ED/Ed(ED refersto thedamageenergy)when theenergy is comparable with Ed.The effective energy ED,eff proposed in the present work allows to describing the number of atomic displacements for most presently studied monatomic materials by the unique variable ED,eff/Ed.Nevertheless,it is noteworthy that the BCA simulation damage energy depends on Ed,whereas the currently used analytical method is independent of Ed.A more accurate analytical damage energy function should be determined by including the dependence on Ed.

Keywords Atomic displacement.Damage energy.Effective energy.SRIM neutron cascade simulations

1 Introduction

Radiation damage is an important issue for materials because it changes the properties of materials(e.g.,radiation-induced segregation,swelling,hardening,and variation of resistivity[1]).Primary radiation damage,widely known as atomic displacement damage,is fundamental for studying the irradiation effect on materials.The number of atomic Displacements per Atom(DPA)was proposed to unify the damage caused by different irradiation sources,such as neutrons,protons,photons,electrons,positrons,and ions.DPA is now used as an essential parameter for studying the secondary(or macroscopic)radiation damage of materials[2].

In the past decades,various methods and computational codes have been developed to calculate the number of stable atomic displacements(i.e.,Frenkel pairs).The Stopping and Rangeof Ionsin Matter(SRIM)Monte Carlo code[3](or the former TRansport of Ions in Matter(TRIM)code[4]included in)iswidely used for computing atomic vacancies because of its ease of use with a userfriendly graphical interface.Presently,it is almost a nonofficial standard step to use SRIM simulations for ion irradiation studies.

Regardless of the discrepancy between the Quick Calculation(QC)and Full-Cascade simulations(FC)[5,6],SRIM ispowerful for modeling ion irradiation in materials.However,it cannot be directly used to compute the displacement damage induced by other radiation sources.To unify the atomic vacancies induced by irradiation with different sources,simulation tools or models using the kinetic energy of the Primary Knock-on Atom(PKA)as a major parameter should be used.Therefore,the so-called neutron cascade simulation has been implemented in the SRIM code[3].With the neutron cascade option,users identify the position and energy of PKA in a separate file(i.e.,TRIM.DAT)to compute the atomic vacancies induced by a PKA rather than an external ion.For irradiation sources other than ions,the PKA spectra can be determined using specific calculation tools and combined with neutron cascade simulations to obtain the number of point defects.Therefore,consistent results can be obtained for the number of atomic vacancies induced by different irradiation sources with SRIM simulations.

However,the numbers of atomic vacancies(from vacancies.txt file)obtained using the two basic options,QC and FC,of SRIM-like codes differ by a factor of about 2[5,6](1.0–2.7 depending on the incident ion and target atom[7]).Stoller et al.[6]recommended the useof the QC option for obtaining results comparable with the Norgett-Robinson-Torrens(NRT)[8]model.Recently,Crocombette and Van Wambeke[9]and Weber and Zhang[10]recommended using the FC option,especially for compound materials,because it is physically more reasonable.Weber and Zhang[10]and Chen and Bernard[11]explained that the discrepancy between QC and FC is due to the displacements induced by low-energy atoms,which are considered unable to produce additional displacements in classical models[12,13](detailsare given in Sect.2.1).Nevertheless,the residual energy transfer effect[11]is not evident for a low initial energy.Consequently,further studies on the atomic displacement model for low PKA energy are crucial for unifying the displacement damage induced by different radiation sources,especially for light particles such as electrons,positrons,and photons.

In addition,because the atomic displacements at low energies have an integrated effect on the model at high energies[10],studies on low-energy cascades can reveal thefeaturesof theprimary radiation damageover theentire energy range.Therefore,the present work investigates the atomic displacements mainly based on SRIM-2013 FC simulations at low PKA energies.The simulation methods are described in Sect.2.2.The simulated results and corresponding discussion are presented in Sect.3.Detailed discussion and comments on the use of damage energy for quantifying displacement damage are presented in Sect.4.The main conclusions of the present study are summarized in Sect.5.

2 Current atomic displacement models and SRIM simulations

2.1 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 Edin the present study,is widely used in analytical atomic displacement models.Using Ed,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 EPKAas[12]

It is noteworthy that the cut-offenergy Ecis 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 KPformula based on several Binary Collision Approximation(BCA)calculations[8,15]:

where EDis 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=ED/EPKA)with Robinson’s analytical fitting[17]for various monatomic materials.

Because some displaced atoms are recombined before reaching thermal equilibrium,the Athermal Recombination-Corrected(ARC)model corrects the NRT model for ED＞2Ed/0.8=2.5Ed[13,18,19].The athermal recombination of displaced atoms cannot be simulated by BCA codes.However,it hasa quitelimited influence for thelow 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 EPKA/Edand ED/Ed,respectively.Therefore,the present study uses the energy normalized by Edas 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,letdenote the number of atomic vacancies using the normalized energy as a unique parameter,i.e.,

2.2 SRIM simulations

In SRIM-likecodes,thereare four methodsto 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 themethod of using damage energy isbased 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 Edbut smaller than～2Ed.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 Ed,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 nm3(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 thegrid of 2nPKAs.One can conclude that 8192(=213)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 Ed≤EPKA＜2Ed,the SRIM FC gives ν＞1 for Ed＜EPKA＜2Ed,which is achievable for atomistic simulations because the TDE is direction-dependent.It ismuch lessevident in SRIM simulationsowing 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.2 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 FCconfirmsthenonnegligible differences between the atomic displacements induced by a PKA and those induced by an incident selfion when the energies are comparable with Ed.Accordingly,the neutron cascade option must be used to study atomic displacements versus PKA energy.

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.

Table 1 Average TDE for monatomic materials

3 Simulation results and discussion

3.1 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 valuesfrom 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.5Ed,a typical discrepancy of a factor of about 2 is found between the SRIM FCand NRT calculation.Such discrepancy hasbeen widely recognized and analyzed[5,6,9–11];therefore,the present work does not emphasize this point.

Fig.4 (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 Ed≤ED＜2.5Ed,SRIM FC show that ν＞1 and is a strictly increasing function of damage energy,whereas the NRT model implies thatνNRT=1.The ratio of SRIM FC to NRT increases from～1 to～2 when EDincreases from～Edto 2.5Ed.This region has received much lessattention becauseit isnot crucial for the displacement damage induced by reactor neutronsand 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.5Edare yet to be studied.

3.2 Atomic displacements for ED＜Ed

As the results shown in Fig.3,atomic displacements are observed when ED＜Edbut EPKA≥Ed.This is a direct consequence of the definition of Ed:a PKA with EPKA≥Edis 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 Ed,it is questionable to use the damage energy as the effective energy for computing the number of atomic displacements.

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

3.3 Atomic displacements between Ed and 2.5Ed

Once EPKAis larger than 2 times the minimum TDE,denoted by Ed，minhereinafter,it is possible to produce two atomic displacements.Because the TDE is direction dependent,it is possible thatν＞1 for Ed＜EPKA＜2Ed.In fact,for some materials,Ed＞ 2Ed，min.For iron,Table 2 reveals that 9 of the 11 interatomic potentials used in Ref.[30]give Ed＞ 2Ed，min.Because the NRT model and NRT-based models use 2.5Edas a demarcation energy to ensure continuity,we investigate the point defects for both PKA energy and damage energy in the range of[Ed,2.5Ed].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 Edand 2.5Ed.For the sake of simplification,they are respectively denoted byKP(left plot)andNRT(right plot)with the variables EPKA/Edand ED/Ed.It is obviously confirmed thatν＞1 for EPKA＞Edfor all monatomic materials.Excluding the material dependence already included in Ed(and ED),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.KPof C and W are quasi-identical but smaller than that of the other five.KPof Si is between those of Al,Fe,Ni,and Cu and the ones of C and W.However,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.of C and W are quite different,whereas theirKPare quasi-identical.

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

The discrepancies shown in Fig.5 can be attributed to the different materials and different TDEs.Therefore,Fig.6 compares the results of Si with two different values of Ed.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 Ed=21 eV obtained by Bourgoin et al.[31].It can bconcluded that the value of Edinfluences bothKPandNRT,even though they are independent of Edin typical models(cf.Section 2.1).In addition,comparing the resultsof C(Fig.5)and Si(Figs.5 and 6),KPandNRTof Si with Ed=24 eV are larger than those of C,of which Ed=25 eV.Consequently,＞1 for EPKA＞Edanddepends on both the material and value of Ed.Therefore,we cannot obtain a unique simple function of EPKA/Edor ED/Edto describe the number of atomic displacements from the SRIM FC in the range of.Nevertheless,as the results shown in Figs.5 and 6,a linear fitting of the number of atomic displacementsvs.PKA energy or damage energy is suitable for each monatomic material.

Fig.5 (Color online)Number of atomic vacanciesversusnormalized PKA(left,noted by KP in thetext)and damage(right,noted by NRT in the text)energies from SRIM FC of neutron cascade.The straight lines are linear fittings

4 Comments on the effective energy in displacement models

Comparing the two plots with the PKA energy shown in Fig.6,Ed=24 eV has large number of atomic displacementsthan Ed=36 eV for agiven EPKA/Ed.For aspecific material with a smaller Ed,the same EPKA/Edimplies a smaller EPKA,so that the partition function is larger,which further implies a larger ED/Ed.Therefore,the two different plots versus PKA energy in Fig.6 confirm that the damage energy better describes the number of atomic displacements than the PKA energy.

However,theresultsversusdamage energy illustrated in Fig.6 show that the damage energy is not necessarily better than the PKA energy for determining a simple unique formula for a specific material.In fact,the inelastic energy loss has little influence on the atomic displacements when the kinetic energy iscomparable with or even smaller than Ed.Therefore,the PKA energy and damageenergy are two extreme energies for computing the atomic displacements. New efficient energy should be determined for more accurate calculations.

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

Table 2 Comparison of Ed and Ed，min for Fe with 11 different potentials [30]

4.1 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.5Eddoes not influence the number of atomic displacements.Thus,the inelastic energy loss when an atom slows down from 2.5Edto 0 should be added to the damage energy for computing theatomic displacements[32].Based on this reasoning,they obtained the effective energy as[32]

It is noteworthy that the demarcation of 2.5Edin the NRT formula is used only to ensure the continuity of the displacement function.According to the reasoning of Kinchin and Pease[12],2Edisaphysically crucial limit.In addition,an atom does not slow down with continuous energy loss.A collision may decrease the energy of an atom from E1＞2Edto E2＜2Ed.Assuming the equiprobableenergy distribution(i.e.,hard-spherecollision[11])for an atom slowed down to E＜2Edfor thefirst 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 2Ed,ED（EPKA）≈ EPKA×P（Ed）when 0≤EPKA≤2Ed.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 Edof 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],theexperimental results of Cu obtained by Averback et al.[34],and many MD simulations for energy around 2Ed.

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 Edand 2.5Ed.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（ED）≈ ED/2Edfor ED＞2Ed.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（ED）= κED/2Edis determined by fitting the BCA calculation results[15].Therefore,if the effective energy ED，eff= ηEDrather than EDis used,the fitted constant(or widely recognized as the correction to the hard-sphere collision cross section) becomes κ′= κ/η ≈ 0.7.

4.2 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 ED（Ed）from SRIM FC and denoted byηSRIM.The effective energy is computed as follows:

Fig.7 (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(R2=0.994)

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

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

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

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

4.3 Further comments on the damage energy

Figure 8 shows that the damage energy from SRIM FC depends on the value of Ed,whereas the Lindhard’s damage energy is independent of Ed.However,these are not physically incompatible.In fact,the original equation governing the damage energy(or atomic vacancies)includes Edas 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 Edin the original equation(i.e.,their approximation(B):Edis 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 Ed:40 eV and 20 eV.The results show that the damage energy depends on Edfor PKA energiesup to 100 keV.For a given PKA energy,the damage energy is larger for a larger Ed.This is a consequence of knocked-on atoms having a smaller kinetic energy for larger Ed.A lower kinetic energy results in lower inelastic energy losses in subsequent collisions.Therefore,once Edis changed,the corresponding number of atomic displacementscannot be directly predicted as the inverse proportion to Ed.Taking the Fe PKA in Fe shown in Table 3 as an example,the number of atomic displacementsisreduced by afactor greater than 2(the last column in Table 3)if Edisdoubled(20 eV→40 eV),whereasthe NRT model predictsareduction of afactor of 2(or smaller than 2 if the slight increase in damage energy is considered).This confirms the conclusion given in Sect.4.2:ED/Edor ED，eff/Edcannot be the unique variable for computing the number of atomic displacements.A variation of Edby a factor of x does not imply a variation of a factor of 1/x for the number of atomic displacements.

Table 3 Damage energy and the number of atomic displacements for Fe PKA in Fe with Ed=40 eV(a)and 20 eV(b)from SRIM FC

5 Conclusion

In SRIM FC,the number of atomic displacementsν＞1 and cannot beuniquely described by ED/Edor EPKA/Edfor PKA energy from Edto a few times of Ed.Because a part of the inelastic energy loss(when the kinetic energy is smaller than～2Ed)does not influence the number of atomic displacements,an effective energy ED，eff= ηEDis proposed.An approximate value ofη≈1.2 is obtained for both the present proposal and that of Robinson and Oen[32].Thisvalueisconsistent with the experimental data for Fe,Ni,and Cu for damage energies of about 2Ed.Using ηSRIM=Ed/ED（Ed）,six of the seven monatomic materials considered in the present work have the same number of atomic displacements as a function of ED，eff/Ed.

However,further investigation shows that the damage energy depends on Ed,whereas the currently used analytical damage energy is independent of Ed.For a given PKA energy,the damage energy is larger for a larger Ed.Consequently,none of ED，eff/Ed,ED/Ed,and EPKA/Edcan be the unique variable for describing the number of atomic displacements.The only general conclusion is that a linear function fitting is suitable for quantifying the number of atomic displacements as a function of ED，eff,ED,and EPKAfor energy comparable with Ed.A more accurate analytical damage energy function should be determined by solving Lindhard’s integro-differential equation with Ed.