Paving a way to suppress hydrogen blistering by investigating the hydrogen–beryllium interaction in tungsten

Abstract

We have investigated the hydrogen–beryllium (H–Be) interaction and the behavior of Be in tungsten (W) in order to explore the influence mechanism of Be on H using a first-principles method. A single Be atom is energetically favorable sitting at the octahedral interstitial site (OIS) instead of the tetrahedral interstitial site (TIS), and prefers to diffuse along the OIS → TIS → OIS path. Interestingly, it has been demonstrated that there is large binding energy between Be atoms in W (>1 eV), leading to them energetically clustering. This can be attributed to the strong Be–W repulsion and the intrinsic electrophobic properties of Be. Further, it is found Be has significant effect on H behavior in W. On the one hand, the interstitial Be atom can enhance the stability of H in W, and thus it can serve as a trapping center for H, due to the redistribution of electron density induced by Be. This will block H diffusing deeper into the bulk, leading to the decrease of H retention in W. On the other hand, the H trapping capability of vacancy will be severely weakened by Be, because Be will provide electrons to vacancy resulting in the increasing of the electron density. Hence, the Be-V complex only can hold 5H atoms, 58.3% less than that of the Be-free vacancy, and there is no H₂ molecule formation. This indicates Be could suppress H bubbles formation in W. Consequently, H retention and blistering in W can be suppressed by doping Be.

---

1. Introduction

Nuclear fusion is considered to be a good way to liberate human beings from fossil fuel sources in future. Magnetic confinement fusion energy is being developed via the International Thermonuclear Experimental Reactor (ITER) project,¹ which is one of the most ambitious energy projects in the world today. It aims to demonstrate the extended burn of the deuterium–tritium (D–T) plasma and thus the scientific and technological feasibility of fusion power. Nowadays, the development of key materials is apparently particularly challenging for the final application of nuclear fusion energy, in which the choice of the plasma facing materials (PFMs) is one of the key issues. Owing to many attractive properties, tungsten (W) has been considered as the most promising candidate for PFMs in fusion reactors. However, as a PFM, W will be exposed to high fluxes of hydrogen (H) isotope ions, which can induce the H bubbles on W surfaces. Hydrogen bubbles can not only degrade the physical and mechanical properties of materials, but also influence the plasma stability and result in great concerns with respect to safety and cost efficiency in fusion reactors.^1,2 Therefore, H isotopes retention and blistering in W are the great challenges for the application of W in fusion reactors.

The formation mechanism and suppression methods for H bubbles in W have been under intensive investigation both experimentally and computationally. In experiments, H isotopes plasma exposure usually was performed in the plasma generators, including Pilot-PSI,³ PISCES-A,⁴ and NAGDIS-I,⁵ etc. Nuclear reaction analysis (NRA) using different ion energies was applied to examine the depth profile of the retained H isotopes. Thermal desorption spectroscopy (TDS) was performed to study the behavior of thermal H isotopes release. Surface blistering was investigated after plasma exposure using scanning electron microscopy (SEM) and transmission electron microscopy (TEM). SEM combined with a focused ion beam (FIB) was used for preparing and imaging cross-sections. Positron Annihilation Doppler Broadening (PADB) was used to monitor formation and clustering of defects. Further, detailed descriptions of H behavior at different scales have become attainable due to the rapid development of advanced algorithms and computing capability. Modeling and simulations provide indeed an indispensable way to explore the underlying physical mechanisms and compare with experiments. Various computational methods, covering different space and time scales, are employed to investigate the behavior of H in W: first principles based on density functional theory (DFT), molecular dynamics (MD), kinetic Monte Carlo (KMC), cluster dynamics (CD) and finite element methods (FEM). It has been found that the defects in W are the origin for H bubble formation, such as vacancies, and grain boundaries.^6–18 An single H prefers to occupy the tetrahedral interstitial site in bulk W,^{10,11,13,15,19,20} but it can't form H₂ molecule at the interstitial cases due to the very weak interaction between H atoms (∼0.1 eV).¹⁰ Further, it has been demonstrated that vacancy and grain boundary defects can serve as trapping centers for H,^{10,12–15,18} which provide large space with low charge density. It is found that an H₂ molecule will form at the vacancy center with as many as 14H atoms in a vacancy¹⁸ at 0 K or 10 in ref. 10. Further, only 6H atoms can be trapped by a vacancy at room temperature,^15,16 and then there is no H₂ molecule formation. Therefore, H₂ molecule is mostly believed to form in vacancy clusters and void when H concentration reaches extreme conditions. The strain induced by H accumulation in defects will enhance H solubility in W and facilitate further growth of H bubble, i.e. the strain-triggered cascading effect on H bubble growth.²⁰

Great efforts have been made to explore the possible approach to suppress the H bubbles formation and reduce the H retention in W. It is interested to find that the existing of other impurity elements have helpful effect on suppressing H bubbles during investigating their synergistic effect on W, especially helium (He).^11,21–25 It is reported that D retention decreases significantly due to the presence of He in comparison with that in pure W.^21,23 This can be attributed to that He is more energetically favorable sitting at defects than H, because of its special closed-shell structure,^11,25 and thus block the formation of H₂ molecules.¹¹ However, it should be noted that He is also harmful to W because of the formation of He bubbles under He irradiation,^26–29 similar to H bubbles. Moreover, other noble gas will induce extra defects when they are implanted into W, due to their large atomic radius and high diffusion energy barrier, and then influence the thermal conductivity and mechanical properties of W. Therefore, we should explore other better approach to suppress H bubbles in W.

Although He cannot be used to suppress H bubbles in actual engineering, it provides a new idea to solve H retention and blistering in fusion reactor. Next, we pay our attention to explore other elements, which can play the similar role in suppressing H bubbles, such as beryllium (Be). The ITER will be operated with Be wall in the main chamber and W surfaces in the divertor. Be will go into plasma under H plasma irradiation, and then has effect on H behavior in W.^30–36 Most importantly, some experimental results show that the existing of Be can really play a role in suppressing H blistering in W.^30,31 It is found that only 0.2% Be impurity in the plasma can significantly eliminate the H bubbles on the W surface at 573 K.³¹ In comparison with He seeding, the Be seeding shows more dominant influence on the D retention in W.³² In addition, the Be/W mixed materials also have been investigated.^37–40 It is found that the solution energy of H in Be₁₂W (1.10–1.41 eV) is similar to that of H in W, but the diffusion barrier of H in Be₁₂W (0.35–0.50 eV) is significant larger than that in W.³⁴

Apparently, the preliminary results show that the existing of Be is useful to suppress H bubbles in W, but the behavior of Be and its effect on H is still unclear. So far little work focuses on this aspect. In the present work, we have systematically investigated the dissolution, diffusion and accumulation of Be and its effect on H in W in order to explore the physical origin of the interplay between H and Be. The understanding of the mechanism will provide a positive contribution to suppress H isotope bubbles in W-PFMs.

2. Computational method

Our first-principles calculations were performed using the pseudopotential plane-wave method implemented in the Vienna Ab initio Simulation Package (VASP) code^41,42 based on density functional theory. We used the generalized gradient approximation of Perdew and Wang⁴³ and projected augmented wave potentials,⁴⁴ with a plane wave energy cutoff of 350 eV. A bcc W supercell of 128 atoms without and with monovacancy was used, and their Brillouin zones were sampled with (3 × 3 × 3) k-points by the Monkhorst–Pack scheme.⁴⁵ The calculated equilibrium lattice constant is 3.17 Å for bcc W, which is in good agreement with the corresponding experimental value of 3.16 Å. A conjugate-gradient algorithm was used to relax the ions into their instantaneous ground state. The relaxations of atomic position and optimizations of the shape and size of the supercell were performed. The structural optimization was truncated when the forces converge to less than 10⁻³ eV Å⁻¹. The version of the VASP code we used is 5.3.3.

The climbing-image nudged elastic band (CI-NEB) method was used to investigate the diffusion energy barrier of Be in W.^46,47 This method is used to find a saddle point and a minimum energy path between the initial and final states of a transition, which constructs a number of intermediate images along the energy path. A spring interaction between adjacent images is added to ensure continuity of the path. The activation energy is obtained by the energy difference between the minimum and the transition states. Seven images were used in all CI-NEB calculations. The forces on all the atoms in each image of the CI-NEB chain converge to 10⁻³ eV Å⁻¹.

The zero-point energy (ZPE) correction has been taken into account for all cases. ZPE of light atoms are calculated by summing up the vibrational energies of the normal modes by the expression: where ℏ and v_i is the Plank's constant and the normal vibration frequencies. Our computation shows the ZPE of Be is ∼0.116 eV and ∼0.095 eV at the TIS and OIS in the intrinsic bulk W, respectively. Therefore, the ZPE of Be has been taken into account in the investigation of Be behavior in W and its effect on H. The ZPE correction of the H solution energy surrounding the interstitial Be (ΔZPE_interstitial) can be calculated by

where ZPE_128W,H,Be is the total ZPE of H and Be atom in W, and ZPE_128W,Be is the ZPE of one Be atom in W. The third term is the ZPE of one H atom in H₂ molecule, which is calculated to be 0.135 eV. It is found that the variation of ΔZPE_interstitial induced by Be is very small (in the range of ∼0.001–0.02 eV), which is similar with that of oxygen.⁴⁸ Further, the ZPE correction of the H solution energy in the Be-vacancy complex (ΔZPE_Be-V) can be defined aswhere ZPE_127W,H,Be-V is the total ZPE of H and Be in W with vacancy, and ZPE_127W,Be-V is the ZPE of one Be atom in vacancy. The ZPE correction also has been taken into account for the H trapping energies surrounding interstitial Be atom and in the Be-vacancy complex.

3. Results and discussion

3.1. The behavior of Be in bulk W

3.1.1. Stability of Be in bulk W. In order to investigate the effect of Be on H behavior in W, we first examine the stability of Be in W. When a Be atom is introduced into a bcc W lattice, it should prefer to occupy an interstitial site for its small size relative to the host atom. There are two kinds of interstitial sites in a bcc lattice, the tetrahedral interstitial site (TIS) and the octahedral interstitial site (OIS). The structural models of TIS and OIS are shown in Fig. 1. The TIS has four nearest neighboring W atoms at 0.559a₀, where a₀ is the lattice constant. The TIS is in the center of the tetrahedron formed by its nearest neighboring W atoms, as shown in Fig. 1(a). The OIS has six nearest neighboring W atoms, with two of them located at 0.5a₀ and four at 0.707a₀. These six W atoms form an octahedron, and the OIS is in the center of the octahedron, as shown in Fig. 1(b). When the Be atom occupies an interstitial site, the solution energy of E^s_Be is defined as

E^s_Be = E^T_128W,Be − E^T_128W − E^ref_Be, (3)

where E^T_128W,Be and E^T_128W are the total energies of the supercell containing 128W atoms with and without one Be atom, respectively. E^ref_Be denotes the reference energy of a Be atom, which is chosen to be the energy per atom of the hexagonal closed packed (hcp) Be in its ground state. Here, we obtained E^ref_Be = −3.73 eV. As shown in Table 1, the solution energy of Be at the OIS is 4.20 eV, which is 1.08 eV lower than that at the TIS, and thus OIS is the most energetically preferable site for Be in bulk W, consistent with the previous study.⁴⁹ Table 1 also shows that the solution energy of Be at both TIS and OIS are positive, indicating that the dissolution of Be in bulk W is an endothermic process, in agreement with that of He in W.^11,19,25

Fig. 1 The interstitial sites in a bcc lattice. (a) TIS, and (b) (OIS). The larger blue spheres represent the W atoms, while the smaller red spheres represent the potential Be atoms sites.

Table 1 The solution energy, mechanical contribution and electronic contribution of a single Be in W at TIS and OIS, and two Be atom at OISs in the most stable configuration. All energies are in units of eV

	1Be		2Be
	TIS	OIS	OIS
Solution energy	5.28	4.20	7.29
Mechanical contribution	3.32	2.58	5.42
Electronic contribution	1.96	1.62	1.87

---

The covalent radii of Be and W atoms are ∼0.995 Å and ∼1.30 Å, respectively; while the initial distances between Be and its first nearest neighboring (1NN) W atoms are 1.77 Å and 1.59 Å at TIS and OIS in bulk W, respectively. This clearly indicates that Be will induce lattice distortion when it occupies interstitial site. Moreover, there is electronic interaction between Be and W atoms. Both of them have contribution to the dissolution of Be in W. To further shed light on the physical mechanism underlying the stability of Be in W, it is helpful to decompose the solution energies into two contributions.⁵⁰ One is the deformation energy induced by the embedded Be atom, defined as the energy release during the course of W relaxation after Be is removed, and called the mechanical contribution. The other is the electronic effect (the direct interaction between Be and W atoms), called the electronic contribution. After optimization, it is found that the mechanical contribution (deformation energy) is 3.32 eV and 2.58 eV for Be at TIS and OIS, respectively, which is accounted for 62.9% and 61.4% of the solution energy. Therefore, the mechanical contribution should be mostly responsible for the poor stability of Be in W. Further, the electronic contribution are also positive values, 1.96 eV and 1.62 eV for Be at TIS and OIS, respectively. This indicates the electronic environment of both TIS and OIS is uncomfortable for Be.

In order to analyze the electronic interaction between interstitial Be and neighboring W atoms, we further investigate the projected electronic densities of states (DOS) for Be in W. It is found that only the p-projected DOS of Be and the d-projected DOS of its 1NN W exhibit significant changes (Fig. 2). There may be hybridization between the p-projected DOS of Be and the d-projected DOS of W, because they are similar in shape. Of course, it should be noted that the Be–W hybridization is passive due to the short distance between Be and W at interstitial cases. The stronger the hybridization, the less is stability of Be, because Be is unwilling to get electron from W. Fig. 2 shows that the p-projected DOS of Be at TIS near the Fermi energy level is higher than that of Be at OIS, indicating the Be–W hybridization at TIS is stronger than that at the OIS, consistent with that the electronic contribution of Be at TIS is larger than that of Be at the OIS (Table 1). This is in agreement with the hybridization between He and metal atoms in interstitial cases.^51–53

Fig. 2 The d-projected DOS for the nearest-neighboring W atom of interstitial Be, the p-projected DOS for Be. The solid and dashed lines show the partial DOS of W and Be with a Be atom at respective TIS and OIS.

3.1.2. Diffusion of Be in bulk W. The diffusion properties are quite important for us to understand the behavior of Be in W. The CI-NEB method has been used to find a minimum energy path. There are two main mechanisms for impurity atoms diffusing in a solid, i.e. the interstitial and substitutional mechanisms, which depend on the occupation behavior of the impurity atoms. As mentioned above, Be is a typical interstitial impurity in W due to its small size, and thus Be will diffuse in W via the interstitial mechanism. It has been demonstrated that Be is energetically favorable sitting at OIS, hence, we consider processes that allowed Be to hop from an OIS to another nearest neighbor OIS. It is found that the OIS → TIS → OIS path is the optimal diffusion path for Be in bulk W with the diffusion barrier energy of 1.08 eV. This indicates that it is difficult for Be entering the inside of W by ion implantation, while it will deposit on the near surface region of W, in agreement with the previous studies observation.^33,37,39,51 It has been reported the diffusion barrier for Be penetrating into the subsurface from W (100) surface is ∼3.60 eV,⁴⁰ which is about 3 times of the diffusion barrier of Be in W bulk.

We further estimate the diffusion coefficient of Be in W in order to explore the effect of temperature on the Be diffusion behavior. Here, we employ the Arrhenius diffusion equation D = D₀ exp(−E_b/kT) to calculate the Be diffusion coefficient, where D₀, E_b, k, and T represents the pre-exponential factor, activation energy (diffusion energy barrier), the Boltzmann constant, and the absolute temperature, respectively. According to Wert and Zener's diffusion theory,⁵⁴ the pre-exponential factor (D₀) can be obtained by where n and a represents the number of equivalent jump paths and the equivalent jump length, respectively. There are four nearest neighbor OISs for a OIS in bcc W, and thus there are four equivalent jump paths of Be. Therefore, n is obtained to be 4. The value of a is equal to the distance between two nearest neighbor OISs in bcc W, i.e. 1.59 × 10⁻¹⁰ m. ν is the vibration frequency of interstitial Be in W, which can be expressed as where m is the mass of Be atom. It is calculated to be 3.03 × 10¹⁴ Hz using the Be mass of 1.49 × 10⁻²⁶ kg. Further, the D₀ is obtained to be 6.38 × 10⁻⁷ m² s⁻¹. The diffusion coefficient of Be as a function of the reciprocal temperature is shown in Fig. 3. It can be found that the Be diffusion coefficient will increase with the increasing of temperature, in agreement with the experimental observation.⁵⁵ This indicates that the diffusion of Be in W will become easier with the increasing of temperature. In addition, it should be noted that the Be diffusion coefficient of the experimental values⁵⁵ (1023 K, 4.33 × 10⁻¹⁵ m² s⁻¹; 1123, 5.83 × 10⁻¹³ m² s⁻¹) is lower than the calculated results (1023 K, 3.08 × 10⁻¹² m² s⁻¹; 1123, 9.17 × 10⁻¹² m² s⁻¹). This may be attributed to the difference between the experimental samples and the calculated model. There is no any defect in the calculated model, and thus the calculated diffusion energy is the ideal diffusion energy. However, there are various defects in experimental samples, which can serve as the trapping center of Be, and make the Be diffusion energy increase. Consequently, the defects will leads to the decrease of the Be diffusion coefficient.

Fig. 3 The diffusion coefficient of interstitial Be in W as a function of the reciprocal temperature.

3.1.3. Interaction between interstitial Be atoms in bulk W. Next, we investigate the interaction between interstitial Be atoms in order to explore the accumulation behavior of Be in W. As mentioned above, Be is energetically favorable sitting at the OIS in bulk W. Therefore, the first Be atom is set at a OIS, and the second one occupies other OIS's surrounding the first one, as shown in Fig. 4. The Be–Be binding energy (E^b_Be–Be) can be obtained by

E^b_Be–Be = 2E^T_{128W,1Be(OIS)} − E^T_{128W,2Be(OIS)} − E^T_128W, (4)

where E^T_{128W,1Be(OIS)} and E^T_{128W,2Be(OIS)} are the total energies of supercell containing 128W atoms with one and two Be atoms at the OIS, respectively. Positive values indicate attraction, while negative ones indicate repulsion.

Fig. 4 The interstitial Be–Be positions in bulk W. The first Be atom is denoted by FB, and the potential site for the second one is denoted by A–E. The larger blue spheres represent the W atoms, while the smaller red spheres represent the potential Be atoms sites.

Table 2 lists the binding energies for different Be–Be configurations as shown in Fig. 4. The FB-B case is the most stable Be–Be configuration with the largest binding energy of 1.11 eV at a distance of 1.93 Å. The Be–Be binding energy should be attributed to the strong Be–W repulsion, which is similar with the case of He in W.^19,56 The Be–W repulsion drives Be atoms to get close to each other, in agreement with the experimental and theoretical studies.^37–39,57 It has observed that Be will form Be-rich alloy phases in W, such as Be₁₂W and B₂₂W. Therefore, we propose that the self-trapping mechanism not only is applicable for He in metals but Be as well. In addition, for the FB-A case, the binding energy is −0.09 eV, suggesting there is repulsive interaction between Be atoms. Correspondingly, the Be–Be distance in the final configuration increases by 0.23 Å in comparison with initial state. This can be attributed to that the Be–Be distance is shorter than the sum of the covalent radii of two Be atoms. The Be–Be interaction will change to be weak with the increasing of the Be–Be distance, such as the FB-E case.

Table 2 The binding energy (eV) and the initial/final distance (Å) between interstitial Be atoms corresponding to the second Be atom situated at A–E position, respectively

Position	Be–Be initial distance	Be–Be final distance	Binding energy
A	1.58	1.81	−0.09
B	2.24	1.93	1.11
C	2.75	2.05	0.88
D	3.54	3.76	0.25
E	4.48	4.46	0.11

---

Next, we investigate the atomic configuration and the electron density in order to explore the origin of the large Be–Be binding energy in W. As a matter of fact, the Be–Be binding energy is equal to the difference between the double of the solution energy of a single Be at OIS and the solution energy of Be–Be pair. Apparently, the solution energy of Be–Be pair (7.29 eV) is much lower than the double of the solution energy of a single Be (4.20 eV) as listed in Table 1, which leads to the large Be–Be binding energy. The solution energy of Be–Be pair also can be decomposed into two contributions as well as that of a single Be mentioned above. It is found that the lattice expansion induced by Be–Be pair larger than that of the single Be atom, indicating the synergistic effect of Be–Be pair will cause more significant damage to the W atoms bonds in comparison with the single Be. Table 1 shows that the mechanical contribution of solution energy for Be–Be pair is 5.42 eV, while it is only 2.58 eV for a single Be. Thus, the mechanical contribution for Be–Be pair is 0.26 eV larger than the two times of that for the single Be. This indicates that the synergistic effect of Be–Be pair will cause more significant distortion to the W lattice in comparison with the single Be. Further, it is interested to find that the electronic contribution of solution energy for the Be–Be pair is 1.87 eV, while it is 1.62 eV for a single Be. Correspondingly, the electronic contribution for Be–Be pair is −1.37 eV lower than the two times of that for the single Be. This should be responsible for the large Be–Be binding energy in W. The electronic contribution is directly related to the electronic environment of Be in W. The electron density at the site of Be–Be pair is ∼0.02 e Å⁻³ lower than that of single Be, due to the significant lattice expansion induced by Be–Be pair.

We further investigate the effect of electron density on Be dissolution by calculating the solution energy of Be in a homogenous gas shown in Fig. 5. The dissolution of Be in a homogeneous electron gas with different electron density is simulated by adding electrons in an empty box with a compensating uniform positive charge background. First, we constructed a vacuum cubic supercell with fixed volume. Then, a given number of electrons is added for a given electron density along with the same density of uniform positive charge background. Next, Be is added in the supercell at different electron density, and the Be solution energy is calculated. It is found that the Be solution energy almost increases linearly with the increasing of electron density, except for a small range of much low electron density (<∼0.015 e Å⁻³), in consistent with previous study.⁵⁸ Since the electron density everywhere in bulk W is much higher than 0.015 e Å⁻³, and thus this small range can be neglected. Here, the slope of the fitting line of the Be solution energy is calculated to be 67.6 eV (e/ų)⁻¹. Therefore, the electron density difference of 0.02 e Å⁻³ will leads to the change of the Be solution energy of 1.35 eV, which is in good agreement with the difference of the electronic contribution between Be–Be pair and the single Be. Consequently, we propose that the decrease of local electron density induced by the Be–Be synergistic effect as well as their electrophobic interaction⁵⁹ drives Be atoms to get together in W.

Fig. 5 The Be solution energy as a function of the electron density for Be in homogenous electronic gas.

3.2. Effect of the interstitial Be on H in bulk W

3.2.1. Stability of a single H surrounding the interstitial Be. As mentioned above, an interstitial Be atom energetically prefers to occupy the OIS instead of the TIS in bulk W. On the contrast, a single H atom has been demonstrated to be energetically favorable sitting at the TIS.^{10,11,13,15,19,20} Next, we investigate the effect of Be on the stability of H in W, one Be atom is set at a OIS while the H atom occupies a series of TIS with the different distance from Be. The H solution energy (E^s_H) surrounding Be can be expressed aswhere E^T_{128W,Be(OIS),H(TIS)} and E^T_128W,Be(OIS) is the total energy of the supercell containing 128W atoms and one Be-OIS atom with and without H-TIS atom. E_H2 is the total energy H₂ molecule.

Fig. 6 shows the H solution energy at different TIS as a function of the H–Be distance. When H occupies A site, the initial distance between H and Be is ∼0.79 Å, which is much shorter than the sum of H and Be covalent radius. This leads to that the H solution energy is much larger than that in bulk W without Be, because of the strong repulsive interaction between H and Be. After A case, the H solution energy surrounding the interstitial Be atom is lower than that in bulk W (Fig. 6), suggesting that the presence of Be can enhance the stability of H in W. The H solution energy increases with the increasing of H–Be distance, and converges to that of H in bulk W at ∼5.09 Å. Further, the H solution energy at the most stable site (B, the H–Be distance of ∼1.69 Å, shown in Fig. S1(a)†) is 0.20 eV lower than that in bulk W, which may be considered as the H–Be binding energy. This indicates there is attractive interaction between H and Be, which is similar with the H–He interaction in W.^11,25 Therefore, the interstitial Be atom can serve as a trapping center of H.

Fig. 6 H solution energy in W with Be as a function of the H–Be distance. A–J represents the potential occupation sites of H.

In order to explore the origin of the H solution energy decreasing due to the existing of Be, we further investigate the distortion of atomic configuration and the redistribution of electron density induced by Be. The doped Be atom leads to significant local lattice expansion. For example, the distance between Be and W (1NN) atom extends to 2.02 Å as compared with initial distance of 1.59 Å. The significant shift of Be neighboring W atoms will cause the movement of other W atoms. These will lead to further expansion of the neighboring TISs surrounding Be, because they share W atoms with the Be-OIS. Most importantly, the volume expansion will further induce the electron density decrease of the TISs. It has been demonstrated that the dissolution of H in metals can be understood via the optimal charge density mechanism.^10,12,58 The optimal charge density for H in homogenous electron gas is ∼0.018 e Å⁻³,⁵⁸ but the electron density in bulk W is much larger than it. Therefore, the lower the electron density in W, the lower is the H solution energy. It is found that the electron density at the normal TIS in W without Be is ∼0.27 e Å⁻³. The presence of Be leads to the significant decrease of electron density at the TISs (B–J), inducing the decrease of H solution energy at these TISs in comparison with that in bulk W. For example, for the TIS of B with lowest solution energy, the electron density decreases to as low as 0.20 e Å⁻³. Therefore, the effect of interstitial Be on H dissolution in W can be attributed to the redistribution of electron density and local lattice distortion induced by Be.

3.2.2. H trapping surrounding the interstitial Be. As mentioned above, the presence of Be can facilitate the dissolution of H in bulk W, and thus it can act as the trapping center for H. Next, in order to investigate the trapping of multi-H atoms surrounding the interstitial Be atom in W, we calculate the trapping energy of additional H atoms segregating to the Be atom and determine the number of H atoms that an interstitial Be atom can trap. As mentioned above (Section 3.1.1.), a single Be atom prefers to occupy the OIS in W bulk. Thus, one Be atom is set at a OIS. After that, we bring the H atoms one by one to the interstitial Be atom and minimize the energy to find the most stable configuration of BeH_n. In each step, we investigate up to 9 possible configurations based on the most stable configuration of BeH_n−1. The H trapping energy can be expressed as

E^trapping_H,Be(OIS) = E^T_{128W,Be(OIS),nH(TIS)} − E^T_{128W,Be(OIS),(n−1)H(TIS)} − E^T_128W,H(TIS) + E^T_128W, (6)

where E^T_{128W,Be(OIS),nH(TIS)} and E^T_{128W,Be(OIS),(n−1)H(TIS)} are the total energies of the supercell containing 128W atoms, one Be-OIS atom, and n/(n − 1)H atoms, respectively. E^T_128W,H(TIS) and E^T_128W are the total energies of the supercell with and without H-TIS atom, respectively. The H trapping energy as a function of H number trapped by the interstitial Be is shown in Fig. 7(a).

Fig. 7 (a) The trapping energy per H as a function of the number of H atoms trapped by the interstitial Be in bulk W, and the zero point is the energy of H at the TIS far away from the interstitial Be. (b) The atomic configuration of the BeH₅ in W.

It can be found that the trapping energy for the first H trapped by Be is about −0.20 eV with H occupying the TIS. Next, the H trapping energy will increase with the increasing of the number of H atoms. This indicates the trapping capability of Be for the subsequent H will change to weaken with the increasing H number. When the fifth H atom is added, the trapping energy increases to −0.063 eV, 0.137 eV larger than that of the first H. Further, the trapping energy becomes to be positive value of 0.11 eV for the sixth H in reference of H at the TIS far away from Be in the bulk W. This suggests that the sixth H energetically prefers to occupy the TIS far away form Be rather than neighboring it. Therefore, a single interstitial Be atom can trap as many as 5H atoms. Fig. 7(b) shows that the atomic configuration of the BeH₅. The shortest distance between these 5H atoms is 2.14 Å, which is much larger than the H–H distance in H₂ molecule (0.75 Å), indicating that there is no H₂ molecule formation surrounding the interstitial Be. The relaxed configuration of BeH₅ with all atoms is shown in Fig. S1(b).† In addition, it can be found that the H trapping capability of a single interstitial Be is much weaker than that of a vacancy. On the one hand, according to the definition of H trapping energy, negative values indicate that H will be trapped by Be, while positive ones indicate H energetically prefers to occupy the TIS far away from Be rather than neighboring it. Therefore, the lower (more negative) the H trapping energy, the stronger is the H trapping capability of Be. Here, the trapping energy is only −0.20 eV for the first H surrounding Be, while it is −1.36 eV for H in vacancy.¹⁵ Thus, the H trapping energy surrounding Be is larger than that in the vacancy, suggesting the H is more energetically trapped by vacancy in comparison with Be. On the other hand, the H number trapped by Be is 7 less than that of vacancy.^13–15

It has been demonstrated that there is large binding energy between Be atoms in W, leading to them clustering. The Be clusters should also have trapping effect on H in W. We further investigate the trapping behavior of H surrounding Be clusters. It is found that the H solution energy at the most stable site surrounding a Be₂ cluster is 0.73 eV, which is 0.27 eV lower than that of H in the bulk W. Therefore, the trapping energy of a single H at the Be₂ cluster is −0.27 eV, 0.07 eV lower than that of H surrounding one Be atom, suggesting that Be₂ cluster has stronger tapping effect on H than a single Be atom. This can be attributed that Be₂ cluster induces more significant deformation to W lattice in comparison with one Be atom, as well as larger decrease of electron density. As listed in Table 1, the deformation energy induced by Be₂ cluster is 5.42 eV, which is larger than two times of that induced by one Be atom. This results in that the electron density decreases to as low as 0.186 e Å⁻³ at the H most stable site surrounding Be₂ cluster, which is 0.014 e Å⁻³ lower than that surrounding Be₁.

3.3. Effect of Be-V complex on H in W

Generally, defects are considered as the main origin for the accumulation of impurities in metals. Vacancy is one of the typical defects in metals, which has been demonstrated to be the trapping center for H in W.^10,13–15 We first investigate the stability of Be in vacancy. It is found that Be energetically prefers to occupy the vacancy center with the solution energy of −2.11 eV, 6.31 eV lower than that of Be at the OIS. This suggests that Be is favorable to segregate into the vacancy with the segregation energy of 6.31 eV. Moreover, the segregation energy for a single H in vacancy is only ∼1.38 eV, which is much lower than that of Be in vacancy. Therefore, Be will be more energetically favorable to occupy the vacancy than H.

Next, we investigate the dissolution behavior of H in the Be-vacancy (Be-V) complex. Fig. 8 shows five potential occupation sites for H surrounding the Be-V complex. The solution energy of H at the Be-V complex can be defined as

where E^T_127W,Be-V,H and E^T_127W,Be-V are the total energies of the Be-V supercell with and without H atom, respectively. The H solution energies are listed in Table 3. It can be found that the solution energy at the site 1 is much lower than that at other sites. Therefore, the site 1 (the TIS at the Be-V lattice surface) is the most stable site at the Be-V complex. Here, the distance between H and Be is ∼1.69 Å, which is similar with that for the interstitial case. With the presence of H, Be atom will slightly deviate from vacancy center shown in Fig. S2(a).† Further, the H solution energy is 0.39 eV, 0.61 eV lower than that of H in the bulk W. This suggests that the Be-V complex can serve as a trapping center for H with the trapping energy for a single H of −0.61 eV, which is 0.67 eV larger than that of H in a Be-free vacancy. Therefore, the presence of Be weakens the trapping capability of vacancy for H. This can be attributed to the redistribution of electron density in vacancy induced by Be. In previous study,¹⁰ it has been demonstrated that the Be-free vacancy can provide an isosurface with low electron density of 0.11 e Å⁻³. However, Be atom is very easy to lose its 2s² to the surrounding environment, leading to the significant increasing of electron density in vacancy. The electron density of the H most stable site in the Be-V complex is ∼0.19 e Å⁻³, much larger than that in the Be-free vacancy. The larger the electron density, the larger is the H trapping energy, in consistent with the H solution energy in the homogenous electron gas as a function of the electron density.^58,59

Fig. 8 Five high-symmetric sites for a single H atom at the Be-V complex.

Table 3 The solution energies (eV) of H at different sites surrounding the Be-V complex

Site	1	2	3	4	5
Solution energy	0.39	0.46	1.07	1.16	0.95

---

We further exam the number of H atoms can be trapped by the Be-V complex. Fig. 9(a) shows the H trapping energy in vacancy with and without Be atom as a function of the number of H atoms, consistent with the previous studies.^10,13–15 It can be found that the trapping energy of per H atom will increase with the increasing H atoms number for both cases. Moreover, for the first–fifth H atom, the trapping energy of per H atom in the Be-V complex is ∼0.58–0.72 eV larger than that in the Be-free vacancy, which further proves the negative effect of Be on the H tapping capability of vacancy. Furthermore, the trapping energy of the sixth H atom is −0.69 eV in the Be-free vacancy, while it becomes positive in the Be-V complex (0.12 eV). Therefore, the Be-V complex can only trap 5H atoms, 7 less than that of the Be-free vacancy.^13–15,48 Further, Fig. 9(b) shows the atomic configuration of the Be-V–H₅ complex (Fig. S2(b)†). The shortest distance between these H atoms trapped by the Be-V complex is 1.89 Å, much larger than that in H₂ molecule (0.75 Å), indicating that there is no H₂ molecule formation. Therefore, the presence of Be can severely weaken the trapping effect of vacancy on H as well as block the formation of H₂ molecule in vacancy.

Fig. 9 (a) The trapping energy per H as a function of the number of H atoms in vacancy with and without Be, and the zero point is the energy of H at the TIS far away from vacancy. (b) The atomic configuration of the Be-V–H₅ in W.

Further, temperature has significant effect on the dissolution and retention of H in W. The number of H atoms trapped by defects will decrease with the increasing of temperature. For example, the maximum number of H that can be trapped by a vacancy is 12 at 0 K,^13–15,48 while it will decrease to 6 at RT.^15,16 In the present work, a Be-V complex can only hold 5H atoms at 0 K. Further, the trapping energy of the first H atom in the Be-V complex is −0.61 eV, which is larger than that of the sixth H in the Be-free vacancy (−0.69 eV). This suggests that the Be-V complex can't contain any H atom at RT, thus the Be-V complex can trap H at low temperature below RT only.

3.4. Discussion

The influence of Be on H behavior in W can be divided into two aspects based on our calculations. On the one hand, the presence of Be will decrease the H solution energy in bulk W, which can be attributed to the local lattice distortion and the redistribution of electron density induced by Be atom. Therefore, the interstitial Be atom can serve as the trapping center for H. Although the trapping effect of a single Be on H is very weaker in comparison with that of vacancy, it will strengthen with Be clustering in W. We have demonstrated that there is large binding energy between Be atoms in W (>1 eV), leading to them energetically clustering. This can be attributed to the strong Be–W repulsion and the electrophobic properties. Further, the presence of Be clusters will trap H atoms, and block H diffusing deeper into the bulk, leading to the decrease of H retention in W as well as the formation of H bubbles. On the other hand, Be can significantly weaken the trapping effect of vacancy on H. It has been demonstrated that Be atom is more energetically favorable occupying vacancy than H. The presence of Be will dramatically decrease the H number trapped by vacancy, and block the formation of H₂ molecule in vacancy. Therefore, the existing of Be can suppress H bubble formation in W. Based on these two aspects, we can well understand the effect of Be on H in W.

Although the effect of Be on H is similar with that of He in some respects, the behavior of Be in W is very different from that of He. (i) He irradiation has a significantly destructive action to the W-PFM. There will form He bubbles and nanostructure fuzz on W surface after He irradiation,^26–29,60 which will severely degrade the thermal and mechanical properties of W, thus reducing the lifetime of the W-PFM. However, these damage phenomenon haven't been observed after Be irradiation on W. (ii) The weakening effects of Be on the H trapping capability of vacancy is more significant than that of He. Although the presence of He in vacancy will leads to the increasing of H trapping energy in the He-V complex in comparison with that in He-free vacancy, the He-V complex still can hold 12H atoms,¹¹ the same with that of the He-free vacancy.^13,14 However, the presence of Be in vacancy can significantly reduce the number of H atoms trapped by vacancy, besides increasing the H trapping energy. The Be-V complex only can contain 5H atoms, 58.3% less than that of the Be-free vacancy. This can be attributed to the difference between H–He interaction and H–Be interaction in vacancy. The H–He repulsive interaction drives H to an isosurface with high electron density, resulting in the increasing of H trapping energy. While Be not only can push H to the high electron density isosurface due to their repulsion, but also increase the electron density in vacancy by providing electrons to vacancy. Both of them leads to the severely negative effect of Be on the H trapping capability of vacancy. (iii) Be can more easily diffuse into W than He. According to the calculation results of solution energy and diffusion energy, we can make a rough relative estimate of Be and He diffusion energy barrier from vacuum into W. It is found that Be atom needs to overcome a diffusion energy barrier about 5.28 eV, 0.92 eV lower than that of He. Consequently, Be is the better choice to reduce H retention and suppress H blistering in W in comparison with He.

4. Conclusion

We have investigated the synergistic behaviors of H and Be in W to explore the effect of Be on H using a first-principles method. It is found that a single Be energetically prefers to occupy the OIS in W with a solution energy of 4.20 eV. The OIS → TIS → OIS path is the optimal diffusion path of Be in W with diffusion energy of 1.11 eV. Most interestingly, there is large binding energy between Be atoms in W (>1 eV), suggesting that Be atoms are energetically favorable getting together in W. It can be attributed to the strong Be–W repulsion and the intrinsic electrophobic properties of Be. Further, we have demonstrated that the presence of Be has significant effect on H behavior in W, which can be divided into two aspects. For the interstitial case, the Be atom can serve as the H trapping center with the trapping energy of 0.20 eV for a single H, and a single Be atom can trap 5H atoms. Moreover, the Be clusters have stronger trapping effect on H than one Be atom. This will block H diffusing deeper into the bulk, leading to the decrease of H retention in W. For the vacancy case, the presence of Be severely weakens the H trapping capability of vacancy. The Be-V complex can only hold 5H atoms without H₂ molecule formation, 58.3% less than that of the Be-free vacancy. Therefore, H retention and blistering in W can be suppressed by doping Be. Our calculations will provide a good reference for developing W materials as a PFM.

Conflict of interest

The authors declare no competing financial interest.

Acknowledgements

This research is supported by the National Magnetic Confinement Fusion Program with Grant No. 2013GB109002, the National Natural Science Foundation of China with Grant No. 11405006 and 11675011, and the China Scholarship Council with Grant No. 201506025054.

References

1. R. Causey, K. Wilson, T. Venhaus and W. R. Wampler, J. Nucl. Mater., 1999, 266–269, 467.
2. J. Roth, et al., Plasma Phys. Controlled Fusion, 2008, 50, 103001.
3. G. De Temmerman, J. J. Zielinski, S. van Diepen, L. Marot and M. Price, Nucl. Fusion, 2011, 51, 073008.
4. D. M. Goebel, et al., J. Nucl. Mater., 1987, 145, 61.
5. S. Masuzaki, N. Ohno and S. Takamura, J. Nucl. Mater., 1995, 223, 286.
6. R. A. Causey, J. Nucl. Mater., 2002, 300, 91.
7. T. Shimada, H. Kurishita, Y. Ueda, A. Sagara and M. Nishikawa, J. Nucl. Mater., 2003, 313–316, 204.
8. T. Funabiki, T. Shimada, Y. Ueda and M. Nishikawa, J. Nucl. Mater., 2004, 329–333, 780.
9. Y. Ueda, T. Funabiki, T. Shimada, K. Fukumoto, H. Kurishita and M. Nishikawa, J. Nucl. Mater., 2005, 337–339, 1010.
10. Y. L. Liu, Y. Zhang, H. B. Zhou, G. H. Lu, F. Liu and G. N. Luo, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 79, 172103.
11. H. B. Zhou, Y. L. Liu, S. Jin, Y. Zhang, G. N. Luo and G. H. Lu, Nucl. Fusion, 2010, 50, 115010.
12. H. B. Zhou, Y. L. Liu, Y. Zhang, S. Jin, G. N. Luo and G. H. Lu, Nucl. Fusion, 2010, 50, 025016.
13. K. Ohsawa, J. Goto, M. Yamakami, M. Yamaguchi and M. Yagi, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 184117.
14. L. Sun, S. Jin, X. C. Li, Y. Zhang and G. H. Lu, J. Nucl. Mater., 2013, 434, 395.
15. Y. W. You, X. S. Kong, X. B. Wu, Y. C. Xu, Q. F. Fang, J. L. Chen, G.-N. Luo, C. S. Liu, B. C. Pan and Z. G. Wang, AIP Adv., 2013, 3, 012118.
16. N. Fernandez, Y. Ferro and D. Kato, Acta Mater., 2015, 94, 307.
17. X. S. Kong, S. Wang, X. Wu, Y. W. You, C. S. Liu, Q. F. Fang, J. L. Chen and G. N. Luo, Acta Mater., 2015, 84, 426.
18. K. Ohsawa, F. Nakamori, Y. Hatano and M. Yamaguchi, J. Nucl. Mater., 2015, 458, 187.
19. K. O. E. Henriksson, K. Nordlund, A. Krasheninnikov and J. Keinonen, Appl. Phys. Lett., 2005, 87, 163113.
20. H. B. Zhou, S. Jin, Y. Zhang, G. H. Lu and F. Liu, Phys. Rev. Lett., 2012, 109, 135502.
21. T. Hino, K. Koyama, Y. Yamauchi and Y. Hirohata, Fusion Eng. Des., 1998, 39–40, 227.
22. H. T. Lee, A. A. Haasz, J. W. Davis, R. G. Macaulay-Newcombe, D. G. Whyte and G. M. Wright, J. Nucl. Mater., 2007, 363–365, 898.
23. H. T. Lee, A. A. Haasz, J. W. Davis and R. G. Macaulay-Newcombe, J. Nucl. Mater., 2007, 360, 196.
24. Y. Ueda, M. Fukumoto, J. Yoshida, Y. Ohtsuka, R. Akiyoshi, H. Iwakiri and N. Yoshida, J. Nucl. Mater., 2009, 386–388, 725.
25. C. S. Becquart and C. Domain, J. Nucl. Mater., 2009, 386, 109.
26. F. Sefta, K. D. Hammond, N. Juslin and B. D. Wirth, Nucl. Fusion, 2013, 53, 073015.
27. C. M. Parish, H. Hijazi, H. M. Meyer and F. W. Meyer, Acta Mater., 2014, 62, 173.
28. F. Hofmann, D. Nguyen-Manh, M. R. Gilbert, C. E. Beck, J. K. Eliason, A. A. Maznev, W. Liu, D. E. J. Armstrong, K. A. Nelson and S. L. Dudarev, Acta Mater., 2015, 89, 352.
29. Q. Yang, H. Fan, W. Ni, L. Liu, T. Berthold, G. Benstetter, D. Liu and Y. Wang, Acta Mater., 2015, 92, 178.
30. M. J. Baldwin, R. P. Doerner, D. Nishijima, K. Tokunaga and Y. Ueda, J. Nucl. Mater., 2009, 390–391, 886.
31. R. P. Doerner, M. J. Baldwin, D. Nishijima, J. Roth and K. Schmid, J. Nucl. Mater., 2011, 415, S717.
32. K. Tokunaga, M. J. Baldwin, R. P. Doerner, D. Nishijima, H. Kurishita, T. Fujiwara, K. Araki, Y. Miyamoto, N. Ohno and Y. Ueda, J. Nucl. Mater., 2011, 417, 528.
33. A. Lasa, K. Heinola and K. Nordlund, Nucl. Fusion, 2014, 54, 123021.
34. A. Allouche, N. Fernandez and Y. Ferro, J. Phys.: Condens. Matter, 2014, 26, 315012.
35. A. Mutzke, G. Bandelow and R. Schneider, J. Nucl. Mater., 2015, 467, 413.
36. H. Bergsåker, I. Bykov, Y. Zhou, P. Petersson, G. Possnert, J. Likonen, J. Pettersson, S. Koivuranta, A. M. Widdowson and JET contributors, Phys. Scr., T, 2016, 167, 014061.
37. A. Wiltner and C. Linsmeier, New J. Phys., 2006, 8, 181.
38. R. P. Doerner, J. Nucl. Mater., 2007, 363, 32.
39. C. Linsmeier, K. Ertl, J. Roth, A. Wiltner, K. Schmid, F. Kost, S. R. Bhattacharyya, M. Baldwin and R. P. Doerner, J. Nucl. Mater., 2007, 363, 1129.
40. A. Allouche, A. Wiltner and C. Linsmeier, J. Phys.: Condens. Matter, 2009, 21, 355011.
41. G. Kresse and J. Hafner, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 47, 558.
42. G. Kresse and J. Furthmuller, Phys. Rev. B: Condens. Matter Mater. Phys., 1996, 54, 11169.
43. J. P. Perdew and Y. Wang, Phys. Rev. B: Condens. Matter Mater. Phys., 1992, 45, 13244.
44. P. E. Blochl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953.
45. H. J. Monkhorst and J. D. Pack, Phys. Rev. B: Condens. Matter Mater. Phys., 1976, 13, 5188.
46. G. Henkelman and H. Jonsson, J. Chem. Phys., 2000, 113, 9978.
47. G. Henkelman, B. P. Uberuaga and H. Jonsson, J. Chem. Phys., 2000, 113, 9901.
48. X. S. Kong, Y. W. You, Q. F. Fang, C. S. Liu, J. L. Chen, G.-N. Luo, B. C. Pan and Z. G. Wang, J. Nucl. Mater., 2013, 433, 357.
49. X. B. Wu, X. S. Kong, Y. W. You, C. S. Liu, Q. F. Fang, J. L. Chen, G. N. Luo and Z. Wang, Nucl. Fusion, 2013, 53, 073049.
50. R. Wu, A. J. Freeman and G. B. Olson, Science, 1994, 265, 376.
51. T. Seletskaia, Y. Osetsky, R. E. Stoller and G. M. Stocks, Phys. Rev. Lett., 2005, 94, 046403.
52. T. Seletskaia, Y. Osetsky, R. E. Stoller and G. M. Stocks, Phys. Rev. B: Condens. Matter Mater. Phys., 2008, 78, 134103.
53. X. T. Zu, L. Yang, F. Gao, S. M. Peng, H. L. Heinisch, X. G. Long and R. J. Kurtz, Phys. Rev. B: Condens. Matter Mater. Phys., 2009, 80, 054104.
54. C. Wert and C. Zener, Phys. Rev., 1949, 76, 1169.
55. M. J. Baldwin, R. P. Doerner, D. Nishijima, D. Buchenauer, W. M. Clift, R. A. Causey and K. Schmid, J. Nucl. Mater., 2007, 363–365, 1179.
56. C. S. Becquart and C. Domain, Phys. Rev. Lett., 2006, 97, 196402.
57. A. Lasa, K. Heinola and K. Nordlund, Nucl. Fusion, 2014, 54, 083001.
58. M. J. Puska, R. M. Nieminen and M. Manninen, Phys. Rev. B: Condens. Matter Mater. Phys., 1981, 24, 3037.
59. H. B. Zhou, J. L. Wang, W. Jiang, G. H. Lu, J. A. Aguiar and F. Liu, Acta Mater., 2016, 119, 1.
60. M. Yajima, et al., Plasma Sci. Technol., 2013, 15, 282.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra20430a

---