Surface stability and H adsorption and diffusion near surfaces of W borides: a first-principles study

Understanding the behavior of tungsten boride (W x B y ) surfaces in a fusion reactor environment is an important topic since boronization is a common wall conditioning method used in fusion Tokamaks. We report the results of density functional theory calculations that investigate the surface stability of W x B y (tetragonal I41/amd-WB, hexagonal P63/mmc-WB2 and tetragonal I4/m-W2B) with low-index orientations, as well as hydrogen (H) energetics near W x B y surfaces. For single element terminated W x B y surfaces, B terminated surfaces are more energetically stable than W terminated as a result of significant reconstruction of B. The H surface adsorption energy and activation energy of H diffusion penetration below W x B y surfaces are mainly related to the outer termination. Specifically, the WB(001) surface terminated with two B layers, referred to as WB(001)-TBB, has higher H adsorption affinity and lower H diffusivity on this surface than other terminations, which is controlled by the significant charge transfer from B to H. However, B atoms on the WB2(0001)-TBB surface decrease both H adsorption and diffusivity on the surface, but enhance H diffusion below the surface in comparison to W terminated WB2(0001) surface. H would be trapped and diffuse within atomic surface gaps on the WB2( 21ˉ1ˉ0 ) surface, while H below the surface layer would jump along the [0001] direction rather than diffuse into bulk. The surface diffusion activation energy of H on the W2B(001) surface slightly varies with terminations. Once H crosses the surface layer of W2B(001) with either termination, it prefers to diffuse into the bulk, or back towards the surface, rather than move parallel to the surface. Interestingly, WB2(0001) and WB2( 21ˉ1ˉ0 ) surfaces will have relatively higher H retention than the other W x B y surfaces evaluated in this work.

Understanding the behavior of tungsten boride (W x B y ) surfaces in a fusion reactor environment is an important topic since boronization is a common wall conditioning method used in fusion Tokamaks. We report the results of density functional theory calculations that investigate the surface stability of W x B y (tetragonal I4 1 /amd-WB, hexagonal P6 3 /mmc-WB 2 and tetragonal I4/m-W 2 B) with low-index orientations, as well as hydrogen (H) energetics near W x B y surfaces. For single element terminated W x B y surfaces, B terminated surfaces are more energetically stable than W terminated as a result of significant reconstruction of B. The H surface adsorption energy and activation energy of H diffusion penetration below W x B y surfaces are mainly related to the outer termination. Specifically, the WB(001) surface terminated with two B layers, referred to as WB(001)-T BB , has higher H adsorption affinity and lower H diffusivity on this surface than other terminations, which is controlled by the significant charge transfer from B to H. However, B atoms on the WB 2 (0001)-T BB surface decrease both H adsorption and diffusivity on the surface, but enhance H diffusion below the surface in comparison to W terminated WB 2 (0001) surface. H would be trapped and diffuse within atomic surface gaps on the WB 2 (2110) surface, while H below the surface layer would jump along the [0001] direction rather than diffuse into bulk. The surface diffusion activation energy of H on the W 2 B(001) surface slightly varies with terminations. Once H crosses the surface layer of W 2 B(001) with either termination, it prefers to diffuse into the bulk, or back towards the surface, rather than move parallel to the surface. Interestingly, WB 2 (0001) and WB 2 (2110) surfaces will have relatively higher H retention than the other W x B y surfaces evaluated in this work.
Keywords: tungsten boride surface, hydrogen, adsorption, diffusion, DFT calculation (Some figures may appear in colour only in the online journal) * Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
To solve future world energy demand, nuclear fusion energy has been investigated for decades. Most fusion research employs a magnetic confinement fusion device known as a tokamak, where the condition of plasma-facing materials (PFMs) is key for the application of fusion energy [1][2][3][4][5][6][7]. Tokamaks containing tungsten (W) plasma facing components have demonstrated better plasma confinement following wall conditioning using boron (B) [8][9][10][11][12][13][14][15][16][17][18][19][20]. During the deposition of B on W surfaces, W borides (W x B y ) may form near the surface with increasing B surface coverage. Experimental studies have indicated the formation of W borides in boronized W [21][22][23][24]. As such, the PFMs that contain W boride surfaces will be exposed to high fluxes of hydrogen (H) isotopes and helium ions, as well as neutron irradiation during a burning fusion plasma operation. Therefore, understanding how W boride surfaces behave in the reactor environment and their interaction with fusion performance is necessary. In addition to neutron irradiation damage and transmutation impurities, sub-surface bubble formation and surface blistering from the plasma surface interaction will result in changes in thermal and physical properties which could significantly degrade the performance of PFMs [25][26][27][28]. In this work, we focus on investigating the stability of W boride surfaces, and the H behavior on, or below, these surfaces.
The phase diagram of W borides is complex with many possible compounds, even for a fixed stoichiometry [29][30][31][32][33][34][35][36][37]. Previous studies have revealed that WB has the lowest average formation enthalpy in W boride compounds [30,31] and single-phase WB was observed to form near W surfaces after pure W was boronized in experiments [21], we therefore focus on the surface stability of low-index WB surfaces, in addition to the H adsorption and diffusion behavior near these WB surfaces. The energetics of intrinsic point defects and H in six ground-state W borides with various stoichiometries, such as W 2 B, WB, WB 2 , W 2 B 5 , WB 3 and WB 4 , have been studied using density function theory (DFT) in our previous work [38] and it was found that the H behavior in bulk P6 3 /mmc-hP6-WB 2 is significantly different from I4 1 /amd-WB. In addition to the second lowest average formation enthalpy and a relatively high formation energy of Frenkel pairs, P6 3 /mmc-hP6-WB 2 has the lowest H formation energy and the lowest H binding energy to single vacancies, as well as nearly isotropic H diffusion among the six evaluated W borides [38]. In consideration of these outstanding properties of WB 2 and the presence of W 2 B as W is coated by low B content [24], we also studied the stability of low-index surfaces of ground-state WB 2 and W 2 B, as well as H behavior near the relatively stable surfaces of these two W borides. Further, we have analyzed the effects of surface termination on surface stability and H behavior for the surfaces of WB, WB 2 and W 2 B. All studies in the present work have been performed by DFT calculations, and the obtained results are helpful in understanding H behavior near the surfaces of the likely formed W borides as a result of fusion tokamak wall conditioning using boron.

Methodology
The Vienna Ab initio Simulation Package (VASP) [39], along with the pseudopotentials of W (6s 2 5d 4 ), B (2s 2 2p 1 ) and H (1s 1 ), has been used in this work. The interaction between ions and electrons was described by the projector-augmented wave method [40], while the generalized gradient approximation in the Perdew-Burke-Ernzerholf form [41] was used for the exchange and correlation interactions. The energy cutoff for the plane-wave basis set is 360 eV. The cutoff energy of 360 eV for PAW-PBE and 350 eV for PAW-PW91 W pseudopotential has been widely used in previous W-related DFT calculations using VASP [42][43][44][45][46][47]. The Gaussian smearing method was applied for the Fermi surface smearing with a width of 0.1 eV. The Gamma centered Monkhorst-Pack scheme was used to sample the Brillouin zone and generate k-point grids. The surfaces evaluated in this work are based on the stable structures of tetragonal WB (I4 1 /amd), hexagonal WB 2 (P6 3 /mmc-hP6) and tetragonal W 2 B (I4/m). The optimized lattice parameters for WB a 1 = b 1 = 3.138 Å and c 1 = 16.957 Å; for WB 2 a 2 = b 2 = 2.926 Å, c 2 = 7.747 Å, and for W 2 B a 3 =b 3 = 5.579 Å and c 3 = 4.753 Å. The structures of these three W borides are available in previous work [31,38].
The stability of low-index W boride surfaces, as well as H adsorption and diffusion near these surfaces has been investigated using a unit slab and supercell slab geometries, respectively, and both consist of multiple atomic layers and a 1.2-2.0 nm vacuum layer.
The surface formation energy (γ) of W x B y is calculated by [48] Here, E total slab is the total energy of the surface slab. N W and N B are the number of W and B atoms. µ W and µ B are the chemical potential of W and B in W x B y , respectively, which are available in [38]. A is the surface area. To gain convergent surface energies, the thickness of the slab and the density of kpoint meshes were increased gradually until the difference of the surface formation energy is within 10 meV. In general, the lower the surface energy, the more stable the surface [49].
The H adsorption and diffusion near the low-index stable surfaces has been further studied using supercell slabs. Based on the thickness convergence of the surface energies, the supercell slab thickness ranges from 1.6 to 2.0 nm. The dimension for each surface slab without counting vacuum above the surface in the H behavior investigation will be provided in section 3. A 4 × 4 × 1 k-point grid was applied to each supercell slab of W boride surfaces except for WB 2 (2110), which used a 4 × 5 × 1 k-point grid. The bottom few layers of the slab (thickness is about 0.4 nm) were fixed, while the remaining layers were relaxed freely. All calculations were performed at constant volume, and using period boundary conditions along the two directions on the lateral plane.
The definition of H adsorption energy on surfaces is available in [50,51]. Negative adsorption energy indicates attraction between surface and H, and a higher adsorption energy corresponds to a lower adsorption ability of surfaces to H. The climbing image nudged elastic band (CI-NEB) method [52] was used to investigate the energy barriers of H diffusion near W boride surfaces. Three or five images were inserted between the initial and final states, and a minimum energy path between given initial and final states of a transition was gained by the CI-NEB method. All calculations in the present work were performed with the convergence criteria for the electronic selfconsistent iteration and the forces on the unfixed atoms of 10 −5 eV and 0.1 eV nm −1 , respectively.
For compound systems, surface characteristics vary with different terminations even for the same index surface [53,54]. We found that the evaluated low-index WB, WB 2 and W 2 B surfaces terminated with single element layers can be identified by the index and elements on the top two surface layers except for W 2 B(001), and thus, W x B y (hkl)-T ij or W x B y (hklm)-T ij is used to represent a surface with terminated elements i and j, where h, k, l and m are Miller-indices, and the i and j correspond to the element on the first and second surface layer, respectively. The W or B terminated W 2 B(001) surface is referred to as W 2 B(001)-T W or W 2 B(001)-T B , respectively, while W x B y (hklm) without extension represents a surface with a termination of both W and B elements on the surface layer. Note that if the distance between atoms normal to a plane is less than 0.05 nm, the atoms are identified to be on the same planar layer.

Surface stability
A WB(001) surface has four termination combinations, including T BW , T WB , T BB and T WW . The simulations demonstrated that most atoms on the surface layer move towards the bulk and some atoms to a less extent move on the lateral plane parallel to the surfaces during structure optimization. The relaxed structures of the WB(001)-T ij surfaces with different terminations and the WB(100) surface are shown in figure 1. The amount of the change in interlayer spacing for the top three layers relative to the value before relaxation (∆d nm /d nm ) and the surface formation energy for each evaluated WB surface are provided in table 1. Negative and positive values of ∆d nm /d nm indicate contraction or expansion of the interlayer spacing, respectively. It is clear that the interlayer spacing between the top two layers of the WB(001)-T ij surfaces, especially WB(001)-T BW , WB(001)-T WB and WB(001)-T BB , is significantly contracted while the spacing between the 2nd and 3rd layers is slightly expanded after relaxation. The surface atoms on WB(100) move slightly on the lateral plan in addition to a weak expansion of the interlayer spacing in comparison to the (001) orientation, as shown in figure 1 and table 1. In general, a large surface modification corresponds to a lower surface formation energy. However, while the magnitude of the change in the interlayer spacing for the WB(001)-T BW is slightly smaller than that of the WB(001)-T BB , the surface energy of WB(001)-T BW is lower than the latter. According to the creation process of surfaces, the surface energy should be attributed to both cleavage and relaxation [53]. To understand the effect of the cleavage and relaxation on the surface formation energy, the projected electronic density of states (PDOS) of atoms in the WB(001)-T BW and WB(001)-T BB surface systems with and without relaxation, as well as in bulk WB, has been calculated. Figure 2 presents the PDOS of unrelaxed surface and bulk B atoms (top panels), as well as the PDOS of a surface B atom with and without relaxation (bottom panels). A comparison of PDOS between the surface (without relaxation) and bulk B atoms is indicative of the contribution of the cleavage to surface formation. As seen in the top panels of figure 2, the PDOS of the s and p states of the surface B atom dramatically shift toward a higher energy level in comparison with the PDOS distribution of a bulk B atom, which is consistent with a significant contribution of the cleavage to the surface formation. A larger cleavage contribution corresponds to a higher surface formation energy [53]. The slight shift of DOS for both s and p states of the surface B atoms toward a lower energy following the relaxation, as shown in the bottom panels in figure 2, indicates that the surface relaxation decreases the surface formation energies. The similar difference in the PDOS distribution of surface B atoms with and without relaxation for the WB(001)-T BB and -T BW surfaces implies that the relaxation effect on the surface stability is similar for the two terminations of WB(001), which is consistent with the interlayer spacing changes. However, in addition to the shift of PDOS, the profile of PDOS of the B atom on the WB(001)-T BW surface is only slightly changed relative to the bulk B atom, while the PDOS profiles are much different for the B atom on the WB(001)-T BB surface and in the bulk, as shown in the top panels in figure 2. Moreover, the DOS of the B-p state on the WB(001)-T BB surface near the Fermi energy level is significantly higher than the bulk B atom, as illustrated in figure 2(b). In short, the cleavage contribution to the surface stability of WB(001)-T BB is stronger than within the WB(001)-T BW surface, resulting in a higher surface formation energy of WB(001)-T BB than WB(001)-T BW . Based on the surface formation energy, the WB(001) surface prefers to be terminated with B atoms.
The ground-state bulk WB 2 with the lowest formation enthalpy has a structure of ReB 2 with a space group of P6 3 /mmc [31,38]. The surface structures of the (0001) orientation with different terminations and the (2110) surface of P6 3 /mmc-hP6-WB 2 are shown in figure 3. The quantity of ∆d nm /d nm and γ for the evaluated WB 2 surfaces are presented in table 2. The WB 2 (0110) surface was not included because of the instability of the surface during structure optimization. Similar to the WB(001) surface terminated with two B layers, the top two B layers of the WB 2 (0001)-T BB surface contract significantly, which to some extent contributes to a lower surface formation energy of T BB than either the T WB or T BW terminated (0001) surface. In addition, it is noticeable that both the WB 2 (0001)-T BB and WB 2 (0001)-T BW surfaces are terminated with a B surface layer, although their formation energies are significantly different. The cleavage of WB 2 (0001)  Table 1. Surface structures and formation energy (γ) of W B surfaces. ∆dnm/dnm (n = 1, 2, m = n + 1) is the change in interlayer spacing between two neighboring atomic layers relative to the spacing before relaxation. surfaces with terminations of T BB or T BW primarily involves breaking W-B and B-B bonds, respectively. The high formation energy of WB 2 (0001)-T BW indicates the B-B bonds are much stronger than W-B bonds in WB 2 , which is proved by the charge distribution in WB 2 [31,38]. Among the four lowindex surfaces of WB 2 , the WB 2 (2110) surface has the lowest formation energy.

WB(001)-T BW WB(001)-T WB WB(001)-T BB WB(001)-T WW WB(100)
In the case of I4/m-tI12 ′ -W 2 B, the B atoms between two adjacent (001) W atomic layers are treated as one B layer, since the normal distance between these B atoms along the [001] direction is less than 0.03 nm. The structures of the (001) and (100) surface orientations for W 2 B with different terminations are shown in figure 4. Table 3 lists the relative change in the interlayer spacing of the top three layers and     Interestingly, for the evaluated WB, WB 2 and W 2 B surfaces terminated with single element layers, the most stable surface is terminated with B. This is likely because B atoms at the surfaces are more prone to reconstruct than W atoms. To some extent, the second layer below the surface also plays an important role in the surface stability. The surface formation energy consists of the cleavage energy and relaxation energy. For similar surface cleavage contributions, the reconstruction of B surface atoms can more significantly decrease the surface energy than W surface atoms.

H adsorption and diffusion near W boride surfaces
3.2.1. WB(001) surface. According to the formation energy of low-index WB surfaces, H adsorption and diffusion near the relatively stable WB(001) surfaces, including WB(001)-T BW , T WB and T BB , has been investigated using a supercell surface slab. The possible adsorption sites on WB(001)-T BW are shown in figure 5(a), where, the TO site is on the top of the surface B atom and the 4F site is above the center of four surface B atoms. Both of the BR and HL sites are located between the B-B centers, but they are above the third B or third W layer, respectively. The H adsorption energies at the 4F, TO, HL and BR sites on the WB(001)-T BW surface are provided in the second column in table 4. These adsorption energies indicate that H prefers to occupy the 4F site, followed by the TO site. Hydrogen adsorption at HL or BR sites is unstable/less stable due to a positive adsorption energy.
The H diffusion behavior near the WB(001)-T BW surface has been investigated using the NEB method implemented in VASP and the results are shown in figure 5. The H diffusion pathways on the surface are indicated by the arrows in figure 5(a), and the corresponding energies are shown in figure 5(b). H at the preferred 4F site can migrate to the TO site with an energy of 0.42 eV, while H at the TO site will move to the nearby 4F site easily due to a low energy barrier of 0.03 eV. At relatively high temperatures, H at 4F sites can jump to the neighboring 4F sites by passing through a HL or BR site, with an energy of 0.91 and 1.07 eV, respectively. The near zero energy barriers from HL or BR to the 4F sites implies that B at these two sites will move to neighboring 4F sites if it is displaced from either the HL or BR sites. In addition to diffusion on the surface, the inward diffusion of H below the surface is important to quantify. H at the preferred 4F sites on the WB(001)-T BW surface can move to the subsurface site S 1 at temperatures higher than 900 K because of an energy barrier of 2.48 eV, as shown in figures 5(c) and (d). The required diffusion temperature is estimated by equation (7) in [38] according to the work of Janotti et al [55]. On the pathway from 4F to S 1 sites, H will pass by HL sites. The next jump for H penetration below the surface (S 1 →S 2 ) requires an energy of 1.98 eV, which is slightly higher than the energy barrier of H diffusion along the same pathway in bulk WB (1.89 eV) [38]. This is similar to the energy map of H perpendicular diffusion below W(100), (110) and (111) surfaces, in which the  [56]. Figure 5(e) presents the energy map for H parallel migration relative to the surface. Note that the energy map of H diffusion between two identical S 1 or S 2 sites along the LW 12 and LW 34 pathways in figure 5(e) corresponds to the preferred pathways, as respectively indicated by the magenta and orange lines in figure 5(c).
Here, LW nm (n = 1 and 2, m = n + 1) represent the diffusion pathways those are located between the nth and mth W layers. The energy barrier of H parallel diffusion along LW 12 below the WB(001)-T BW surface is only 0.14 eV, which is significantly lower than the energy barrier for H jumping from S 1 to S 2 (1.98 eV) and also lower than from S 1 back to the 4F surface sites (0.44 eV). This implies that after penetration the top two atomic layers of the surface, H at the S 1 site prefers to diffuse between the first and second W layers below the surface. In addition to the WB(001) surface with terminations of T BW or T WB , the H adsorption and diffusion near the WB(001)-T BB surface has also been investigated. The four possible adsorption sites on the WB(001)-T BB surface are presented in figure 7(a). As shown in the fourth column in table 4, H prefers to adsorb at TO sites above the surface B atoms with an adsorption energy of −1.33 eV, while the adsorption energies at the other three sites (4F, HL and BR) are positive. Comparing the H adsorption energies on the WB(001) surfaces with different terminations, it is clearly seen that the WB(001)-T BB surface has stronger attraction for H atoms than the other two terminations, i.e. WB(001)-T BW and WB(001)-T WB . H can migrate between two neighboring TO sites along two pathways indicated by the red solid or dashed lines in figure 7(a). The migration energies in figure 7(b) demonstrate that H jumping along the pathway TO→HL→TO or TO→BR→TO has an energy barrier of 1.39 eV and 1.73 eV, respectively. The energy profile of H diffusion on the WB(001)-T BB surface demonstrates that H cannot stay at BR or HL sites due to a zero energy barrier of jumping to the nearby TO sites.
Although the surface formation energy of WB(001)-T BB is slightly larger than WB(001)-T WB , the WB(001)-T BB surface is stable as H penetrates below the surface. Like on the WB(001)-T BW surface, H direct migration from the preferred adsorption TO sites on the WB(001)-T BB surface to the sublayer site S 1 has an energy barrier of 2.40 eV, which can occur at temperatures higher than 900 K. The two steps of the H motion from S 1 →S 2 →S 3 have energy barriers of 0.41 and 0.23 eV. It is noticeable that the jump from S 3 to S 4 requires an energy of 1.74 eV, which is slightly lower than the bulk value (1.89 eV). The highest energy barrier for the penetration diffusion of H along the preferred pathways TO→S 1 →S 2 →S 3 near the WB(001)-T BB surface is 2.40 eV. As seen in figures 7(c) and (d), the energies barriers for the outward diffusion from S 4 to the surface (S 4 →S 3 →S 2 →S 1 →TO) are 1.38, 0.03, 0.20 and 0.19 eV, respectively.
The parallel H diffusion between LW 12 along the [100] direction below the WB(001)-T BB surface, specifically through the site below the B atoms on the second layer, has the lowest  Table 5 indicates that the H diffusion activation energies on the three different terminated WB(001) surface (T BW , T WB and T BB ) with a 2D diffusion mechanism are 0.42, 0.37 and 1.73 eV, respectively. Obviously, the WB(001)-T BB surface has a much stronger adsorption ability and a significantly higher diffusion activation energy for H than the other two terminated WB(001) surfaces. This demonstrates that boron-rich environment of the WB(001)-T BB surface enhance H trapping and impede H diffusion, which may be due to the formation of strong H-B bonds on the surface. Considering the instability of WB(001)-T WB as H penetrates below the surface, the remaining discussion will focus on the WB(001) surface with either T BB or T BW terminations. Figure 8 shows the electronic charge distribution near H and its neighboring atoms for H adsorbed at the site of 4F or TO on the WB(001)-T BW surface and at the TO site on the WB(001)-T BB surface. To quantitatively show the charge transfer due to the H adsorption, the Bader charge analysis [57] has been performed, and the effective charge of H and host atoms on the top two layers are provided in units of |e| near atoms. Note that the charge transfer is not the only criterion for understanding the strength of a bond; however, charge transfer, as indicated by the charge distribution, can be used to understand the interaction difference between atoms in similar systems. As seen in figure 8(a), H adsorption at the 4F site impacts the charge distribution of its neighboring atoms, especially the first nearest neighboring (1NN) W, and thus H mainly interacts with the W atom on the second layer, forming a W-H bond with a length of 1.71 Å. Figure 8(b) demonstrates that H adsorption at the TO site on the WB(001)-T BW surface is mostly attributed to the interaction between H and its 1NN B, and most of the charge of the 1NN B transfers to H, resulting in the effective charge change of the 1NN B from −0.80 to −0.29 |e|. It is surprising that the adsorption energy of H at the TO site is −0.11 eV, which is larger than at the 4F site (−0.50 eV) on the same WB(001)-T BW surface. This is likely due to the fact that H adsorption at the TO site weakens the interactions between the neighboring B and W atoms (see figure 8(b)), leading to an energy increase of the surface slab. Likewise, H at the TO site on the WB(001)-T BB surface strongly interacts with the neighboring surface B atom, as indicated by the significant charge accumulation near H in figure 8(c). In general, B atoms in bulk WB [38]    the WB(001)-T BB surface is more significant than H at the 4F or TO site on the WB(001)-T BW surface. A comparison of figures 5 and 7 reveals that the energy maps for H perpendicular and parallel diffusion below the WB(001) surface are to some extent different for the two terminations of T BW and T BB . In particular, the energy barrier of H parallel diffusion between the first and second W layers (LW 12 ) below the WB(001)-T BW is clearly lower than below the surface terminated with two B layers, as shown in table 5. Moreover, there are more metastable sublayer sites for H solution below the WB(001)-T BB surface than below the WB(001)-T BW surface. However, the activation energies of H perpendicular diffusion below the two different terminated WB(001) surface are similar (2.48 eV for T BW vs 2.40 eV for T BB ). In addition, the accumulation tendency at LW 12 or jumping back to the surface for H atoms below the WB(001) surface is similar for the two different terminations. Thus, the H behavior on the WB(001) surface is sensitive to the surface termination of two layers, but the H diffusion below the WB(001) surface only slightly depends on the two terminations of T BW and T BB due to the same outer-most termination of B.

WB 2 (0001) and (2110) surfaces
In consideration of the relative surface stability of low-index WB 2 surfaces, we are interested in H behavior near the WB 2 (0001) surface with terminations of T BB or T WB and the WB 2 (2110) surface. Four adsorption sites, including TO, FCC, HCP and NBC, were considered for H adsorption on the WB 2 (0001) surface, as shown in figures 9(a) and 10(a), where the NBC site is near the bond center of two neighboring surface atoms and the . The higher energy barrier along the pathway TO→FCC→TO indicates lower possibility for H diffusion than along TO→NBC→TO. The zero energy barriers from FCC or NBC to TO sites reveal the instability of H at FCC and NBC sites, similar to the HCP site. In addition to diffusion on the surface, H at the TO sites can jump to the sublayer site S 1 below the WB 2 (0001)-T BB surface with an energy barrier of 1.78 eV, as shown in figures 9(c) and (d). Similar to H diffusion in bulk WB 2 [38], H can easily jump between S 1 and S 2 sites by crossing a W layer due to an energy barrier of less than 0.05 eV. However, the H parallel diffusion between the second B layer and the first W layer (S 1 →S 1 ′ in figures 9(c) and (e)) requires an energy of 0.94 eV. As expected, the activation energy for H diffusion from S 2 to S 2 ′ increases to 1.26 eV, approaching the bulk value of 1.22 eV [38]. As illustrated in figures 9(c) and (d), the energy barrier for H indirect jumping from S 2 to S 3 is 1.11 eV, while the direct jumping from S 2 to S 3 requires an energy of 1.43 eV. Same as the situation in bulk WB 2 , H near the (0001) surface is inferred to avoid the For the WB 2 (0001) surface terminated with W at the top layer and B at the second layer, i.e. WB 2 (0001)-T WB , the FCC adsorption site is preferred for H, closely followed by the HCP site, as seen from table 6. H at the NBC site moves to a nearby FCC site after relaxation. In contrast to the WB 2 (0001) surface terminated with two B layers, the TO site on the WB 2 (0001)-T WB is energetically unfavorable for H adsorption. Figures 10(a) and (b) demonstrate that H can diffuse easily between neighboring FCC sites by passing through HCP sites on the WB 2 (0001)-T WB surface at room temperature due to an activation energy of 0.29 eV. While, H diffusion below the WB 2 (0001)-T WB surface requires an energy of at least 2.77 eV from FCC to S 1 sites. Also, H diffusion along the [0001] direction prefers the indirect pathway, as indicated by the blue arrows in figure 10(c). With increasing depth, the energy barrier from S 2 to S 3 decreases to 0.99 eV. Normally, the activation energy for H parallel diffusion increases with depth below the surface, as seen in the WB 2 (0001)-T BB , WB(001)-T BB and WB(001)-T BW surfaces. Unexpectedly, the energy barrier of H parallel jumping between S1 and S 1 ′ below the WB 2 (0001)-T WB surface is higher than between S 2 and S 2 ′ , as illustrated in figures 10(c) and (e). In addition, it is interesting that the energy barrier from S 2 to S 3 (0.99 eV) is lower than the parallel diffusion activation energies (1.31 and 1.21 eV for S 1 →S 1 ′ and S 2 →S 2 ′ , respectively), and it is also lower than the energy barrier for H jumping back to the surface from S 1 to FCC sites (1.24 eV), which suggest that once H crosses through the top three atomic layers of the WB 2 (0001)-T WB surface, it prefers to reside at S 1 and S 2 sites, and then diffuse deeper below the surface towards the WB 2 bulk. Further, the high energy barrier from S 1 to FCC sites will block the H diffusion from WB 2 bulk to surface and subsequent desorption, resulting in relatively high H retention below the WB 2 (0001)-T WB surface. The H behavior near the WB 2 (2110) surface has also been investigated in this work, based on the relatively low surface formation energy of this surface. Many possible sites for H adsorption on the WB 2 (2110) surface have been considered, and the adsorption energies for stable or metastable sites are displayed in table 6, where the NBL 2 , NWL 2 , BBL 2 and HL sites are indicated in figure 11(a) by stars, squares, triangle and crosses, respectively. Clearly, H on the WB 2 (2110) surface prefers to occupy the NBL 2 sites, then followed by the NWL 2 and BBL 2 sites. The energetics of H diffusion on the surface have been calculated using the NEB method in VASP, and the energy profiles along five H jumping pathways between two NBL 2 sites, as shown in figure 11(a), are provided in figure 11(b). The energy barrier of H diffusion along the P1 path is 0.67 eV, which is similar to the barrier of 0.69 eV along P2. The energy barriers for P3, P4 and P5 pathways are 0.87, 1.29 and 1.62 eV, respectively. It is noticeable that there are atomic gaps on the WB 2 (2110) surface, as highlighted by two black dashed lines in figure 11(a). According to the energy map, H prefers to reside near the gap edges and diffuse within the surface gap. An energy of 1.62 eV is required for H jumping between gaps. For H diffusion penetration below the surface, the energy barriers vary with different jumping steps. Figures 11(d) and (e) show the minimum energy profiles of some considered jumping paths. H at a preferred site (NBL 2 ) on the surface can jump into the sub-surface by either path P1 sub (NBL 2 →S 1 →S 3 ) or P2 sub (NBL 2 →S 2 → S 4 ) since the first step of the two pathways has a similar energy barrier of 2.74 eV and 2.76 eV, respectively. The NEB calculations for H jumping from S 1 to S 2 sites did not converge. H at the S 1 site would diffuse along the [0001] direction with an energy of 0.64 eV rather than jump to S 3 with an energy of 0.84 eV, as shown in figures 11(c)-(e). Alternatively, on the pathway P2 sub , H at the NBL 2 site might indirectly jump to S2 by passing through a NWL 2 site with a same energy saddle as the direct jump. H migration from S 2 to S 4 requires an energy of 1.21 eV, which approaches the bulk value of 1.22 eV along the same pathway [38]. In addition, the activation energy of H diffusion along the [0001] direction at the depth of S 2 is 1.08 eV, which is only 0.03 eV higher than the bulk value. Note that the diffusion along the [0001] direction is indirect, similar to the perpendicular diffusion below the WB 2 (0001) surface. In addition, it is of interest to see from figure 11(d) that the energy barrier of a H jump from the S 1 to NBL 2 sites is 0.89 eV, which is higher than the barrier of 0.84 eV from S 1 to S 3 . Likewise, on the P2 sub path, the energy barrier of 1.43 eV for H jumping from S 2 to NBL 2 sites is higher than the barrier of 1.21 eV from S 2 to S 4 . This implies that H atoms tend to accumulate at the S 1 and S 2 sites below the WB 2 (2110) surface, similar to below the WB 2 (0001)-T WB surface.
The activation energies of H diffusion near the three WB 2 surfaces are summarized in table 7. Comparing to the WB 2 (0001)-T BB surface, the top W layer of the WB 2 (0001)-T WB surface decreases both the adsorption energy and surface diffusion activation energy of H, but increases the energy barriers of H perpendicular and parallel diffusion below the surface, as shown in tables 6 and 7. In other words, the boronrich environment on the WB 2 (0001)-T BB surface decreases the H adsorption and diffusion on the surface, but enhances the H diffusion right below the surface in comparison to the WB 2 (0001)-T BW surface. Because of the surface atomic gap, the H diffusion on the WB 2 (2110) surface is more difficult than on the WB 2 (0001) surface. In addition, the activation surface is similar to below the WB 2 (0001)-T WB surface, while the activation energy of H parallel diffusion below the WB 2 (2110) surface is significantly lower than below the WB 2 (0001) surfaces, as shown in table 7.

W 2 B(001) surface.
As described in the introduction, W 2 B might form as W is coated by boron at a low B content. The H behavior near the B terminated W 2 B(001) surface, i.e. W 2 B(001)-T B , has been studied in this work since it has the lowest surface formation energy among the five evaluated W 2 B surfaces. Five possible sites on the W 2 B(001)-T B surface, as shown in figure 12(a), were considered for H adsorption, including the TB, TW, HL, BR and 3F sites. The TB and TW sites are above the first B layer and the first W layer, respectively. The HL, BR and 3F sites are indicated by the stars, square and triangles in figure 12(a), respectively. The H adsorption energies at these sites are listed in the second column in table 8. According to these adsorption energies, H prefers to occupy HL sites on the W 2 B(001)-T B surface, closely followed by BR sites. Comparing to the WB and WB 2 surfaces evaluated above, the W 2 B(001)-T B surface demonstrates a noticeably weak H adsorption ability. The NEB calculations reveal that H can easily jump between two nearest neighboring HL sites (HL→1NN HL) by passing      To assess the effect of surface termination on H behavior, the H energetics near the W 2 B(001) surface terminated with W has also been investigated by DFT calculations. The H adsorption sites and energies on the W 2 B(001)-T W surface are shown in figure 13(a) and table 8. As seen in table 8, the four adsorption sites (HL, BR, 3F and TB) on the W 2 B(001)-T W surface all demonstrate an attraction to H, and H prefers to adsorb at the HL sites. The TW site on this surface is not stable for H following relaxation. The H migration pathways between two preferred sites on the W 2 B(001)-T W surface and the corresponding energies are shown in figures 13(a) and (b). The energy barriers for H jump between two 1NN HL sites by passing through a BR site is 0.35 eV, which is higher than the barrier of 0.23 eV for jumping between two 2NN HL sites by passing two 3F sites. Despite a negative adsorption energy, the zero-energy barrier from BR or 3F sites to adjacent HL sites indicates the instability of H at the BR or 3F sites. A comparison of H adsorption behavior on the W 2 B(001) surface with a termination of W or B indicates the H preferred adsorption site is independent on the termination although the adsorption energy does vary with termination. In contrast to the B terminated W 2 B(001) surface, H at HL sites on the W 2 B(001)-T W surface would jump to the 2NN HL rather than the 1NN HL. However, the H activation energies for 2D diffusion on the W 2 B(001) surface are similar for B or W termination, as shown in table 9. H at HL sites on the W 2 B(001)-T W surface can jump to S 1 sites below the surface with an energy barrier of 1.78 eV, which is clearly higher than the barrier of 1.19 eV below the W 2 B(001)-T B surface. For depths deeper than S 1 sites, the activation energies of H perpendicular diffusion into the W 2 B bulk are similar for the two terminations. The NEB calculations for H jumping between the first and second W layers below the W 2 B(001)-T W did not converge. However, the energy barriers of perpendicular diffusion are evidently lower than parallel diffusion between adjacent W layers, which is independent on the surface termination of W 2 B(001). In summary, H adsorption energies and the critical energy barrier of H diffusion penetration below the surface are related to the termination of the W 2 B(001) surface, while the H diffusion activation energies on the surface and below the first stable subsurface sites only slightly vary with termination.

Conclusion
The surface stability of ground-state WB, WB 2 and W 2 B with low-index orientations has been evaluated by DFT calculations. For the evaluated low-index WB, WB 2 and W 2 B surfaces terminated with single elements, B terminated surfaces are more stable than W terminated. The reconstruction of B surface atoms plays a key role in the surface formation energy.
In consideration of the surface formation energy of lowindex WB, WB 2 and W 2 B surfaces, the H adsorption and diffusion near the relatively stable WB(001), WB 2 (0001) and WB 2 (2110), as well as W 2 B(001) surfaces has been investigated. The DFT calculations demonstrate that the H adsorption energy and the activation energy of H diffusion penetration below the W boride surfaces are sensitive to the outer-most termination. Moreover, the energy barrier of H parallel diffusion below the surfaces of W borides is sensitive to the top few surface layers.
Comparing to the WB(001) surface terminated with either T BW or T WB , the boron-rich environment on the WB(001)-T BB surface enhances H adsorption but impedes H diffusion on the surface with lower adsorption energies and higher diffusion activation energy barriers. The significant charge transfer from B to H can account for the strong adsorption ability of the WB(001)-T BB surface to H atoms. The parallel diffusion tendency near surfaces or jumping back to the surface for H atoms below the WB(001) surface is similar for the termination of T BW or T BB . B atoms on the WB 2 (0001)-T BB surface decrease both H adsorption and diffusion on the surface, but enhance H diffusion below the surface by decreasing the activation energy in comparison to the WB 2 (0001)-T WB surface. H on the WB 2 (0001) surface would rather diffuse into the bulk than move parallel to the surface. While, H on the WB 2 (2110) surface prefers to reside near the surface atomic gap edges and diffuse within the gap. H below the surface layer of WB 2 (2110) would jump along the [0001] direction rather than diffuse into bulk. The high energy barrier from the sublayer to surface of both WB 2 (0001) and WB 2 (2110) block the H desorption from bulk to surfaces, resulting in relatively high H retention below the surface, especially for WB 2 (0001)-T WB and WB 2 (2110).
The preferred adsorption site and the surface diffusion activation energy of H on the W 2 B(001) surface do not, or at most, only slightly depend on the surface termination, although the H adsorption energy does depend on surface termination. The activation energy of H diffusion penetration below the B terminated W 2 B(001) surface is 1.19 eV, which is lower than the barrier of 1.78 eV below the W terminated surface. Once H cross the top surface of W 2 B(001), H would rather diffuse into bulk than move parallel to the surface, which is independent on the surface termination.