Defect production in tungsten under sub-threshold energy irradiation: role of hydrogen and surface effects

Despite the low solubility of hydrogen isotopes (HIs) in tungsten (W), their concentration can reach up to ∼10 at.% after low-energy plasma irradiation. This is generally attributed to the vacancies that may accommodate excessive HIs. However, the kinetic energy of incident HIs transferred to W is far below the energy threshold to create a Frenkel pair, the underlying mechanism of defect production is still unclear. Here, we investigate the influence of H on the defect production in W using the molecular dynamic (MD) simulations. It is found that the threshold displacement energy (TDE) in bulk W slight decreases with the increasing of H concentration. This is due to the formation of H-vacancy complexes, which prevents the vacancy-interstitial recombination. More importantly, the H effects are significantly magnified in the surface region. On the one hand, the maximum kinetic energy transferred from 400 eV H to W can reach up to ∼21 eV due to the double-hit process, which is two times higher than that predicted by elastic collision model. On the other hand, the momentum transferred to W is completely random, including both the recoil direction upward and downward from the surface. Accordingly, the lowest TDE in W surface is only 15–21 eV at sub-surface layers with the depth of 6.7–11.1 Å, which is lower than the maximum kinetic energy transferred to W. Therefore, the low-energy HIs irradiation can create the defects in W surface directly. Our findings provide deep insight into defect production in W at sub-threshold energy and have wider implications for materials performance under low-energy ions irradiation.

Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Tungsten (W) is considered as the most promising candidate for plasma-facing materials (PFMs) in future fusion reactors due to its excellent thermal properties, high sputtering threshold and low hydrogen isotopes (HIs) retention [1].However, W-PFMs will be irradiated by high flux of HIs (e.g.deuterium and tritium), helium plasma, high-energy neutrons, and high heat load [2,3].Such extremely severe conditions will induce the significant performance degradation of materials and raise a great concern for the stable operation of fusion devices [4].Especially considering the radioactively and self-sufficiency of fusion fuel tritium, the interaction of HIs with W-PFMs is of great importance and triggered enormous research effort [5][6][7].
Interestingly, although one of the favorable properties of W-PFMs is the extremely low solubility of HIs (∼10 −18 atomic fraction at room temperature) [8], their concentration in W can reach up to ∼1 at.% after high-flux and low-energy plasma irradiation, especially in the surface region (∼10 at.% [9,10]).Such remarkable increase of HIs retention, called 'HIs supersaturation surface layer (HSSL)', is generally attributed to the irradiation-induced vacancies in W that serve as the trapping center for HIs and accommodate additional HIs atoms [11,12].It is important to note, however, that the kinetic energy of incident HIs transferred to the W atom (based on the binary collision model) is far below the threshold displacement energy (∼41 eV [13]) in W to create stable Frenkel pairs [9].Therefore, the vacancies induced by conventional elastic collision alone cannot account for the supersaturated HIs retention in low-energy plasma-loaded W, and the intrinsic mechanism of defect production at sub-threshold energy has yet to be emerged.
First-principles calculations and molecular dynamic (MD) simulations have been used to explore the influence of H on the defect production in W [14][15][16][17][18][19][20][21][22], and several mechanisms have been proposed accordingly, including the Hreduced vacancy formation energy [14,15], the self-trapping of planar H clusters [16][17][18][19], the suppressing effect of H on defect annihilation [20,21], and the H-assisted self-interstitial structure formation [22].Although these mechanisms provide the thermodynamic driving force for superabundant vacancies, they require an extremely high concentration of interstitial H in W and the kinetic processes of H-induced defect production are still unclear.More recently, Gao et al [23] clearly observed the flat nanocavities with high aspect ratios at the HSSL in plasma-loaded W, indicating the crucial role of open surface on vacancy formation.It is, therefore, important to clarify how defects are created at sub-threshold energy in W from dynamic analysis, especially for the synergistic effect of H and surface.
Herein, employing the MD simulations, we explicitly elucidated the critical role of H and surface on the defect production in W at sub-threshold energy.Surprisingly, a 400 eV H atom can transfer 21 eV of kinetic energy to a W atom through double-hit process.The momentum transferred to W can be upward from the surface, and the corresponding threshold displacement energy (TDE) are lower than the maximum kinetic energy transferred to W. These suggest that the low-energy HIs irradiation will induce the defect production in the surface region directly, which should be responsible for the formation of HSSL and has significant effect on the evolution of HIs in W-PFMs.

Computational methods
The MD simulations were performed in the large scale atomic molecular massively parallel simulator code.An embeddedatom method (EAM) [24] potential developed by Wang et al [25] was used for recoil simulations, which shows a good performance on H-defect interactions.Similar to previous studies [26][27][28], a simulation box of 16a 0 × 16a 0 × 24a 0 (a 0 ∼ lattice constant of W) was used for all calculations.The temperature was set to be 300 K, consistent with the experimental conditions [9,10].The timestep was adapted automatically to ensure that the maximum displacement of atoms per step is less than 0.01 Å.In bulk system, as seen in figure 1(a), the periodic boundary conditions (PBCs) were applied in all three directions.The systems were thermally equilibrated at 300 K for 5 ps under an NPT ensemble (pressure is 0 pa).After that, an additional kinetic energy was given to a specific atom at the center of simulation box, which is regarded as the primary knock-on atom (PKA).During the recoil simulation, the NVE ensemble is performed for all atoms, except for three atomic layers at the boundaries with NVT ensemble.
As for surface simulation, as seen in figure 1(b), the top and bottom faces in the x-direction were non-periodic, while the PBCs were used on the lateral faces along y-and zdirections.Note that the W(110) surface was selected in present work, because it has lowest surface energy among various orientations [29].Similar to previous studies [30,31], the bottom two lattice atomic layers (x direction) were fixed.The thermal equilibrium was reached firstly under the NVT system.Then, the recoil process was simulated under NVE ensemble except for two atomic layers along x direction above the fixed layers and one layer at the periodic boundaries (y and z directions) with NVT ensemble.In the simulation of H atom incidence, a H atom was randomly placed at the position upon the W surface (∼5 Å to surface plane).Then the H atom was given a perpendicular downward surface velocity.In particular, direct reflection of H atom by surface W atoms (probability is about 8.8%) had been excluded.
For TDE calculations, each recoil process was simulated up to 5 ps.The Wigner-Seitz analysis [32] was used to determine the defect formation.For each specific direction, the recoil energy was set to be 1 eV initially, and then increases gradually with a very small step size (∼1 eV) until a Frenkel pair was generated.Accordingly, the minimum recoil energy for defect production is regarded as TDE in this direction.Note that the simulations were repeated four times for each recoil direction and the average values were presented in the following sections.

Threshold displacement energy in bulk system
In order to determine the influence of H on the defect production, we firstly examined the threshold displacement energies (TDEs) in bulk W with and without H.Considering the anisotropy of atomic structure, the TDEs are intimately linked with the recoil direction.Accordingly, we examined the stereographic projection of the crystallographic direction, as illustrated in figure 2(a).Each point on the spherical surface is equated to the projection of a recoil direction, and all of them are included in a triangle enclosed by <100>, <111> and <110> because of the symmetry of bcc lattice.Here, a set of random recoil directions (including fifty arbitrary points in the triangles) and seven typical directions (i.e.<100>, <111>, <110>, <122>, <133>, <135> and <235>), are considered to clarify the direction dependence of TDE in W. Note that these specific recoil directions contain both lowindex and high-index directions, and have been widely used in previous studies [33][34][35][36].
Figure 2(b) shows the TDE profile of different recoil directions in pure W. It is found that the <100> recoil holds the minimum value (∼30 eV) of TDE, which is in good agreement with previous experiments (∼41 eV in [13] and ∼38-42 eV in [37]) and simulations (∼31 eV in [38], ∼34 eV in [39]).Besides <100> recoil, there are several typical directions that have relatively small TDE, such as 39 eV for <111> direction and 46 eV for <122> direction.As for other typical directions, they hold the relatively large TDE, including 73 eV for <235> direction and 84 eV for <135> direction.Such high threshold energy is due to the scatter of PKA after recoil, which dissipates the kinetic energy and leads to the increase of TDE.Generally, the obtained TDE profile in pure W is qualitatively similar to that obtained by Banisalman et al [36], but there are also apparent discrepancies between our calculations and previous results.For example, the average value of TDE for 50 random recoil direction is calculated to be 54 ± 2 eV in present work.This value is reasonably consistent with those obtained by Vajda et al (∼50.5 eV) [40] and by Maury et al (∼55.3 eV) [13], but somewhat lower than that calculated by Banisalman et al (∼85 eV) [36] and Setyawan et al (∼98 eV) [28].The difference can be mainly attributed to the different interatomic potentials.We use EAM2 potential developed by Marinica et al [41], while Banisalman et al and Setyawan et al select the potential reported by Bjorkas et al [42].Physically, the value of TDE mainly depends on the interaction between W atoms at short distance, because the sufficient kinetic energy is required during collision process to displace W from its lattice site, creating vacancy and interstitial W atoms. Therefore, we examined the interaction of W-W pair in vacuum using both the potential reported by Bjorkas et al [42] and the EAM2 potential developed by Marinica et al [41].It is found that when the distance between W atoms is lower than 1.60 Å, their repulsive interaction using the potential developed by Bjorkas et al is much stronger than that using the EAM2.Namely, the former potential is 'harder' and the latter one is 'softer'.The strong repulsive interaction between W atoms at short dis-tance is expected to increase the energy barrier required to generate the displacement defect, leading to the increase of TDE.
Next, we explored the effect of interstitial H on the TDEs in W. Here, two different cases are considered.For the first one, all H atoms are randomly introduced at the tetrahedral interstitial sites (TISs), and the systems are fully optimized before the recoil simulation.As seen in figures 2(c) and (d), the TDE profiles in this case are similar to that in pure W, implying that the H addition only has a slight effect on the its distribution at different directions.However, the value of TDE slightly decreases with the increasing of H concentration.For example, the TDE value of <135> (or <100>) recoil decreases from 84 eV (or 30 eV) in pure W to 80/73 eV (or 29/29 eV) with 5/10 at.% H.Moreover, the average value of TDE in W-H system is only 54 ± 2 eV and 51 ± 2 eV with H concentration of 5 at.% and 10 at.%, respectively, which is lower than that in pure W.
The second case corresponds to an extreme situation, in which the multiple H atoms are set at the different TISs around the PKA.Similarly, the systems are fully optimized before the recoil event, and the atomic configurations are illustrated in figure 3(a).Since the average value of TDE for fifty random directions is in good agreement with that for seven typical directions, as displayed in figure 2, only the latter ones are considered in the following part.Obviously, the presence of interstitial H atoms in the vicinity of PKA reduces the TDE in bulk W, and their influences become stronger with the increasing of H numbers, from 30-84 eV for pure W to 23-61 eV with six H atoms.This should be attributed to the suppressing effect of H on the recombination of Frenkel pairs.For example, although the Frenkel pairs can be temporarily created by 45 eV PKA in pure W along <122> direction, the ejected interstitial W atom is close to the vacant site, leading to the vacancy-interstitial recombination spontaneously.Nevertheless, if there are multiple H in the vicinity of vacant site, these H atoms can segregate into the vacancy and prevent the annihilation of Frenkel pair, resulting in reduction of TDE.This may be attributed to the variation of stress field induced by H segregation, from tensile stress for a single vacancy to compressive strain for H-vacancy complexes [20].These results suggest that the addition of interstitial H atoms may inhibit the recombination of temporary Frenkel pairs, but it only has a relative low impact on the TDE (reduction of ∼6.3%) in bulk W even with 10 at.%H concentration.
As a non-hydride forming metal, the formation of hydride in W is very difficult and only occurs at extremely high pressures [43].Although there is no conclusive evidence of hydride formation in plasma-loaded W, recent computational studies [17,18,44] suggest that interstitial H atoms may aggregate spontaneously (known as H self-clustering) in W, leading to the formation of H platelets.Within the platelets, the H-to-W atomic ratio is close to 1:1, corresponding to the stoichiometry of the most stable hydride phase.Also, by using the molecular dynamic simulations and grand canonical Monte Carlo, Hou et al [17] found that the nanohydride may be formed at dislocations in W when the H concentration or chemical potential exceeds a critical value, and the growth of the nanohydride is in the form of H platelets that are enhanced by anisotropic strain.Therefore, we also investigate the influence of interstitial H clusters and hydride phases on the TDE in W. Here, as illustrated in figures 3(b) and (c), two different sites surrounding H platelet are considered, including the edge site and center site.Generally, as displayed in figure 3 3( f ).It can be found that the TDE in hydride phase is significantly lower than that in pure W, regardless of the recoil directions.For example, the TDE value along <100>/<135> direction is 21/30 eV and 17/17 eV for WH and WH 2 respectively, while that in pure W can reach up to 30/84 eV.These suggest that the hydride formation substantially reduces the TDE in W, which can be rationalized by the low cohesive energy of hydride phases in comparison with pure W.
Since severe lattice distortion is clearly observed in plasmaloaded W [9,10], it may play an important role on the defect production and hydrogen isotopes supersaturation.This is confirmed by theoretical studies [7,16,18,45], in which the behaviors of H (e.g.dissolution, migration and aggregation) and displacement cascade in W also depend on the strain filed.For example, the total number of surviving defects, the fraction of defect in cluster and the average size of defect cluster increase with the increasing of hydrostatic tensile strain and decrease at compressive state [45].Therefore, we further investigate the influence of hydrostatic tensile strain on the TDE in W. As shown in figure 3(g), the TDE monotonically decreases with the increasing of tensile strain, from 30-84 eV in strain-free W to 24-54 eV under 4% strain, which is in good agreement with previous studies [46].This is originated from the reduction of cohesive energy, which weakens the interaction between W atoms and facilitates the defect production in W.

Bombardment of H on W surface
As demonstrated in recent experiments [23], the vacancy formation at sub-threshold energy and thus HSSL in W are closely related to the surface.Therefore, in order to realize the lowenergy H irradiation in a straightforward manner, we further investigated the bombardments of H atoms on W surface.The incident energy of H is selected to be 400 eV, which is consistent with the experimental conditions [10]. Figure 4(a) shows the frequency distribution of kinetic energy (>1 eV) transferred from incident H atom to W in 1000 times simulations.For each simulation, the kinetic energy of H is dissipated gradually by extensive binary collisions and the average number of collisions is about 33 (excluding the case of transferred energy lower than 1 eV).Interestingly, although the frequency decreases with the increasing of kinetic energy transferred to W, the maximum value of transferred energy can reach up 21 eV.This is much higher than that predicted by the elastic collision model (∼8.6 eV) based on the conservation of energy and momentum [10], suggesting that the conventional model significantly underestimates the kinetic energy of W during low-energy H irradiation.Such remarkable discrepancy should be attributed to the double-hit process observed in MD simulation, in which the incident H hits a W atom twice in a very short time.One typical example is illustrated in figures 4(b)-(e).Once the incident H atom collides with a W atom at lattice site, its moving direction is changed due to H-W interaction and the kinetic energy will partly transfer to W (∼6.7 eV for this case).After that, the H atom strikes another W via central impact, and it will be bounced back quickly due to the large difference in mass between H and W. Hence, the scattered H atom will collide with the first W atom again, leading to the cumulation of kinetic energy.Despite low probability of double-hit process, its total number of times is higher than 144 in whole simulations (corresponding to the cases that energy transferred to W higher than 8.6 eV), and should be important for the low-energy H irradiation.Different from the bulk system, the TDE in the surface region is not only related to the kinetic energy transferred to the W atoms, but also dependent on the relative angle to the surface and the depth of PKA.Here, to quantitatively analysis the directional distribution of the velocity (or momentum) transferred from incident H to W, we examined the relative angle between surface and velocity of W within 1000 times incidence simulations.Note that a large number of collisions occur at each simulation and only the cases with high energy transferred to W (>1 eV) are counted.As illustrated in figure 5(a), the relative angle between surface and W velocity shows a uniform distribution from 0 to π, suggesting the random directional distribution of W velocity.Namely, the probability of W gaining a velocity upward from the surface (corresponding to the relative angle from 0 to π/2) is about 0.5, which is as same as that downward from the surface (π/2 ∼ π).This unexpected result is mainly attributed to the significant mass difference between W and H.For example, if we consider the central collision of incident H and W using the elastic collision model, the velocity of H will be changed to the opposite direction and the energy loss of H (corresponding to the energy transferred to W) is only 0.02E 0 (where E 0 is the kinetic energy of H before central collision).Therefore, the direction of H velocity is completely random (except for the first collision) and the kinetic energy dissipation is very slow, which induce the variable direction of velocity transferred to W during H incidence.Besides the directional distribution, we also examined the depth (or sub-surface layers) distribution of the collision process, as displayed in figure 5(b).It is found that the elastic collisions spread in a wide range of depth, from 1st sub-surface layer (corresponding to the depth of 2.2 Å) to the 29th layer (i.e.64.4 Å), but their maximum frequency occurs at 11th layer (i.e.24.4 Å).Such depth distribution is qualitatively consistent with the SDTrimSP calculations, in which the peak collision frequency emerges at 20 Å [10].Generally, if the velocity of W at sub-surface layers is upward from the surface, the difficulty of defect production should be dramatically reduced in comparison with the velocity downward from the surface.To verify this speculation, we further investigated the TDEs in the surface region (from 3rd to 6th sub-surface layers), and both the recoil directions upward and downward from the surface are considered.As expected, the TDEs for the downward cases are basically consistent with that in bulk W, indicating the negligible effect of surface on defect production in W.However, the results for the upward from the surface are completely different.As seen in figure 6(a), if the PKA is located at the 3th sub-surface layer (depth of 6.7 Å) and its velocity is upward from the surface, the corresponding TDE in pure W ranges from 15 to 55 eV at different recoil directions.These values are much lower than that in bulk W, which is mainly attributed to the strong absorption capability of surface to self-interstitial atom.This obviously promotes the defect production and reduces the energy required to create the defect.More importantly, the minimum value of TDE for PKA along [111] direction is only 15 eV at 3rd sub-surface layer, which is much lower than the maximum kinetic energy transferred to W (∼21 eV, see figure 4).Besides, if there are interstitial H atoms in the vicinity of PKA at sub-surface layers, the corresponding TDE decreases with H numbers for high-scattering-rate directions (i.e.[135] and [110]) but is almost invariable for lowscattering-rate directions (i.e.[100] and [111]).This is consistent with H effects in bulk W and can be attributed to the suppressing effect of H segregation on the vacancy-interstitial recombination.Furthermore, as illustrated in figures 6(b)-(d), the influences of surface on the defect production have waned with the increasing of PKA depth, leading to the slight increase of TDE.However, even at 5th sub-surface layer (i.e.depth of 11.1 Å), the lowest TDE is only 21 eV in pure W, which is also no higher than the maximum kinetic energy transferred from H to W. These results suggest that the low-energy H irradiation may generates the Frenkel pairs directly in the surface region.
As demonstrated in recent experiments [9,10,23], there is a supersaturated HIs retention (∼10 at.%) in W surface exposed to low-energy HIs plasma (215 eV for deuterium or 415 eV for H), which is generally ascribed to the creation of vacancy by incident HI ions.It is interesting to note that the Hinduced superabundant vacancy formation (SAV) was widely observed in both bcc and fcc metals [8,15,[47][48][49][50].In previous studies [8,15,[47][48][49][50], this effect is mainly attributed to the strong attractive interaction of H with vacancies.Namely, the presence of H can significantly reduce the vacancy formation energy, which increases the vacancy concentration in thermal equilibrium.Therefore, based on the energetics of H-vacancy complexes and thermodynamic models, the concentration of SAV in materials can be predicted as a function of temperature and H chemical potential (related to pressure and/or flux) [47].However, although these studies provide a plausible interpretation for H-induced SAV in materials on the basis of thermodynamics, the kinetic processes of defect production at subthreshold energy are still unclear (since the kinetic energy of HIs transferred to W is far below the TDE).
Our calculations suggest that, despite the slight decrease of TDE (due to the stabilization of temporary defect by H segregation), the influence of interstitial H on defect production in bulk W is very weak.For example, even when the H concentration reaches ∼10 at.%, the average and lowest TDE in bulk W is about 50.8 eV and 29 eV, respectively.These values are much higher than the kinetic energy transferred from incident H to W and should not be responsible for the vacancy formation at sub-threshold energy.The formation of hydride phases can significantly reduce the lowest TDE (∼21 eV in WH and ∼17 eV in WH 2 ), but it only occurs at extremely high H concentration (>10 at.% [18]).Surprisingly, the H effects are significantly magnified in the surface region.On the one hand, the incident H may hit a W atom twice in a very short time during the H bombardment simulation.This double-hit process increases the maximum kinetic energy (∼21 eV) transferred to W, which is two times higher than that predicted by elastic collision model (∼8.6 eV).On the other hand, the velocity (or momentum) transferred to W is completely random, including both the recoil direction upward and downward from the surface.More importantly, once the velocity of W at sub-surface layers is upward from the surface, the corresponding TDE is reduced significantly.The lowest value of TDE is only ∼15/18/21 eV at 3rd/4th/5th sub-surface layer (corresponding to the depth of 6.7-11.1 Å), which is lower than the maximum kinetic energy (∼21 eV) transferred to W. Therefore, the low-energy HIs irradiation can create the defects directly in W surface, which provides a reasonable interpretation for the experimental observations in W and may also be applied in plasma-loaded and/or ion-irradiated fcc metals.

Conclusion
By using the MD simulations, we systematically investigated the effect of H on the defect production in W. For the bulk system, the H addition slightly decreases the average/lowest TDE, from 54/30 eV in pure W to 51/29 eV with H concentration of 10 at.%, which is due to the inhibiting effect of H segregation on the recombination of temporary Frenkel pairs.Nevertheless, these values are still much higher than the kinetic energy transferred to W. With the further increase of H concentration, the hydride phases (i.e.WH and WH 2 ) can be formed, leading to the significant reduction of TDE.More importantly, different from the negligible effect of lowconcentrated H on defect production in bulk W, its influence becomes stronger in the surface region.This can be attributed to the double-hit process and outward collision, which increases the kinetic energy transferred to W (∼21 eV) and significantly reduces the TDE (lowest value ∼ 15-21 eV in the depth of 6.7-11.1 Å).These results explicitly clarified the kinetic process and underlying mechanism of defect production at sub-threshold energy, and provide an important reference for understanding the performance of W-PFMs under low-energy HIs irradiation.

Figure 1 .
Figure 1.Schematic representation of simulation cell, (a) bulk system and (b) surface region.

Figure 3 .
Figure 3. (a) TDE of seven typical directions as a function of the number of interstitial H atoms.The atomic configurations show the structures around the PKA before recoil.The atomic structures of (b) center site and (c) edge site of H platelet, (d) WH and (e) WH 2 .The H platelet is placed on the xy plane (below the PKA), with the size of 8 × 8 a 0 .The blue, golden and red balls indicate W atoms, the PKA and H atoms, respectively.( f ) TDE along the seven typical directions of pure tungsten, H platelet and hydrides.(g) TDE and cohesive energy as a function of hydrostatic strain.

Figure 4 .
Figure 4. (a) The frequency distribution of energies transferred (>1 eV).The number is the sum of 1000 simulations.(b) The H atom moves towards the target W atom. (c) The H atom collides with the target W atom.(d) The H atom collides with another W atom and bounces back.(e) The H atom collides again with the target W atom.The green sphere and blue spheres represent the H atom and W atoms, respectively, while the gold sphere denote the target W atom.

Figure 5 .
Figure 5. (a) Directional distribution of W velocity (>1 eV).φxy and φxz are the angles of the velocity vector to the x-axis after projection to the xy-plane and xz-plane, respectively.The light blue square indicates the surface, perpendicular to the x-axis.The angle between 0 and π/2 means that W gaining a velocity upward from the surface, while that between π/2 and π a velocity downward from the surface.(b) The depth (or sub-surface layers) distribution of the collision process.

Figure 6 .
Figure 6.TDE values for the (a) third, (b) fourth, (c) fifth and (d) sixth sub-surface layer with different direction and interstitial H atoms.The black dotted line indicates the maximum kinetic energy (21 eV) that can be transferred to W by a 400 eV H atom.