Interaction of 1/2〈111〉 interstitial dislocation loop with hydrogen and helium in tungsten: molecular dynamics simulation

The interaction of hydrogen and helium atoms with 1/2 〈111〉 interstitial dislocation loop in tungsten is investigated by molecular dynamics simulation. The binding energies of hydrogen and helium atoms around dislocation loop are calculated by molecular statics method. The results show that the outer region of the loop is attractive to the two atoms and the inner region is repulsive. Notably, the maximum binding energies are located in the core region of the dislocation loop. We have also studied the influence factors of the interaction between the dislocation loop and two atoms: free volume, lattice distortion degree, the radius and shape of the dislocation loop. The results show that large free volume benefits the retention of hydrogen and helium atoms, especially for helium. The less lattice distortion caused by the impurity atom, the more favorable for the dislocation loop to trap it. In addition, the larger dislocation loop with higher defect concentration results in stronger capture ability for the hydrogen and helium atoms. The different dislocation loop shapes lead to different binding energy distribution patterns. And the hydrogen and helium atoms tend to occupy the groove region of the concave dislocation loop. Finally, we employ the nudged elastic band theory and dynamics method to investigate the diffusion pattern of the hydrogen atom in the dislocation loop and find that the hydrogen atom tends to migrate spirally around dislocation line. Based on the obtained results, a reasonable interpretation of the interaction behaviors between the dislocation loop with hydrogen and helium atoms are discussed, which can provide essential parameters for mesoscopic scale simulations.


Introduction
Tungsten (W) and its alloys are considered the most promising candidate material for plasma-facing components in fusion reactors (such like ITER and DEMO) [1][2][3][4][5][6], due to their high melting point, high conductivity and low sputtering rate [7][8][9][10][11][12][13][14][15]. Under neutron irradiation, there will be collisions between energetic particles and lattice atoms, leading to a large number of displacement cascades, ultimately creating plenty of point defects (vacancies and self-interstitial atoms) in plasma-facing materials [16][17][18]. These defects may accumulate to form large-scale defects, which will induce the degradation of thermal and mechanical properties [8,10,[19][20][21][22]. Due to the high binding energies and low migration energy barriers, self-interstitial atoms (SIA) tend to cluster rapidly and eventually accumulate to form the interstitial dislocation loop (IDL) [23][24][25], which can be observed by transmission electron microscopy (TEM) [26,27]. In BCC tungsten, two types of IDLs are confirmed along the Burgers vectors: 1/2〈111〉 and 〈100〉, with different properties [24,[28][29][30][31]. Notably, the former types of IDLs can diffuse quickly at room temperature, while the latter is almost immobile [24,32,33]. These dislocation loops result in the hardening, embrittlement and swelling of tungsten, which has been extensively studied [23,34,35]. The high flux irradiation may lead to high permeation and retention of the hydrogen (H) atoms, which is harmful to the metallic materials [36][37][38]. Moreover, H can be easily trapped in dislocation loops, which could lead to hydrogen embrittlement and hydrogen-assisted cracking [39][40][41][42]. Helium (He) is produced in fusion reactions and transmutation reactions and can also be firmly trapped by radiation defects [43,44]. But unlike H, He tends to trap other He atoms to form the He bubbles, which can also evolve to form the nanostructured 'fuzz' on the surface observed in experiments [45][46][47]. Furthermore, 'fuzz' can increase the sputtering rate of tungsten atoms, release dust to pollute the plasma, change the hydrogen isotope retention rate, and influence the material's thermal conductivity and microscopic properties [48][49][50][51][52]. Therefore, understanding the interaction between the interstitial dislocation loop with hydrogen and helium atoms is essential to explore their microscopic mechanisms.
Ab initio calculation and elasticity theory indicated that the 1/2〈111〉 IDL presents the most stable configuration in W because the formation energy is proportional to the square of the Burgers vector [28,53]. For 1/2〈111〉 IDL, there exist three common habit planes as {111}{110}{211}. All these three planes are easily transformed into other because of the low rotational energy barrier [16]. Moreover, the IDL can act as the biased sink to attract vacancies and SIAs [54]. The IDL prefer to grow up by absorbing SIAs and can also undergo annihilation reactions with vacancy clusters [55][56][57]. It was found that IDL can act as substantial obstacles to other dislocation loops and dislocation glides, which may induce the hardening and embrittlement of the W materials [58][59][60]. First-principles calculations reveal that a single hydrogen atom tends to occupy zones with depleted (not zero) charge density, and MD studies demonstrate that it prefers to occupy the tetrahedral interstitial site (TIS) in BCC metal [61,62]. According to the kinetic Monte Carlo (KMC) calculations, it was found thatH atoms do not form bound states with other H atoms in W, which excluded the possibility of H selfclustering without the defects [63]. Another impurity atom He can modify the surface properties and affect the behaviors of hydrogen isotopes in the materials [64,65]. Furthermore, strong binding interactions between He atoms in W lead to the formation of the He clusters and blisters [66]. As lightweight atoms, the diffusion rates of H and He in bulk Ware high even at room temperature [67,68]. However, the IDL could trap H and He to limit their migration abilities [69]. Although several experimental and computational works have focused on the interaction between H/He and dislocation loops, there are lack of systematic studies of the factors influencing the interaction. Therefore, it is crucial to understand how the IDL to capture H and He atoms and explore the tritium's permeation and retention patterns around the IDL. These in turn provide theoretical and data support for the experiments.
In this study, we present molecular dynamics (MD) simulations to study the interaction between 1/2〈111〉 IDL with H and He in bcc W. The binding energy distributions are obtained by the molecular statics (MS) method to determine the attraction and repulsion regions for the IDL on impurity atoms. Furthermore, we investigate the effects of free volume, lattice distortion energy, radius and shapes of dislocation loops on the interactions deeply. Ultimately, the nudged energy band (NEB) and dynamics methods are applied to explore the diffusion pattern of H around the IDL.

Simulation details
Molecular static (MS) and molecular dynamics (MD) calculations are employed with the Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) [70] to investigate the interactions. The Dislocation Extraction Algorithm(DXA) and Voronoi analysis used for the post-processing atomic simulations are embedded in OVITO [71]. The interatomic potentials developed by Wang et al [72] and Bonny et al [73] are used to describe the W-H and W-He interactions, which is applicable to investigate the capture properties of IDL.
The principal axes x, y and z of the simulated bcc crystal box are oriented along the [111],  and [−110] directions, and there are more than 60000 atoms in the computational boxes. Periodic boundary conditions are applied in all three directions. The initial 1/2〈111〉 IDL is constructed by replacing three atomic layers with SIAs along the 〈111〉 direction on the {111} atomic plane, and the Burgers vector is normal to the habit plane [74]. Remarkably, the different insertion distances of the three layers of interstitial atoms ensure that the model is stable and to prevent changes in the habit plane or Burgers vector after optimization. Relaxation optimization is performed using a conjugate gradient algorithm embedded in LAMMPS with relative energy change between relaxation steps of 10 −20 as a stopping criterion.
We obtain the atomic stress distribution of the IDL model with the virial calculation procedure implemented in LAMMPS, and the equations of the hydrostatic and shear stresses are given [75]: Here s h represents the sum of the normal stress components and t x means the shear stress along the direction perpendicular to the habit plane. But note that the calculated stress components need to be divided by the volume Ω of each atom. Furthermore, we evaluate the trapping capability of the interstitial dislocation loop to impurity atoms by binding energy values, and the formula is as follows [76]: Here E f denotes the formation energy in this system and assesses the degree of formability of the system. N denotes the number of atoms and E c is the cohesive energy, and the E p denotes the potential energy of system. Therefore, the positive value of binding energy indicates attraction and the negative value indicates repulsion. We apply the NEB method to calculate the migration energy barrier to evaluate the diffusion ability of H and He atoms around the dislocation loop. This method is based on calculating the energy of the transition state between the initial and end states to obtain energy barriers.

The model and stress distribution
The IDL model is constructed as the figures 1 and (b) reveals that IDL is the loop made up of a dislocation line essentially. After relaxation optimization, IDL, as a line defect, can generate local lattice stresses that lead to the displacement of the surrounding atoms, as shown in figure 1(c). Figure 2 shows the stress distribution of the IDL, it is found that the hydrostatic and shear stresses are both symmetric about the habit plane's center in IDL. The hydrostatic stress is positive in the outer region of the IDL, which is expressed as the tensile stress. But it is negative in the habit plane, where is compressive stress. Notably, the hydrostatic stress leads to the alteration of lattice volume, so we denote the outside of the IDL as the stretch region and the inside as the compression region. This is because that the inner part of the IDL has three layers of atoms along the Burgers vector, while the outer part only has one layer of atoms. The presence of interstitial atoms on the habit plane leads to the local lattice swelling that the atoms in the outer region around the dislocation line are under tensile stress. The region where the hydrostatic stresses reach a maximum around the dislocation line is called the core region of the IDL. According to figures 2(c), (d), the maximum shear stress (19.05 GPa) is close to the hydrostatic stress (14.00 GPa). The result indicates that the IDL has both normal and shear stress components, which is consistent with the elasticity theory. Furthermore, the maximum shear stress is concentrated around the dislocation line. Moreover, we have calculated the stress distributions for 1/2〈111〉{112} edge dislocation (ED) with the same coordinate system, as illustrated in figure 3. The hydrostatic stresses in the ED are divided into stretch and compressive regions by the slip plane and the maximum shear stresses are concentrated around the dislocation line, which exhibits the similar stress distribution pattern to the IDL. Therefore, we assume that the 1/2〈111〉 IDL could be obtained by rotating the 1/2〈111〉{112} ED around the green line (as shown in figure 3(a)). The hypothesis helps to understand the interaction of H/He atoms with IDL in the following.

The binding energy calculation
We have calculated the binding energies of a single H or He atom around 1/2〈111〉 IDL, as illustrated in figure 4. Normally, H and He tend to occupy the TIS in bulk crystal. However, due to the displacement of atoms around IDL, the most stable sites for H and He are altered and the alteration extent depends on the distance to the core region. Figures 4(a), (c) shows the binding energies of H and He are positive in the outer region of IDL indicating that they are attracted. However, the binding energies are negative in the inner region indicating that they are repelled. Moreover, as seen in figures 4(b), (d), the attraction or repulsion of the IDL to H and He atoms becomes weaker with increasing distance from the habit plane. It can eventually converge to the binding energies of the bulk W to the impurity atoms. Notably the core region with the maximum binding energy is the optimal capture site. The binding energy distribution is symmetrical around the center of the habit plane, which is consistent with the hydrostatic stress distribution in the IDL. The hydrostatic stress induces variations in the lattice and interstitial volumes, which may affect the retention of H and He atoms in IDL. Due to the differences in the two potentials (the W-He potential does not consider the contribution of He to the electron density), the values of the binding energy of H and He in IDL are not directly comparable. Nevertheless, the difference between the capture patterns of two impurity atoms can be found in the maps. When the surrounding atoms are moved away from the original lattice site by IDL's stress, and H tends to occupy the new TIS formed by the displaced atoms. As a noble gas element, the He atom has a stronger elastic effect on W atoms, and therefore it tends to occupy the region with large interstitial space. Furthermore, we have counted the optimal trapped sites of H and He around the IDL, as shown in figure 5. It can be found that the trapped region is a hollow cylindrical geometry. The axis of the cylinder coincides with the center of the habit plane, while the radius of the cylinder is slightly larger than the IDL.

Influencing factors
The interaction between IDL and H/He atoms is investigated from two perspectives: the effect of the IDL stress field and the influence of the residual stress caused by the impurity atoms. Firstly, we have investigated the influence of the IDL stress field. As hydrostatic stress can cause the expansion and contraction of the lattice system, we have calculated the free volume distribution to represent the effect of the hydrostatic stress field. The free volume means the Thiessen Polygon volume of the atomic site, which can also represent the spatial dimension of the interstitial site. Figure 6 illustrates the shear strain field, the free volume at the interstitial sites, and the binding energy distribution of H/He atoms in the IDL. The maximum shear strain of the IDL is distributed around the dislocation line. The H atom's maximum binding energy site is located at the midpoint between the maximum shear strain and the maximum free volume sites. Therefore, it is predicted that the choice of the capture site for the H atom is synergistically affected by the two stress fields. But as the He is a noble gas element, it shows a different trapping pattern in the IDL with the H. The maximum binding energy sites for He almost coincide with the maximum free volume sites. It implies that the choice of the capture site for He atom is just influenced by the hydrostatic stress field but not the shear stress field. This phenomenon is due to the strong elastic interaction between the He and W atoms. Moreover, the interstitial He atom will repel the neighboring atoms and cause atomic displacement. Thus, the He tends to occupy the large space interstitial site to reduce the local lattice distortion of the system. Next, we investigate the influence of introducing the H and He atoms on the IDL system. We calculate the distortion energy distribution around the IDL system, as shown in figure 7. The distortion energy means the IDL's lattice distortion increment degree caused by the H and He atoms. It is found that the minimum distortion energy increment of the system is obtained when the two impurity atoms occupy the region outside the IDL around the dislocation line. According to the results of the fitted curve, these sites may coincide with the maximum binding energy sites. However, it is worth noting that the distortion energy caused by the H atom around the IDL core is negative, while the energy of the He atom is positive. It means that the H atom may cause the IDL structure to transfer closer to the bulk. To explore the reason, we have investigated the neighboring atomic configuration of the H and He atoms that are trapped in the IDL core, as illustrated in figure 8. The He atom tends to occupy the maximum free volume interstitial sites. Moreover, the He atom will repel the neighboring atoms along the Burger vector and cause atomic displacement. The repulsion of the He atom and the displaced W atoms will lead to the displacement of the other neighboring atoms in other directions, eventually increasing the degree of the IDL lattice distortion. The H atom tends to occupy the interstitial site at the midpoint of the maximum and minimum hydrostatic stress atomic sites. Figure 8(a) demonstrates that the H atom trapped in the core will cause the W atoms in loop and out of loop to displace in reverse directions, resulting the slip of the IDL along the Burger vector direction. Therefore, the capture sites of the H atoms in the IDL will also displace, resulting in the actual capture range being larger than the theoretical range. Furthermore, this is the reason that the capture range of H along the x-direction is longer than the He atom. We explored the capturing capability of the IDL for H/He atoms and find it is related to its size, as illustrated in figure 9. It can be seen from figures 9(a), (b) that the binding energy of H and He atoms rises with increasing distance from the habit plane along 〈111〉 direction from the IDL's center. It is mainly because the stress is negative in this region, and the absolute value decreases with increasing distance from the habit plane. Thus, the repulsion of hydrogen and helium by IDL is weaker along 〈111〉 direction. Figures 9(c), (d) displays the variation of the binding energy of H and He on the dislocation line along the direction of the Burgers vector with distance from the habit plane. From the diagram, we find that the binding energy reaches its peak around the dislocation line and drops dramatically as the impurity atoms leave the capture region along the P2 direction. Moreover, it indicates a well-defined capture range for H/He atoms in IDL, and the shape is consistent with figure 5. Therefore, it is concluded that the binding energy of the H/He atoms in the capture region increases with the radius of IDL. The possible reason is that a larger radius means a higher concentration of SIAs in the IDL, and it leads to an increase in the local lattice distortion local around the core region. Thus, the higher defect concentration prompts IDL to have a stronger capability to capture H and He. Nevertheless, it cannot be determined how the radius of IDL affects the range of the capture region yet.
The most common IDL shapes in bcc metal are circle and hexagon. However, we have constructed other kinds of IDL to investigate the effect of shape on the capture ability. The different shapes of IDL are constructed by modifying the distribution of SIA atoms on the habit plane, and we have created three types of IDLs: rectangular, concave and butterfly, as shown in figures 10(a), (g), (m). Figure 10(b) is the hydrostatic stress distribution on the rectangular IDL. We find that the tensile stresses are distributed in the outer region of IDL and the compressive stresses are distributed in the inner part. The maximum stress is concentrated around the dislocation line, consistent with the circle IDL discussed above. Additionally, the binding energy distributions of H and He around the IDL are shown in figures 10(c), (d). The trapped sites for both atoms are at the outer region of the IDL close to the dislocation line. From investigating H's free volume and distortion energy distribution around the IDL, we have concluded the same pattern with circle IDL: impurity atoms tend to occupy sites with large free volumes and a low degree of lattice distortion caused. Therefore, it is found that the capture patterns for both H and He atoms in different convex IDLs are similar. Furthermore, to investigate whether IDLs with grooves have the same capture pattern as the convex IDLs, they are divided into concave IDL and butterfly IDL based on the depth of the grooves, as illustrated in figures 10(g), (m). It is worth noting that the butterfly IDL corresponds to a concave IDL with a deeper groove instead of two IDLs being joined together. Figures 10(h), (n) displays the hydrostatic stress distributions for the two types of IDLs, it is found that the outer region of the IDLs exhibits tensile stresses and the inner part shows compressive stresses. However, it is notable that the maximum stresses are concentrated in the groove area of the IDL rather than being evenly distributed around the dislocation line. The reason is that the outer atoms in the non-groove region are joined to the inner atoms in only    figure 1(a)).  energies are distributed inside the grooves. It demonstrates that the optimal trapping sites for impurity atoms in both IDLs are transformed to the groove region. Subsequently, we get the free volume distribution of H in both IDLs, as shown in figures 10(k), (q). Based on the effect of hydrostatic stress, the maximum free volume sites also become the innermost zone of the groove region close to the dislocation line. The distribution of the distortion energy exhibits the same pattern, indicating that these two influences also apply to different shapes of IDLs. However, due to the difficulty in unifying the dimensions of IDLs of different shapes, we are unable to determine the effect of the IDL's shapes on the binding energy. Overall, the IDLs with grooves have different capture patterns for H and He from the convex IDLs: the impurity atoms tend to be retained in the innermost part of the groove in the former rather than the whole area around the dislocation line.

Migration patterns
The IDL has stronger mobile ability than the H atom [69]. Hence, when the H atom remains in the original position, the IDL may have migrated. The migration trajectory of the IDL at 600 K is shown in figure 11(a). We found that the IDL tends to diffuse in the direction of the Burgers vector. However, the different migration speeds of SIAs at different locations lead to the alterations in the habit plane throughout the displacement process. Considering the strong binding capacity and diffusion capacities' differences, H atoms can act as pinning sites for IDL. Therefore, we have placed H atoms evenly in the capture region of the IDL and obtain the migration trajectory of the IDL in this case, as seen in figure 11(b). At this time, the IDL's migration capacity is significantly limited. To further determine the pinning effect of H on the IDL, we have placed H atoms evenly around the half dislocation line and obtained the migration trajectory of the IDL, as shown in figure 11(c). It is interesting to note that the half IDL without H pinning can still swing back and forth driven by the temperature, but the IDL cannot displace as a whole. Overall, the IDL can act as the trapping site for H atoms and limit its diffusion. Conversely, H atoms can also act as pinning sites for the IDL and limit its migration. This phenomenon is attributed to the strong binding between IDL and H atoms, and the huge differences between their diffusion capacities. However, it is worth noting that the H atoms do not remain stationary during the pinning of the IDL driven by temperature. Therefore, we have explored the actual diffusion path of a single H atom, as illustrated in figure 12. It is clear to discover that the H atom tends to diffuse along the circle channel around the dislocation line on the outside of the IDL. In order to confirm the theory further, we use the nudged elastic band (NEB) method to calculate the migration energy landscapes of trapped H along different paths, as illustrated in figure 13. Based on the binding energy distribution around the IDL, we have predicted two possible diffusion paths: circling round the outside of the IDL and along the Burgers vector direction. The comparison reveals that the migration pattern with H circling the dislocation line has the lowest energy barrier, which is the most likely diffusion mode, and this conclusion agrees with [77]. From the results in figure 6, we note that the channel with the largest free volume around the dislocation line at the outside of the IDL can be used for the fast diffusion of interstitial atoms. Furthermore, another explanation for this phenomenon is that the binding energy of the impurity atom around the dislocation line is the highest, where the potential energy of the system is the lowest. Therefore, the H atom tends to migrate from one low-energy site to another in the core region instead of climbing across the energy barrier from low-energy sites to other high-energy sites. However, it is essential to note that this pattern only applies at a low temperature. The H can escape from the binding effect of the IDL and the IDL will escape from the pinning effect of the H atoms at a high temperature. In general, the H atom tends to diffuse along the dislocation line under the trapping effect of the IDL, and a certain number of H atoms can pin down the IDL and inhibit its free migration.

Conclusion
The interaction of H and He atoms and 1/2〈111〉 IDL is studied using molecular dynamics. Static calculation results reveal that the H and He tend to occupy the outer sites of the IDL around the dislocation line similar to the 1/2〈111〉{112} edge dislocation. Additionally, the IDL's capture region of the impurity atoms is the cylinder with a radius slightly larger than the IDL. Based on the elasticity theory and the properties of the IDL, we predict the factors influencing the elastic interaction: free volume, distortion energy, size and shape of the IDL. It is found that the H tends to occupy the sites between the maximum free volume positions and the dislocation line and the He tends to occupy the maximum free volume sites. The reason for this distinction is: as a noble gas element, He has a stronger elastic effect and is minimally affected by shear stress. He atoms will repel its neighboring atoms and cause the atomic displacement. H atoms will cause the reverse displacement of atoms on either side of the dislocation line. The size of the IDL is positively related to its ability to attract H and He because larger IDL sizes correspond to higher defect concentrations with stronger trapping effects. The shapes of the IDL will affect the trapped sites of H and He, and the impurity atoms tend to retain in the grooves of the concave and butterfly IDLs. Finally, we have explored the H atoms' migration paths and pinning effect by NEB and dynamics methods. H tends to diffuse in the circle along the dislocation line under the trapping effect of the IDL, while the pining effect of H atoms inhibits the migration of the IDL. The IDL increases the retention of H in the material, limiting the diffusion capacity of both. Moreover, because He atoms tend to aggregate and form clusters in W, the retention of them is synergistically influenced by the IDL and He clustering. Therefore, the IDL around the surface promote the hydrogen embrittlement and 'fuzz'. In turn, the pinning effect of H atoms limits the migration of the IDL and influences annihilation reaction between the IDL with the vacancy and their clusters, which result the pore swelling in the material. The results of this work help understand the interactions between the IDL with H and He, and they can be employed as input data for larger-scale simulations.