Theoretical evaluation of the role of crystal defects on local equilibrium and effective diffusivity of hydrogen in iron

ABSTRACT Hydrogen diffusion and trapping in ferrite is evaluated by quantum mechanically informed kinetic Monte Carlo simulations in defective microstructures. We find that the lattice diffusivity is attenuated by two to four orders of magnitude due to the presence of dislocations. We also find that pipe diffusivity is vanishingly small along screw dislocations and demonstrate that dislocations do not provide fast diffusion pathways for hydrogen as is sometimes supposed. We make contact between our simulations and the predictions of Oriani's theory of ‘effective diffusivity’, and find that local equilibrium is maintained between lattice and trap sites. We also find that the predicted effective diffusivity is in agreement with our simulated results in cases where the distribution of traps is spatially homogeneous; in the trapping of hydrogen by dislocations where this condition is not met, the Oriani effective diffusivity is in agreement with the simulations to within a factor of two. This paper is part of a thematic issue on Hydrogen in Metallic Alloys


Introduction
Hydrogen can be dissolved in most metals and alloys [1], although it has remarkably low solubility in body-centred cubic (bcc) iron. The interaction of hydrogen with the microstructure is one reason for embrittlement and premature failure in high-strength steels and other engineering metals including nickel, titanium and zirconium alloys [2]. The deleterious effects of hydrogen are frequently connected to its rapid diffusion through the crystal lattice, in particular in bcc and martensitic steels. The interactions of hydrogen with crystal defects and their consequences are fundamentally less well understood than diffusion of hydrogen in the perfect crystal lattice, despite generally dominating the influence of hydrogen in metals.
Hydrogen atoms in the lattice are frequently divided into two categories: diffusing hydrogen, which can move freely through the normal interstitial sites, and trapped or non-diffusing hydrogen, residing around various crystal imperfections (vacancies, dislocations, grain boundaries and so on) [3][4][5]. The driving force for diffusion is the chemical potential gradient generated by hydrogen gradients combined with hydrostatic stress gradients. If fast trapping and detrapping kinetics are assumed, the freely diffusing and trapped hydrogen concentrations need to be in equilibrium. Based on the McLean isotherm, Oriani's [6] fractional occupancy θ t of a trap site is connected to the occupancy of normal interstitial sites θ L by where E B is the corresponding trap binding energy, k is the Boltzmann constant and T is the absolute temperature. If N L and N t represent the number of normal interstitial sites and number of traps per unit volume, then the volume concentrations of hydrogen in normal interstitial sites are c L = N L θ L and in traps c t = N t θ t . Using Equation (1), Oriani [6] finds an effective diffusivity of hydrogen given by (2) in which D L is the hydrogen diffusivity in the perfect bcc lattice. In this work the validity of Equations (1) and (2) is investigated by comparison to quantum mechanically informed kinetic Monte Carlo (kMC) simulations. Using the kMC, we may calculate the diffusivity of hydrogen in the perfect lattice, D kMC L . We can then create blocks of crystal containing defects which act as trap sites, and from the evolution of the simulation we can extract an effective diffusivity, D kMC eff . From the kMC simulations, the evolution of the occupancy of normal interstitial and trapping sites in pure bcc Fe is extracted and used as inputs to Equations (1) and (2) to calculate corresponding trap barriers and effective diffusivity, D Or. eff , in steady state, which are in turn compared for consistency with the input trap depths and output effective diffusivity, D kMC eff , to the kMC. Simulations of a perfect lattice are used to benchmark the resulting effective diffusivity in the lattice with defects (vacancy, 1 2 111 edge or 1 2 111 screw dislocation).

Diffusion model and parametrisation
Diffusion of hydrogen atoms through normal interstitial sites is governed by a random walk between occupied and empty sites. Under the framework of a kinetic Monte Carlo (kMC) method (cf. [7]), the rate probability for a jump from site A to B is given as where ν AB is the attempt frequency and E AB is the energy barrier between metastable position A through the saddle point to the position B, dependent on the local atomic configuration. For a given local configuration, a process k is selected according to where ρ 1 is a uniform random number between 0 and 1, and n is the number of all possible events, i, during one transition. The kMC simulates Poisson processes and the simulation time t kMC of each event is evaluated as [8][9][10] where ρ 2 is a second random variable in (0, 1], and directly corresponds to physical time t when only interstitial diffusion is considered. The time-dependent diffusion coefficient is given by the Einstein expression, [11,12] as where d is the system dimensionality (d = 1,2,3) and |r(t) − r(0)| 2 is the average squared displacement of particles in time t, [13]. In this study, all simulations were done at a constant temperature of 300 K, using a box with 20 3 bcc Fe unit cells, unless stated otherwise. Simulations were performed with full periodic boundary conditions in all directions. Fe atoms were not moved throughout the simulations and the permitted events were jumps between tetrahedral (T) sites, from T sites to trapping sites, vice versa and between trapping sites. The parameters used in the kMC model were for a vacancy obtained by using the quantum mechanical tight-binding approximation and molecular dynamics methods based around the Feynman path integral formulation of the quantum partition function [14,15]; or taken from the literature for the hydrogen interaction with 1 2 111 edge and screw [16] dislocations on {110} slip planes, also determined within the framework of path integral molecular dynamics at 300 K. The use of 20 3 bcc Fe unit cells results in very large defect concentrations in the simulation box; therefore, hydrogen concentrations were increased to obtain a realistic ratio between defect and interstitial concentrations. Simulations were performed with 5-30 hydrogen atoms and a vacancy, edge or screw dislocation, and a combination of edge and screw dislocations in the simulation box. The corresponding concentrations of interstitial hydrogen in atomic and weight parts per million are given in Table 1. To demonstrate the effect of simulation box size, test runs were also performed using 25 3 , 30 3 and 50 3 bcc Fe unit cells and compared with results obtained with 20 3 bcc Fe unit cells.

Parameters for diffusion around vacancy
For simulations on the movement of interstitial hydrogen, the energy landscape around a defect is needed. Energies are then used as inputs in Equation (3) to calculate the rate of all possible events. In Table 2, energies and calculated rate constants needed for interstitial movement around a vacancy at 300 K are given. Several events are considered: TT is a jump between adjacent tetrahedral sites, TTV1 is a jump from tetrahedral site to an unoccupied vacancy trap, TTV2 is a jump from a tetrahedral site to the vacancy trap if there is  already a hydrogen atom occupying one trapping site around the vacancy, TVV is a jump between adjacent trap positions around the vacancy and TVT is used for detrapping from the vacancy. According to quantum calculations up to six hydrogen atoms can be trapped around a vacancy; however, only two have high probability of being found there [14,15,17]. Therefore, in the model only up to two hydrogen atoms are allowed to be trapped around the vacancy to increase the computational efficiency of the model without sacrificing significant reality. Figure 1 shows a schematic representation of the events we consider in the kMC simulations of the vacancy microstructures. In all the calculations performed at 300 K, we use an activation barrier for lattice diffusion of 0.045 eV, calculated using zero-point energy-corrected density functional theory [18]; and an Arhenius prefactor of 5.12 × 10 12 Hz.

Parameters for diffusion around edge and screw dislocation
The energetic landscape of an interstitial hydrogen around edge or screw dislocations has many more local minima in comparison to the vacancy. As a remote hydrogen atom approaches the core, it encounters several metastable binding sites distributed between two parallel adjacent crystal planes which are perpendicular to the dislocation line sense [16]. In particular, the cores of edge and screw dislocations are reached through seven and six deep metastable binding sites, respectively. Movements from one metastable binding site to another and energy barriers for that transition were determined in [16] at 300 K using path integral molecular dynamics. However, these authors did not address transitions from bulk tetrahedral to metastable binding sites. The saddle point energy between a tetrahedral interstitial (T) site and a metastable binding site is assumed here to be similar to values obtained for transitions from a tetrahedral interstitial site to metastable binding site next to a vacancy [19]. In the case of edge dislocations, the present kMC model does not account for the change of binding energies and energy barriers between the metastable sites in compressive and tensile areas on both sides of the planar core. Instead, we assume that each metastable binding site is surrounded by four tetrahedral sites with energies for transition between T site to metastable binding points given in Tables 3 and 4 for 1 2 111 edge and 1 2 111 screw dislocation, respectively. TB and TZ represent transitions from four nearest T sites to metastable binding sites close to the dislocation, while BT and ZT represent the reverse. Transitions are grouped into TB/BT and TZ/ZT based on the distance between the T sites and metastable binding sites close to the dislocation (cf. Figure 2). A larger number labelling of a binding site corresponds to a position closer to the dislocation core, Table 3. Energies for transitions to and from the metastable binding site used in constructed kMC model for edge dislocation at 300 K.   where, for an edge dislocation, metastable point 7 is the closest to the dislocation line and, in the case of a screw dislocation, metastable point 6 is the closest to the core. As in the case of vacancy, at 300 K we use an activation barrier for lattice diffusion of 0.045 eV and an Arhenius prefactor of 5.12 × 10 12 Hz. Schematic representations of edge and screw dislocation energetic landscapes are shown in Figure 3. In the case of an edge dislocation, the core is depicted as an inverted T, while numbers represent metastable binding points for a hydrogen atom diffusing along the (110) slip plane of the edge dislocation (Figure 3(a)). Figure 3(b) shows schematically energetic landscape and metastable points around a screw dislocation. The dotted lines serve as a guide to view threefold symmetry of the (110) slip planes, and solid squares indicate the three deepest traps closest to the dislocation core. The arrows indicate those transitions of hydrogen atoms along the slip plane in the direction of the core that are included in our kMC simulations.
The activation energies for a proton to jump from a metastable binding site close to the core to another metastable binding site close to the core in the direction of the dislocation line are given in Table 5.

Diffusivity in ideal lattice
To begin with, we show in Figure 4 the calculated lattice diffusivity, D kMC L , in a perfect lattice at 300 K as a function of simulation time. (In this and subsequent figures we use 'diffusivity' loosely to mean the simulated diffusivity as a function of simulation time. But our calculated diffusivities are of course the steady-state limit of this quantity.) It is vital to note that simulation times of at least one microsecond are required before a convergent result is reached. This is an important observation since it emphasises the need to go beyond timescales accessible to molecular dynamics. We have used an activation barrier for lattice diffusion of 0.045 eV and an Arhenius prefactor of 5.12 × 10 12 Hz at 300 K, which reproduces the diffusivity measured by Nagano et al. [20], namely 8.98 × 10 −9 m 2 s −1 . As shown, there is negligible effect of hydrogen concentration on the interactions between interstitial atoms despite the small simulation box, as our results are in the range from 8.978 × 10 −9 m 2 s −1 to 8.981 × 10 −9 m 2 s −1 .

Diffusivity around point defect
In order to determine appropriate simulation box size in the presence of a trap, we performed simulations with different box sizes of 20 3    5 hydrogen atoms in the simulation box. Our results depicted in Figure 5(a) show that box size only has an effect on diffusivity at the outset of the simulation. Since very similar final results for diffusivity was obtained for all tested box sizes (cf. Figure 5(b)), the decision was made to use 20 3 bcc Fe unit cells (5.74 nm with 16000 Fe atoms and 100800 tetrahedral sites) to keep the simulation runtime achievable and reasonable.
The effects of different initial hydrogen positions were evaluated. It was found that different initial positions of interstitial hydrogen have negligible effects on the diffusivity and time when the steady state is reached. Figure 6 shows evolution of effective diffusivity from kMC simulations for a single vacancy and different hydrogen concentrations. As shown, the time needed to reach steady state for a vacancy is about 1 × 10 −6 s, indicating that hydrogen diffusivity around point defects achieves steady state almost instantly. For point defects this is expected since the activation barrier for trapping is lower than that for lattice diffusion, while the barrier for escape from the trap is so large that escape is a very rare event on the timescale of our simulations.
A comparison between calculated diffusivity in the perfect lattice and calculated effective diffusivity with different vacancy or hydrogen concentrations determined from kMC simulations is shown in Figure 7.  It can be seen that similar ratios between defect and hydrogen concentrations lead to comparable effective diffusivities.

Diffusivity around dislocations
Simulations with dislocations were performed mainly using a box size of 20 3 Fe bcc unit cells and different hydrogen concentrations. Results of effective diffusivity for a single screw, a single edge and a combination of screw and edge dislocations are shown in Figure 8. Depending on the defect type and hydrogen concentration, the time needed to reach steady state is for both dislocation types, around 1 × 10 −4 s. Results shown in Figure 8 indicate that hydrogen diffusivity around traps, simulated with kMC achieves steady state almost instantly. This discovery is important for mesoscopic and continuum simulations of crack propagation where the redistribution of hydrogen during crack growth plays an important role.
During the simulations, we also followed diffusion along dislocation cores (pipe diffusion) and results show that pipe diffusivity is several orders of magnitude lower compared to the effective lattice diffusivity. The resulting diffusivity in the case when hydrogen atoms are trapped inside the most stable points close to the dislocation core and move along the [111] direction of the dislocation line, clearly indicate that the diffusivity resulting from hydrogen transport along the dislocation line is several orders of magnitude lower compared to the effective lattice diffusivity. It follows that dislocation lines do not act as fast pathways for hydrogen pipe diffusion in bcc iron. In Table 6 are given results obtained for various hydrogen concentrations in a simulation box with a single edge or screw dislocations, where D kMC eff is the effective diffusivity and D kMC pipe is the diffusivity along the [111] direction within the dislocation core. Figure 9 shows the evolution of the diffusivity along   a screw dislocation core, showing that the diffusivity along the screw dislocation line is more than order of magnitude smaller than the effective diffusivity of the dislocated lattice. Pipe diffusion along the screw dislocation [111] direction is several orders of magnitude smaller than the lattice diffusivity.
From this we can conclude that the screw dislocation core is not a pathway for fast diffusion of hydrogen.
The origin of this is twofold: first the activation barrier is large, Table 5, and second the trap occupancy is relatively large leading to blocking of the diffusion path.
In Figure 10 a comparison between calculated effective diffusivity for several cases with different dislocation or hydrogen concentrations is shown. It can be seen that similar ratios between defects and hydrogen concentration lead to comparable effective diffusivities.

Comparison between kMC and Oriani's trapping theory
To test the validity of Oriani's assumption (cf. Equation (1), [6]) of equilibrium between freely diffusing and trapped hydrogen concentrations, the kMC model was developed in a way that trap and lattice occupancies could be followed. The total time hydrogen atoms spent at each site can be calculated at any given time during the simulation as the time spent in the trap or the lattice. Fractional occupancy of sites is obtained by dividing the time hydrogen atoms spent in a site by the total time. Since a kMC simulation aims to reach thermodynamic steady state, values at the end of simulations were used to determine trap binding energies of defects from resulting data for occupancies between lattice hydrogen and trapped hydrogen using equation (1) and to test the validity of Oriani's theory of local equilibrium by comparing these binding energies with those used as input. After suitable averaging of the trap occupancies, the diffusivity obtained from kMC is compared to the effective diffusivity calculated with Equation (2) using occupancies from kMC simulations. Detailed results and calculated trap binding energies for different defect and hydrogen concentrations are presented in Tables 7-9. Results for the vacancy in Table 7 show a minimal effect of concentration in the cases where the number of hydrogen atoms is smaller than the number of trapping sites. The average value for vacancy trap binding energy calculated using equation (1) is determined to be 0.22 ± 0.03 eV. In the case of a single type of dislocation (cf. Table 8), the number of possible trapping sites far exceeds the number of hydrogen atoms in the box. The trap binding energy determined for an edge dislocation is 0.21 ± 0.01 eV and is independent of the number of dislocations in the box. The trap binding energy of the screw dislocation is slightly dependent on the number of dislocations in the box and is 0.28 ± 0.002 eV for a single screw dislocation in the box and 0.34 ± 0.01 eV in the case of simulations with two screw dislocations. This increase occurs due to overlap between both dislocations. Oriani's effective diffusivity, D Or. eff , calculated using equation (2) for vacancies using occupancies from kMC simulations is presented in Table 7. Comparison between D Or. eff and D kMC eff shown in Figure 11 reveals that the diffusivity obtained from kMC simulations is slightly lower than that obtained using equation (2). The largest differences are in the case of lowest hydrogen concentrations, and the least differences when the hydrogen concentration is highest. This can be explained by the small difference between the number of hydrogen atoms that can freely move compared to possible trap sites. In the cases where the number of traps exceeds the number of hydrogen atoms, we get very similar results when comparing both approaches. From our results, we may show that Oriani's effective diffusivity is a very good Table 7. Results of kMC simulations for various vacancy and hydrogen concentrations. Note: Vac. is the number of vacancies, c H is hydrogen concentration, D kMC eff is the effective diffusivity, θ L is the lattice occupancy, θ T is the trap occupancy and E B is the trap binding energy determined from assumption of local equilibrium using Equation (1) and D Or. eff is the effective diffusivity calculated from Equation (2).  approximation for the description of effective diffusivity of hydrogen in the case of microstructures with a single type of trap uniformly distributed, for example, the case of vacancies. Results from kMC simulations with edge and screw dislocations in the box are presented in Table 9. The trap binding energy for combination of both types of dislocations in kMC simulation is 0.22 ± 0.003 eV and is close to the simulations where only an edge dislocation is present. Comparison of effective diffusivity shows that though vacancies can bind hydrogen to their nearest tetrahedral sites, the number of these sites is too small to have any large effect unless the number of possible trapping sites exceeds the number of hydrogen atoms. In the case of dislocations effective diffusivity is significantly lower than that in the ideal lattice or even in a lattice with only vacancies. Moreover, results show that the effective diffusivity in the presence of screw dislocations is a couple of decades smaller compared to simulations with edge dislocations or combinations of defects. This is due to deeper trap sites around screw dislocations and interactions between them if several screw dislocations are present in the box. Resulting values for the binding energies for vacancy and edge or screw dislocations are in good agreement with published results from the literature [21,22]. Figure 12 shows comparison between calculated diffusivity from kMC simulations and determined effective diffusivity using Oriani's theory for dislocations. In the case of a screw dislocation, Oriani's effective diffusivity calculated with inputs for occupancies at the end of kMC simulations is slightly higher compared to diffusivity  from kMC simulation (cf. Figure 12(a)). In the case of an edge dislocation and combination of edge and screw dislocations depicted in Figure 12(b), diffusivity obtained from kMC simulations is higher than calculated using equation (2). Larger discrepancies between diffusivities obtained from kMC simulations and calculated Oriani effective diffusivity, D Or. eff , using data for occupancies from kMC simulations can be explained by the spatial inhomogeneity of the energy landscape due to trap sites.

Conclusions
A kMC model for diffusion of hydrogen has been developed. This is able to simulate diffusion in ideal lattices and with vacancies and edge or screw dislocations with 1 2 111 Burgers vector. The following conclusions can be reached from the present work.
• Diffusivity in the ideal lattice is comparable to experimental results. • Direct comparison can be made of explicit simulations and the theory of Oriani which assumes local equilibrium and one type of trap. • The kMC confirms the establishment of local equilibrium in all the cases we have looked at.
• The Oriani effective diffusivity, D Or. eff , forms a good approximation to that calculated by kMC, D kMC eff . The best agreement is in the case of a homogeneous distribution of vacancies. For the case of dislocations, because these are spatially inhomogeneously distributed and there is a range of trap types with different binding energies, the Oriani effective diffusivity may differ by a factor of two from the kMC value.
As a result of kMC simulations of hydrogen transport in the proximity of a screw dislocation, we can conclude as follows: • Hydrogen atoms are confined to the region around the dislocation line, where they jump between the adjacent binding and low-energy metastable sites. • Hydrogen atoms occupy predominantly the deepest trap sites close to the dislocation line. • Hydrogen diffusivity significantly decreases with increasing dislocation density. • Dislocation pipe diffusion of hydrogen results in significantly lower diffusivity compared to lattice diffusion. We calculated the jump diffusion coefficient for hydrogen diffusion in the cases when hydrogen atoms are trapped inside the core and move in the [111] direction along the dislocation line. Our kMC simulations clearly indicate that the diffusivity resulting from hydrogen transport along the dislocation line is several orders of magnitude lower than that resulting from the lattice diffusion. Thus, the dislocation lines do not act as fast pathways for hydrogen diffusion in bcc Fe. The two reasons for this are the large activation barrier and the near saturation of trap sites leading to a blocking of the diffusion path.
Similar simulations of the hydrogen diffusivity in proximity of 1 2 111 edge dislocations show the following.
• The diffusivity near edge dislocations is higher compared to that in the vicinity of screw dislocations. Apart from the higher number of trap sites, diffusivity near a screw dislocation is also affected by its non-planar core structure. The planar core of the edge dislocation allows hydrogen atoms to move in directions parallel to the lines of trap sites, which increases their pipe diffusivity. We repeat that, in these preliminary simulations, we have neglected possible attraction and repulsion of hydrogen by the tensile and compressive strain fields above and below the slip plane of an edge dislocation. • There are no dominant trap sites. Hydrogen atoms are confined more or less evenly in all binding sites and jump between trap sites and the adjacent lowenergy metastable sites.
• There is no low-energy pathway for hydrogen diffusion in the [112] direction along the dislocation line. We used our kMC model for screw and edge dislocations to simulate diffusion of ten hydrogen atoms. Hydrogen diffusivity in this case is higher than the diffusivity around the same number of screw dislocations and lower than the diffusivity around edge dislocations. The reason for this is the higher hydrogen diffusivity in the proximity of edge dislocations.

Disclosure statement
No potential conflict of interest was reported by the author(s).

Funding
We