Stability and mobility of vacancy–H complexes in Al

The effect of hydrogen loading on the stability and mobility of vacancy–H complexes in aluminum is determined by applying DFT and the minimum-mode-following method. The binding energy per H-atom within a complex is found to range from −0.36 eV/atom to −0.34 eV/atom for an occupancy of, respectively, a single and eight H-atoms. When eight H-atoms are neighboring the vacancy the total binding energy becomes −2.72 eV. However, already at a load level of two H-atoms the total binding energy reaches −0.70 eV, which fully compensates the vacancy creation energy. It is observed that for complexes with four or more H-atoms the vacancy gets pinned, as the diffusion barrier increases by a factor of two, reaching a value of 1.03 eV or more. The explanation for the increased energy barrier is that at the higher hydrogen load levels the system must traverse an energetically unfavorable configuration where two or more H-atoms are separated from the vacancy. As a possible consequence of the decreased mobility and increased stability, highly loaded vacancy–H complexes are likely to act as nucleation sites for extended defects.


Introduction
The stability and mobility of vacancy-hydrogen (vac-H) complexes in metals are topics of both technological and academic importance. The wide usage of aluminum and its alloys in aircraft and other high-performance structures makes any knowledge regarding the dynamics of defects and ongoing embrittlement processes, which might lead to fracture, crucial technological information. From an academic viewpoint hydrogen in aluminum is an important system since the results and findings describing the dynamical behavior of H within the Al matrix qualitatively are transferable to Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. other metals where H also behaves as a contaminant of low solubility.
So far, several computational studies have been conducted addressing both vacancies, hydrogen and vac-1H complexes in Al. By applying density functional theory (DFT), Ho et al [1] have determined the vacancy formation energy to be 0.73 eV. Furthermore, they have computed the energy barrier for vacancy migration within a pure Al FCC crystal to be 0.63 eV-in good agreement with the experimentally determined value, which ranges from 0.61 to 0.67 eV. From a series of calculations addressing H in Al, which also applied DFT, Wolverton et al [2] have determined that within a perfect FCC crystal a H-atom has a preference for the interstitial tetrahedral holes (T). However, to embed the H-atom in the crystal an energy of 0.69 eV/atom must be supplied. From experimental observations the reported embedding energy at a T-site is within the range from 0.6 to 0.7 eV and the simulations have yielded values lying at the upper end of this range. For the octahedral holes (O) the embedding energy is computed to be 0.82 eV. From a T-site within the crystal, the barrier for the H-atom to diffuse was determined to be 0.18 eV. If a vacancy and a H-atom are present in the sample these might form a vac-1H complex and gain −0.33 eV, when compared to an isolated vacancy and an H-atom at a T-site in a crystalline environment. It should be noted that the definition-for the energy values reported here-is that a negative value represents a thermodynamically more stable configuration. The structure of the formed complex is that the H-atom resides in a T-like site, which has the vacancy as a nearest neighbor. Finally, Wolverton et al have also determined that for the H-atom to escape the vacancy and break the vac-1H complex a barrier of 0.54 eV must be surpassed. In the computational study by Lu et al [3], a vac-H complex was exposed to an increasing H-load, which yielded the surprising result that a total of twelve H-atoms could be embedded into a single vacancy. It should be noted that the number of T-sites neighboring the vacancy equals eight and the observed phenomenon was named 'superabundancy' for obvious reasons. However, at the higher load levels the H-atoms pair up and arrange themselves in the O-like sites in the molecular form rather than being located in the T-like sites as individual atoms. In an earlier work by one of the authors, the behavior of hydrogen at extended planar defects was studied by applying an empirical potential for the atomistic interactions [4]. From the obtained results it was realized that a twist grain boundary blocks perpendicular diffusion. Furthermore, it was observed that a boundary being both twisted and tilted provides pockets of low-energy sites in between the grains. The presence of these closed pockets greatly reduces the parallel diffusivity in sharp contrast to the commonplace assumption that such grain boundaries furnish a fast diffusion network for hydrogen.
The present work addresses how the stability and mobility of vac-H complexes are affected as the hydrogen occupancy of the vacancy increases. The study is restricted to address hydrogen in its atomic form and the maximal occupancy of the vac-H complex is therefore limited to eight. The simulations are conducted by applying DFT for the atomistic interactions, which is combined with the minimum-mode-following (MMF) method to determine the possible restructuring mechanisms.

Models and methods
To form a structure containing a single vacancy a cubic supercell of 32 atoms in an FCC ordered lattice is generated-using the DFT lattice constant-from which a single atom is removed. The resulting 31 Al atoms are then subject to periodic boundary conditions in the XYZ directions to represent an infinite bulk environment. Translational degrees of freedom are eliminated by keeping a single atom at a fixed position with all other atoms free to move during the simulations. By adding an additional one, two, four or eight H-atoms at interstitial sites neighboring the vacancy, vac-H complexes at increasing H-load levels are constructed. The rather limited system size is a consequence of the cubic shape of the system and the applied DFT method for the atomic interaction.
The specific implementation of DFT used to compute the forces and energy is VASP [5]. An ultrasoft Perdew-Wang 91 GGA functional is applied as the approximation for the inner electrons and exchange/correlation effects. When determining possible restructuring mechanisms the kpoint sampling is 3 × 3 × 3 on a Monkhorst-Pack grid and the plane waves are expanded up to an energy of 150 eV. Of the determined mechanisms, the ones involving the lowest barrier are refined using a cutoff energy of 300 eV and 6 × 6 × 6 kpoints. The determination of the electronic degrees of freedom is considered to be converged when the energy difference of two consecutive iterations drops below 10 −7 eV. When comparing energies obtained with these two parameter sets they are usually converged to within 5 , but it should be kept in mind that, during the searching phase, the priority is to sample as many mechanisms as possible, and performance is desired rather than accuracy.
A scheme for an automatic determination of the position of the vacancy within the structure has been formulated under the assumption that the Al atoms somewhat retain their crystalline structure. With this assumption the vacancy can be considered as the center of mass of the twelve neighboring Al atoms. A simple way to determine these neighbors is through an inspection of the coordination numbers and, by extracting atoms having only eleven nearest neighbors, the required subset results, see e.g. figure 1.

Minimum-mode following method
Navigating on the potential energy surface (PES) and locating saddle points (SPs) by utilizing eigenmodes of the Hessian matrix was initially suggested by Cerjan and Miller [6]. By inverting the force component parallel to the minimum eigenmode-corresponding to the lowest eigenvalue-of the Hessian matrix, an SP is mapped onto a local minimum on the transformed surface, which can be located by applying an ordinary minimization method. In its original formulation the scheme relies on the full Hessian matrix, i.e. the square matrix of the second-order derivative of the energy, which for larger systems quickly becomes an impossibly large task to compute. However, in later works it was proposed to utilize only the minimum-mode. It should be noted that the lowest eigenmode can be estimated efficiently using schemes that only require the first-order derivative of the energy. One such scheme by Henkelman and Jónsson is known as the dimer method [7] and Malek and Mousseau proposed to use the Lanczos method for the mode estimation [8]. Both of the two mentioned schemes greatly reduce the computational demand and enable mode guided SP searches both for systems of larger sizes and for systems where the interactions are determined within a plane wave basis-set.
To initiate a search for an SP a local displacement must be imposed on the system, and it has been determined that elements drawn from a Gaussian distribution result in a better sampling than random elements from a uniform distribution [9]. From the displaced configuration an uphill walk on the PES is conducted. This walk is guided by the minimum-mode, which yields an effective force given by the expressions: where F i,eff and F i are the effective and the true force acting on the system, v i,min is the minimum-mode, which corresponds to the lowest eigenvalue of the Hessian matrix and λ i,min is the lowest eigenvalue in iteration i. In the present work the initial displacement is constructed using elements from a Gaussian distribution with a standard deviation of 0.01 Å. The displacement is always centered on a H-atom and involves all atoms within a 3.0 Å radius. To estimate the minimum-mode, Lanczos method is applied, which is considered to be converged when the change in direction of the estimated mode gets below 0.001 • . The convergence criteria when locating SPs is that the maximal force component on any of the atoms must be smaller than 0.01 eV Å −1 . Following a search it is determined which minima the SP connects and two local minimization runs are conducted. By displacing the system slightly forward and backward along the direction of the minimum-mode at the SP, the two obtained copies of the system are placed within the attraction basins for both the connected minima and will result from ordinary minimization runs. With the minimizations converged the combined outcome of an SP-search is a set of three configurations-two minima and the SP in-between-and the corresponding energy barrier to overcome at 0 K.

Cloud computing, Amazon EC2
Despite the fact that the computational demand is greatly reduced by applying MMF it still imposes a significant computational challenge to conduct SP searches when DFT is applied to account for the atomic interactions. It is, however, important to note that the task of conducting many such searches is an embarrassingly parallel problem that is highly suited for both parallel, distributed and cloud computing. By applying distributed computing Xu et al became able to address the oxidation of Ca on MgO [10] and their later paper on the diffusive behavior of Li on MgO surfaces [11] is another work where DFT and MMF were combined to describe the atomistic dynamics.
The calculations presented here were possible to complete within a relatively short time period by utilizing Amazon's cloud resources. To gain access to their elastic cloud the existing EON software package 5 [12] was extended with a new communicator implemented utilizing boto 6 , which is a module providing interfaces to Amazon's cloud services for the Python programming language. The computations were done on first generation m1.large and c1.xlarge compute instances 7 that are equipped with either four EC2 compute units and 7.5 GB memory or twenty EC2 compute units and 7 GB memory. Using a m1.large compute instance a single SP-search and the following determination of the two connected minima could on average be completed within two days. The flexible communication framework implemented in EON enable the execution on cloud facilities through a local master node. This master node managed a pool of slaves-the Amazon compute instances-by requesting and collecting SPs. Besides the bookkeeping of the located SPs the master also progressed the system into a connected minimum when a sufficient sampling of the local environment was obtained. A total of 278 unique SPs were determined (120 for 1H, 86 for 2H, 36 for 4H and 36 for 8H) within a total of 4295 m1.large wall hours and 3881 c1.xlarge wall hours, which had been granted by Amazon.

Results and discussion
To inspect and confirm the accuracy of the applied DFT parameters and evaluate eventual size effects, three series of simulations are conducted. The first series is to determine the formation energy and migration barrier for a vacancy in Al  of eight such sites surround a vacancy at a distance of 1.66 Å from its center. Furthermore, a total of six O-like sites are also present in the first shell of sites, at a distance of 1.62 Å from the vacancy. However, these sites are only stable for the calculations using the 3×3×3 kpoints where a barrier of a few meV results. As the O-like sites are completely destabilized by the presence of the vacancy when applying 6 × 6 × 6 kpoints, these are considered as quasi-SPs connecting the four T-like sites within the plane spanned by these. The barrier for the H-atom to reach an O-like configuration from a T-like site is 0.24 eV. Besides connecting T-like sites neighboring the vacancy, the O-like site also provides a path for the H-atom to leave the first shell and break the complex. To reach a T-site in the second shell, which is located at a distance of 3.40 Å from the vacancy, an additional barrier of 0.23 eV must be overcome, which yields an effective barrier of 0.47 eV in order to break the complex. This effective barrier is in good agreement with the barrier determined by Wolverton et al. In the broken configuration the energy is increased by 0.36 eV and a reverse barrier of 0.11 eV prevents the complex from reforming. To reach a T-site within the third shell a barrier of 0.12 eV must be surpassed from a site within the second shell; however, the effective barrier from the first shell remains basically unchanged at 0.48 eV. A surprising feature of the third shell is that the binding energy is 80 meV stronger than for sites belonging to the second shell. To determine if this is an artifact of the limited system size the energy difference is also determined using the larger system, as described above. In this case the same qualitative behavior is observed but the difference is reduced to 40 meV. Furthermore, both LDA and GGA give this result, which could, however, an effect related to the surface effect observed by Carlin et al. The T-sites within the third shell are located at a distance of 4.38 Å from the vacancy. The relative energies for the three shells and effective barriers for the H-atom to move between these are listed in table 1.

Vacancy and hydrogen complex
Another mechanism that breaks the complex is if the vacancy diffuses. For such a transition to take place, a 0.54 eV barrier must be surpassed and the vacancy migrates to the lattice site diagonally opposite to the H-atom, which therefore recedes into the third shell of T-sites. Of the describer mechanisms, the most relevant energies and barriers are listed in table 3 and shown in figure 2. Finally, concerted events-where both the vacancy and the H-atom move-are also observed. Of the events that preserve the strongest bound vac-1H complex, the lowest barrier is 1.02 eV, see figure 3 for a schematic illustration. At the SP the H-atom is located at an O-like site, while the Al-atom passes below, causing the distance between these two atoms to decrease to 1.71 Å, which is a slight compression compared to the 2.04 Å when the H-atom resides in an O-like site. The observed compression is considered to explain the relatively high energy barrier for this mechanism.

Complex with two hydrogen atoms
For a complex with two H-atoms the strongest bonding arises when both H-atoms neighbor the vacancy; however, they should be placed at diagonally opposite T-like sites. For this configuration the total gain in energy is −0.70 eV or −0.35 eV per H-atom. It is noteworthy that the total binding energy fully compensates the vacancy formation energy. The mutual interaction between the H-atoms appears to be insignificant, as the energy for a configuration where the H-atoms are placed next to each other-within a plane, see figure 1-results in an energy change of only less than 10 meV. For this transformation to take place a barrier of 0.27 eV must be traversed.
As for the complex with a single H-atom the O-like configurations neighboring the vacancy are destabilized and the barrier to be surpassed for a H-atom to enter the second shell of T-sites remains 0.47 eV. For this configuration an increase of 0.34 eV results as compared to the strongest bound complex. When comparing this change to the value for the similar configuration of the vac-1H complex these are very similar, and it seems reasonable to consider the complex as broken when a H-atom leaves the first shell. To reform the broken complex a barrier of 0.13 eV must be surpassed. To reach the third shell the effective barrier to overcome is similar to the barrier to enter the second shell. Again the third shell offers more favorable configurations as compared to the second shell and a configuration with a single H-atom within this shell has an energy of 0.30 eV.
For both H-atoms to leave the vacancy and enter the second shell an effective barrier of 0.86 eV needs to be Circles mark H-atoms within the first shell and crosses mark H-atoms in the second shell. Leftmost column, a pinned vacancy neighboring four H-atoms, where the pinning occurs because a configuration with two or more H-atoms within the second shell-sites of higher energy-must be traversed for the vacancy to move. Middle column, the vacancy remains mobile when two H-atoms are embedded. Rightmost column, concerted mechanism where both the vacancy and the H-atom move; however, the barrier for this event is more than twice the value for an event only displacing either the H-atom or the vacancy. Table 2. Relative energies and effective barriers for vac-2H complexes. The first three columns describe the shell configuration. The latter two columns are the effective barrier for the transformation to take place and the energy relative to the strongest binding complex, respectively. All energies and barriers are reported in eV. From a configuration where the H-atoms reside in the first and second shell it is possible for the vacancy to diffuse via low barrier mechanisms. The mechanism with the lowest barrier is when the vacancy moves and thereby shifts the shell configuration of the H-atoms. The occurring shift is that the H-atom that was located within the second shell enters the first and vice versa, see figure 3 for a schematic representation. The barrier for this event is 0.65 eV. It should be noted that an analogous mechanism exists for configurations with a H-atom in the first and third shell involving a similar barrier.

Complexes with four or eight hydrogen atoms
The final two series of simulations address vac-H complexes with either four or eight H-atoms. When four H-atoms are next to the vacancy they arrange themselves in a diagonal plane. In this configuration the binding energy per H-atom remains −0.35 eV and the total binding energy at this coverage reaches a value of −1.40 eV. For a H-atom to rearrange within the first shell an energy barrier of 0.21 eV must be traversed, while for a H-atom to escape the complex the barrier to be surpassed has decreased to 0.47 eV. When a total of eight H-atoms are within the complex, all T-like sites neighboring the vacancy are occupied and no processes for rearrangements within the first shell are determined. In this complex the energy per H-atom has decreased slightly to −0.34 eV, while the total binding energy reaches a value of −2.72 eV, which is the maximum for all the addressed vac-H complexes. The energy barrier for a H-atom to escape this complex is 0.38 eV, which is the smallest determined value, but still more than twice the value for H diffusion in bulk Al.
The most striking consequence of the increased load level is that a pinning of the vacancy occurs. For the configuration with four H-atoms the vacancy diffusion barrier increases to 1.03 eV, which is approximately a factor of two greater than the value for diffusion within crystalline Al. To understand this sudden increase one should note that only for systems containing one or two H-atoms is it possible to displace the vacancy without traversing a configuration where more than a single H-atom is separated from the vacancy. However, when four or more H-atoms are embedded, the system is required to enter a configuration where two or more H-atoms are separated from the vacancy. If one subtracts the determined barrier for vacancy diffusion within crystalline Al (0.49 eV) then an additional barrier of 0.54 eV results, which must be surpassed for the vacancy to move at this load level. It should be noted that this additional energy barrier is within a 0.1 eV energy window of the penalty resulting when a vac-4H complex is broken into a vac-2H complex and two H-atoms in a crystalline environment. A schematic example of this is shown in figure 3.

Conclusions
The lowest energy of a vac-H complex containing a single H-atom is determined to be a configuration where the H-atom resides in the tetrahedral-like environment next to the vacancy. For this configuration a binding energy of −0.36 eV is obtained. In order to break the complex the H-atom must leave the first shell of sites neighboring the vacancy and enter the second shell, which involves a barrier of 0.53 eV. When two H-atoms are neighboring the vacancy the binding energy per atom decreases slightly, reaching a value of −0.35 eV/atom; however, the resulting total binding energy exceeds the vacancy formation energy. For a complex with a H-atom at T-sites in the first and second shell an increase of 0.34 eV results, and the barrier for the transformation breaking the complex is 0.47 eV. For a configuration with eight H-atoms the binding energy per H-atom further decreases to −0.34 eV/atom, whereas the total gain in energy reaches a maximum of −2.72 eV. The barrier for vacancy diffusion for a complex with a single H-atom is 0.54 eV, which is similar to the conditions in bulk Al. When two H-atoms are present the barrier increases slightly to 0.65 eV. However, when the H-loading of a complex reaches four H-atoms then pinning of the vacancy is observed to occur. For the vac-4H complex the vacancy diffusion barrier is increased by a factor of two, reaching 1.03 eV. This increased barrier is due to the fact that for the vacancy to move the system is required to traverse through energetically unfavorable configurations where two or more H-atoms are outside the first shell of T-like sites neighboring the vacancy. A possible consequence of the observed pinning is that a highly loaded vac-H complex might act as a nucleation site for the formation of an extended defect by capturing the more mobile vacancies, H-atoms, and vac-H complexes.