Energetics of silicon in the bulk and near surfaces of tungsten: a first-principles study

Siliconization of the tokamak walls is a candidate method to improve plasma confinement in fusion tokamaks containing tungsten plasma facing components (W PFCs). To understand the interactions of silicon (Si) with W, the Si behavior in bulk W, and near three low-index W surfaces ((100), (110) and (111)) has been investigated using first-principles density functional theory. In bulk W, Si interstitial atoms have a low solution ability and high mobility, and Si atoms can be strongly trapped by W vacancies. The interaction between two Si adatoms is responsible for the stability of adatom superstructures on W surfaces, consistent with previous experimental observation (Tsong and Casanova 1981 Phys. Rev. Lett. 47 113). Although the coverage dependence of Si adsorption and diffusion energetics on surfaces is related to surface orientation, the W(110) surface has lower Si adsorption affinity and higher Si diffusivity than either the W(111) or W(100) surfaces. The most stable Si adatom superstructure on W surfaces is: square c(2 × 2) pattern on W(100) covered with 0.5 ML Si; rectangular c(4 × 2) pattern on W(110) with 0.25 ML Si; and rhombus p(1 × 1) pattern on W(111) with 1 ML Si. The coverage dependence of Si mobility on/toward W surfaces is generally related to the stability of the Si superstructures as a function of coverage on each surface. Interestingly, Si adatoms prefer to transport below the surface and into W subsurface by an exchange mechanism with W atoms, indicating the likelihood of epitaxial growth of W silicide layers on W surfaces during the operation of W PFCs.

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 meet increased energy demands through the world, fusion energy has been investigated for decades; however, the development of energy-producing fusion reactor still faces many issues, such as the material challenge associated with the interaction between the plasma and the plasma-facing components (PFCs) in the divertor.Particularly, arcing is a persistent plasma-surface interaction issue in nuclear fusion devices, which can cause surface erosion and contamination of the core plasma by ejecting dusts and droplets into the plasmas [1][2][3].Arcing has been observed and investigated in several tokamaks [3][4][5][6][7][8].To improve tokamak performance, silicon (Si) has been considered as a wall conditioning element [9][10][11][12][13][14][15][16][17][18].Siliconization of the tokamak wall has been experimentally observed to suppress impurities and reduce particle recycling [10][11][12][13][14]. Tungsten (W) is a promising divertor candidate in the future fusion reactors due to its good mechanical and thermal properties [15][16][17][18][19]. Tang et al [10] investigated the influence of Si coatings on W erosion in arc plasma and confirmed the W substrates can be effectively protected by Si coatings from arc erosion.Therefore, understanding the fundamental materials science of the interaction of Si with W surfaces is helpful for improving, or optimizing, the Si coating on the W divertor.
The adsorption of semiconductors on metallic substrates, including Si adsorption on tungsten surfaces, has been investigated to explore the formation mechanism of metal silicides and field emission [20][21][22][23][24][25][26][27][28][29][30].Tsong and Casanona [22,23] studied the behavior of Si adatoms on the W(110) surface using the field ion microscope (FIM), estimated the binding energy of a pair of Si adatoms on the W(110) surface from the direct FIM observations, and predicted the minimum-energy superstructure from the binding energy, which was subsequently confirmed by their thermal equilibration experiments with deposited Si adatoms at a low degree of coverage at about 300 K.In addition, Sault and Goodman [24,25] investigated the Si superstructure formed on the W(110) surface as a function of Si coverage through reactions of silane with the surface using temperature programmed desorption, Auger electron spectroscopy and low energy electron diffraction (LEED).The adsorption of Si on polycrystalline tungsten surfaces has also been studied by thermal desorption [29].These studies provided some information about Si behavior on W surfaces; however, most were focused on Si behavior on the W(110) surface.Moreover, the Si LEED patterns as a function of Si coverage reported in [24,25] might be impacted by the hydrogen atoms in silane.To assess the surface orientation and coverage dependence of Si behavior, the W(100), W(110) and W(111) surfaces covered with different Si contents, have been studied in this work.
Previous experiments have demonstrated that Si adatoms at high coverages on W surfaces can dissolve into W bulk, forming silicide WSi 2 across a temperature range of 1000-1400 K [26], which corresponds to an activation energy of about 3 eV.Moreover, repeated cycles of silane adsorption followed by annealing to 1050 K resulted in the epitaxial growth of stoichiometric WSi 2 layers on the W(110) surface [24].Therefore, to reveal the formation mechanism of W silicide as Si atoms are adsorbed on W surfaces, we have modeled the Si diffusional penetration below W surfaces with either an interstitial or an exchange mechanism using density functional theory (DFT), in addition to evaluating the adsorption and diffusion of Si adatoms on W surfaces.The reported results in this work will be useful in understanding Si behavior near the W surfaces as a result of fusion tokamak wall conditioning using Si.

Methodology
In this work, the Vienna Ab initio Simulation Package (VASP) [31] has been used to investigate the Si behavior in bulk W and near W surfaces, in which the pseudopotential of W contains six valence electrons (6s 2 5d 4 ) and Si contains four valence electrons (3s 2 3p 2 ).The projector-augmented wave method was used to describe the interaction between ions and electrons [32], while the generalized gradient approximation in the Perdew-Burke-Ernzerholf form [33] was employed for the exchange and correlation interactions.The Gaussian smearing method was applied for the Fermi surface smearing with a width of 0.1 eV.The energy cutoff for the plane-wave basis set is 360 eV, which corresponds to the optimized lattice of 0.3171 nm (a) for body-centered cubic (bcc) W.
A cubic bcc supercell of 4a × 4a × 4a or 6a × 6a × 6a was used to obtain the energetics of Si in bulk W. A slab geometry with a 2 nm vacuum was used to model W surfaces.Considering accuracy and the computational costs, the size of surface slabs varies with different goals, as shown in table 1.Except the atoms at the bottom three, two and four layers for W (100), (110) and (111) surfaces, respectively, the remaining atoms have been allowed to relax into their minimum energy positions.A surface reconstruction model with zigzag chains, in which the top surface layer atoms shift 0.027 nm and the second layer atoms shift 0.004 nm along the <011> directions, was considered for the W(100) surface [34].The Brillouin zone was sampled using the Monkhorst-Pack scheme [35], and the k-point grid for all surface slabs is 3 × 3 × 1 and the number of grids for bulk calculations is provided in table 2.
The average adsorption energy of Si on W surfaces relative to a Si crystal (E ads ) or the average desorption energy (E des ) required for adsorbate Si far away from the surface is calculated as [36]: where n is the number of Si atoms, E tot W and E tot W+nSi is the total energy of the slab before and after Si adsorption.If µ Si is the average energy of diamond cubic Si, equation ( 1) is E ads [20]; if µ Si is the energy of an isolated Si atom at the ground electronic state, equation ( 1) is E des .A larger adsorption energy means that Si atoms are more energetically adsorbed on the surface.
To assess the adsorption ability of surfaces to additional Si atoms, the sequential adsorption energy is calculated by  direction.E m represents the migration energy barrier of a Si interstitial atom between two adjacent OISs.E b is the binding energy of Si to a mono-vacancy.ISIF2 indicates that the atomic positions are relaxed, and ISIF3 means the atomic positions, cell volume and shape are allowed to change.where m−n is the number of additional adatoms.

E tot
W+nSi and E tot W+mSi is the total energy of surfaces covered with n or m Si adatoms, respectively.
The energy barriers of Si migration in bulk W and near W surfaces have been obtained by the climbing image nudged elastic band (CI-NEB) method [37] implemented in VASP.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 determined by the CI-NEB calculations.The calculations related to surfaces were performed at constant volume and with a periodic boundary condition along the two lateral directions.All calculations were converged by the criteria for the electronic self-consistent iteration and the forces on the unfixed atoms of 10 −5 eV and 0.1 eV nm −1 , respectively.

Si in bulk W
First, the behavior of Si in bulk bcc W has been investigated using a 128-atom and 432-atom supercell, respectively.The tetrahedral and octahedral interstitial sites (TISs and OISs), as well as dumbbell configurations were considered.Note that the dumbbell configuration consists of an interstitial Si atom and a nearby W atom.The [100] dumbbell configuration transformed into an OIS after relaxation.The energetics of Si in bulk W is presented in table 2. The solution energy of an interstitial Si obtained using the 4 × 4 × 4 supercell only slightly depends on the use of a 3 × 3 × 3 or a 6 × 6 × 6 k-point grid for the same relaxation method.The following discussion will be based on the results of the 4 × 4 × 4 and 6 × 6 × 6 supercell with a k-point grid of 6 × 6 × 6 or 4 × 4 × 4, respectively.The difference in the solution energy of a Si interstitial atom at a given site is less than 2% for the two supercells with atomic position relaxation (ISIF2).For the full relaxation (ISIF3), the solution energy of Si also slightly depends on the supercell size used in this work, except for the [111] mixed W-Si dumbbell configuration.In general, a smaller available volume change for an interstitial defect in a larger supercell after the relaxation with ISIF3 results in a smaller decrease of the solution energy relative to the relaxation with ISIF2.The DFT calculations demonstrate that the OIS is the most stable site for Si interstitials, closely followed by the TIS.The Si-W dumbbell complex along <111> or <110> directions may appear because of the only slightly larger (∼8%) solution energy compared to the Si OIS.The formation of Si-W dumbbell complexes is likely due to the comparable radius of Si (1.10 Å) with W (1.35 Å) 3 and nearly equal electronegativities of Si and W [38].The migration energy barrier between two first nearest neighboring (1NN) OISs approaches 0.2 eV.The high binding energy of Si to a mono W vacancy (∼10 eV) indicates that Si will be strongly trapped by W vacancies.It is noticeable that the interaction between Si and W is different than between carbon (C) and W in bulk W, although Si and C atoms have the same number of valence electrons [39,40].The solution energy of Si at OISs in W is significantly higher than C in W [39] (6.96 eV vs 0.78 eV).The solution energy of interstitial impurities can be divided into a mechanical and electronic contribution [41,42].The mechanical contribution to the solution energy for Si at OIS is calculated to be 74% by the method proposed in [42].The significantly higher solution energy of Si than C in W is thereby due to the fact that the radius of Si is larger than C and the solution energy of interstitial defects is primarily attributed to the mechanical interaction between Si and W atoms.The adsorption of a Si adatom at different sites on the W (100), ( 110) and (111) surfaces has been studied, and the adsorption energies calculated by equation ( 1) are provided in table 3 (recall equation (1) defines stronger adsorption with a larger adsorption energy).On the reconstructed W(100) surface, the 4F sites are above the second layer, while LB and SB indicate the long and short bridge sites, respectively.For the W(110) surface, TOP and LB indicate the sites above the first and second atomic layers, respectively, while the SB site is above the mid-point of two 1NN atomic sites on the surface layer.TOP, HCP and FCC on W(111) are the adsorption sites above the first, second and third atomic layers, respectively.A detailed representation of these adsorption sites can be found in [36].As shown in table 3, Si on W(100) prefers to occupy the 4F site with an adsorption energy of 2.15 eV.Si on W(110) resides at the LB site rather than the SB and TOP sites.For the case of W(111), the FCC site is the preferred occupation site for Si adatoms.Note that these preferred Si adatom sites are the same as W adatoms on the three W surfaces [43,44].The maximum adsorption energy of Si on W(100) is larger than W(111), which has a larger Si adsorption energy than W(110).The stability order of the three surfaces is W(110) > W(111) > W(100) [45].This indicates the more stable surface has a lower adsorption ability to Si atoms relative to the other two surface orientations.
To assess the thermal desorption activation energy of Si adatoms away from W surfaces, the desorption energy relative to the isolated Si in vacuum was calculated by equation ( 1) and the data are presented in table 3.Although the adsorption energies of Si at the TOP site on W(110), as well as at the HCP and TOP sites on W(111) have negative values, the desorption energies are still larger than 2.8 eV, which indicates that electronic interactions between Si adatoms and W surfaces exist for the adsorption sites evaluated in this work.

A pair of Si adatoms.
In addition to a single Si adatom, the behavior of Si with different coverages adsorbed on W surfaces has been investigated.Previous experiments have shown that Si superstructures with various patterns were formed on the W(110) surface covered with different densities of Si [22][23][24].It was believed that adatom-adatom interaction was responsible for the adatom layer (adlayer) superstructure formation [22,23].To explore the surface orientation dependence of Si superstructures, the binding energy between two Si adatoms on the three W surfaces as a function of distance was calculated first.The configuration of a pair of Si adatoms (Si 2 -xNN, x = 1−6, ranges from the first to six nearest neighbor) and the corresponding binding energy between two Si adatoms are shown in figure 1.Note that except the second Si atom in Si 2 -1NN on W(111), the two Si atoms both are located at the preferred 4F, LB and FCC sites on the W(100), W(110) and W(111) surfaces, respectively.On the W(100) surface, the maximum binding occurred at a distance of √ 2a in the [110] direction (Si 2 -2NN) and the maximum anti-binding was at a distance of a in the [100] direction (Si 2 -1NN).On the W(110) surface, the binding energy at a distance of √ 2a in the [110] direction (Si 2 -3NN) is 0.01 eV larger than at √ 3a/2 in the [111] direction (Si 2 -1NN); and the maximum anti-binding on W(110) is at a in the [001] direction (Si 2 -2NN), whose binding energy is 0.41 eV lower than Si 2 -1NN.Thus, Si 2 -3NN and Si 2 -2NN on W(110) has the highest and lowest binding energy, respectively, which is consistent with the derived values according to experimental observation frequencies of Si adatom pairs on W(110) in [22].The binding energy of a pair of Si adatoms on the W(110) surface as a function of distance is slightly different than that of a pair of W adatoms on W(110) [44], in which W 2 -1NN is more stable than W 2 -3NN.For the W(111) surface, Si 2 -2NN along the [ ] direction weakly binds together with a binding energy of 0.01 eV and Si 2 -1NN along the direction strongly anti-binding with a binding energy less than −1 eV.It is interesting to see that the maximum binding for Si 2 on the three W surfaces all occurred at a distance of √ 2a in the <110> direction, although the bond strength of Si 2 on the W(100) surface is stronger than the W(110) surface, followed by the W(111) surface.

Si coverage effect.
Based on the binding energy of two Si adatoms, the preferred superstructures on the W (100), ( 110) and (111) surfaces as a function of Si coverage were determined by exploring many possible configurations, and the results are presented in figure 2. At low coverages (θ < 0.5 ML), Si atoms on W(100) prefer to occupy 4F sites at a distance of √ 2a in the <110> directions instead of occupying separated 4F sites, as shown in figure 2(a) for the coverage of 0.25 ML.At a coverage of 0.5 ML, the minimum-energy configuration of Si adatoms on W(100) has a square c(2 × 2) pattern.Additional Si atoms would reside at the unoccupied 4F sites, forming a square p(2 × 2) pattern at a coverage of 0.75 ML.It is expected that 1 ML Si atoms will occupy all 4F sites on W(100), as shown in figure 1S in supplemental material; surprisingly, this configuration is not the most energetically favorable.Some Si adatoms among the 1 ML coverage prefer to reside on a second adlayer, rather than occupy the 4F sites on the first adlayer, which is likely due to the strong anti-binding interaction between 1NN Si adatoms on W(100).Based on more than ten evaluated configurations at a Si coverage of 1 ML, the configuration, as shown in the bottom panel in figure 2(a), has a relatively low energy.The Si patterns on the W(100) surface for coverages larger than 1 ML are shown in figure 1S in supplemental material.
Low coverages of Si atoms on the W(110) surface tend to form chains lining up along the [110] direction with a bond distance of √ 2a.The preferred Si configuration at a coverage of 0.25 ML has a rectangular pattern of c(4 × 2), which is the same as the LEED pattern observed in experiments [22][23][24][25].However, the present DFT calculations show that the rhombus p(2 × 2) pattern for a coverage of 0.5 ML is more energetically favorable than either the rectangular c(2 × 2) or c(2 × 4) LEED pattern following adsorption of silane at 350 K and annealing to 700 K reported in [24].Note that all Si atoms in the p(2 × 2) pattern for the coverage of 0.5 ML align along the ] direction on the W(110) surface, as shown in figure 2(b).This is not surprising since the binding energy of the Si 2 -1NN configuration along the <111> direction is only slightly lower than Si 2 -3NN on W(110) (see figure 1).In addition, the preferred pattern predicted in our DFT calculations at a coverage of 0.75 ML is c(4 × 2), consisting of a Si trimer, which is slightly different than the experimental LEED c(4 × 2) pattern [24].For the 1 ML coverage, the p(2 × 2) pattern predicted by DFT calculations is similar to the LEED pattern reported in [24], although the four Si atoms in one pattern observed in the latter is more compact than the former.The difference between the DFT calculations and previous experimental observations may be due to the influence of hydrogen introduced by the silane exposure and/or finite temperature.
On the surface of W(111), the Si patterns with increasing Si coverage is different than on the W(100) and W(110) surfaces.At coverages less than 1 ML, all Si atoms on W(111) prefer to occupy the adjacent FCC adsorption sites, forming adatom clusters, as illustrated in figure 2(c).However, the binding among Si adatoms in these clusters is very weak; in fact, the pattern of Si on W(111) for θ < 1 ML is protean at temperatures higher than 0 K.At the coverage of 1 ML, all FCC sites are occupied by Si.As the coverage is increased from 1 ML to 2 ML, Si adatoms located on the second adlayer deviate from the HCP sites and approach the sites above the bond center between two adjacent top W atoms, forming dimers with the pre-existing Si at FCC sites on the first adlayer.At surface coverages larger than 2 ML, Si adatoms on the third layer occupy the TOP sites.The Si adatom patterns on W(111) for coverages larger than 1 ML can be found in figure 1S in the supplementary material.
The average adsorption and desorption energy as a function of Si coverage, corresponding to the configurations in figures 2 and 1S, are plotted in figure 3. It is seen that the average adsorption energy of Si on W(100) increases with coverage until 0.5 ML.For coverages larger than 0.5 ML, the average adsorption energy on W(100) decreases steeply with increasing coverage.On the W(110) surface, the average adsorption energy of Si increases and reaches the maximum at the coverage of 0.25 ML.Then, the average adsorption energy gradually decreases with additional Si adatoms.However, the maximum average adsorption energy of Si on W(111) occurs at the coverage of 1 ML.Before reaching the maximum value, the average adsorption energy of Si increases slightly with coverage, but then the average absorption energy deceases for coverages beyond 1 ML.Among the three surfaces, the average adsorption energy of Si on the W(110) surface is significantly lower than either the W(100) or W(111) surface.At low coverages, the average Si adsorption energy on W(100) is larger than on W(111); while at coverages larger than 0.75 ML, the W(111) surface has stronger adsorption ability to Si than the W(100) surface for a given coverage.The most stable superstructure of Si adatoms on W(100), ( 110) and (111) occurs at the coverage of 0.5 ML, 0.25 ML, and 1 ML, respectively, which we refer to as W(100)-0.5ML, W(110)-0.25 ML, and W(111)-1 ML Si, respectively.It is noticeable that the patterns at these three coverages consist of the strongest Si-Si pair bonds on each surface, as seen from figures 1 and 2. The present DFT calculations confirm the previous proposal that adatom-adatom interaction is responsible for the adlayer superstructure formation on W surfaces [22], and is independent of the surface orientation.However, the surface orientation dependence of the adsorption energy for the most stable superstructure is not consistent with the order of the maximum binding energy of a pair of Si adatoms on W surfaces.Specifically, the maximum adsorption energy per Si adatom in the most stable superstructure on W(110), 1.18 eV, is clearly lower than the value of 1.83 eV and 2.46 eV on W(111) and W(100), respectively, while the maximum binding energy of a pair of Si adatoms on W(110) is in between that of W(111) and W(100).
The projected electronic density of states (PDOS) (states/eV) of a Si adatom in the most stable superstructure on each of these three W surfaces was calculated to reveal the underlying effect of surface orientation on Si adsorption.As expected, the majority PDOS of the s state of the Si adatom (Si-s) is distributed at a low energy region, while most of the Si-p state PDOS is distributed near the Fermi level, as shown in figure 4. The position of the peaks in the DOS relative to the Fermi level and the height of DOS near the Fermi level determine the occupation of the states and the bonding nature between adatoms and surface atoms [46].Figure 4 indicates that the PDOS peak of the Si-s state on W(100) is located at a lower energy in comparison to the W(111) and W(110) surfaces.Moreover, the PDOS of the Si-p state near the Fermi level for W(100) and W(110) is the lowest and highest, respectively, as seen from figure 4(a).The interaction strength of Si adatoms with the three W surfaces is also indicated by the change in the PDOS distribution of W surface atoms resulting from Si adsorption [47].In general, a stronger interaction corresponds to a higher adsorption affinity of the surface to adatoms. Figure 4(b) demonstrates the change in the PDOS of the W(100) surface atom, especially near the Fermi level, is the most significant, followed by W(111) and then W(110) with and without Si adatoms.The PDOS distributions are indicative of the highest and the lowest adsorption ability of W(100) and W(110) for Si, respectively.In addition, a comparison of the PDOS of the Si-p state with the W-d state indicates a similar electronic character of the Si adatoms to W atoms near the Fermi level, i.e. some states appear near the Fermi level.
To assess the Si coverage effect on desorption of Si from surfaces, the average desorption energy as a function of coverage was calculated.As expected, the average desorption energy as a function of coverage is the same as the average adsorption energy as a function of coverage.According to the average desorption energy shown by the right scale in figure 3, it is inferred that Si adatoms on the W surfaces can be more easily desorbed off the W(110) surface than the other two surface orientations.The average desorption activation energy of Si from the three W surfaces is in the range of 5.4-7.0 eV/atom for coverages lower than 1.0 ML, and decreases with increasing coverage.Our DFT values are close to the previously reported values of 4.8-6.0eV/atom in [26], but larger than the 3.82 eV/atom for the coverage range from 0.1 to 0.3 ML in [27] and 3.90-4.77eV/atom on polycrystalline W surfaces [29].As indicated by the DFT calculations, the desorption energy dependence on the surface orientations and Si coverages is consistent with the average adsorption energy.The present DFT calculations were based on three single surface orientations of a single W crystal at 0 K, which is different than the experimental environment in [26,27,29].Therefore, the discrepancy in the desorption energy predicted by the DFT calculations at 0 K and experiments at higher than room temperature may be caused by the impact of surface orientations and other factors.
To explore the adsorption in a 'sequential' way [48], the sequential adsorption energy is evaluated using equation ( 2) and the results are presented in figure 5.It should be noted that the Si adatoms are not always added one by one, and thus, the sequential energy at a coverage is the average value  between the current and prior coverage.The trend of the sequential adsorption energy with coverage before the adsorption energy reaches the maximum value is similar to that of the average adsorption energy (see figure 3); however, it is dependent on the surface orientation following additional Si atom adsorption beyond the maximum value.Particularly, the sequential adsorption energy of Si on W(100) drops significantly as θ > 0.5 ML, and then generally decreases with some fluctuations.It is surprising that for the W(110) surface, the significant drop in the sequential adsorption energy occurs at the coverage beyond 1 ML instead of 0.25 ML.The big turn at the 1 ML coverage indicates the weak stability of the p(2 × 2) Si pattern on W(110) (see the bottom panel in figure 2(b)), which is consistent with the reconstruction of the pattern following additional Si adsorption.Similar to W(110), the major turn in the sequential energy of Si adsorbed on W(111) occurs at a surface coverage of 1 ML.For coverages in the range 1 < θ < 2.5, the sequential adsorption energy of Si on W(111) slightly changes with coverage.According to the sequential adsorption energies and the adatom configurations, the three W surfaces can be covered by approximately one monolayer Si, and additional Si atoms beyond 1 ML reside on higher adlayers with lower energies.Thus multiple Si adlayers are predicted to form on W surfaces at high Si coverages, particularly on W(100) and W(111), by epitaxial growth above the initial Si adlayer.Previous experiment results have demonstrated the localized growth around the W(100) planes for average Si coverages larger than 2 ML [27].In addition, the saturation Si coverage of one monolayer on W(110) predicted by the present DFT calculations is in good agreement with the experimental observation in [24].

Si diffusion on W surfaces
To reveal the formation process and stability of the superstructures of Si adatoms on W the migration energy of a Si adatom on W surfaces with different Si coverages (θ ⩽ 1 ML) has been calculated using the NEB method.Since the initial configuration of the superstructure is more stable than after Si jumping, the migration energy barrier during the first jump, referred to as the threshold migration energy barrier, is generally larger than the following jumps [44].Thus, only the migration energy for the first jump is presented in figure 6.The threshold migration energy barrier and the required temperature are provided in table 4, in which the temperature was estimated by equation (7) in [49] according to the work of Janotti and van de Walle [50].As shown in figure 6(a) and table 4, the energy barrier for Si jumping on W(100) increases with Si coverage and reaches the maximum value of 2.78 eV at a coverage of 0.5 ML, while it decreases to 0.72 eV at a coverage of 1.0 ML.Note that in the superstructure of 1 ML, the migration of the adatom on the first adlayer requires 0.2 eV more than the adatom on the 2nd adlayer.
In addition, the adatom on the 2nd adlayer prefers jumping around the vacant 4F site with an energy of 0.62 eV; while the Si jumping into the vacant 4F site requires 0.82 eV, as indicated by the dashed and dotted lines in the bottom panel in figure 2(a), respectively.A single Si adatom (θ = 0.03 ML) can diffuse on W(100) at temperature higher than 750 K, and then react with each other to form an adatom cluster.Once adatom clusters are formed, the diffusion of a Si adatom away from the cluster requires a higher energy than a single Si adatom it has to overcome the dissociation energy from the cluster.The Si diffusion in the c(2 × 2) superstructure at the coverage of 0.5 ML on W(100) is unlikely to happen at temperatures lower than 1077 K.With further increasing surface coverage, the Si adatom on W(100) becomes more mobile.
The DFT calculated diffusion activation energy of a single Si adatom on W(110) (θ = 0.02 ML) is 0.74 eV, which is in good agreement with the experimental vale of 0.70 ± 0.07 eV [23].This indicates single Si adatoms at low coverages can diffuse on W(110) at room temperature.The energy barrier for Si adatom jumping at a coverage of 0.25 ML reaches a maximum value of 1.16 eV.The energy required for Si jumping at coverages between 0.5 and 0.75 ML is slightly lower than a single Si adatom diffusion on W(110).The NEB calculation of the Si migration at a coverage of 1 ML is not presented because of many possible migrations in the weakly stable superstructure, but the energy barrier is roughly estimated to be less than 0.6 eV.Thus, the required temperature for Si diffusion on W(110) is estimated to be less than 450 K.This implies that the patterns at various coverages, as shown in figure 2(b), may vary due to the migration of Si adatoms, even for the most stable superstructure at a coverage of 0.25 ML.These DFT calculations are consistent with the experimental observation, which demonstrated that although chain-form atomic clusters were often formed when more than one adatom was present on the W(110) plane, the clusters were not stable during diffusion [23].
Interestingly, the energy barrier for Si migration on W(111) at coverages less than 1 ML is around 1.7 eV, nearly independent of the coverage, as shown in figure 6(c) and table 4.Although Si adatoms on W(111) prefer to cluster, the bond strength among these adatoms in the clusters is very weak, as indicated by the near zero binding energy of two Si adatoms on W(111) in figure 1(d).The saddle point on the migration pathway of Si is close to the HCP sites.As such, provided there is an unoccupied HCP site on the pathway, the migration energy of Si between two 1NN FCC sites only varies slightly with surface coverage for θ < 1 ML.The required temperature for Si adatom diffusion on the W(111) surface for surface coverages less than 1 ML is approximately 660 K. Once all available FCC sites are occupied by 1 ML Si, the Si adatom mobility decreases.The Si adatom at FCC sites prefers jumping to its 1NN HCP site with an energy of 1.96 eV, and then the Si can jump to nearby HCP sites by an exchange mechanism with an energy of 0.51 eV, as indicated by the dashed arrow in the bottom panel in figure 2(c).
The threshold migration energy barrier of Si adatoms as a function of coverage on each surface is similar to the stability (adsorption/desorption energy) of Si superstructures as a function of coverage, except for W(110) covered with 0.75 and 1 ML Si.The most stable superstructure of Si adatoms in W(100)-0.5ML, W(110)-0.25 ML and W(111)-1 ML Si has the highest threshold migration energy barrier for each surface, as expected.At surface coverages less than 0.75 ML, the mobility of Si adatoms on the three W surfaces follows this order: (110) > (111) > (100), which is consistent with the surface stability order.In other words, Si adatoms on a more stable surface have a higher mobility at relatively low Si coverages.Although only three low-index surfaces were evaluated in this work, the W (320), (310), ( 211) and (210) surfaces are also more stable than W(100) at temperatures lower than 1500 K [45].Therefore, it is anticipated that Si adatoms at relatively low coverages are mobile on W surfaces at the operational temperature of W PFCs in tokamaks.

Si diffusion below W surfaces 3.4.1. A single Si atom.
To explore the formation mechanism of W silicide as Si atoms are adsorbed on W surfaces, the Si diffusional penetration below W surfaces through the interstitial or exchange mechanism was investigated by NEB calculations.Figure 7 shows the energy map of a Si atom diffusion from the surface adsorption 4F site to the W(100) subsurface (θ = 0.03).The interstitial sites above the third W layers are not stable for Si following relaxation.For the interstitial diffusion mechanism, as indicated by the red lines in figures 7(a) and (b), the Si adatom at the 4F site first jumps to the LB site on the W(100) surface with an energy of 1.94 eV.Then, the Si atom at the LB site perpendicularly migrates to the OIS at the third W layer (O L3 ) with an energy of 7.24 eV.The following jump of the interstitial Si atom from O L3 to O L5 by passing through O L4 only requires an energy of 0.72 eV.The energy map shows the activation energy for Si diffusion penetration below W(100) with an interstitial mechanism is approximately 9.2 eV.Without a doubt, this diffusion cannot happen at the operational temperature of W PFCs in a fusion environment [51].However, the Si adatom at 4F sites on W(100) can replace a W surface atom by kicking out the surface atom along the <111> direction with an energy of 2.21 eV (4F → Ex L1 ).Meanwhile, this kicked-out surface atom becomes a W adatom, as shown in figure 7(c).This process is referred to as Si diffusion toward/below W surfaces with an exchange mechanism.The W adatom may move to the 1NN site with an energy of 1.61 eV since the W adatom at the 1NN site has a 0.79 eV lower energy than close to Si.Likewise, Si at the surface layer can move to the second layer (Ex L1 → Ex L2 ) with an energy of 2.77 eV.The movement of Si from the second to third layer (Ex L2 → Ex L3 ) requires an energy of 4.82 eV.Obviously, the required energy for Si diffusion into W sublayers increases with depth.The required temperature for the first and second jumps is estimated to be 856 and 1073 K, respectively, and thus, Si can diffuse into the W sublayers by an exchange mechanism during operation of a W divertor in a fusion tokamak reactor.This diffusion process results in introducing Si in solution below the W(100) surface and the segregation of W adatoms, as shown in figures 7(b)-(e).With additional Si diffusion into the W sub-surface region  [36].
by an exchange mechanism, W silicide is anticipated to form at the W surface.Our DFT calculations indicate the activation energy of Si absorption into W from the (100) facet is 2.77 eV, which is close to the experimental value of ∼3 eV reported in [26].
In addition to W(100), Si diffusion below W(110) has also been evaluated.The energy map for Si diffusion below the W(110) surface and the configurations for the exchange diffusion mechanism are shown in figure 8. Si at LB sites on the W(110) surface can replace the surface atom by kicking out a W atom along the <111> direction with an energy of 2.41 eV, as indicated by the arrow in figure 8(b).Similar to the W(100) surface, the substitutional Si at the W(110) surface layer may diffusively penetrate below the surface by an interstitial or exchange mechanism.For the interstitial mechanism, the Si atom moves slightly toward the subsurface and the W adatom moves to the top of the Si atom simultaneously, forming a [110] Si-W dumbbell nearly centered at the first layer, (Ex L1 → d[110] L1 ), which requires an energy of 2.03 eV.Then, the dumbbell Si atom moves to the second layer with an energy of 4.05 eV.With increasing depth, Si will diffuse deeper to OISs at, and below, the third layer.Note that the migration energy barrier from O L3 to O L3−4 approaches the bulk value of ∼0.2 eV.For the interstitial diffusion mechanism, Si dissolution into the W bulk has to overcome an energy barrier of 6.53 eV.However, the W adatom in the configuration of Ex L1 prefers jumping to the 1NN site with an energy of 0.92 eV (Ex L1 → W ad -1NN), as shown in figures 8(c) and (d).The subsequent jumping from the first layer to the second layer by the exchange mechanism (W ad -1NN → Ex L2 ) requires an energy of 4.08 eV.Thus, the energy map in figure 8(a) indicates that Si adatoms can move to the W surface layer with an exchange mechanism at temperatures higher than 934 K (2.41 eV); while, it cannot diffuse deeper at temperatures lower than 1581 K (4.08 eV).Note that for the exchange mechanism, the W atoms move outward from the surface as Si adatoms move inward below the W surface.Our results are consistent the experimental observations in [24], which reported epitaxial growth of stoichiometric WSi 2 layers on the W(110) surface at 1050 K.
Different than the W(100) and W(110) surfaces, Si adatoms on W(111) cannot diffuse below the W surface as an interstitial atom regardless of energy levels.Our simulations demonstrate that a Si atom initially placed at interstitial sites above the 8th W layer below W(111) always replaced a nearby W atom and produced a W adatom during the relaxation.A Si atom can occupy the OIS on the 8th layer, but its solution energy is much higher than the substitutional configuration with a W adatom, as shown in figure 9(a).Thus, only the exchange mechanism was evaluated for Si adatom diffusion penetration below the W(111) surface.According to the NEB calculations, the energy barrier for the Si to jump from an FCC site to the first layer (FCC → Ex L1 ) is 2.20 eV.The Si atom at the first layer can jump to the second layer, and then arrive at the third layer (Ex L1 → Ex L2 → Ex L3 ) gradually with an energy of 4.58 eV.Alternatively, the Si atom at L1 can directly jump to L3 with an energy of 5.94 eV.Both of these are unlikely to be observed at temperature lower than 1775 K.
On the three low-index W surfaces studied in this work, Si adatoms prefers diffusing into W subsurface by an exchange mechanism, in which a W lattice atom is kicked out by a Si atom along the <111> direction, forming a substitutional Si  atom and a W adatom simultaneously.The required energy for the first jump of Si below W(100) (2.21 eV) is similar to W(111) (2.20 eV), and both are slightly lower than W(110) (2.41 eV).According to the estimated migration temperatures for continued direct exchange mechanism of Si and W atoms, Si adatoms on the W(100) surface can reach the 2nd W layer, while they only can replace the surface layer atoms of W(110) and W(111) during the operation of W PCFs.

Si coverage effect.
To understand the effect of Si coverage on the Si adatom diffusion toward W surfaces, the migration energy of a Si atom from a surface adsorption site toward the W surface with different Si coverages (θ ⩽ 1) was calculated and the energy map is plotted in figure 10.The migration energy barrier and the estimated temperature for significant diffusion are provided in table 5. Note that only the energetics of the Si adatom jumping to the first W layer  by an exchange mechanism was calculated here.As shown in figure 10(a) and table 5, the energy barrier for the coverage of 0.25 ML on W(100) is similar to the coverage of 0.03.However, the Si diffusion energy required for a coverage of 0.5 ML is evidently larger than the other coverages.Beyond 0.5 ML, the energy barrier for the coverage of 0.75 ML is similar to the 1 ML coverage, and both are slightly lower than those at relatively low Si coverages on the W(100) surface.Figure 10(b) and table 5 demonstrate that the energy barrier for Si adatom moving toward the W(110) surface layer with a coverage of 0.25 ML is clearly higher than the other four evaluated coverages.Surprisingly, the energy barrier for the other four coverages depends only slightly on the Si coverage.For the W(111) surface, the energy barrier for Si jumping toward the surface is not sensitive to the Si coverage with θ < 1 ML, similar to the coverage dependence of threshold energy barriers for Si migration on this surface.Note that the Si adatom migration to the W(111) surface layer at θ < 1 ML is involved the Si adatoms at the edge of surface adatom clusters (see figure 2(c)).The energy barrier for the coverage of 1 ML is approximately 0.8 eV higher than for coverages less than 1 ML.As expected, the coverage dependence of Si migration toward the surface is related to the stability of the Si superstructures as a function of coverage, except for the W(110) at a Si coverage of 1 ML, similar to the coverage dependence of Si on the surfaces.However, the maximum energy barrier of Si moving toward the W(110) subsurface, which occurs at the coverage of 0.25 ML, is higher than the maximum values for W(111)-1 ML, closely followed by W(100)-0.5ML.
A comparison of figure 10 (table 5) with figure 6 (table 4) shows that although the coverage dependence for the energy map of Si migration toward W surfaces is different than Si migration on these W surfaces, the coverage dependence for the migration energy barrier is similar for both cases near each surface.Si adatoms diffuse on the W surfaces rather than move toward the surface.Particularly, the energy barriers of Si adatom diffusion on the W(110) surface is significantly lower than penetration below this surface.Although the coverage dependence of Si diffusion energetics is related to surface orientation, for the most stable Si superstructure on each surface, the surface orientation dependence of the Si diffusion penetration below the surface is contrary to the Si diffusion on the surface, i.e. when the surface is more stable, Si diffusion on the surface is enhanced with a lower activation energy, while Si diffusion below the surface is reduced due to a larger activation energy.
The present DFT calculations indicate W silicide can form through an atomic exchange mechanism of Si migration toward W surfaces.The coverage and surface orientation dependence of Si moving to the first layer of the three surfaces implies that the W silicide formation condition is sensitive to the surface orientation and surface coverage, which will be studied in future work.W silicide formation on W surfaces may affect the advantage of a W divertor in fusion reactors, such as a low hydrogen retention.According to the present DFT and previous experiment results [52,53], the temperature and surface orientation are two key parameters for W silicide formation on W surfaces.To avoid or diminish the formation of W silicide on W divertors during siliconization, the temperature should be lower than 830 K and the exposed surface of W divertors should not be the (100) surface.

Conclusions
The solution and diffusion of Si in bulk W, as well as the adsorption and diffusion of Si near the W(100), W(110) and W(111) surfaces as a function of Si coverage have been simulated by first-principles DFT using VASP.In bulk W, Si prefers to reside at OISs and diffuse with an energy of 0.2 eV.Si atoms can be strongly trapped by W vacancies.Single Si adatoms prefer to occupy 4F, LB and FCC sites on the W (100), (110) and (111) surfaces, respectively.A pair of Si adatoms tend to combine at a distance of √ 2a along the <110> direction on each of these three W surfaces, which establishes the preferred superstructure of Si adatoms on the surface.The most stable Si adatom superstructure on these three W surfaces is a square c(2 × 2) pattern on W(100)-0.5ML Si, rectangular c(4 × 2) pattern on W(110)-0.25 ML Si, and rhombus p(1 × 1) pattern on W(111)-1ML Si, respectively.
The three evaluated W surfaces can absorb approximately one monolayer Si.In general, the additional Si adatoms beyond one monolayer reside on higher adatom layers with relatively low energies.The desorption activation energy of Si from the three W surfaces is in the range of 5.4-7.0 eV/atom for coverages less than 1.0 ML.The adsorption or desorption energy of Si adatoms on W(110) is significantly lower than either W(100) or W(111).At surface coverages less than 1 ML, the mobility of Si adatoms on W(110) is significantly higher than W(111) and W(100) at a given coverage.
Si adatoms diffusion into the W subsurface is energetically favored with an exchange mechanism, in which a W lattice atom is kicked out by a Si atom along the <111> direction, forming a substitutional Si atom and a W adatom simultaneously.The exchange diffusion mechanism indicates the possibility of epitaxial growth of W silicide layers on W surfaces during the operation of W PCFs in fusion reactors.
The coverage dependence of Si migration on/toward W surfaces is generally related to the stability of the Si superstructures as a function of coverage on each surface.For the most stable Si superstructure on each surface, the mobility of Si on W(110) is higher than W(111), then followed by W(100); however, the Si mobility toward the W(100) surface is similar to W(111), and both are higher than W(110).

Figure 1 .
Figure 1.Configuration of two Si adatoms on surfaces of (a) W(100), (b) W(110) and (c) W(111), as well as (d) binding energy of two Si adatoms as a function of distance, where the green and black spheres represent W atoms and one Si atom, the number range from 1 to 6 represents the second Si atom at the xNN (x = 1−6) site, respectively.

Figure 2 .
Figure 2. Si surface pattern as a function of Si coverage on surfaces of (a) W(100), (b) W(110) and (c) W(111).Green spheres are W atoms; black and magenta sphere are Si adatoms on the first and second layer, respectively.Red lines highlight the patterns.Blue arrows indicate the pathways of Si migration on W surfaces.

Figure 3 .
Figure 3. Average adsorption energy (left scale) and average desorption energy (right scale) of Si on W surfaces as a function of coverage.

Figure 4 .
Figure 4. Projected electronic density of states (PDOS) (states/eV) of a Si adatom (a) and a nearby W surface atom (b) in the most stable Si adatom superstructure on W surfaces.The top and bottom panel in (a) are the s state and p state of a Si adatom, and black, red and green lines represent PDOS of Si adatom on the W surface of (100), (110) and (111), respectively.The blue and orange lines in (b) represent the d state of a W surface atom with and without (clean) Si adatoms.The W surface atom on W(110) is the 2NN W to Si, and the others are the 1NN surface W to Si.The cyan lines indicate the Fermi level at 0.0 eV.

Figure 5 .
Figure 5. Sequential adsorption energy (eV) for additional Si atoms adsorbed on W surfaces as a function of Si coverage.

Figure 6 .
Figure 6.Migration energy of Si between two preferred adsorption sites on surfaces of (a) W(100), (b) W(110) and (c) W(111).The migration pathways are illustrated by blue arrows in figure 2.

Figure 7 .
Figure 7. (a) Migration energy of Si as a function of depth below W(100) at a coverage of 0.028, (b) configuration for Si adsorbed at a 4F site, configuration following the first (c), second (d) and third (e) jump with an exchange mechanism.In (a) dots and squares indicate the energy map for the Si migration with an interstitial and exchange mechanism, respectively.In (b)-(e), the cyan and green spheres represent the top surface and subsurface W atoms, black and red spheres represent Si and W adatoms.The red solid arrow in (b) indicates the first jump pathway for the interstitial mechanism.The blue arrows in (b)-(d) indicate the subsequent migration pathways by the exchange mechanism.The number 1 in (c) denotes the 1NN site for the W adatom.The depth of 0.0 nm in (a) is the W surface. O L3 , O L4 , and O L5 are the OISs right below the LB adsorption site and located at the third, fourth and fifth W layers, respectively, see figure 7(a) in [36].

Figure 8 .
Figure 8.(a) Migration energy of Si below W(110) at a coverage of 0.021, (b) configuration for Si adsorbed at a LB site, configuration following the first (c), second (d) and third (e) jump of Si with an exchange mechanism.In (a) blue dots indicate the energies for the first jump, red dots and squares represent the energy map for the Si migration with an interstitial and exchange mechanism, respectively.d[110] L1 and d[110] L2 represent the [110] Si-W dumbbell nearly centered at the first and second W layer, respectively.O L3 indicates OIS on the third W layer and O L3-4 indicates OIS between the third and fourth W layer, see figure 8(a) in [36].The representations of spheres and blue arrows in (b)-(e) are the same as those in figure 7.

Figure 9 .
Figure 9. (a) Energy of Si below W(111) at a coverage of 0.06 relative to Si adsorbed at FCC sites, (b) configuration for Si at a FCC site, and configuration following the first (c), second (d) and third (e) jump of Si with an exchange mechanism.In (a) black lines with squares indicate Si from a FCC site migration to the 3rd W layer by an exchange mechanism.The magenta lines with squares indicate substitutional Si migration from the 1st to the 3rd W layer directly.The blue squares denote the energy of a substitutional Si with a W adatom relative to Si adsorbed at FCC sites.The red dot indicates the relative energy of an interstitial Si atom at the 8th layer.In (b)-(e), the cyan and purple spheres represent the 1st and 2nd W layer, respectively.The green spheres represent W atoms below the 2nd layer.In (e), the Si atom is right below the W adatom and located at the 3rd W layer.The representations of the blue arrows in (b)-(d) are the same as those in figure 7.

Figure 10 .
Figure 10.Migration energy of Si adatom toward the top surface layer of (a) W(100), (b) W(110) and (c) W(111) by an exchange mechanism.

Table 1 .
Size of surface slabs used in each section, where nx and ny are the number of unit cell along the x and y directions on lateral planes, respectively.

Table 2 .
Energetics (eV) of Si in bulk bcc W, where E s OIS and E s TIS represent the solution energy of a Si atom at OIS and TIS, respectively; E s d<111> and E s d<110> represent the solution energy of a Si-W dumbbell in the [111] or [110]

Table 3 .
Adsorption energy (E ads ) and desorption energy (E des ) for a Si atom adsorbed at different sites on W surfaces.

Table 4 .
Threshold migration energy (E m /eV) and required temperature (T/K) of Si on W surfaces covered with different surface densities of Si.

Table 5 .
Migration energy barrier (E m /eV) and required temperature (T/K) of Si toward W surfaces covered with different surface densities of Si by an exchange mechanism.