First-principles energy and vibration spectrum simulations of Cr/V interacting with H in W-based alloy in a fusion reactor

ABSTRACT

We perform first-principles total energy and vibration spectrum calculations to study the effect of Cr/V on the H formation and diffusion at temperature range 300–2100 K in dilute W–Cr/W–V alloy. Temperature and H chemical potential are two important factors to affect the H formation energy and migrating energy. The H formation energy referring to the static/temperature-dependent H chemical potential decreases/increases with the temperature. At each given temperature, the presence of Cr in W reduces the H formation energy, while the existence of V in W has little effect on the H formation energy. The diffusion energy barrier and pre-exponential factor strongly depend on the temperature and increase with the temperature from 300 to 2100 K. Both Cr and V additions in W has a large consequences on H migrating energy. The energy barriers at any given temperature can be reduced by ∼0.05 and 0.10 eV in W–Cr and W–V alloys, respectively. The current study reveals that vibration free energy plays a decisive role in the formation energy and migrating energy of H with the temperature, while the thermal expansion energy has little influence on the H formation energy and migrating energy with the temperature.

KEYWORDS:

• Tungsten-based alloys
• fusion reactor
• hydrogen formation and diffusion
• temperature and H chemical potential
• first-principles calculations

1. Introduction

Tungsten (W) is the potential candidates for the first-wall in fusion reactor because W exhibits the low hydrogen (H) sputtering erosion, the high thermal conductivity and high melting temperature [1–4]. Many efforts have been performed to study the behaviors of H isotopes in W [5–22]. In the role of the first-wall, W has to be subjected to high flux H isotopes, which can adsorb on the W surface and further permeate into the bulk and migrating in the bulk. The retention of H isotopes may result in the modification of physical and mechanical behaviors of W. For example, H in W can cause the embrittlement, decrease of mechanical strength and bubble formations.

Adding alloying elements (AEs) has been considered as an important method to improve the basic physical and mechanical properties of W. However, the AEs will inevitably interact with the impurity H from edge plasma so as to possibly result in the change of microstructure and mechanical properties of W material [23]. Correspondingly, the understanding the interaction between impurity H and AE in W is an extreme important work in the current fusion study.

Up to now, several studies including experiments and computational simulations have addressed some of the relevant issues concerning the interactions of H with AEs (solutes) in W. Experimentally, a few works have been carried out to study the effect of the AEs on the H retention in W, which have demonstrated that the AEs can alter obviously the H retention [24–30]. This suggests that the investigation of interaction between H and AE can play a decisive role in the reliable predication of H recycling and retention in W. Theoretically, Kong et al. recently gave the relevant studies of this subject. They have carried out a series of first-principles simulations to study the interaction of H with interstitial AEs in W [31,32]. Later, they extend their study to the investigation of the interaction between H and substitution AEs in W [33]. Among them, the structures and thermodynamics between AEs and H at the interstitial sites and vacancy are discussed in detail using first-principles static calculations. However, there are still many open questions, particularly those which concern the effect of temperature on the interaction between AE and H in the interstitial lattice sites in W. In other words, the temperature effect has been missed on H thermodynamic properties in these theoretical studies. For instance, how is the H formation energy affected by the temperature, i.e. how does the H formation energy changes with the changing temperature. Furthermore, how are the H diffusion energy barrier and pre-exponential factor influenced by the temperature? In previous most calculations all atomic nuclei are treated as point-like classical particles which can lead to inaccuracies in the description of H diffusion including migrating energy barrier and pre-exponential factor where the temperature effect has been negligible, similar to H thermodynamics in W as well as W-based alloys. Here, we will touch these topics from temperature angle, paying close attention to the influence of temperature on H formation and diffusion in W-based alloys.

In this study, we systematically carry out first-principles total energy and vibration spectrum calculations to investigate the influence of temperature on the formation and diffusion behaviors of H in the interstitial positions in W-based alloys. For W-based alloys, we choose vanadium (V) and chromium (Cr) as the model AEs to study the interaction between AE and H in the interstitial lattice by taking into account the temperature effect. Our choice of choosing V and Cr as the AEs is partly because V has been proposed as a possible AE in W with special research emphasis for the applications in a fusion reactor [34], while Cr is generally regarded as the main alloying component in most transition metals. Our findings demonstrate the significant influences of temperature and alloying addition on the formation and diffusion behaviors of H in W, and provide a reasonable interpretation for the larger scatter of the above-predicted results on H formation and diffusion in W-based alloys.

2. Methodology

2.1. Total energy computations of H–W (W-alloy) supercell

For all the energy calculations, we carried out first-principles simulations using the Vienna Ab-initio Simulation Package (VASP) code [35,36] based on density functional theory (DFT). The electron exchange and correlation were treated within the generalized gradient approximation (GGA) using the Pardew–Burke–Ernzerh form [37–39]. The ionic cores were represented by the projector augmented wave potentials [40,41]. During the geometry optimization, we used the 128-atom supercell containing (4 × 4 × 4) unit cells with the lengths of 12.68 Å in the [100], [010], and [001] directions, respectively. The plane wave energy cutoff and k-point density are both tested for each system to reach an energy convergence of 10⁻³ eV per atom for the electronic minimization. According to a series of convergence test calculations, the energy cutoff of plane wave for W is chosen as 500 eV, while the Brillouin zone of supercell was sampled with the Monkhorst–Pack scheme [42] using a (3 × 3 × 3) k-points mesh. The calculated equilibrium lattice parameter (a) of W lattice is 3.17 Å, which is in good agreement with the corresponding experimental value of 3.16 Å [43]. The Methfessel–Paxton [44] smearing method was used to integrate the Brillouin zone and account for partial occupancies of the metals near the Fermi level with a smearing width of 0.1 eV. During the calculations, the supercell size and atomic positions are relaxed to equilibrium, and energy minimization is converged until the forces on all the atoms are less than 10⁻³ eVÅ⁻¹.

2.2. Computations of formation energy and concentration of H in W and W-alloy: thermodynamic and statistical model

In W or W-based alloys, the formation energy of H in an interstitial position at the finite-temperature can be attributed to the variation of the Gibbs free energy at the temperature T and pressure P [45], (1) $$\begin{array}{c}\hfill G_{\mathrm{H}}^{\mathrm{f}}(T,P)=G_{\mathrm{ Total }}^{\mathrm{H}}(T,P)-G_{\mathrm{ Total }}^{\mathrm{ bulk }}(T,P)-{\mu }_{\mathrm{H}}(T,P),\end{array}$$(1) where G^H_Total(T, P) and G^bulk_Total(T, P) are, respectively, the total Gibbs free energies of the supercell system (W or W-based alloy) with an interstitial H atom and perfect supercell system (W or W-based alloy) at the temperature T and pressure P, and they can be given as (2) $$\begin{array}{c}\hfill G(T,P)=E_{a_T}+F_{a_T}^{vib}\left(T\right)+pV(T,P),\end{array}$$(2) where $E_{a_T}$ is the DFT static energy, which is obtained from the current first-principles calculation as the ground-state energy at a given lattice parameter a_T. As to the lattice parameter a_T, it changes with the increase of temperature and will be clearly explained in Section 3. The term PV can be neglected for solid state systems [45], $F_{a_T}^{vib}\left(T\right)$ is the vibration Helmholtz free energy [45] and in the quasi-harmonic approximation it is written as (3) $$F_{a_T}^{vib}\left(T\right)=\sum _{i=1}^{3N}\frac{\hslash {\omega }_i}{2}+kT\sum _{i=1}^{3N}ln[1-exp(-\hslash {\omega }_i/kT)],$$(3) where ℏω_i is the DFT-obtained vibration energy of normal modes, ℏ is Planck's constant divided by 2π, k is the Boltzmann's constant, and 3 N is the total number of vibration modes.

In Equation (1)(1) $$\begin{array}{c}\hfill G_{\mathrm{H}}^{\mathrm{f}}(T,P)=G_{\mathrm{ Total }}^{\mathrm{H}}(T,P)-G_{\mathrm{ Total }}^{\mathrm{ bulk }}(T,P)-{\mu }_{\mathrm{H}}(T,P),\end{array}$$(1) , the μ_H(T, P) is the chemical potential of H, which is related to temperature T and pressure P. The H chemical potential μ_H(T, P) is generally written as [46–48] (4) $$\begin{array}{c}\hfill {\mu }_H(T,P)={\mu }_H(T=0K)+{\mu }_H(T\ne 0K,P),\end{array}$$(4) where μ_H(T = 0 K) is half of the energy of an H₂ molecule at T = 0 K. The total energy of an H₂ molecule is calculated by putting an H₂ molecule in a cubic vacuum box with 12.68 Å sides and simultaneously carrying out a zero-point energy (ZPE) calculation at T = 0 K. According to the current calculation, the energies of an H₂ molecule without and with ZPE-correction are, respectively, −6.76 eV and −6.45 eV with the H-H bond length of 0.75 Å. Thus, the μ_H(T = 0 K) without and with ZPE-correction are −3.38 eV and −3.255 eV. As to the μ_H(T ≠ 0 K, P), it rests with the temperature and H₂ gas pressure. Experimentally, to study the retention of H isotopes under low energy high flux plasma irradiation at the different temperature in W [49], the background and working procedure of H₂ gas is less than or close to one atmosphere pressure (i.e. P = ∼10⁵ Pa) around the sample holder. In order to well compare with experiments, here we will directly use the temperature-dependent μ_H(T ≠ 0 K, P_1atm), which is given as μ_H(T ≠ 0 K, P_1atm) = −1.75 kTln (T/7.55) provided by Ji et al. in the previous study [50] because this result has been derived in detail by Sugimoto and Fukai [47]. Thus, according to Equation (4)(4) $$\begin{array}{c}\hfill {\mu }_H(T,P)={\mu }_H(T=0K)+{\mu }_H(T\ne 0K,P),\end{array}$$(4) , the total H chemical potential at one atmosphere pressure can be written as (5) $${\mu }_H\left(T,P_{1\mathrm{ atm }}\right)=-1.75kT\mathrm{ ln }\left(T/7.55\right)-3.255\mathrm{ eV }.$$(5)

At the low concentration case, the equilibrium concentration of H dissolving in W and W-based alloys can be established according to the Sievert's law and it is given as (6) $$\begin{array}{c}\hfill S_H=\sqrt{(P/P_0)}exp(\Delta S/k)exp(-G_H^f/kT),\end{array}$$(6) where P and P₀ represent the background and reference pressures, respectively. The reference pressure is generally chosen as one atmosphere pressure. The formation entropy ΔS is equal to −4.7k for H in W and W-base alloys [51].

2.3. Computations of H diffusion in W and W-alloy

According to Arrhenius diffusion equation, the diffusivity of an H atom random migrating between two interstitial sites in metal [52] is given by (7) $$\begin{array}{c}\hfill D=D_{0}exp(-G_H^m/kT),\end{array}$$(7) where D₀ is diffusion pre-exponential factor, G^m_H is the H-migrating energy barrier and it is the difference between the energies of H at the transition state (G^ts_H) and the initial state (G^is_H), i.e. (8) $$\begin{array}{c}\hfill G_H^m=G_H^{ts}-G_H^{is}.\end{array}$$(8)

Wert and Zener’ theory [53] gives that D₀ is defined as $\frac{n}{6}{\left(\lambda \right)}^2\sqrt{\frac{2G_H^m}{m{\lambda }^2}}$ for the impurity atom along the diffusion route in a cubic metal, and here n, λ and m are the number of the nearest neighboring interstitial position, the jumping length and the mass of H atom, respectively.

According to the transition state theory [54], the diffusivity can be also given as (9) $$\begin{array}{c}\hfill D=\frac{n}{6}{\left(\lambda \right)}^2\frac{kT}{\hslash }\frac{Z_H^{ts}}{Z_H^{is}}exp(-G_H^m/kT).\end{array}$$(9)

Equation (9)(9) $$\begin{array}{c}\hfill D=\frac{n}{6}{\left(\lambda \right)}^2\frac{kT}{\hslash }\frac{Z_H^{ts}}{Z_H^{is}}exp(-G_H^m/kT).\end{array}$$(9) is a general expression for the impurity diffusion in solid, where Z^ts_H and Z^is_H are the total partition functions of the transition state and the initial state, respectively. Also, this formula can be solved for the diffusivity employing either classical or quantum mechanical solution for the vibration partition function. Employing the quantum mechanical solution of the vibration partition function, within the harmonic oscillation approximation, the diffusion rate is written as. (10) $$\begin{array}{c}\hfill D=\frac{n}{6}{\left(\lambda \right)}^2\frac{kT}{\hslash }\frac{{\prod }_i^{3N-1}\frac{exp\left(\frac{-\hslash {\omega }_i^{+}}{2kT}\right)}{1-exp\left(\frac{-\hslash {\omega }_i^{+}}{kT}\right)}}{{\prod }_i^{3N}\frac{exp\left(\frac{-\hslash {\omega }_i}{2kT}\right)}{1-exp\left(\frac{-\hslash {\omega }_i}{kT}\right)}}exp(-G_H^m/kT),\end{array}$$(10) where ω⁺_i and ω_i are, respectively, two real normal vibration modes for the transition state and the initial state. N represents the number of the vibrating atoms. Note that along the reaction coordinate, the potential energy exhibits a negative curvature so as to yield an imaginary vibration mode (ω*) at the transition state (saddle point). Therefore, there is one real normal mode less at the transition state than at the initial state. The climbing image nudged elastic band (CI-NEB) method [55] was employed to determine minimum energy migration route for the diffusion of H at the interstitial position.

To specifically illustrate the H vibration modes in W and W-alloys, we evaluate the vibration modes of H at the tetrahedral site (TS) and the octahedral site (OS) as well as transition state in the harmonic approximation. According to the quasi-harmonic approximation method of the lattice oscillation, there are 3 N vibration frequencies mentioned in Equation (3)(3) $$F_{a_T}^{vib}\left(T\right)=\sum _{i=1}^{3N}\frac{\hslash {\omega }_i}{2}+kT\sum _{i=1}^{3N}ln[1-exp(-\hslash {\omega }_i/kT)],$$(3) in a supercell system including N atoms so that each atom exhibits three vibration frequencies. The vibration frequencies of H at ground state are listed in Table 1. In the TS, H has three real vibration modes with two degenerate frequencies and one smaller frequency. This indicates that the TS for H is one stable ground state. For the transition state of H diffusion, there are two real vibration frequencies (normal vibration modes) and one imaginary frequency (imaginary vibration mode) in reference to the negative curvature of the saddle point in the direction of the reaction path, based on CI-NEB method [55]. In the OS, however, H has one real vibration frequency and two imaginary frequencies, meaning that the OS for H is only one saddle-point position (unstable state).

Table 1. The vibration frequencies (ω in unit of THz) of H at the TS (ω_TS) and OS (ω_OS) as well as transition states (ω_Tran) in W and W-based alloy at 0 K

	ωTS	ωTran	ωOS
W	46.1, 46.1, 34.9	63.0, 46.1, i25.1	75.6, i26.1, i26.1
W–Cr alloy	47.8, 39.0, 30.9	64.5, 38.5, i22.2	77.4, i21.8, i23.1
W–V alloy	48.1, 43.1, 33.2	64.9, 42.5, i23.9	77.8, i24.1, i24.8

3. Results and discussion

3.1. Establishment of the temperature-dependent lattice parameter

For the temperature effect in this work, we consider two-part contributions, i.e. the thermal expansion of lattice and vibration free energies of the system. The thermal expansion of lattice rests with the increase of the temperature-dependent lattice parameter. Previously, there were many studies on the thermal expansion behaviors of different structural metals by experimental and theoretical methods. Recently, Lu et al. [56] have given the thermal expansion coefficients of body center cubic (bcc) W using the PARROT module [57], and the obtained data are consistent with the reported experimental results. Although the 0-K lattice parameter of 3.16 Å obtained by Lu et al. is in agreement with that (3.17 Å) calculated in this work, these temperature-dependent lattice parameters given by Lu et al. cannot be directly applicable for our current DFT calculations because of the small difference between Lu et al. and this work. However, we may use the similar analytical method proposed by Lu et al. to correlate the temperature with the lattice parameter. Our temperature-dependent lattice parameter (a_T) can be calculated to use in applying our DFT data by multiplying the ratio a'_T/a'₀ given by Lu et al. with the current optimized 0-K lattice parameter (a₀), i.e. a_T = a₀ · (a'_T/a'₀). Employing the above method, each $E_{a_T}$ will be obtained from first-principles calculations as the ground-state energy. According to Equation (3)(3) $$F_{a_T}^{vib}\left(T\right)=\sum _{i=1}^{3N}\frac{\hslash {\omega }_i}{2}+kT\sum _{i=1}^{3N}ln[1-exp(-\hslash {\omega }_i/kT)],$$(3) , F_vib(T) is related to the vibration modes (ω_i) of the system at a given temperature. The vibration modes and a_T are one-to-one relationship, thus the vibration modes of the system will be re-calculated at each given a_T. After establishing the $E_{a_T}$ and F_vib(T), the formation energy and diffusion energy barrier of H at a given temperature in W and W-based alloys can be calculated by substituting the corresponding ‘$G(T,P)=E_{a_T}+F_{\mathrm{ vib }}\left(T\right)$’ into Equations (1)(1) $$\begin{array}{c}\hfill G_{\mathrm{H}}^{\mathrm{f}}(T,P)=G_{\mathrm{ Total }}^{\mathrm{H}}(T,P)-G_{\mathrm{ Total }}^{\mathrm{ bulk }}(T,P)-{\mu }_{\mathrm{H}}(T,P),\end{array}$$(1) and (8(8) $$\begin{array}{c}\hfill G_H^m=G_H^{ts}-G_H^{is}.\end{array}$$(8) ). Finally, we can get how the formation energies and energy barriers of H in W and W-based alloys alter with the increasing temperature.

To investigate the effect of thermal expansion on H formation and diffusion in W-based alloy, we first calculate the temperature-dependent lattice parameter of W-based alloy. Here, we need point out that only one substitute Cr or V atom is added into W supercell system to form a much dilute W_0.9922Cr_0.0078 or W_0.9922V_0.0078 alloy. Therefore, we believe that one AE Cr or V will have little influence on the total W supercell changing with the temperature so that we can obtain the temperature-dependent lattice parameters of W as those of the dilute W–Cr or W–V alloy. As mentioned above, the method was reported by Lu et al. [56] who characterized the molar volumes of W, which is given as (11) $$\begin{array}{ccc}\hfill V_T& =& V_{0}exp(1.3624\times {10}^{-5}T+4.932\times {10}^{-13}{T}^3\hfill \\ & & +5.5762\times {10}^{-2}{T}^{-1}),\hfill \end{array}$$(11) where V_T and V₀ are the molar volumes at the temperature T and 0 K, respectively. Considering the ratio (a_T/a₀) of lattice parameter, Equation (11)(11) $$\begin{array}{ccc}\hfill V_T& =& V_{0}exp(1.3624\times {10}^{-5}T+4.932\times {10}^{-13}{T}^3\hfill \\ & & +5.5762\times {10}^{-2}{T}^{-1}),\hfill \end{array}$$(11) will further expressed as (12) $$\begin{array}{ccc}& & a_T/a_0=\hfill \\ & & \sqrt[3]{exp(1.3624\times {10}^{-5}T+4.932\times {10}^{-13}{T}^3+5.5762\times {10}^{-2}{T}^{-1})}.\hfill \end{array}$$(12)

According to Equation (12)(12) $$\begin{array}{ccc}& & a_T/a_0=\hfill \\ & & \sqrt[3]{exp(1.3624\times {10}^{-5}T+4.932\times {10}^{-13}{T}^3+5.5762\times {10}^{-2}{T}^{-1})}.\hfill \end{array}$$(12) , we can obtain the temperature-dependent lattice parameters of W and W-based alloys, as listed in Table 2.

Table 2. The ratio (a_T/a₀) of lattice parameter, the lattice parameter (a_T) and the lattice expansion strain (%) of W-based alloy at the relevant temperature

Temperature (K)	aT/a0	aT (Å)	Strain (%)
0	1.0000	3.1715	0
300	1.0014	3.1745	0.14
600	1.0028	3.1789	0.28
900	1.0042	3.1834	0.42
1200	1.0058	3.1883	0.58
1500	1.0074	3.1934	0.74
1800	1.0092	3.1991	0.92
2100	1.0111	3.2053	1.11

3.2. Formation of H in interstitial lattice at ground state in W-based alloys

Due to their large atom size, the AEs are generally believed to substitute W atoms in the lattice site. This has been really confirmed by our current calculations, where we find that the formation energies of the AEs Cr and V in substitute positions in W are 0.43 eV and −0.51 eV, respectively, consistent with recent results [58].

For H atom, there are two possible interstitial positions in perfect W, i.e. TS and OS. The preferable site for one H atom is the TS, which has been confirmed by our previous work [11]. However, it is important to check whether the site preference of H atom alters or not when H atom is located at the vicinity of the substitute Cr or V element in W; due to that this is the cornerstone of the following study. Around the substitute Cr/V in W atomic structure plane supercell (top view from the [001] direction), the different neighbor TS and OS are shown in Figures 1(a) and (b), respectively. Here, we point out that the AE Cr/V is exactly placed into the substitute center of 128-atom supercell system since the supercell is made up of 4 × 4 × 4 unit cells, as given in Section 2.2. Accordingly, the distances of the first nearest neighbor (1nn), 2nn, 3nn, 4nn TS (OS) are $\sqrt{5}a/4$(a/2), $\sqrt{13}a/4$($\sqrt{2}a/2$), $\sqrt{29}a/4$($\sqrt{5}a/2$), and $\sqrt{37}a/4$(3a/2), respectively. We first calculate the formation energies of H initially placed in the 1nn TS/OS of the Cr/V. After structural optimization, we found that the position of the H initially occupying the OS remains nearly unchanged around the Cr/V, whereas the H atom initially placed in the TS has been off the perfect TS close to an OS in the vicinity of Cr/V, and the deviation distance is about 0.23 Å from the original TS. However, it can be noted that the deviation distance is not very large, and thus we still call this site as the quasi-TS in dilute W–Cr and W–V alloys. We further find that the formation energies of H in the quasi-TS are still lower than those in the OS in both W–Cr and W–V alloys. This result fully demonstrates that the H atom is still preferable to occupy the TS over the OS in the vicinity of both AE Cr/V in W. Additionally, we calculated the formation energies of H atom located at the 2nn, 3nn and 4nn TS and OS surrounding the AE Cr/V. The calculated results indicate that the TS is always more energetically favorable than the OS in dilute W–Cr or W–V alloy. Thus, in the following analysis we will only give the corresponding discussions on the situation of H at the 1nn TS in dilute W–Cr and W–V alloys. Our choice of choosing 1nn TS as a research object is mainly because our purpose of present study is to investigate the effect of the AE on the H dissolution in W.

Figure 1. (Top view from the [001] direction) The different neighbor TS (a) and OS (b) around the substitute Cr/V (gray ball) in W (red ball) atomic structure plane supercell. The distances of 1nn, 2nn, 3nn, 4nn TS (OS) are $\sqrt{5}a/4$(a/2), $\sqrt{13}a/4$($\sqrt{2}a/2$), $\sqrt{29}a/4$($\sqrt{5}a/2$), and$\sqrt{37}a/4$(3a/2), respectively. The TS and OS are denoted by triangle and circle, respectively.

3.3. Finite-temperature total H formation energy in reference to the temperature-dependent H chemical potential

Employing Equation (1)(1) $$\begin{array}{c}\hfill G_{\mathrm{H}}^{\mathrm{f}}(T,P)=G_{\mathrm{ Total }}^{\mathrm{H}}(T,P)-G_{\mathrm{ Total }}^{\mathrm{ bulk }}(T,P)-{\mu }_{\mathrm{H}}(T,P),\end{array}$$(1) , we calculate the formation energies of H in the TS with the increasing temperature in W–Cr and W–V alloys. For comparison, we also show the corresponding results in pure W, as plotted in Figure 2. The H formation energies in the TS increase with the increasing temperature in W and W-based alloys (Figure 2(a)). This demonstrates that the heat of solution of H dissolving in W and W-based alloys will increase with the increase of temperature. In other words, the dissolution of H becomes more and more difficult with the increasing temperature. The definition of formation energy states that the negative formation energy should be an exothermic process, while the positive formation energy represents an endothermic process. According to the current results from Figure 2(a), one sees that all the formation energies are positive over the whole temperature range from 300 to 2100 K, meaning that the solution process of H in W and W-based alloys is always endothermic with the considered temperature range. Furthermore, we note that with the increasing temperature, the H formation energy is always lower in dilute W–Cr alloy than in dilute W–V alloy as well as pure W, and the energy difference between H in W–Cr and W–V alloys is close to 0.04 eV. This means that the dissolution of H becomes easier in the vicinity of Cr than in the vicinity of V in W. While the formation energy difference of H in W–V alloy and pure W is very small and only ∼0.015 eV, suggesting that the presence of AE V can hardly enhance the H dissolution over the whole temperature regime. Therefore, the appearance of AE Cr in W should be favorable for the H dissolution, while V has little effect on the dissolution of H in W.

Figure 2. The formation energies of H at the TS with the increasing temperature in W and W-based alloys. (a) The total formation energy in reference to the total H chemical potential μ_H(T, P_1atm) = μ_H(T = 0 K) + μ_H(T ≠ 0 K, P_1atm). (b) The corresponding static formation energy in reference to the static H chemical potential μ_H(T = 0 K) based on Equation (13)(13) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ strain }\right)}^f=E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}-{\mu }_H(T=0K),\end{array}$$(13) and it is related to the lattice expansion. (c) Phonon vibration formation energy referring to the temperature-dependent H chemical potential μ_H(T ≠ 0 K, P_1atm) based on Equation (14)(14) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ vibration }\right)}^f& =& F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)\hfill \\ & & -{\mu }_H(T\ne 0K,P_{1\mathrm{ atm }}),\hfill \end{array}$$(14) .

In order to understand specifically the effect of temperature on the H formation energy in dilute W-based alloys, we decompose the total formation energy into two part energies, i.e. the thermal expansion energy and phonon vibration energy. The thermal expansion contribution is the static energy from the DFT calculation for each strained lattice constant a_T at each given temperature, then the H formation energy G^f_H(strain) can be given as. (13) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ strain }\right)}^f=E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}-{\mu }_H(T=0K),\end{array}$$(13) where ${E_{a_T}}_{(\mathrm{H}-\mathrm{W}-\mathrm{ alloy })}$ and ${E_{a_T}}_{(\mathrm{W}-\mathrm{ alloy })}$, respectively, denote the DFT static energies of H–W-alloy and pure W-alloy system for each strained lattice constant a_T at a given temperature. While the H formation energy of phonon vibration contribution can be written as (14) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ vibration }\right)}^f& =& F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)\hfill \\ & & -{\mu }_H(T\ne 0K,P_{1\mathrm{ atm }}),\hfill \end{array}$$(14) where $F_{a_T(\mathrm{H}-\mathrm{W}-\mathrm{ alloy })}^{vib}\left(T\right)$ and $F_{a_T(\mathrm{W}-\mathrm{ alloy })}^{\mathrm{ vib }}\left(T\right)$, respectively, are phonon vibration energies of H–W-alloy and pure W-alloy system for each strained lattice constant a_T at a given temperature. μ_H(T ≠ 0 K, P_1atm) is equal to −1.75 kTln(T/7.55). Actually, here the H chemical potential has been divided into two parts, i.e. the static μ_H(T = 0 K) and the temperature-dependent μ_H(T ≠ 0 K, P_1atm).

Figure 2(b) shows the changing curves of the H formation energy from the thermal expansion contribution with the increasing temperature in W and W-based alloys, respectively. For all the cases, the H formation energy decreases in a close-linear form with the increasing temperature (strain). This demonstrates that the thermal expansion (note that here the thermal expansion corresponds to the isotropic tensile strain) can reduce the H formation energy in W and W-based alloys. The reduced energies from 300 to 2100 K is about 0.02 eV for H in both W and W-based alloys. Also, it can be clearly seen that with the increasing temperature, the H formation energy is always lower in the dilute W–Cr alloy than in the dilute W–V alloy as well as pure W. The energy difference of the thermal expansion contribution between W–Cr and W–V alloy is always ∼0.04 eV, which is not dependent of the temperature. Such a behavior could be understood according to the linear elasticity theory provided by Zhou et al. [16]. According to their study, the H formation energy under tensile strain (thermal expansion) can be obtained as E^f_ϵ = E_{ϵ = 0}^f + (σ_[100]ϵ_[100] + σ_[010]ϵ_[010] + σ_[001]ϵ_[001])V_{ϵ = 0}, where E^f_{ϵ = 0} is the H formation energy without strain, σ_[mln ] is the thermal expansion-induced lattice stress at each crystal-direction in the equilibrium system, ϵ is the strain, and V_{ϵ = 0} is the volume of the system without strain (here, the system is W/W-based alloy supercell) at equilibrium state. For the present thermal expansion (isotropic strain), ϵ_[100]=ϵ_[010] = ϵ_[001] due to that W and W-based alloys belong to the cubic crystal system, while σ_[100] + σ_[010] + σ_[001] can be approximately equal to $3\bar{\sigma }$ ($\bar{\sigma }$is the average stress and it is given as $\bar{\sigma }=\frac{1}{3}({\sigma }_{\left[100\right]}+{\sigma }_{\left[010\right]}+{\sigma }_{\left[001\right]})$); thus, the H formation energy can be further expressed as $E_{\epsilon }^f=E_{\epsilon =0}^f+3\bar{\sigma }{V}_{\epsilon =0}\epsilon$. It can be seen that there is a linear dependence of H formation energy E^f_ϵ on strain ϵ with a slope $3\bar{\sigma }{V}_{\epsilon =0}$. Due to that the average stress $\bar{\sigma }$ is compressive (negative by convention when the lattice is stretched in our calculations), the H formation energy negatively changes in a close-linear form with the increasing tensile strain in both W and W-based alloys. Therefore, we can predict that the H formation energy from the lattice thermal expansion contribution in W and W-based alloys shows nearly the same behavior under the isotropic strain and therefore to the increase of temperature.

Different from the lattice thermal expansion contribution, the H formation energy from phonon vibration contribution shows an obvious temperature-dependent behavior in W and W-based alloys, as shown in Figure 2(c). Based on Equation (14)(14) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ vibration }\right)}^f& =& F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)\hfill \\ & & -{\mu }_H(T\ne 0K,P_{1\mathrm{ atm }}),\hfill \end{array}$$(14) , the energy reference point is not the H chemical potential μ_H(T = 0 K), but the temperature-dependent μ_H(T ≠ 0 K, P_1atm) = −1.75 kTln (T/7.55). It is very clear to see that the H formation energies in W and W-based alloys increase in quasi-linear form with the increasing temperature, although the phonon vibration energy is relevant with the term of the ∼exp ( − ℏω_i/kT) form given in Equation (3)(3) $$F_{a_T}^{vib}\left(T\right)=\sum _{i=1}^{3N}\frac{\hslash {\omega }_i}{2}+kT\sum _{i=1}^{3N}ln[1-exp(-\hslash {\omega }_i/kT)],$$(3) . This is because the absolute value of phonon vibration energy is far lower than that of H chemical potential μ_H(T ≠ 0 K, P_1atm) at the corresponding temperature. For example, the phonon vibration energies of H at T = 1500 K in W, W–Cr alloy and W–V alloy are only 0.043 eV, 0.009 eV and 0.053 eV, respectively, while the μ_H(T = 1500 K, P_1atm) is close to −1.19 eV. Therefore, according to Equation (14)(14) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ vibration }\right)}^f& =& F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{vib}\left(T\right)\hfill \\ & & -{\mu }_H(T\ne 0K,P_{1\mathrm{ atm }}),\hfill \end{array}$$(14) , the temperature-dependent H chemical potential directly dominates the temperature-dependent H formation energy so that the energy-temperature curve exhibits a close-linear changing trend (Figure 2(c)). Here, it should be pointed out that we have checked vibration entropy that does not affect the temperature-dependent H formation energy in both W and W-based alloy, since the contribution of vibration entropy is nearly negligible for the H formation energy with the rising temperature. For instance, we have estimated the free energy contribution from the vibration entropy of H at T = 900 and 1800 K in W, which are about two orders of magnitude smaller than the H formation energy and thus can be neglected. On the other hand, unlike the situation of the thermal expansion energy contribution, the vibration energy contribution in both W and W-based alloys exhibits nearly the same temperature-dependent behavior (see Figure 2(c)). The energy differences among them can be negligible with the increasing temperature.

3.4. Finite-temperature total H formation energy referring to the static H chemical potential

The above results demonstrate that the total formation energies (two parts: the lattice thermal expansion and phonon vibration) of H increase with the increasing temperature at one atmosphere pressure in W and W-based alloys. It is interesting to take into account whether the changing trend of the total H formation energy with the temperature alters or not when the static H chemical potential μ_H(T = 0 K) is only used as the energy reference point. To answer this question, we have calculated the total H formation energy merely referring to μ_H(T = 0 K) with the rising temperature. Note that although such an analysis might be not in accordance with conventional physical logic, our purpose is here to compare with the results of the above part (Section 3.3). Surprisingly, there will be an opposite conclusion on the total H formation energy with the increasing temperature in both W and W-based alloys, i.e. the total H formation energy decreases with the increasing temperature, as shown in Figure 3(a). Similar to the above analysis method, the total H formation energy can be also decomposed into two part contributions, i.e. thermal expansion (G^f_H(strain)) and phonon vibration (G^f_H(vibration)), which are, respectively, given as. (15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) and (16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16)

Figure 3. The formation energies of H at the TS with the increase of temperature in W and W-based alloys. (a) Here, the total formation energy only referring to the static H chemical potential μ_H(T = 0 K). (b) The corresponding static formation energy in reference to the static H chemical potential μ_H(T = 0 K) based on Equation (15)(15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) and it is related to the lattice expansion. (c) Phonon vibration formation energy. Here, it does not refer to any H chemical potential since it is only the difference between H–W-alloy supercell and pure W-alloy supercell, according to Equation (16)(16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16) .

Equation (15)(15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) is equal to Equation (13)(13) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ strain }\right)}^f=E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}-{\mu }_H(T=0K),\end{array}$$(13) , whereas Equation (16)(16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16) denotes phonon vibration energy difference between H–W-alloy supercell and pure W-alloy supercell, i.e. without referring to any chemical potential.

According to Equations (15)(15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) and (16(16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16) ), we recalculate the formation energy of H at the TS with the increase of temperature in W and W-based alloys. Figure 3(b) gives the formation energy of H in reference to the static H chemical potential μ_H(T = 0 K)with the increasing temperature. It is expected that the H formation energy in both W and W-based alloys from thermal expansion contribution displays the same changing trend as the above results, due to that Equation (15)(15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) is equal to Equation (13)(13) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ strain }\right)}^f=E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}-{\mu }_H(T=0K),\end{array}$$(13) . Employing Equation (16)(16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16) , we further analyze the H formation energy from phonon vibration energy contribution. As clearly seen from Figure 3(c), due to that the vibrational Helmholtz free energy involves the exponential term (∼exp ( − ℏω_i/kT)), the H vibration formation energies in W and W-based alloys exponentially decrease with the increasing temperature, and all the values reduce about 0.40 eV when the temperature changes from 300 to 2100 K. In order to demonstrate the validity of these results, we need to discuss the anharmonic effect on the H formation energy from vibration contribution when the temperature alters from 300 to 2100 K. This is because phonon frequencies are the volume dependence. Physically, Grüneisen constants (γ) can effectively describe the vibration frequency (ω) with the change of lattice volume (V), i.e. $\gamma =-\frac{dln\omega }{dlnV}$ [59], this formula clearly indicates that the vibration frequency decreases with the increase of the lattice volume. Therefore, an expansion of the lattice will result in the softening of vibration modes of all the atoms in material system. Taking H atom in W as an example, the calculated vibration modes of H at the TS have two degenerate frequencies and one smaller frequency, respectively, as listed in Table 1. The first degenerate frequency can be found to shift from 46.12 to 43.83 THz if the volume is expanded by 3.232% (corresponding to the temperature increasing from 0 to 2100 K), and then the H vibration formation energy G^f_H(vibration) = ΔF_aT^vib(T) can increase slightly, ∼0.02 eV, when the vibration frequencies of the system with the lattice parameter of 3.205 Å (2100 K) are employed instead of the system with the lattice parameter of 3.171 Å (0 K). Therefore, according to our calculations, it is reasonable to speculate a relatively weak anharmonic effect for H in W even at high temperature 2100 K. Moreover, to some extent, it is also reasonable to employ the vibration frequencies of the system with the optimized lattice parameter at all temperatures (see Table 1). Similar results have been found for impurity carbon in Pd and Pd-based alloys [60]. At the same time, it can be noted that the H vibration formation energy in W and W–V alloy exhibits nearly the same the value at each given temperature, while the H vibration formation energy in W–Cr alloy is lower than those in W and W–V alloy. The difference of the vibration formation energy of H in W–Cr alloy and W–V alloy (as well as pure W) is always ∼0.03 eV, independent of the temperature. This should be understood by the different vibration modes of H in W–Cr and W–V alloys, as listed in Table 1.

In the above discussions, we gave the change of the vibration formation energy of H at the TS with the temperature in W or W-based alloy. However, a crucial question arises of whether there is any possibility to see that the OS-H becomes more stable than the TS-H with the temperature, especially at high temperature, since the stability of H at the interstitial site can be determined by some kind of balance between amplitude of H vibration and thermal expansion of metal matrix. Therefore, we examined this situation and found that the TS for H atom is always more stable than the OS over the whole considered temperature range in W and W-based alloys. For example, using Equations (15)(15) $$\begin{array}{ccc}\hfill G_{H\left(\mathrm{ strain }\right)}^f& =& E_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}-E_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}\hfill \\ & & -{\mu }_H(T=0\mathrm{K}),\hfill \end{array}$$(15) and (16(16) $$\begin{array}{c}\hfill G_{H\left(\mathrm{ vibration }\right)}^f=F_{a_T\left(\mathrm{H}-\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right)-F_{a_T\left(\mathrm{W}-\mathrm{ alloy }\right)}^{\mathrm{ vib }}\left(T\right).\end{array}$$(16) ), the total H formation energy at the OS in W/W–Cr/W–V are calculated to be 0.88/0.81/0.88 eV and 0.83/0.76/0.83 eV at temperatures T = 1800 and 2100 K, respectively. As potted in Figure 3(a), at temperatures T = 1800 and 2100 K, the total formation energy of H at the TS in W/W–Cr/W–V are 0.44/0.37/0.44 eV and 0.33/0.26/0.33 eV, respectively, which are significantly lower than those at the OS. Therefore, the TS-H is invariably more stable than the OS-H over the whole temperature region 300–2100 K in both W and W-based alloy.

3.5. Solubility of H in W and W-based alloys

The above study employs two different calculated methods to determine the formation energy of H at the TS in W and W-based alloys, correspondingly resulting in different changing trends of the H formation energy with the temperature. This should be attributed to two different choices of H chemical potential. If we choose the static H chemical potential μ_H(T = 0 K)as the energy reference point, the H formation energy will gradually reduce with the increase of temperature. If the energy reference point is chosen as the temperature-dependent H chemical potential μ_H(T, P_1atm), the H formation energy will increase significantly with the increasing temperature. Objectively speaking, the temperature-dependent H formation energy in reference to the temperature-dependent H chemical potential should be more in line with the actual case. In other words, the H formation energy referring to μ_H(T, P_1atm) should increase with the temperature in the actual materials. Physically, it is well known that the total energy of the whole material system will increase with the increasing temperature because the electrons will produce the violent motion with the temperature [61,62]. If H or other impurity wants to live in a higher energy system, it will have to get a much higher energy from the external environment to overcome such a higher energy system. This will directly result in an increase of the H formation energy in the higher energy material system. In view of this, in the following solubility calculations, we only take into account the H formation energy obtained by referring to the temperature-dependent H chemical potential μ_H(T, P_1atm).

Using Equation (6)(6) $$\begin{array}{c}\hfill S_H=\sqrt{(P/P_0)}exp(\Delta S/k)exp(-G_H^f/kT),\end{array}$$(6) , we calculate the H solubility as a function of the reciprocal of temperature from 300 to 2100 K at one atmosphere pressure in W and W-based alloys, as shown in Figure 4. Interstitial H concentrations in both W and W-based alloys increase with the temperature as expected. Generally, the lower the H formation energy is, the higher the dissolving concentration of H in the interstitial site is in material. According to the above calculation, the H formation energy in W–Cr alloy is always lower than that in W–V alloy as well as pure W over the whole temperature range. Thus, interstitial H concentration in W–Cr alloy is always higher than that in W–V alloy as well as pure W, as plotted in Figure 4. The agreement with experiment is an important criterion to validate the theoretical calculations. In pure W, the typical H solubility equation at one atmosphere pressure is experimentally given as S(H/W) = 9.30 × 10⁻³exp(−1.04 eV/kT) at temperatures from 1100 to 2400 K, which was reported by Frauenfelder [63]. This experimental data were usually believed to be the most reliable data and have been widely applicable for many experimental and theoretical researches. It can be noted from Figure 4 that at temperature regime from 1200 to 2100 K, our predicted H solubility exhibits a significant deviation from the experimental data and is about two orders of magnitude lower than Frauenfelder's data. This suggests that in W and W-based alloys, the H solubility should be strongly influenced by some factors, which are not taken into account in the current study. First, it is well known that defects such as vacancies and dislocations as well as stacking faults must appear in material. Both vacancies and stacking faults have been demonstrated to serve as the capturing center of H atoms in many transition metals [11,17,64]. Moreover, these captured H atoms can significantly weaken and break the intrinsic metallic bonds in the vicinity of vacancy and stacking fault, which in turn leads to further accumulation of more H atoms in the vicinity of vacancies and stacking faults. For example, Ding et al. [64] predicted a great accommodating ability of the stacking fault for H atoms to be ∼3.40 × 10¹⁵ cm⁻², agreeing with the experimentally observed H bubble in W-based fusion material. Next, high energy neutron irradiated W as well as W-based alloys leads inevitably to a local strain so as to cause the lattice expansion or compression. Accordingly, strain may change the solubility of H in W and W-based alloys. Recently, Zhou et al. [16] discovered that the H solubility can be enhanced by anisotropic strain in some transition metals. Their results demonstrated a cascading effect of H blistering formation in transition metals, i.e. the H dissolution can result in H blistering formation that causes the local strain that by turns promotes H solubility to further enhance blistering growth. In addition, some transition solutes may still remain in the manufacturing process of W material, accordingly they can strongly interact with H to further affect the H dissolution in W. For instance, Kong et al. [33] carried out first-principles ground-state (0 K) simulations to study the interaction of transition solutes with H in the interstitial lattice as well as vacancy in W and found that most interactions between interstitial H and solutes are attractive, which is consistent with our current result that the temperature-dependent H formation energy can be reduced due to the presence of Cr/V in W, as shown in Figure 2(a). This result also means that some transition solutes may capture multiple interstitial H atoms to enhance the H solubility in W. Moreover, Kong et al. [33] also demonstrated that there are the large positive binding energies among transition solute, H and defect-vacancy, which further suggests that the strong interaction between solute and defect vacancy can significantly affect the H solubility in W. Thus, these above factors can fully interpret why the experimental data of H solubility is higher than the theoretical results from current calculations in defect-free W.

Figure 4. Solubility of H as a function of the reciprocal of temperature at one atmosphere pressure in W and W-based alloys. The experimental solubility of H in W is given from Frauenfelder [63].

3.6. Finite-temperature H diffusion behaviors in W and W-based alloys

Now, we investigate the diffusion behaviors of H at finite temperature in W and W-based alloys. In pure W, the TS for H is the most stable location. Thus, there are two possible migrating routes. In the first, H atom moves to one 1nn TS passing through a ‘saddle-point’ position. Then we name this route as t→t route. In the second path, H atom migrates to a 2nn TS, passing through an OS, correspondingly we call this route as t→o→t route. The detailed definition of migrating route can refer to Ref. 4. For the t→t route, it is well known that the transition state of H is located in a position which is off-middle-position between two 1nn TSs close to an OS. As to the t→o→t route, the OS is just the saddle point so that the diffusion energy barrier is the energy difference of H at the TS and OS. At the ground state, the favorable diffusion path for H atom in W is the t→t route, where the migrating energy barrier is about 0.20 eV, less than the value for the t→o→t route, ∼0.38 eV. Thus, the optimal diffusion path for H atom in W is the t→t route. This has been indeed confirmed by our or other previous works [4,11,13,17]. It is important to take into account whether the optimal diffusion path of H atom alters or not when there is a substitute AE atom nearby in W. To solve this problem, we examine the H energy barrier in the vicinity of Cr or V along the t→t and t→o→t routes, respectively. Figure 5 shows the different neighbor (1nn, 2nn and 3nn) t→t and t→o→t diffusion routes around the substitute Cr/V in W atomic structure plane supercell (top view from the [001] direction). The calculated results indicate that around the AE Cr and V, the H energy barriers for the 1nn t→t route are about 0.10 and 0.15 eV, respectively, less than the values for the 1nn t→o→t route, ∼0.27 and ∼0.32 eV, respectively. We further investigated the energy barriers of H along the 2nn and 3nn t→t and t→o→t routes of the AE Cr or V, and found that the t→t route for H migration is still more favorable than the t→o→t route. Moreover, the diffusion routes are much far away from the AE Cr or V, the H-migrating energy barriers are much close to those in pure W. Thus, similar to the above H dissolution analysis, here we still discuss only the 1nn t→t and t→o→t routes of the AE Cr or V so that we can effectively obtain the influence of the AE on the H diffusion in W.

Figure 5. (Top view from the [001] direction) The different neighbor (1nn, 2nn and 3nn) t→t and t→o→t diffusion routes around the substitute Cr/V (gray ball) in W (red ball) atomic structure plane supercell. The red and black arrows represent t→t and t→o→t routes, respectively. The TS and OS are denoted by triangle and circle, respectively.

Figure 6(a) gives the predicted temperature-dependent energy barriers of H along the t→t route in W and W-based alloys. It is important to note that the energy barrier of H migrating along the t→t route in W increases with the increasing temperature, and the value alters from 0.17 to 0.32 eV with the temperature from 300 to 2100 K, consistent with recent calculated results [17]. With the AE additions, we find that the H-migrating energy barriers in both W–Cr and W–V alloys exhibit the same changing trend with the temperature as in pure W. The H energy barriers increase from 0.06 to 0.22 eV and from 0.12 to 0.27 eV when the temperature changes from 300 to 2100 K in W–Cr and W–V alloys, respectively. These results mean that more and more energy should be required for H migrating between two nearest neighbor TSs in W and W-based alloys. On the other hand, we can note that the alloying Cr and V additions have no effect on the changing trend of H-migrating energy barrier with the temperature, but they have large consequences on the absolute values of H-migrating energy barrier at each given temperature. As shown in Figure 6(a), the presence of V can decrease the H energy barrier to some extent, while the existence of Cr can pronouncedly reduce the H energy barrier. At each given temperature point from 300 to 2100 K, the absolute values can be reduced by ∼0.05 and 0.10 eV in W–V and W–Cr alloys, respectively. Therefore, we deduce that the AE Cr or V addition can be more favorable for H diffusion in W.

Figure 6. Along the t→t diffusion route, the H-migrating energy barriers with the temperature in W and W-based alloys. (a) The total energy barriers. (b) The energy barriers from the lattice expansion contribution. (c) The energy barriers from phonon vibration contribution.

To gain some physical insights on why the H energy barriers increase with the increasing temperature in both W and W-based alloys, we also decompose the total energy barrier into two part, i.e. the thermal expansion and phonon vibration contributions. Based on Equation (8)(8) $$\begin{array}{c}\hfill G_H^m=G_H^{ts}-G_H^{is}.\end{array}$$(8) , the energy barrier from thermal expansion contribution can be written as $G_{H\left(a_T\right)}^m=G_{H\left(a_T\right)}^{ts}-G_{H\left(a_T\right)}^{is}$, while the energy barrier from phonon vibration contribution is given as G^m_H(vibration) = G_H(vibration)^ts − G^is_H(vibration). Note that the migrating energy barrier is different from the formation energy and it is independent of the H chemical potential, owing to that it is the energy difference between H at two nearest neighbor TSs. As plotted in Figure 6(b), the H energy barriers from thermal expansion contribution in W, W–Cr alloy and W–V alloy are invariably 0.20, 0.09 and 0.14 eV, respectively, almost independent of the temperature. However, as seen from Figure 6(c), the H energy barriers from phonon vibration contribution significantly depend on the temperature and increase with the increasing temperature. When the temperature alters from 300 to 2100 K, the H energy barriers in W, W–Cr alloy and W–V alloy increase from −0.03 eV to 0.12 eV, from −0.025 eV to 0.135 eV and from −0.025 eV to 0.13 eV, respectively. These temperature-dependent properties of H diffusion can be explained by a similar analysis to that applicable for the H formation energy above. Also, the obvious relationship between the total migrating energy barrier and phonon vibration energy contribution can clearly demonstrate that phonon vibration energy contribution dominates the increase of total migrating energy barrier of H along the t→t route with the temperature. By a similar validation to that applied above to the H dissolution, we also found a relatively weak anharmonic effect for H energy barrier even at high temperature.

Figure 7 shows the temperature-dependent energy barriers of H along the t→o→t route in W and W-based alloys. Similar to the properties of H along the t→t route, the total energy barrier (Figure 7(a)) and phonon vibration energy barrier (Figure 7(c)) of H also increase with the rising temperature, while the H energy barrier from thermal expansion contribution (Figure 7(b)) hardly depends on the temperature and almost remains constant. However, as observed from Figures 6(a) and 7(a) in both W and W-based alloys, the total energy barriers of t→t route are always lower than those of t→o→t route. Therefore, we can infer that the diffusion route of H in the interstitial lattice chiefly advances through the t→t route in W and W–Cr /V alloy. In addition, the changes of phonon vibration energy barriers for t→t and t→o→t routes are 0.15 and 0.24 eV, respectively, when the temperature rises from 300 to 2100 K, as shown in Figures 6(c) and 7(c). Such an obvious difference is mainly originated from that the transition states of t→t and t→o→t routes exhibit two and one real vibration modes, respectively, and thus the real vibration modes of transition state of t→t route are one more than that of t→o→t route, as listed in Table 1. Note that the OS site is exactly the transition state for H diffusing along the t→o→t.

Figure 7. Along the t→o→t diffusion route, the H-migrating energy barriers with the temperature in W and W-based alloys. (a) The total energy barriers. (b) The energy barriers from the lattice expansion contribution. (c) The energy barriers from phonon vibration contribution.

Along the t→t route in W and W-based alloys, the H diffusivity (D) as a function of reciprocal of temperature is presented in Figure 8. Although the diffusion energy barrier and pre-exponential factor are relevant with temperature, the H diffusivities in both W and W-based alloys still present the nearly quasi-Arrhenius changing trend with the reciprocal of temperature. This suggests that both pre-exponential factor D₀ and energy barrier G^m_H should exhibit the temperature-independent fitted values, which corresponds to exhibit the fitted Arrhenius diffusion equation D = D₀exp(−G^m_H/kT) m²/s, as potted in Figure 7. Experimentally, many researchers have reported the H diffusivity in W. Among them, Frauenfelder gave the H diffusivity as D = 4.10 × 10⁻⁷exp(−0.39/kT) m²/s at temperature range from 1100 to 2400 K [63]. This result was usually believed to be the most credible data, since it was provided at a broad temperature regime and was almost less affected by both surface and defect-capturing influence in W. In the current study, our predicted H diffusivity changes from 3.46 × 10⁻⁸ to 6.25 × 10⁻⁸ m²/s with the temperature from 1200 to 2100 K, while the experimental H diffusivity alters from 0.95 × 10⁻⁸ to 5.02 × 10⁻⁸ m²/s at the same temperature regime. Therefore, the order of magnitude (∼10⁻⁸ m²/s) of our predicted H diffusivities are in quantitative agreement with experiment. Such a consistency in order of magnitude with experiment for interstitial H diffusion in W suggests that the temperature effect by considering the thermal expansion and phonon vibration contributions is an extremely important factor to accurately interpret experiments and further predict H diffusion behaviors in other metals and metal-alloys. However, for the front coefficients, there are still certain differences between experiment (3.46 and 6.25) and calculation (0.95 and 5.02). This implies that the interstitial H diffusivity in W should be also affected by some factors such as vacancy, dislocation, local strain and possible transition solute, similar to the interstitial H solubility discussed above. Of course, it can be pointed out that the effect of defect distribution on interstitial H diffusivity in metal is still the hardest problem of theoretical estimation. In addition, due to that the experimental data of H diffusivity in W–Cr and W–V alloys are not available up to now, we thus provide a fairly good reference in the future study of interaction between H and W-based alloys.

Figure 8. Along the t→t route, the diffusivity (D) of H as a function of reciprocal of temperature in W and W-based alloys.

4. Conclusions

Based on first-principles total energy and vibrational spectrum calculations, we investigate systematically the finite-temperature H formation and diffusion behaviors in W and W-based alloys. For the temperature effect involving in the H formation energy and diffusion energy barrier, we take into account the contributions of the two parts, i.e. the thermal expansion of lattice and phonon vibration energies. The specific conclusions are as follows.

As to the H formation energy calculations, the tetrahedral interstitial site is only considered in both W and W-based alloys, due to that the tetrahedral interstitial site is always more favorable for H occupying than the octahedral interstitial site. We show that the H formation energies in W and W-based alloys either increase or decrease with the increasing temperature from 300 to 2100 K. There are two cases occurring for the H formation energy with the temperature: one through choosing a static H chemical potential as energy reference point and the other through choosing a temperature-dependent H chemical potential as energy reference point. If we choose the static H chemical potential as the energy reference point, the H formation energy will decrease with the increase of temperature. If the temperature-dependent H chemical potential is chosen as the energy reference point, the H formation energy will increase significantly with the increasing temperature. In both cases, phonon vibration free energy contribution plays a key role in the H formation energy with the temperature, while the thermal expansion energy alteration has little influence on the H formation energy with the temperature. At each given temperature point, the H formation energy is always lower in dilute W–Cr alloy than in dilute W–V alloy as well as pure W, while the H formation energy in dilute W–V alloy is nearly the same as in pure W. These suggest that the presence of AE Cr in W can promote the H dissolution, while the existence of V in W has little effect on the H dissolution. Using the obtained H formation energy referring to temperature-dependent H chemical potential, we further predict the interstitial H solubility at one atmosphere pressure in W, W–Cr alloy and W–V alloy. Our predicted H solubility in W is about two orders of magnitude lower than the reported experimental data at the temperature range from 1200 to 2100 K. Such discrepancy between theoretical results and experimental data is mainly originated from that we do not take into account practical factors such as vacancies, dislocations and local strain, which can effectively capture large amount of H atoms so as to enhance the H solubility in W.

For H diffusion in W and W-based alloys, we only take into account the optimal migrating route t→t route. The H-migrating energy barriers in both W and W–Cr/V alloys significantly increase with the increasing temperature. When the temperature alters from 300 to 2100 K, the H energy barriers in W, W–Cr alloy and W–V alloy increase from 0.17 to 0.32 eV, from 0.06 to 0.22 eV and from 0.12 to 0.27 eV, respectively. This indicates that the H diffusion in both W and W-based alloys becomes more and more difficult with the increasing temperature. Similar to the H formation, phonon vibration free energy contribution plays a crucial role in the H-migrating energy barrier with the temperature, while the thermal expansion energy contribution has little influence on the H-migrating energy barrier with the temperature. At each given temperature point, the alloying Cr or V addition in W has a large consequences on H-migrating energy barrier. The energy barriers can be reduced by ∼0.05 and 0.10 eV in W–Cr and W–V alloys, respectively. Based on the obtained H-migrating energy barriers, the H diffusivities are predicted in both W and W-based alloys. With the temperature from 1200 to 2100 K, the calculated H diffusivity in W alters from 3.46 × 10⁻⁸ to 6.25 × 10⁻⁸ m²/s. The order of magnitude is in quantitative agreement with the reliable experimental data from 0.95 × 10⁻⁸ to 5.02 × 10⁻⁸ m²/s at the same temperature regime. Due to that the experimental data of H diffusivity in W–Cr and W–V alloys are not available up to now, thus the current predicted H diffusivity in W–Cr and W–V alloys might provide a fairly good reference in the future study of interaction between H and W-based alloys.

Disclosure statement

No potential conflict of interest was reported by the authors.

Additional information

Funding

This work was supported by the National Magnetic Confinement Fusion Program [grant number 2015GB109001] and National Natural Science Foundation of China [grant number 11575153].