Enhanced H-H binding and consequent H-aggregation around dislocation in α-Fe lattice

H-H binding energies in tetragonal and octahedral sites (TS & OS) of α-Fe lattice were calculated by density functional theory (DFT) under correction of elastic energy. Strong attractive interactions were identified as H atoms were incorporated in 3rd and 4th nearest OS neighbors. OS-type H-aggregated clusters with binding energies exceeding 200 meV were identified. Monte Carlo simulation of H-loading indicates abnormal H-aggregation behavior around a 1/2[111] ( 10 1 ¯ ) edge dislocation in relation with the enhanced H-H binding specifically in OS of α-Fe lattice.


Introduction
Enormous efforts have been put into understanding and prevention of hydrogen embrittlement (HE) in metals for over a century [1]. Understanding the HE mechanisms is a prerequisite for HE prevention. The key scientific issue that how the H atoms interact with lattice defects like solute atoms [2], dislocations [3,4], voids [5,6] etc. Some issues such as H interaction with metal surfaces is of fundamental significance not only for HE but also for understanding of effective catalysis [7,8]. Generally, atomistic details of H is tough to be revealed by present experimental methods. Ab initio calculation based on density functional theory (DFT) has been an indispensable implementation to problems mentioned above. Thanks to the increasing computational power of supercomputers, DFT calculations can handle models containing atoms upto thousands, which lead to the recent progress on understandings of H behaviors in large nanovoids [9] and H-bubble formation [10].
H-H interactions in lattice are the key parameters to understanding of related phenomena such as H aggregation and H-bubble formation in metals [10]. However, uncertainty of the H-H interactions in α-Fe exists among different works, leading to controversial predictions about H behavior. For example, as in reference [2], H-H binding energy of 4th nearest-neighbor configuration is about 0.03 eV, while in reference [6] the value is −0.008 eV. Besides, H-H binding energy in OS has never been reported, even though it is well known that shear strain can change H preference from tetragonal interstitial sites (TS) to octahedral sites(OS) [11,12]

Binding energies of H-pairs and H-clusters
The H-H interaction of H-pairs was investigated for H atoms located in different neighboring sites, which is evaluated by binding energy (E b HH ) according to equation (1): Similarly, average binding energy of single H atom in a H-cluster containing N hydrogen atoms can be calculated according to equation (2):  [13] with the projector augmented wave method(PAW) and ultrasoft pseudopotentials [14]. The exchange correlation were described by the generalized gradient approximation (GGA) [15] with the Perdew-Burke-Ernzerhof functional (PBE) [16]. Spin-polarized calculations were employed in all cases. The cutoff energy for the plane-wave basis set is 450 eV, at which the convergence of calculation had been achieved. Γ-centered´3 3 3kpoints was generated by Monkhorst Pack method [17] and used in structural relaxation and the Methfessel-Paxton smearing method with a width of 0.1 eV was used [18]. Structural relaxation is terminated when the maximum force acting on the movable degrees of freedom becomes less than 20 meVÅ −1 and energy change is within 10 -4 eV. In self-consistent calculations Γ-centered´5 5 5kpoints was employed and tetrahedron method with Blöchl corrections was applied to obtain accurate total energy of the system [19]. Convergence of self-consistent calculations were acquired with a convergence threshold of 10 −6 eV. DFT calculations was carried out under periodical boundary condition [20]. Unfortunately, this approach inevitably includes artificial interactions among defects and their periodical images, jeopardizing the accuracy of DFT prediction. Different researchers adopt different size of supercells or different constraints (zero-stress or zero-strain) while doing DFT calculations, this may explain the inconsistency of these results. Accuracy of these DFT results is of crucial importance because they are usually adopted as basic inputs in developing interatomic potentials used in molecular dynamics (MD) or Metropolis Monte Carlo (MC) simulation. In order to address the problem mentioned above, Varvenne et al proposed an approach coupling DFT calculations with linear elasticity theory to obtain the properties of the isolated point defect with reduced supercell sizes, which is also an effective tool adapt to different choices of constraints [21]. In their work, energy of isolated defect ( ¥ E D ) was obtained by energy of defect in finite size ( e= E D 0 ) minus one half of the sum of its interactions with periodical images ( / E 1 2 p int ), described by equation (3)- (6).
As in this work, defects were interstitial H atoms or H-H pairs. P ij is the elastic dipole tensor, C ijkl is modulus tensor of α-Fe, e kl is the homogeneous strain applied on the supercell, s ij is the homogeneous stress of the periodic simulation cell of volume V containing a defect.
is the second derivative of the anisotropic elastic Green's function with respect to the coordinates R .
3 corresponds to the position of the defect periodic images, with A , 1 A 2 and A 3 being periodic vector of the supercell. The term = = = n m p 0 has been excluded from the sum in equation (5). Elastic modulus used in elastic correction were quoted from the DFT study by Psiachos et al (C 11 =284 GPa, C 12 =149 GPa, C 44 =105 GPa) [22].
The zero-point energy correction (ZPE) correction of the H binding energy within α-Fe was calculated from the Hessian matrix eigenvalues. We assumed that Fe atoms were heavy enough comparing to H atom so that the ZPE can be approximately calculated by the sole displacement of H atoms. The corresponding Hessian matrix was calculated by both shifting the H atom in each of the ±X, ±Y and ±Z directions with 0.015 Å and observing the forces acting on them. The ZPE was then calculated according to the following equation: where k i are the eigenvalues of the3 3 or6 6Hessian matrix, m H is the mass of an H atom and h is the Planck's constant. 3. Results and discussion 3 [6]. Different sizes of supercells (froḿ3 3 3to´10 10 10were adopted to test the convergence of the E HH calculated by the stand-alone MS and MS corrected by elastic approach. Cubic cell with 2 atoms per cell was used to construct supercells. The MS calculation was carried out via LAMMPS [26] under constant volume condition and internal stresses of different configurations were obtained simultaneously to serve as the input of the elastic correction. Results are demonstrated in figure 1. In figure 1, black lines represent H-H binding energies without elastic correction and red lines with elastic correction, denoted as 'HH' and 'HHels' respectively.
Generally, E HH converges rapidly as the size of supercell increases. With or without elastic correction, convergence within 1 meV can be reached when the size of supercell is as large as´10 10 10.As shown in figure 1, E HH with elastic correction exhibits better convergence. For instance, even the supercell size is as small as´3 3 3, E HH with correction is only 2∼3 meV higher than the converged value. On the contrary, the deviation of the uncorrected one is approaching 10 meV. It is also shown in figure 1 that in small supercells the uncorrected E HH calculated under constant volume condition is always lower than the converged values, indicating overestimation of repulsive interactions of H-H pairs calculated in small supercells. Thus, it is confident that accuracy of H-H binding energies within 2 meV can be attained with elastic correction in supercells larger than´4 4 4.

Binding energies of H-pairs and H-clusters
Calculation of H-H binding energies in nearest neighboring TS(NNTS) and OS(NNOS) were carried out by DFT methods in´4 4 4supercells under constant volume condition. Convergence with´4 4 4supercells was assured by the convergence test.
H-H binding energies of TS and OS calculated by DFT method with elastic correction are schematically shown in figure 2(b). The calculation covers TS configurations from 3rdNNTS to 12thNNTS, and OS from 3rdNNOS to 4thNNOS. There is no stable configuration of either 1stNNTS, 2ndNNTS, 1stNNOS or 2ndNNOS because repulsive interactions of these configurations are very strong. For H-H pairs with multiple variants such as 8thNNTS, 9thNNTS, 10thNNTS and 4thNNOS exhibited in, only the most stable one is shown in figure 2 As is shown in figure 2(b), DFT prediction of H-H interactions is attractive beyond 3rdNNTS, stronger than those calculated by Song and Curtin [6] but weaker than those calculated by Counts [2]. For the first time, H-H binding energies of 3rdNNOS and 4thNNOS are reported as 60.6 meV and 60.5 meV respectively, indicating much stronger attractive interaction of H-H pairs in OS than those in TS. It can be predicted that H behavior within strain concentrated area are more dependent on H-H interactions in OS rather than TS.
Generally, binding energy of H atoms increases with the size of H-cluster with proper siting of H atoms. Known the attractive configurations of H-pairs, it is possible to construct H-clusters with H atoms mutually attractive to their nearest and second-nearest neighbors. A H-cluster with every single H located in TS was constructed by repeatedly incorporating H-atoms in 4th NNTS and 8th NNTS till a (110) plane was filled, denoted as TS-cluster. Analogously, a monolayer (110) OS-cluster and a bilayer (100) OS-cluster with H-atoms incorporated in 3rd NNOS and 4th NNOS was constructed, denoted as OS-cluster1 and OS-cluster2 respectively. Configurations of these H-clusters after structural relaxation are shown in figure 3 Average binging energy of single H atom within H-clusters mentioned above were calculated to be 32 meV, 266 meV and 236 meV respectively for TS-cluster, OS-cluster 1 and OS cluster 2 according to equation (2). It is obviously shown in figure 3 In the equation above, T is the temperature, k B is the Boltzmann constant andĒ b H N is average binding energy of single H atom in a H-cluster. At temperature of 300 K, it can be calculated that C H * for OS-cluster 1 is 34 appm, for OS-cluster 2 is 110 appm, and for TS-cluster is 22240 appm. In other words, H-cluster will not form at 300 K unless bulk concentration of H is over 34 appm.

H-loading around a 1/2[111](¯) 101 edge dislocation
Results of H-loading around a 1/2[111](¯) 101 edge dislocation are shown in figure 4. there is almost no difference in H distribution around the dislocations when the bulk concentration of H is 1 appm and the potentials used for H-loading are not concerned. And this is consistent with the prediction based on DFT that critical bulk-concentration C H * below a certain value will not form H-cluster. However, difference among these potentials begins to emerge as the bulk concentration reaches 10 appm. For Song and Wen19, H concentrations around the dislocation cores are higher than that of Wen16. Specially, H concentration around the dislocation core charged by Wen19 exhibits a double-humped profile. As the bulk concentration further rise to 100 appm, H-cluster growth phenomenon can be observed in H-loaded structures using Song or Wen19. On the contrary, H concentration in charged structure by Wen16 only exhibit double-humped profile.
H occupation site around the dislocation was investigated by direct observation from different directions. Data visualization and structure analysis was carried out using OVITO [27]. Take H-loaded structure by Potential Song for example, projection of atomistic structure along [100], [010] and [001] zone axes were taken and shown in figure 5. For clarity, Fe atoms were set as small yellow spheres and H atoms were big ones. The whole area was factitiously separated into four quadrants centered on the dislocation core as is shown in figure 5(a), and H atoms located in different quadrants were colored blue or pink accordingly. H atoms right at  the dislocation core is not shown. It is clearly shown in figure 5 that H atoms around the dislocation always in line with a 〈100〉column, indicating they are occupying OS rather than TS around the dislocation. Secondly, H atoms located in different quadrants are in different types of OS. It is obviously due to the opposite sign of shear strain within different quadrant around the dislocation [28]. Thirdly, projected images indicate these atoms incorporate 3rd NNOS and interatomic distance of the H atoms confirmed this observation, which strongly indicates the growing mechanism of H clusters is via occupying 3rd NNOS (probably also 4th NNOS) repeatedly. It is worth to mention that there are 8 equivalent 3rd NNOS for a H atom in OS pointing to 8 equivalent 〈111〉 directions, allowing the H-cluster to grow in 3 dimensions (3D). If 3rd NNOS is unfavorable  such as the case of Wen16, the only source of driving force for cluster growth is the attractive interaction of H in 4th NNOS, allowing expansion of H-cluster within a (100) plane along [010] and [001] direction, which is only a 2D growth. As a result, H-loading with Wen16 didn't cause H-aggregation or hydride-growth as that with Song and Wen19. Clearly, the double-humped hydrogen distribution profile around the dislocation in this work is inconsistent with the binding energy mapping of single H atom around the dislocation [28]. Such anomaly should only attribute to the strong attractive H-H interactions in OS demonstrated in this work. With bulk H concentration up to 100 appm, the 'H hump' could grow into iron-hydride plates which even extended far from the stress center of the dislocation, indicating a self-clustering of H, and within the cluster OS occupation is strictly preserved. It is clearly shown in figure 5(a) that plane of iron atoms shuffled a bit within the H-aggregated region, forming hydride of FCC structure which was discussed in detail in our earlier work [29], being consistent with the DFT prediction of structure of OS-cluster 1 in this work. Dissociation of dislocation is one possible mechanism for cross-slip of dislocation. Following our earlier discussion, formation of iron hydride (H-clusters) lower the critical stress for dissociation of dislocation and facilitates plasticity [29]. With spontaneous growth of H-cluster, it can be expected that mobility of dislocations can be enhanced. Increased defects accompanying dislocation movement and cross-slipping can also be realized. A strain relaxation test on H-loaded pure iron reveal steep reduction of activation volume and rise of dislocation velocity at a critical H concentration of approximately 20 appm [30]. Moreover, a '30 appm' anomaly was reported recently by TEM observation of H-loaded pure iron, demonstrating increased dislocation cross-slip or climb as well as dislocation debris of dislocation slipping [31]. These experimental evidence of a critical H concentration agrees quite well with C H * for spontaneous growth of OS-cluster 1 (27 appm at 293 K) determined in this work, indicating a possible connection between H-clusters and enhanced plasticity in α-Fe. H-clusters enhanced dislocation dissociation might be an important supplement to mechanisms of HE.

Conclusions
A comprehensive study of H-H interactions in α-Fe lattice was performed with correction of elasticity and zeropoint energy, including interactions of H in octahedral sites. Conclusions are as follows: (1) Strong H-H attraction was revealed when H atoms were incorporated in 3rd and 4th nearest neighbour octahedral sites, with binding energies of 60.6 meV and 60.5 meV respectively.