Anisotropic hydrogen diffusion in α-Zr and Zircaloy predicted by accelerated kinetic Monte Carlo simulations

This report presents an accelerated kinetic Monte Carlo (KMC) method to compute the diffusivity of hydrogen in hcp metals and alloys, considering both thermally activated hopping and quantum tunneling. The acceleration is achieved by replacing regular KMC jumps in trapping energy basins formed by neighboring tetrahedral interstitial sites, with analytical solutions for basin exiting time and probability. Parameterized by density functional theory (DFT) calculations, the accelerated KMC method is shown to be capable of efficiently calculating hydrogen diffusivity in α-Zr and Zircaloy, without altering the kinetics of long-range diffusion. Above room temperature, hydrogen diffusion in α-Zr and Zircaloy is dominated by thermal hopping, with negligible contribution from quantum tunneling. The diffusivity predicted by this DFT + KMC approach agrees well with that from previous independent experiments and theories, without using any data fitting. The diffusivity along  is found to be slightly higher than that along , with the anisotropy saturated at about 1.20 at high temperatures, resolving contradictory results in previous experiments. Demonstrated using hydrogen diffusion in α-Zr, the same method can be extended for on-lattice diffusion in hcp metals, or systems with similar trapping basins.

. (a) The hcp unit cell of Zr (gold) with tetrahedral (green) and octahedral (red) interstitial sites for hydrogen. (b) Schematic of an energy basin with two transition states (T 1 and T 2 ). p 1 (p 3 ) represents the probability for the system to transition from T 1 to T 2 (T 2 to T 1 ) in next event. p 2 (p 4 ) is the probability for the system to exit the basin from the left (right) side in next event. By definition p 1 + p 2 = p 3 + p 4 = 1. (c) 1NN interaction between possible impurity trapping sites (blue) with hydrogen interstitials at tetrahedral and octahedral sites.
hopping paths, namely first nearest neighbor (1NN) TT, TO, OT and OO jumps, are found responsible for the three-dimensional diffusion of H. The activation barriers and vibrational frequencies needed to calculate the hopping rates are predicted by DFT calculations, as listed in Table 1. The activation barriers calculated here agree well with those computed by Domain et al. 17 and Christensen et al. 18 . In addition to these 1NN jumps, second nearest neighbor (2NN) TT and OO jumps (denoted as TT2 and OO2, respectively) were also suggested in Domain et al. 17 . These 2NN jumps are however found to be unstable in this work. Our calculations show that TT2 and OO2 paths spontaneously relax into TO + OT and OT + TO, respectively. Therefore, in our KMC simulations, only 1NN TT, TO, OT and OO jumps are considered. The TT jump has a much lower barrier than others, forming a trapping energy basin in the PES ( Fig. 1(b)), which significantly reduces the efficiency of KMC modeling.
Although the absolute solution energy is not needed in KMC simulations, it is calculated in this work for comparisons with the literature. Using Eq. (21) in the Method Section, H solution energy is given as − 0.429 eV and − 0.366 eV at the tetrahedral and the octahedral site, respectively. Therefore, the tetrahedral site is more stable than the octahedral site for H, with an energetic preference of − 0.063 eV, which compares well with the − 0.059 eV reported in Domain et al. 17 and − 0.086 eV in Burr et al. 29 . For a more general comparison, most previous DFT calculations 17,18,29,30,33 have predicted that H prefers the tetrahedral site over the octahedral site based on the solution energy as defined in Eq. (21), which was reported to be − 0.60 eV in Domain et al. 17 and Lumley et al. 33 , − 0.46 eV in Burr et al. 29 , and − 0.52 eV in Udagawa et al. 30 . Experimentally, a value of − 0.66 eV was suggested 32 .
It has been pointed out that while evaluating the relative stability of various H interstitial configurations, the contribution of zero-point-energy (ZPE) in the solution energy should be included due to the small mass of H 18 . At 0 K, the ZPE of a H atom at interstitial site i is given by ∑ hv /2 Here h is Planck's constant, and v i,j the j th vibration frequency at site i. The summation runs over all three vibrational frequencies at energy minima and only the two real frequencies at TSs. Using the results listed in Table 1, the ZPE within the harmonic approximation is found to be 0.134 eV for the octahedral site, and 0.221 eV for the tetrahedral site, respectively. For the ZPE of H at the O site, our result agrees well with the value (0.123 eV) predicted by Christensen et al. 18 . With ZPE included, the octahedral site becomes more stable than the tetrahedral site with an energetic preference of 0.026 eV at 0 K. Such a small difference is actually within the uncertainty caused by the harmonic assumption in DFT calculations. It has been shown by Christensen 18 that by fully taking into account the anharmonic effect, the ZPE at the T site could be reduced from 0.18 eV to 0.166 eV if the neighboring Zr atoms were allowed to relax. A detailed consideration of this anharmonic effect is beyond the scope of this work. Nevertheless, the results indicate that T and O sites are both occupied by H atoms at finite temperatures due to their similar energetic levels.
The binding energies between H and alloying elements such as Sn, Fe, Cr, and Ni are calculated using Eq. (23). The results are listed in Table 2. The concentrations of the alloying elements used in the KMC simulations are also listed. Here, only the 1NN interactions are considered. The 2NN interactions are found to be much weaker than 1NN ones and are therefore neglected. For each type of impurity, three types of 1NN interactions with H exist, depending on the site of H and the relative positions of impurity atoms in reference to the H atom. As shown in Fig. 1(c), H at an octahedral site can be trapped symmetrically by six impurity sites within 1NN distance (number of trapping sites N t = 6) with a binding energy of E b O . Tetrahedral H interstitials can be trapped by one impurity site along < c> (N t = 1) with a binding energy of E b T c , and three others symmetrically (N t = 3) with E b T a . As shown in Table 2, Fe, Cr and Ni are attractive to H. The presence of these alloying elements will increase the overall H solubility and reduce H diffusivity by trapping H atoms locally. Sn is repulsive to H. In particular, H at the tetrahedral sites surrounding a Sn atom is not stable and automatically relaxes into a nearby interstitial site, meaning that these sites are blocked for H diffusion. Accordingly, the rate for H to jump to these sites is set to be zero in KMC simulations. Except for Ni, the trapping energies calculated here are much lower than those given in Christensen et al. 18 , which are also listed in Table 2. The comparison is made by comparing with the highest Path E m (eV) ν i,1 (THz) ν i,2 (THz) ν i,3 (THz) ν ts,1 (THz) ν ts,1 (THz) ν ± (THz)  (16) and (17) are exact, the acceleration is expected to result in no change in the H diffusion kinetics. To confirm this, KMC simulations are carried out every 100 K from 300 K to 1100 K. In each simulation, 100,000 basin exiting events are used to obtain good statistics. As shown in Fig. 2(a,b), for both exiting time and probability, the results from accelerated KMC perfectly match the analytical solutions and the results from regular KMC as well, indicating the correctness of the analytical solutions and the preservation of diffusion kinetics with acceleration. At 300 K and 400 K, regular KMC is too inefficient to get sufficient sampling of KMC events. Still, the results from accelerated KMC follow the analytical solutions. It is also interesting to note that at low temperature, the hopping rate ω 1 ≫ ω 2 (thus p 1 ≫ p 2 ), so that P L ≅ P R as given by Eq. (17). This implies that after a long trapping time in a basin, the memory of the basin-entering site is lost, leading to nearly equal exiting probability to the left (from the path H enters basin) and right existing path. However, at high temperatures, the trapping time is short and the escaping probability is biased favoring the basin-entering site. Capturing the escaping probabilities is critical for accurate H diffusivity in hcp crystals. Since exiting to the left results in zero net displacement along < c> , overestimate (underestimate) of its probability effectively suppresses (accelerates) H diffusion along < c> , leading to error in both overall diffusivity and diffusion anisotropy.
The average basin exiting time depends on temperature, as does the rate of acceleration relative to regular KMC. At low temperatures, p 1 approaches unity and the acceleration rate n acc , as defined in Eq. (19) in the Method Section, approaches 1/p 2 . It should be noted that the acceleration in Eq. (19) applies only to basin exiting. The computing time for jumps outside the basin remains constant. Therefore, the overall acceleration depends on how often the system visits basins and how long it takes to exit a basin by average. For H in α -Zr, the tetrahedral site is frequently visited since its population is twice of that of the octahedral site, with similar energies for both sites. As shown in Fig. 2(c), the acceleration is significant at low temperatures and decreases upon increasing temperature, mainly because of decreasing basin exiting time as given in Eq. (16). At room temperature, over 2,000 times overall acceleration can be obtained. Even at high temperature such as 800 K, the KMC simulations can still be accelerated by a rate of about 10. Since DHC is usually active in the temperature range of about 150 to 300 °C in Zircaloy claddings 27 , acceleration is important for efficient calculations of H diffusivity at low temperatures.
Hydrogen diffusion in pure Zr. The three-dimensional H diffusivity in α -Zr calculated by accelerated KMC is plotted in Fig. 3. The calculations are performed from 300 K to 1100 K, one data point every 100 K, covering the temperature range used in previous experiments. To minimize stochastic scattering, at each temperature 100,000 KMC simulations are carried out, each lasting until 100,000 basin exits are detected for sufficient statistics. The averaged mean-square-displacements (MSD) display the expected linear relationship with time ( Fig. 3(a)), with the slope proportional to diffusivity. At high temperatures (> 600 K), the KMC results are at the upper bound of the scatter in the experimental data, within a factor of 2 different from most experimental results [10][11][12][13][14] , as shown in Fig. 3(b). It should be emphasized that the modeling prediction is based on first principles and is completely independent of the experiments (i.e., no data fitting is used). Thus, this work demonstrates the predictive power of computer modeling. The minor discrepancy between KMC and experiments may come from uncertainty in DFT calculations due to underbinding of GGA or that in experiments due to the differences in measuring method and sample preparation. In particular, in Kearns et al. polycrystalline samples with various types of impurities were used 14,34 . Depending on their interactions with H, impurities and grain boundaries may As will be shown later, better agreement can be achieved by including impurity effects. At low temperatures (< 500 K), good agreement between KMC results and experiments 15 is also observed, with the former being higher but no more than a factor of 2 than the latter. It needs to be pointed out that in ref. 15, H diffusivity was measured separately along < a> and < c> , and the overall diffusivity plotted in Fig. 3(b) is converted using Eq. (13). As will be elaborated later, the results on D c reported by that work was unreasonably low compared to D a , possibly leading to a low overall diffusivity as well.
Interstitial diffusion in hcp metals involving octahedral and tetrahedral sites has been treated theoretically by different groups. Using an on-lattice random walk model, Ishioka and Koiwa 22 proposed an algorithm to derive the diffusivities of impurity atoms on a crystal lattice containing multiple sublattices, such as the octahedral and tetrahedral sites in hcp crystals. For H in hcp metals, the diffusivities along < c> and < a> directions are given by: TO  OO TT  TT OT   TT  TO  TO  OT   a  TO OT   TO  OT   2   2 where the hopping rate ω ij can be calculated using Eq. (5). More recently, starting from Fick's law and considering the balanced flux between equilibrated tetrahedral and octahedral sites, in Klyukin et al. an expression was derived for H diffusivity in hcp metal as 23 : Here Δ E TO is the difference between the solution energies at the tetrahedral and the octahedral sites, including ZPE. K B is the Boltzmann constant. Note that in Eq. (2) there is no differentiation of D c and D a . The parameters listed in Table 1 are used while applying the above theories. As shown in Fig. 3(b), our KMC results agree perfectly with the Ishioka model, which involves all relevant jumps. The Klyukin model is found to overestimate H diffusivity. It ignores the TT and OO jumps and essentially considers only H diffusion along < a> . The overestimation probably comes from the fact that it ignores the time spent for TT and OO jumps. The above results are obtained without considering quantum tunneling 19,20 , which could be important at low temperatures. Given the relatively high barriers for OO and OT jumps, for H diffusion in α -Zr, tunneling is expected to be less important than in some other metals 16 . According to the SC-HTST theory, the contribution of tunneling becomes significant below a critical temperature, T c = (hv ± E ZP /K B )/(2πE ZP − hv ± ln2). Here ν ± is the imaginary vibrational frequency at the TS (Table 1), and E ZP the migration barrier with ZPE correction. For TT, TO, OT and OO jumps, the critical temperature is estimated as 72 K, 68 K, 68 K and 38 K, respectively, indicating that tunneling is not significant above room temperature. To better estimate the tunneling effect, Eqs (8) to (10) are applied to the Ishioka theory and in KMC from 300 K to 1100 K. The ratio of tunneling corrected diffusivity (D T ) over that without correction (D) is plotted in Fig. 4. Since the correction factor in Eq. (9) is always larger than 1, so is the ratio of D T /D. Also, D T /D decreases with increasing temperature, upon which the tunneling effect diminishes. Above room temperature, the tunneling correction increases H diffusivity by less than 10%, with a  Table 2 in unit of atomic percent. As shown in Fig. 5, with all four alloying elements present, the overall hydrogen diffusivity is reduced, with greater reduction at lower temperatures. Fitting of the KMC results using D = D 0 exp(− E m /K B T) gives 1.08 × 10 −6 m 2 /s and 0.46 eV for D 0 and E m , agreeing nicely with the averaged values of 7.0 × 10 −7 m 2 /s and 0.46 eV for α -Zr and Zircaloy from experiments 14 . The higher fitted migration barrier for Zircaloy in reference to that of pure Zr (0.41 eV) is mainly due to the reduction in diffusivity at low temperatures. In the high temperature region (≥ 600 K), the reduction as shown by the KMC results is less than 15% (mostly within 10%). Such a small reduction is within the experimental scatter, in good agreement with Kearns et al. that H diffusivities in Zr, Zircaloy2 and Zircaloy4 are hardly distinguishable from each other 14 . In the previous experiments, the pure Zr samples contained impurities with very dilute concentrations 34 , which may have brought the H diffusion lower and closer to that in Zircaloy. Again, the minor discrepancy between KMC and experiments for Zircaloy could be caused by errors in DFT calculations, or possible presence of other impurities and grain boundaries in experiments.
As shown in Table 2, the binding energies calculated here are quite different from that given in Christensen et al. 18 . For a comparison, KMC simulations were also performed with the binding energies taken from ref. 18. As no information about the trapping site was given in ref. 18, in KMC only the trapping of a H atom at a tetrahedral site by a nearby solute atom along < c> direction (with N t = 1) is considered, so that the least trapping (i.e., minimum reduction of H diffusivity) is expected. Still, significant reduction in H diffusivity is observed using the binding energies in ref. 18 (Fig. 5), in contrast to previous experiments that showed similar H diffusivity in α -Zr  and Zircaloy 14 . This suggests that the H-solute binding might be overestimated in Christensen et al. 18 , which is possibly due to the neglect of local magnetic moments of Fe and Cr atoms in DFT calculations.
Theoretically, the trapping of impurities can be estimated using the Oriani model 25,35 when the concentrations of the diffuser and impurities are not very high. Specifically, the diffusivity with impurities D im is given by: . D L is the diffusivity without impurities. The physical meanings of other symbols are the same as in Eq. (20) in the Method Section. As shown in Fig. 5, the results from KMC simulations agree well with the Oriani model, while applying which the impurity-free diffusivity D L is obtained using the Ishioka model and the binding energies in Table 2 are used. According to Eq. (3), alloying elements that are attractive to H, such as Fe, Cr and Ni, reduce H diffusivity. In contrast, solutes such as Sn that are repulsive to H slightly increase H diffusivity. The results from KMC on the separate effect of each element are also in good agreement with the Oriani model, which adopts the same assumptions as used in the mean-field KMC method.
Hydrogen diffusion anisotropy. Due to the hexagonal symmetry of α -Zr, H diffuses anisotropically along < c> and < a> directions, with the latter representing the isotropic 2D diffusion in the basal plane. Since multiple jumps are involved, the anisotropy is not readily evident by just examining the individual hopping rates. KMC simulations allow for the calculation of diffusivity along both < c> and < a> and thus the anisotropy. Moreover, the relative importance of the jump paths can be evaluated at given temperatures to fully elucidate the controversy as reported in previous experiments: D c /D a > 1 but not exceeding 2 at temperatures above 600 K in Kearns et al. 14 , and D c /D a < 0.1 below 500 K in Zhang et al. 15 . Even though these experiments were done in different temperature ranges, the data reported were sufficient to establish contradicting trends over a wide range of temperatures.
In KMC simulations, D c and D a can be calculated using Eq. (12) by decomposing the MSD into components along < c> and < a> . The results obtained from KMC simulations follow the Ishioka model well, as shown in Fig. 6(b). Some scatter in the KMC results is observed due to the stochastic nature of the KMC method. The D c /D a ratio is found to increase with temperature and saturates to about 1.2 at high temperatures ( Fig. 6(b)), again agreeing well with the Ishioka model, where the anisotropy is given by: The temperature dependence in D c /D a comes from the temperature dependent hopping rates. At very low temperatures, ω OO ≪ ω OT and ω TO ≪ ω TT , so that D c /D a approaches 3c 2 /8a 2 , which is 1.0 for hcp metals with the ideal c/a ratio and about 0.96 for α -Zr using the lattice constants from our DFT calculations. This means that at very low temperatures, D a > D c for H in α -Zr. Upon increasing temperature, OO jumps become more important and the first term in the parentheses of Eq. (4) increases. Since the migration barrier for TT jump is significantly lower than that for TO jump (see Table 1), the second term in the parenthesis of Eq. (4) remains close to 1 and only slightly decreases with increasing temperature. The overall trend is that D c /D a continuously increases with temperature. As a result, at high temperatures D c becomes higher than D a . The transition between D c /D a < 1 and D c /D a > 1 occurs at about 270 K. In the above analysis the anisotropic thermal expansion of α -Zr is ignored since its effect is orders of magnitude lower.
Both KMC and the Ishioka model agree with Kearns et al., which predicted D c /D a > 1 with a ratio no more than 2.0 at temperatures above 600 K 14 . The KMC results on D a agree perfectly with those measured by Zhang et al.  Fig. 6(a)), where pure single-crystal α -Zr samples were used 15 . However, substantial discrepancy on D c is noticed. In KMC simulations, D c is not distinguishable from D a from 300 K to 500 K, while in the experiments the former was reported to be over an order of magnitude lower than the latter 15 (see Fig. 6(a)). Such a low D c /D a ratio is very unlikely if the 3-D, on lattice random walk of H is not altered. H diffusion along < a> is via TO and OT jumps, both with a component along < c> and thus contributing to diffusion along < c> . TT and OO jumps contribute only to < c> diffusion but via two different ways: i) the net displacement induced by these two jumps, and ii) providing a path for H to jump from one < a> plane to another, i.e., bridging TO (OT) jumps in neighboring < a> planes, so that H can perform 3-D random walk. Without the bridging effect, H diffusion will be confined in < a> planes, and the < c> component of TO and OT jumps will cancel each other. During 3-D random walk, the geometry of TO and OT jumps sets a lower bound of D c /D a , determined by the ratio of l l / OT TO is the total displacement along < c> (< a> ) summed over all OT and TO jumps. The KMC analysis also proves our assumption of the 3-D, on-lattice random walk of H, because otherwise OT TO a OT TO 2 2 cannot converge to 0.13. Note that this ratio is obtained without including the contribution from TT and OO. Therefore, D c /D a should be no less than 0.26 as long as H performs 3-D random walk, suggesting that in Zhang et al., D c was likely be measured at a situation where H diffusion along < c> was suppressed, deviating the random walk behavior. Given the small barrier of TT jump, confined diffusion along < a> is very unlikely in α -Zr or Zircaloy in the dilute concentration regime. Actually, our KMC simulations show that regardless of temperature, TT jumps that bridge two neighboring < a> planes accounts for about 20% of all KMC jumps ( Fig. 7(b)), and provides an effective path bridging OT and TO jumps in neighboring < a> planes. Here, each TT jump represents a basin exit event involving a net TT displacement, not an actual TT move as in regular KMC. With increasing temperature, OO jump becomes more important in contributing to < c> diffusion, for its increasing fraction of jumps and long absolute length.
The unreasonably low D c /D a ratio measured in Zhang et al. 15 may be explained by the so-called blocking layer effect observed in experiments on H in hcp Mg 36 . In Zhang et al., H diffusivity was measured using the surface segregation approach 15 . It is possible that although the matrix H concentration was kept below the solubility limit, the local H concentration near surface could have exceeded that limit when surface segregation occurred. Consequently, precipitation of hydrides such as coherent hydride clusters 37 may have occurred. As the hydrides habit on basal planes with small thickness along < c> , they block H diffusion along < c> , reducing D c without affecting D a much, similar to the situation for H in hcp Mg 36 .

Discussion
In this section, some discussion is given regarding the comparison between the present results with previous experiments and theories, as well as the applicability of the accelerated KMC method in other material systems.
The diffusivity of H in α -Zr predicted by KMC is at the upper bound of the scatter in the experimental data. The minor discrepancy, which is within the experimental error limit, may come from uncertainties in DFT calculations or in previous experiments. While calculating H vibrational frequencies using DFT, it is assumed that the vibration of H can be decoupled from that of the α -Zr lattice due to the large difference in their masses. Specifically, Zr atoms are fixed in these calculations since their vibration frequencies are orders lower. Minor error might be induced by such a treatment. For instance, an error of about 0.30 THz (~10 cm −1 ) was identified for the vibrational frequencies of H (at the order of ~30 THz) on a Ni (111) surface by ignoring lattice atom vibration 38 . Similar anharmonic effect has been noticed for H in α -Zr by Christensen et al. 18 . We also note that in the diffusivity calculations (Figs 3 and 5), the effect of quantum tunneling is neglected because its contribution is negligible above room temperature. It should be pointed out that the previous experiments are not exempted from uncertainties either due to the differences in measurement method and sample preparation, as indicated by the scatter in results.
For Zircaloy, again the KMC results are higher but within a factor of 2 than those from experiments. Compared to that in α -Zr, the addition of alloying elements slightly reduces H diffusivity, with negligible effect at high temperatures. This is in line with previous experiments, which found similar H diffusivity in α -Zr, Zircaloy2, and Zircaloy4 14 . The mean-field KMC approach for impurity trapping has two major assumptions. First, Eq. (20) is used assuming that the impurity atom located in the proximity of an H atom modifies only the energy of the initial state without altering the TS. This is inaccurate for H hopping around impurity atoms. However, it still captures the effective migration barrier for an H atom to move away from an impurity atom. To confirm this, the diffusion barrier for an H atom to diffuse away from a Ni atom, from a 1NN tetrahedral site to a 2NN tetrahedral site along < c> , is calculated using DFT. The directly calculated barrier is 0.350 eV, very close to that estimated using Eq. (20), 0.341 eV (E b T c , 0.212 eV, plus E m for T-T jump, 0.129 eV). Second, in this work only interaction between H and impurity atoms is considered, with no interactions between H atoms and between impurity atoms. Such an assumption applies when the concentrations of H and impurities are low (i.e., dilute solution). For high H concentration or high impurity concentrations, the interactions between H atoms and between impurity atoms need to be included. The mean-field approximation in Eq. (20) needs to be replaced by more accurate calculations of rate parameters considering local atomic configurations. The same assumptions are shared by the Oriani model 35 , which essentially predicts the same results as those from the mean-field KMC.
By comparing KMC with previous theories, this work demonstrates that for hcp metals the Ishioka model is accurate for H diffusivity, and a correction using the Oriani model is sufficient for dilute alloys, provided that all hopping rates are available. However, KMC simulations will be more useful in cases involving spatial heterogeneities, e.g., stress fields induced by cracks and hydrides. The presence of stress fields alters the local solubility and hopping rates as well at each atomic site, making the mean-field approach not applicable.
While the current work focuses on α -Zr and Zircaloy, the same method can be directly applied for H diffusion in other hcp metals such as Mg and Ti and their alloys, where similar diffusion mechanisms and energy basin hold. In Mg, the barrier for T-T jump is about 0.1 eV, much lower than that for T-O (~0.23 eV) and O-O (~0.21 eV) jumps 23 . In Ti, the barrier of T-T jump (0.061 eV) is nearly an order lower than that of T-O (0.424 eV) and O-O (0.625 eV) jumps 2 . H diffusion is critical for H storage in Mg and its alloys 1 , and for hydride embrittlement in Ti-based alloys 2 . In a more general sense, the acceleration method applies to interstitial diffusion in hcp metals, or to problems with similar energy basins 39 as shown in Fig. 1(b).

Conclusion
In this work, an accelerated KMC method is developed for efficient calculations of H diffusivity in hcp metals, and demonstrated using H diffusion in α -Zr. Using the hopping rates predicted by DFT, the method accurately predicts H diffusivity in α -Zr and Zircaloy, providing reliable data that could be used for upper scale modeling 5,6 . The results from KMC agree very nicely with previous independent experiments at various temperature ranges [10][11][12][13][14][15] and the Ishioka theory 22 . The perfect agreement between KMC and Ishioka theory validates the correctness of the analytical equations in this theory. The microscopic diffusion mechanisms obtained from DFT and KMC provide a systematic understanding of H diffusion in α -Zr, which helps to resolve the controversy in previous experiments regarding H diffusion anisotropy. Above room temperature, H diffuses by thermal hopping involving 1NN OT, TT, TO and OO jumps. At low temperatures, an effective diffusion path is via OT-> TT-> TO moves, with the contribution of OO increasing with increasing temperature. The diffusion of H is anisotropic, with the D c /D a ratio increasing from below 1.0 at very low temperatures to about 1.20 at high temperatures. In addition to H diffusion in hcp metals, the accelerated KMC method may be applied in other material systems with similar energy basins.

Methods
Residence time lattice kinetic Monte Carlo. In this work we use lattice KMC to calculate the diffusivity of H in Zr. Due to its small atomic size, H stays in Zr as an interstitial, residing in either a tetrahedral or an octahedral site, as shown in Fig. 1(a). A H atom at a tetrahedral site, T 1 for instance, can jump to the neighboring T 2 site (TT jump) or one of the three neighboring O sites (TO). A H atom at an octahedral site, say O 2 , can jump to one of the two neighboring O sites (OO), or one of the six neighboring T sites (OT). All these jumps are involved in three-dimensional diffusion. According to the quantum corrected harmonic transition state theory 40 , the hopping rate from site i to j is given by: Here T is temperature; and h is the Planck constant. E m is the classical migration barrier associated with the i-> j path. Z TS and Z i are the partition functions for the transition state (TS) and the initial state. Their ratio accounts for the zero-point energy correction and it is given by: where f(x) = sinh (x)/x. Note that Eq. (6) takes into account the ZPE correction at low temperatures and also reproduces the classical hopping rate at high temperatures. In Eq. (6) v i,j (v TS,k ) is the j th (k th ) vibrational frequency of H at site i (TS). For each jump, the migration barriers and vibrating frequencies are calculated using DFT calculations at 0 K. The results are listed in Table 1. At each KMC step, the hopping rates are calculated using Eq. (5) for all possible moves. A random number is drawn to pick one jump from the list. The time is advanced following the residence time algorithm by: with R being another random number between 0 and 1. The use of a random number in Eq. (7) better mimics the stochastic nature of kinetic events. With sufficient sampling, the two expressions in Eq. (7) converge to each other. The above formulation describes diffusion via thermally activated hopping without considering quantum mechanical tunneling, which has been shown to have substantial contribution to H diffusion in some metals at low temperatures 16 . Quantitative prediction of tunneling requires time-consuming calculations of the PES of H interstitial for each jump. Alternatively, assuming that the PES is harmonic at TSs and energy minima, the effect of tunneling can be estimated by applying a semiclassical correction to the hopping rates following the SC-HTST 20,21 , given by: Here ω SC ij is the corrected hopping rate with quantum tunneling, accounted for by the coefficient Γ ij SC . At a given temperature T, Γ ij SC is given by: In Eq. (9), E ZP is the ZPE corrected migration barrier, by and v ± is the imaginary vibration frequency at the transition state, given as negative in Table 1. The upper limit of the integration is given by To obtain diffusivity, the mean-square-displacement (MSD) of the migrating H atoms at a time t is calculated by: where N is the number of total H atoms in a KMC simulation, and r i (t) and r i (0) are the atomic positions of the i th H atom at time t and time 0, respectively. When multiple H atoms are used in a simulation, they are treated as non-interacting with each other and allowed to overlap, to give the diffusivity in the dilute concentration regime. Decomposing the total MSD into that along < c> and < a> , the diffusivities along < c> , D c , and < a> , D a , can be calculated by: The different coefficients in the denominator for D c and D a are due to the fact that D a is two-dimensional in < a> plane (normal to < c> ), while D c is one-dimensional along < c> .
The three-dimensional diffusivity is given by: By plotting the MSD as a function of time, diffusivities at various temperatures can be obtained by linearly fitting the MSD curves. Usually, averaging over a large number of independent KMC runs is needed to minimize the stochastic effect.

Acceleration of KMC.
For H in hcp metals, the migration barrier of TT jump is usually much smaller than that of TO, OT and OO, forming an energy basin ( Fig. 1(b)). At low temperatures, in regular KMC simulations a large fraction of the KMC moves are the back and forth moves between 1NN tetrahedral sites, with zero contribution to the long-range diffusion. This drastically affects the simulation efficiency, assuming that each KMC step costs about the same CPU time. (This assumption holds well for the case here. In general it depends on the number of events at each KMC step and the complexity to calculate the rates of all events). In general, the energy basin problem can be overcome by solving the master equations of the absorbing Markov chain for the occupation probability of the system as a function of time for each transition state 28 . The solutions can be used to replace regular KMC events to improve the efficiency. Following the same method, the expected exiting time and exiting probability for energy basin containing two energy minima ( Fig. 1(b)) are derived as: Here P L is the probability for the system to exit to the left, i.e., the site from where it enters the basin, and P R that to the right. By definition P L + P R = 1 once the system exits the basin. t 1 (t 2 ) are the residence time at the T 1 (T 2 ) states before next event occurs, which is the time step in a KMC simulation as given in Eq. (7). p 1 (p 3 ) represents the probability for the system to transition from T 1 to T 2 (T 2 to T 1 ) in next event. And p 2 (p 4 ) is the probability for the system to exit the basin from the left (right) side in next event. By definition p 1 + p 2 = 1 and p 3 + p 4 = 1.
In pure hcp crystals, the two neighboring T sites (T 1 and T 2 in Fig. 1(b)) are symmetrical and so are the O sites (for each exit in Fig. 1(b) there are three symmetrical O sites) for the H atom to exit to. Therefore, we have p 1 = p 3 and p 2 = p 4 . Similarly, t 1 = t 2 = t 0 . Eqs (14) and (15) can thus be reduced to: with ω i being the hopping rate of the ith jump, as given by Eq. (5). As can be seen from Eq. (17), a H atom always exits with a higher probability from the tetrahedral site where it enters the basin, since P L > P R always holds. In cases of ω 1 ≫ ω 2 , e.g., at low temperatures, P R approaches P L , giving nearly equal probabilities to exit from both T sites. To assure that the diffusion kinetics is preserved during acceleration, it is critical to capture the probability of each exiting path in addition to the basin exiting time. It is also interesting to note from Eq. (16) that for this case the expected exiting time is independent of the barrier between the two minima in the basin. For acceleration, the KMC events in basins are replaced using the solutions in Eqs (16) and (17). In a KMC simulation with only one H atom, once it enters a basin, it will be taken out in next event, with the exiting path chosen according to Eq. (17). The time advancement is: where R is another random number. Statistically, a H atom in a basin will perform p 1/ 2 moves before exiting. Assuming that each KMC step consumes the same CPU time, the acceleration rate of basing exiting time can estimated by the below equation: The above approach is for KMC simulations with one H atom in each. In simulations with multiple H atoms, the same acceleration can be achieved by combining neighboring tetrahedral sites (T 1 and T 2 in Fig. 1(b) for example) into a super site (T * ), and modifying the hopping rates towards neighboring octahedral sites based on the path along which a H atom enters the combined T * site. We note that since the time step is inversely proportional to the number of H atoms used in the simulation, the use of multiple non-interacting H atoms results in no change in the diffusion kinetics and theoretically no change in the efficiency either.
Mean-field impurity trapping. To consider impurity trapping, the mean-field KMC approach 25,26 developed for cubic lattice is extended to hcp, where each impurity atom can trap H in three different ways with different binding energies (see Fig. 1(c)). At each KMC step, the migration barrier of each jump is modified by: Here E m 0 is the barrier without trapping. E t i is the binding energy of H at a trapping site t with a concentration c t i , induced by impurity i, whose concentration is c i . R is a random number to be drawn each time when the jump rate is evaluated. The concentration of trapping site t is given by = c N c t i t i , with N t being the number of trapping sites induced by an impurity atom.
DFT calculations. Ab initio calculations are performed using the all-electron projector augmented wave method within the generalized gradient approximation of Perdew, Burke, and Ernzerhof (PBE-GGA) 41 , as implemented in VASP 42 . Large 96-atom supercells, which can be constructed from a 4 × 4 × 3 extension of the 2-atom hcp Zr unit cell, are used throughout our calculations in order to minimize the unphysical interactions between a H atom and its periodic images. A high plane-wave cut-off energy of 500 eV and a dense 5 × 5 × 5 Monkhorst-Pack k-point mesh are used to ensure high numerical accuracy for total energy calculations. All internal atomic positions are fully optimized using a conjugate gradient method until forces are less than 0.01 eV/Å. Further relaxations of supercell volumes have been found to have negligible effect on the final results.
Following Domain et al. 17  The TSs for all diffusion paths of H are obtained using the climbing image nudged elastic band (CI-NEB) technique 43 with 3 intermediate images.
Here we only consider 1NN TT, TO, OT and OO jumps since the 2NN jumps are found to be unstable and they spontaneously relax into 1NN jumps. The normal-mode vibrational frequencies of H are obtained from the eigenvalues of the Hessian matrix constructed using finite differences with a small displacement of 0.05 Å. All metal atoms are rigidly constrained during such calculations.
The binding energy between a substitutional solute element X (X = Sn, Fe, Cr, Ni) and an interstitial H atom in hcp Zr can be calculated using the following equation: where E(Zr N−1 X 1 H 1 ) is the total energy of the supercell containing one impurity atom X and one H interstitial in close proximity of each other. E(Zr N ) and E(Zr N−1 X 1 ) are the total energy of the perfect Zr supercell with N Zr and the Zr supercell with a substitutional X atom, respectively. E(Zr N H 1 ) is the energy of a Zr supercell with one H occupying the same type of interstitial site as that in the supercell used for E(Zr N−1 X 1 H 1 ). As shown in Fig. 1(c), there exist three different 1NN interactions between H and an impurity within 1NN distance. All these interactions are considered in our calculations. For binding between H and Fe, Cr and Ni, spin-polarized calculations have been performed. According to our calculations, Fe and Cr develop large local magnetic moments (> 2 μ B ), while Ni is essentially non-magnetic. The final results are summarized in Table 2.