Hydrogen in α-iron: role of phonons in the diffusion of interstitials at high temperature

Ab-initio Density Functional Theory has been used to compute phonons for interstitial hydrogen in α-iron. In the harmonic approximation, these phonons yield Helmholtz’s free energy as a function of temperature, which can be used to obtain diffusion barriers from an Arrhenius plot. By comparing with the experimental database compiled by Kiuchi and McLellan, we show that the role of phonons is crucial to understand the diffusion of interstitial hydrogen at T > 300 K. The computed specific heat for Fe16H and Fe behaves quite differently due to the appearance of optical modes and could be used to calibrate the amount of interstitials in the iron matrix.

Civil engineering structures like buildings, bridges, dams, and others, rely on cores of ferritic steels made of body-centred-cubic (bcc) iron (α-Fe). To understand structural performances, it is of paramount importance to study the conditions that may degrade their properties, like mechanical strength, ductility, toughness, etc. Hydrogen embrittlement is arguably the most frequent cause for catastrophic failure in structures subjected to tensile loads and consequently, one of the most active research areas in metallurgy [1][2][3][4][5][6][7][8] . High strength ferritic steels are particularly vulnerable to hydrogen embrittlement, which despite having been first reported as early as in 1875 9 , is a complex phenomenon that has so far evaded a fully satisfactory explanation. There are currently several promising strategies trying to describe, understand and link the intervening mechanisms at several scales (multiscale approach) 3,10 . Likewise, many experimental studies have adopted the same approach and characterised these mechanisms from the atomic scale to the macroscopic one 6,[11][12][13] .
Among the relevant phenomena, hydrogen transport in the bcc iron lattice plays a prominent role in hydrogen embrittlement. This process has been addressed both from the theoretical and experimental viewpoints. First-principles calculations show the stability of the cubic lattice high symmetry points (tetrahedral and octahedral, respectively) to host interstitial hydrogen and the associated diffusion barriers [14][15][16][17][18] . There are many experimental results to determine the diffusion coefficient of hydrogen in α-Fe 19,20 . However, lattice defects trap hydrogen and affect these measurements, which can only be interpreted as an effective diffusion coefficient. Analysis of thermal desorption of hydrogen is the method of choice to characterise the influence of lattice defects [21][22][23] . Since this technique also depends on the value of the diffusion coefficient of H in the bcc iron lattice, a realistic measurement of the diffusion coefficient can be obtained through successive iterations. In this regard, a reasonable initial theoretical estimate of the defect-free diffusion coefficient over an extended temperature range is of great importance to decouple observed diffusion from the effect of H traps. Such an estimation will prove valuable for the interpretation of desorption experiments and to feed computational models for diffusion (such as standard or Kinetic Monte Carlo models).
The goal of this paper is to study from first-principles the contribution of phonons to the diffusion of hydrogen in α-Fe. Recent work by different groups has shown that ab-initio spin-polarised density functional theory (DFT) 24,25 based on pseudopotentials 26-28 is a reliable tool to describe the physical properties of iron, including structural parameters, magnetic properties, elastic constants, and different phases [14][15][16]18,29 . Phonons for α-Fe with interstitial H yield the vibrational internal free energy, allowing to extend previous DFT studies to the T ≠ 0 K domain and to probe the influence of temperature on the diffusion of hydrogen. As shown, phonon contributions are a crucial factor to explain the diffusion of interstitial hydrogen at high temperature.

Results
To study the effect of phonons on the diffusion of interstitial hydrogen we set up a model based on a 2 × 2 × 2 body-centred cubic unit cell (Fe16), which is large enough to minimise the interaction of the interstitial with its periodic-generated images, and at the same time, it is amenable for the task of computing phonons. According to the precision of our DFT calculations (c.f. Methods: Density Functional Theory) the optimum cubic cell has a = b = c = 5.668 Å, α = β = γ = 90°, symmetry group IM-3M (IT # 229), and a magnetic moment per atom of 2.18 μ B . Accepted experimental values for these magnitudes are 5.734 Å and 2.22 μ B 30 ; therefore, the fractional errors in these calculations are 1.2% and 1.8% respectively, which show the adequacy of the theoretical setup for our purposes. tetrahedral and octahedral high-symmetry sites. Optimization of forces and stresses for a hydrogen concentration Fe 16 H results in a global minimum where the interstitial is located at the tetrahedral high-symmetry site. We recall that for large concentrations of interstitials, the internal stresses originated by the lattice distortions needed to accommodate these interstitials alter the relative equilibrium between tetrahedral and octahedral sites, making the octahedral site the global minimum 16 . Therefore, it is clear that for the task of finding a global energy minimum, both forces and stresses have to be minimized in a consistent way. We find that the internal stress related to the tetrahedral interstitial is accommodated in the cubic unit cell by the following internal deformation: phonons. The main result of our calculations are phonon spectra for Fe 16 H with (i) the interstitial located in a tetrahedral site ( Fig. 1), or (ii) in an octahedral site (Fig. 2). Figure 3 gives the corresponding density of states, ρ(w), which allows us to obtain the thermodynamical magnitudes listed in the section of Methods (Thermophysics). When comparing with pristine bcc-Fe, the significant feature of these spectra is the appearance of high-frequency optical modes. Regarding the internal electronic energy, the tetrahedral site corresponds to a stable stationary point; therefore, three vibrational modes appear approximately between 900 and 1200 cm −1 . On the other hand, the octahedral site represents an unstable point (local maximum) with two imaginary modes, and one real with w ≈ 2045 cm −1 . These are related to Fe-H vibrations, as a gaussian calculation confirms for FeH: this simple model has a single stretching frequency at 1780 cm −1 for an equilibrium distance of 1.57 Å between H and Fe (to be compared with 1.47 Å, the distance between H in the octahedral site and the nearest Fe) 31 .
Physical insight on the role of these new optical modes as a function of T can be obtained by computing the specific heat at constant volume (C V (T), Equation (9). By comparing with pristine bcc-Fe, we see how the main effect of high-frequency optical modes is to shift C V towards higher values of temperature, e.g. compare black dashed (Fe 16 H) and black continuous (Fe) in Fig. 4. This result constitutes a precise prediction that can both be tested with experiments, and can also be used to calibrate the number of interstitials due to its sensitivity to concentration. While a Debye's model represents well pristine Fe with an acceptable value T D = 476 K 32 , we observe that Fe 16 H demands a higher value, T D = 900 K, which usually is interpreted as a hardening of the material. However, we notice that in this case the higher value for T D is more related to the appearance of optical modes than to the hardening of the acoustic ones; in practice, such a large value for T D shows the difficulty to populate high-frequency modes at low temperatures. This interpretation is supported by a calculation of the elastic constants for Fe 16 H (for reference, ab-initio values for Fe computed with the same formalism are given in parantheses): C 11 = 177(187), C 12 = 87(76), C 44 = 52(62) GPa. These values can be transformed into an averaged sound speed by using Debye's theory, which could be useful to test the theory with experiments. We get 2.7 km/s for Fe 16 H, to be compared with the value 2.9 km/s for pure α-Fe (an experimental value for Fe is 3.6 km/s 32 ). Therefore, we conclude that despite the need for a larger best-fit Debye temperature, a comparison with the elastic constants for the pristine material indicates that the presence of interstitials tends to make it softer.
For the process of diffusion, which is our primary concern, we turn now to the vibrational properties of Fe 16 H while H is near the transition state. The calculation of phonons properly converged right at the transition state was a challenging task due to the unstable nature of the transition state and its low symmetry. Therefore, we investigate the phonons at the nearby octahedral site, which is also unstable but belongs to a specific group symmetry that can be used to perform a constrained calculation. We have remarked on the importance of eliminating internal stresses to find equilibrium configurations. However, for the barriers in the diffusion process, we find that this is not a crucial point since the time spent by the interstitial near the barrier is not enough to allow a full relaxation of the whole lattice. Therefore, we compute energy configurations near the octahedral site using parameters for the periodic lattice which correspond to the global equilibrium configuration for the system, i.e., the tetrahedral site. In such a configuration the material is under the following internal stress: σ xx = σ yy = 1.3, σ zz = −2.7 GPa. As expected, for such a non-stable configuration, we find that a doubly degenerated optical mode becomes soft, while a third optical branch hardens considerably, cf. compare Figs 1 and 2 (top panels). Such a high-frequency optical mode is responsible for a significant increased zero-point energy correction between the tetrahedral and the octahedral sites, which becomes a relevant feature for the diffusion of the interstitial, as we discuss in the next section. Figure 5 shows the difference between Helmholtz's free energies for tetrahedral and octahedral sites related only to phonons, obtained by using Equation (4), which gives us access to the variation with temperature of diffusion barriers (the electronic contribution to the free energy, which is temperature independent, is considered separately). The value ΔF(t − o, T = 0) gives the difference between the zero-point energy correction for both sites, +70 meV, which reduces the difference in the internal electronic energy between these two sites from meV, which can be compared to the zero-point correction (T = 0 K) obtained from the difference in free energies: F t − F o = 70 meV. Therefore, both by the quantitative determination of the difference in the zero-point correction from only the optical phonons, or by the negligible contribution of the rest of the modes (ΔF < ±1 meV), it is convincing to assign the main effect in the relative stabilization/destabilization of tetrahedral vs. octahedral site to the high-frequency end of the spectrum Fig. 5 shows the evolution of this difference with T when all modes are included in Equation (4). We notice that for high temperatures the tetrahedral site is further stabilized compared to the octahedral one.

Discussion
We benchmark our model for the free energy by computing the diffusion coefficient for an interstitial H jumping where ΔF(T) > 0 is the total free energy barrier between the transition state and the absolute minimum (tetragonal site), a is the length of the jump, and f is a phenomenological factor representing the reduced mobility due to traps and other defects.
To evaluate Equation (1), and to produce an usable expression for the diffusion of interstitial hydrogen in α-iron, we take for a a typical value for the distance between tetragonal sites, a ≈ 1 Å. Regarding f, it depends in general on the number of traps and their trapping energy 19,23,34 ; here we use a constant value, = f 1 2 , which yields an overall best agreement in the range of temperatures of interest between experiments and DFT ab-initio computed values. ΔF(T) is obtained using Equation (4) for the tetragonal and octahedral sites. As commented above, the calculation of converged phonons at the transition state is quite involved; therefore, we approximate the phonons at the transition state by those corresponding to the neighbouring symmetrical octahedral site. Since phonons involve smaller energies than typical energies in chemical bonds, their calculation using DFT implies more stringent conditions on the convergence thresholds. In that regard, conditions derived from imposing a symmetry group facilitate the process of finding self-consistency and convergence of energies and forces, because symmetry helps to compensate small round errors for slightly different but equivalent geometrical configurations. The resulting variation with temperature of Helmholtz's free energy, Equation (4), is shown in Fig. 5. As the temperature increases, the effective barrier increases by about 15% because of the differences in the phonon spectra between the tetrahedral and octahedral sites. Therefore, we obtain the following parametrization to our full ab-initio calculations, valid over the temperature range   T 300 850 K (at Curie's temperature, T C = 1043 K, iron adopts a paramagnetic phase different from the ferromagnetic one assumed in these calculations; furthermore, at T * = 1185 K, iron undergoes a phase transition to a face-centered cubic structure: as these temperatures are approached, the assumptions underlying the model become less valid). where T is in Kelvin.  www.nature.com/scientificreports www.nature.com/scientificreports/ The behaviour of this expression is displayed as the continous line in the Arrhenius plot in Fig. 6. The dashed line represents an extensive compilation of experimental results, as selected by Kiuchi and MacLellan 20 ; these authors have carefully analysed an extensive experimental dataset to report a representative value for the diffusion barrier of ≈60 meV for T < 300 K and ≈70 meV for T > 300 K. In addition, for the sake of comparison, we include a similar theoretical result by Jiang and Carter, where the total electronic energy has been corrected by the zero-point energy, yielding a barrier independent of temperature of ≈40 meV (dotted line) 18 . Comparing with such a DFT approach, the effect of phonons introduces a correction on the diffusion coefficient of about 67% towards the high temperature region, improving the agreement with the experimental dataset because the phonon-corrected barrier goes approximately from 45 meV at T = 300 K to 60 meV at T = 600 K. To quantify the improvement in the theor y-experiment agreement, we compute a simple figure of merit, T thr e xp 9 1 2 , which behaves like a metric and goes to zero when the pair of functions converge uniformly to each other. In the interval of temperatures of interest, i.e. between 300 and 850 K, we find for the conventional DFT approach (dotted line) R = 12.2, and for the present approach that includes the phonons in the free energy (continuous line) R = 3.1. Therefore, these R values confirm the better tendency of theory towards high temperatures. We remark that we restrict ourselves to the region 300 < T < 850 K; for higher temperatures iron starts to approch a different phase, while for lower temperatures other effects should be taken into account, like tunnelling-aided diffusion. In our calculations, we do not interpret the change in the diffusion barrier between high and low temperature as a transition between two different regimes, as it is assumed in ref. 20 Instead, we associate the difference in the diffusion-barrier with a continuous variation due to the population of vibrational modes with increasing T. We notice that this is a result that we expect to be independent of the approximation we have made to substitute the phonons at the transition state by the ones at the local maximum at the octahedral site. From this result, we conclude that the contribution of phonons is almost entirely responsible for the former discrepancy between experiment and theory observed in the region of high-temperatures. The discrepancy between theory and experiment for low temperatures is likely related to H tunnelling assisted diffusion. For such a mechanism we refer to the path integral calculation by Kimizuka et al., that is better suited to describe the region T < 300 K 35 .

Conclusion
In conclusion, these results help to identify a physical mechanism that is relevant to describe the diffusion of interstitials; the contribution of phonons to the free energy. These calculations prove that the contribution to Helmholtz's free energy of phonons is enough to improve the description of interstitial diffusion in a range of temperatures between 300 and 850 K, which are relevant for industrial and civil engineering applications. Furthermore, these results suggest an experimental route to confirm the specified role of phonons via the computed C V (T), or the predicted softening of the speed of sound.

Methods
Density functional theory. We utilize spin-polarized ab-initio Density Functional Theory, as implemented in CASTEP (plane-waves basis) 36 37 . At a microscopic level, magnetism derives from spin, a quantum operator that can only be properly justified in a relativistic theory, i.e. using Dirac's Equation. Already for H, a simple estimate based in Heisenberg's Uncertainty Principle results in an effective velocity for the electron of about one-tenth of the speed of light. As the nucleus grows, the electrostatic interaction with the inner electron shells grows accordingly and makes relativistic corrections more important. This physical effect is well known, although usually it is ignored in favour of faster standard non-relativistic methods. However, Koelling and Hamon devised a technique easy to implement in standard band-structure calculations, and at the same time able to include part of the desired relativistic corrections. A 1300 eV energy cutoff and a 5 × 5 × 5 Monkhorst-Pack (MP) grid has been used to sample wavefunctions in the Brillouin zone 38 . The exchange and correlation energy has been computed within the Perdew-Burke-Ernzerhof generalized-gradient approximation (PBE) 39 . www.nature.com/scientificreports www.nature.com/scientificreports/ Thresholds for converged geometry configurations are: 10 −8 eV for internal electronic total energies, 10 −3 eV/Å for the maximum residual force, 0.01 GPa for the residual stress and 10 −4 Å for the maximum displacement of atoms in the last three iterations (the window for convergence). By using these parameters, we find that total energy for optimal geometries can be given with a precision of ≈±5 meV. For the calculation of phonons, the internal electronic total energy has been converged further to 10 −10 eV. Geometry optimisations have been performed using the BFGS minimiser 40 . Phonons have been computed from a finite displacement method based on a supercell defined by a cutoff radius of 4.5 Å, and a smearing width of 0.25 meV has been used to obtain the phonon density of states. The transition state has been computed using a quadratic synchronous transit method 41 . Such a formalism yields a reasonable approach to the phonons of pristine α-Fe 42,43 , as it is well documented in the literature, and as it is visualized e.g. in Fig. 7. The phonon density of states in this figure has been used to compute thermodynamical magnitudes for clean Fe.
thermophysics. Vibrations of atoms in the solid have been modelled as a set of independent harmonic oscillators with a characteristic dispersion relation, → w k , where → k takes values in the first Brillouin zone. These phonons have been obtained from the ab-initio DFT calculation, as explained above. In the canonical ensemble (NVT), the partition function of this system is 44 : where the number of particles (N) is fixed, and the volume (V) corresponds to a zero external stress condition.
All the thermodynamic magnitudes of interest can be obtained from the partition function. In particular, our primary interest is to get Helmholtz's free energy: ln( ) ln ( , ) l n 2 sinh 2 ln 2 sinh 2 ( ) ,   the sum is substituted by an integral in the first Brillouin zone using the phonon density of states, ρ(w k ), that obeys the sum rule: Since F w is the characteristic function for the independent variables N, V, T we have: