Ab initio threshold displacement energies in iron

ABSTRACT The threshold displacement energy in iron is determined using ab initio molecular dynamics. This is the most fundamental input parameter for radiation damage assessments. The predictions agree well with the available experiments and provide a significantly lower average value for iron than the standard one. This result impacts radiation damage assessments in iron alloys and steels and especially so for dose estimations and conditions close to the threshold. The importance of using an appropriate description of the core and valence electrons is highlighted. Energy loss simulations provide important fitting parameters for improved interatomic potentials. IMPACT STATEMENT Ground-breaking ab initio calculations of the threshold displacement energies in iron show significant differences in the angular anisotropy and predicted average value with respect to previous literature. GRAPHICAL ABSTRACT

The interaction of energetic particles, such as electrons, ions and neutrons, with matter can damage the crystal structure, if the energy of the particle transmitted to the lattice atoms exceeds a threshold. The most fundamental parameter defining the primary state of radiation damage for a given material is this threshold displacement energy (TDE). It is defined as the minimum kinetic energy needed to displace an atom in a material in order to create a stable Frenkel pair (FP), that is, a vacancy and a self-interstitial (SIA). The TDE has played a key role in radiation damage theory since its first explicit appearance in the literature. [1] The Kinchin-Pease [2] and Norgett-Torrens-Robinson (NRT) [3]  given amount of deposited energy will create in a material. It is the key material parameter used to convert the absorbed irradiation dose, or damage energy, [3] into the commonly used unit of displacements per atom (dpa). The estimated number of dpa provides the most convenient correlation parameter for comparing irradiation effects obtained in irradiation environments that produce very different recoil spectra such as fast compared to thermal fission reactors or neutron compared to charged particle irradiation. Therefore, it is of high importance to obtain the best possible description of the TDE to provide a suitable dpa calculation. The TDE has been studied experimentally and computationally for a wide range of materials [4][5][6][7] but to the authors' best knowledge never before with an ab initio method in a metal. The TDE has been recently investigated ab initio for a few selected semiconductors, such as silicon, [8] SiC, [9] GaN [10] and for some oxides. [11,12] Radiation damage in iron is an intense area of study [13][14][15][16][17][18][19] due to the technological interest from the nuclear fission and fusion communities and from steel manufacturing industries. The reference average TDE in iron according to the ASTM standard is 40 eV. [20] This value is, however, not based on experiments, but on computer simulations using a rather simple short-ranged interaction potential as physical basis, [21] as has been noted in a comprehensive review [22] and as is stated in the ASTM standard. [20] Experiments measuring the TDE for the high-symmetry lattice directions in bcc Fe ( 100 , 110 and 111 [23,24]) have been performed, but these data are not sufficient in order to determine the average TDE over all directions, which is the parameter of technological interest. Lucasson estimated an average TDE of 24 eV in iron [25] but noted the strong dependence on the choice of damage model and the assumed resistivity of a single FP. This average TDE is conspicuously close to the 23 eV one obtains if applying the analytical model by Jan and Seeger, based only on the measured TDE in the high-symmetry directions. [26] However, these TDE values are not used in a later review by Lucasson, [7] casting doubts on the validity of the estimates which depend intrinsically on the applied damage model. Furthermore, the analytical model is bounded by the energy in the three high-symmetry directions, in disagreement with molecular dynamics simulations. [22] Modeling short-range interactions between atoms in solids is at the stage that extreme compression is well understood in terms of the Coulomb and the screened Coulomb models, even though the screened Coulomb interaction (ZBL) is empirically fitted. [27] Near equilibrium, where the electronic and chemical interactions dominate the bonding, many, if not most, materials are very well understood. In the range between high compression (either local or global) and near equilibrium, the field is more open. In the construction of semi-empirical interatomic potentials that should be capable of treating high kinetic energy trajectories through a solid, this intermediate range is most often treated by a simple analytical interpolation with four well-defined parameters in order to observe continuity and first-order differentiability at the two knot points. Through full-potential ab initio calculations, one can approach this interaction range more explicitly and provide an insight into more reasonable choices for potential construction.
Density functional theory molecular dynamics (DFT-MD) simulations in the Born-Oppenheimer approximation have here been performed over a large range of angles and knock-on energies in the micro-canonical (NVE) ensemble. The projector augmented wave formalism [28,29] has been applied in the generalized gradient approximation. [30] Both the standard potential for Fe with electronic configuration Ar4s 1 3d 7 (from here on denominated Fe sd ) and the harder semi-core potential with configuration Mg3p 6 4s 1 3d 7 have been used (from here on denominated Fe psd ). These methods are implemented in the Vienna Ab initio Simulation Package (VASP). [31] In order to minimize the image interactions and possible glide recombination, non-cubic supercells of at least 504 bcc sites (7 × 6 × 6) have been used for the dynamic simulations. For verification purposes, 6 different supercell geometries with up to 672 atoms have been used for certain displacement simulations. Furthermore, classical MD simulations of the TDE have been performed using a state-of-the-art embedded atom method (EAM) potential for iron, from here on denoted AM04, [32] with cell sizes up to 431,880 bcc sites (59 × 60 × 61). The short-range interaction in AM04 is determined by ZBL up to 1.0 Å and then by a C1 fourthorder exponential polynomial up to 2.05 Å where the equilibrium part takes over. The time step used in the DFT-MD simulations was set to 3 fs after careful tests were performed with time steps ranging from 0.5 to 5 fs. The trajectories were seen to not change appreciably for time step values up to 3 fs. All results presented were obtained using the gamma point in the Brillouin zone and a cutoff of 250 eV. The convergence of these parameters on the dynamic results was rigorously checked. The initial temperature of the simulation cell was 0 K. Lowtemperature thermalization (30 K) was also performed before the impact event, but this had no discernible effect on the results. At higher temperatures, thermal disorder will affect the response and is expected to smear out the directional anisotropy slightly, as has been recently shown using classical MD simulations in iron. [33] The discretization of the knock-on energy sampling was 1 eV. This sampling resolution is sufficient to encompass the statistical fluctuations of multiple simulations in the same direction. The threshold energy for a given direction is determined by confirming the creation of a stable FP over a sufficient amount of time. In order to determine if a stable FP has been created, a dynamic in the NVE ensemble of 2 ps was simulated, tracing the atomic trajectories after giving an atom an impulse. In the cases where the FP ended up with a possibly athermally activated recombination distance after this time, the simulation was continued up to a maximum of 5 ps.
Due to the image interaction through the periodic boundary conditions, the heat dissipation in the relatively small simulation cells used in the DFT-MD study can bias the results if they are continued for time scales of, typically, more than 5 ps. The shock wave that is transmitted through the cell then has time to recouple unphysically with the defected region and affect the dynamics. The average simulation time needed to stabilize the FP is 1.7 ps while the minimum is about 0.9 ps and the maximum seen here is 4.4 ps. In the classical MD simulations with large enough simulation boxes and proper heat dissipation, we note that the stability of the created FP is determined by the separation distance of the dynamically created 110 dumbbell and vacancy for all but the close-packed glide directions 111 and 110 . This was also confirmed by DFT and MD simulations of the athermal dynamical instability of an FP. The ensuing recombination volume of the standard 110 FP is a quasi-spherical ∼ 10 a 3 0 for DFT (both for Fe sd and Fe psd ) and 11 a 3 0 for AM04, while for the unstable 111 FP it is an anisotropic cigar-shaped volume of ∼ 22 a 3 0 for DFT and 24 a 3 0 for AM04, in good agreement with previous MD simulations. [34] The recombination volumes calculated by Nakashima and Stoller at elevated temperatures [35] are larger than this, but they are not directly comparable due to the significant thermal component.
The dissipation of the magnetic degrees of freedom, or spin-lattice coupling, [36,37] could have an impact on the evolution of a damage event. However, since a spin flip event is typically associated with a 0.1 eV change in energy, see the supplementary material, we estimate that the impact is very small in the threshold energy scenario where dissipation phenomena play a smaller role than it does in high-energy cascade events. Even if a few atoms may have magnetic moment inversions during the threshold energy trajectory, most of the compressed atoms simply experience quenching of their magnetic moments. [38] The average TDE, integrated over all angles, is given by where E d (θ , ϕ) is the TDE for a given direction (θ , ϕ). The directions are here presented in the Miller index convention hkl . The DFT-MD simulations using the Fe sd and Fe psd potentials provide estimates of the TDE for the highsymmetry directions in good agreement with the available experiments, see Table 1. In addition to the highsymmetry directions, the 135 direction is also included, since it is the most common choice used in the literature for cascade simulations [39][40][41][42][43] and is a reasonable choice of a low-symmetry representative direction. The stiffer Fe psd generally predicts higher threshold values. Between Fe sd and Fe psd , the largest discrepancy is found for the 110 direction, where the experiments are more in line with Fe sd than with Fe psd .  [24]. b Direction specific experiments [23]. c Global minimum experiment [25]. The largest discrepancy between DFT and EAM is found for the close-packed 111 direction. There, the DFT-MD results are significantly lower than those of the state-of-the-art EAM potential AM04, which is currently by far the most widely used in defect and cascade simulations in iron. [44] The softness of the close-packed directions is further emphasized in Figure 1 where the full angular distribution of the calculated TDE are shown. The semi-empirical potential agrees rather well with DFT-MD for the lower part of the angular space, notably so for the 100 and 110 directions, but not at all for the close-packed direction 111 , where DFT-MD predicts a much lower TDE than AM04 does.
For the purpose of determining an average TDE, no significant difference from the varying box size and geometry using DFT-MD was observed. From the classical MD simulations, the difference in TDE for the selected angles was no larger than 2 eV between the box sizes and for most angles less than 1 eV. The AM04 potential conforms very well to the ASTM standard, by providing an average value close to 40 eV, in agreement with many other semi-empirical potentials. [22] Here it is shown, however, that the DFT-MD calculations provide significantly lower average values, closer to 30 eV, both for Fe sd and Fe psd potentials. The average TDE for Fe psd (Fe sd ) is calculated by Equation (1) using 27 (45) different angles uniformly selected in the irreducible solid angle, see Figure 1. A full tabulation is present in the supplementary material.
The difference between Fe psd and Fe sd is generally a shift of a few eV in TDE, see Figure 2. The angular space near the 110 direction is hyperbolic according to DFT but a local minimum according to AM04, see Figure 2. The experimental measure [24] is close to AM04 and to Fe sd . However, it is not trivial to determine the TDE experimentally for directions where there is no local minimum, since the detection can sample scattering events in adjacent angles. Considering that DFT predicts a hyperbolic TDE surface around 110 , this experimental measure may well be an under-estimation with respect to the exact direction, which may explain the difficulty Maury had on determining the TDE in the 110 direction. [24] For the other high-symmetry directions, this is not an issue since they are local minima. It is clear from Figure 2 that the two DFT potentials describe the same physics and anisotropic character, with Fe psd being stiffer, while AM04 exhibit quite different anisotropy except for some general features.
The DFT-EAM difference is due to a combination of effects. The recombination of a defect that is gliding along a high-symmetry direction occurs either if the crowdion reverses its glide direction, or if the vacancy tracks the crowdion. The degree of localization of the point defects is the major reason for the discrepancy. Both the 111 crowdion and the vacancy are more localized as predicted by DFT than using AM04, see the supplementary material.
The stronger localization of the defects causes the vacancy to remain close to its site of creation while Table 2. Defect formation energies in iron, using 300 eV cutoff, in 250-atom supercells with 3 × 3 × 3 kpoints and 1024-atom supercells with the gamma point only.

Fe sd
Fe psd the crowdion glides through the lattice, losing energy through thermal dispersion to its neighbors. In the classical MD simulations the vacancy instead starts to track the crowdion in a recoil wake and thus the FP can recombine at substantially higher energies than in DFT. The dynamic defect creation process is governed by a springlike force connecting the vacancy and the SIA. If low impact energy is used, then the vacancy is pulled along after the SIA. If enough energy is used, then the spring breaks. A secondary reason for the discrepancy is the difference in instability of the 111 defect. According to DFT the 111 defect is about 0.7 eV less stable than the ground state SIA 110 , see Table 2, while the AM04 only has a 0.4 eV difference. The larger instability of the 111 as predicted quantum mechanically increases the probability of transforming to a 110 defect and thus, the recombination through glide vector reversion is more efficiently hindered, leading to a lower TDE. The issue of possible bias in the results for the closepacked lattice directions due to long ranged 1D defect transport is addressed by performing replacement collision sequence (RCS) simulations. These simulations provide information on how much kinetic energy is lost during defect glide. The energy loss data are used to estimate the supercell size needed in order to avoid image interactions due to glide for TDE simulations in the high-symmetry directions.
RCS are simulated by giving an impulse in a highsymmetry direction to an atom and measuring the kinetic energy loss per collision ( E l ), by averaging over a statistically representative number of collisions (typically [8][9][10][11][12]. The energy has to be high enough to create that many collisions (typically well above the TDE) and low enough not to introduce a significant image interaction in the plane orthogonal to the impulse vector. In order to use the first peak in the energy loss spectrum one has to subtract the FP formation energy. The energy loss was measured by giving an Fe atom a 50 eV impulse in the two most close-packed directions in bcc, 111 and 100 . In the former case, a supercell of 576 atoms was used, elongated along the impulse vector (2 × 3 × 16 cells of  Table 3. Comparing energy loss results using semiempirical interatomic potentials, it is seen that the more recent Fe potential compares rather well with the DFT-MD results. The bias in TDE from the close-packed directions was estimated by performing TDE simulations on the RCS boxes. The results on 111 and 100 TDE were within the 1 eV energy discretization limit and thus indistinguishable from those presented in Table 1. The charge density difference distributions have been analyzed during the simulations; see Figure 3. The vacancy is inside the left-most isosurface cage (a) and the SIA is located to the far right (d). This snapshot is taken 0.45 ps after the initial impulse. The wake of the passage of the defect can be clearly seen (c). The typical vacancy type isosurface cage is delocalized in the 100 direction (b), while the SIA is continually decelerated by the very short-ranged frontal charge localization. Most of the surrounding lattice is, during the kinetic phase, completely unperturbed. Later, during the thermalization stage, all atoms in this rather small simulation cell will begin to vibrate.
In order to better understand the difference observed between the DFT potentials, we have performed quasistatic drag (QSD) simulations of the initial stage in the defect creation process. The QSD is used to analyze how the different interaction models respond to local compression. An atom is dragged step-wise along a displacement vector and the energy response of the system is recorded, see Figure 4. The two DFT potentials are used, as well as AM04 and ZBL. It is clear that the standard Fe sd potential is consistently too soft and fails to properly capture the important short-range interactions. However, for energies up to about 10 eV, there is essentially no difference between Fe sd and Fe psd and thus near-equilibrium properties and thermal vibrations are expected to be very similar for the two electronic descriptions. For the close-packed directions -100 and 111 -AM04, ZBL and Fe psd are in good agreement. For the more crystallographically open directions -110 and 135 -the divergence is larger between all methods, which can be understood by the significant difference in binding that they exhibit. The ZBL, for instance, is purely repulsive. We propose that QSD simulations can provide crucial, and computationally inexpensive, information for fitting improved interaction potentials, for instance by adjusting the many-body part (electronic density and embedding function) on a QSD isoenergy angular surface.
In conclusion, the DFT-MD simulations of the TDE in bcc Fe performed here agree very well with the available experimental data, but have significant differences with respect to the previous state-of-the-art interatomic potentials. The average value, which is of highest technological interest, is about 20% lower than the accepted standard value from the literature. The reason for this  discrepancy is dominated by the approximation of the interatomic forces that have been previously applied and the importance of an accurate and detailed description of the crowdion structure.