Thermodynamic analysis of the interaction between metal vacancies and hydrogen in bulk Cu

Using grand canonical thermodynamic analysis with inputs from DFT calculations we calculated equilibrium molar fractions of copper vacancies (V Cu ), H interstitials (H i ) and their complexes in bulk Cu in a wide range of temperature and hydrogen pressure values. The results show that the equilibrium molar fractions of both V Cu and H i are low in most conditions of interest, in good agreement with available experimental data. Although H i –V Cu complexes have significantly lower formation energies than the isolated defects, the low molar fraction of H is predicted to have little impact on the rise in vacancy molar fraction for external hydrogen pressures below 100 bar. Only at relatively high hydrogen pressures exceeding 10 kbar in the presence of Cu vacancies, the H molar fraction was found to reach the same order of magnitude as the molar fraction of vacancies. These results put thermodynamic limits on the hydrogen-induced vacancy clustering and void formation in bulk Cu.


Introduction
7][8] SAV formation can be induced because H atoms trapped in vacancies decrease the formation energy of larger vacancy clusters. 9,10However, theoretical studies of H trapping in fcc metals have yielded inconclusive results 10 regarding the role of vacancies.In particular, Lu and Kaxiras, 11 who employed density functional theory (DFT) in the local density approximation, demonstrated that up to 12 H atoms may fit into a single vacancy of aluminum (Al).Ismer et al. 12 have shown that up to 10 H atoms can be trapped in an Al monovacancy using the generalized gradient approximation (GGA) of DFT.Similarly, it has been reported that several H interstitials can be trapped in Pd vacancies. 10,13In Cu, DFT calculations have shown that a single Cu vacancy can trap up to 6 H atoms, 14 whereas other theoretical studies have suggested that it is possible to trap up to 8 H atoms in a Cu vacancy. 15though these theoretical studies demonstrate that one metal vacancy can accommodate a significant number of H atoms, the amount of H that can be incorporated into a metal sample under different conditions remains difficult to predict because it is strongly affected by both thermodynamic and kinetic factors.In metal films, H can be introduced as a result of the use of water-based solutions during the deposition process. 16,17In such a case, H can be kinetically trapped inside the sample due to defects like dislocations and grain boundaries, which act as efficient sinks. 18,19On the other hand, H can be incorporated into a metal under H 2 gas pressure at high temperatures. 20In particular, high molar fractions of H have been reported in Ni, Fe, and Pd even at temperatures of 600 K and partial H 2 gas pressure of 1 bar.Experimental studies in electrodeposited (ECD) Ni have measured H molar fractions of 7 Â 10 À2 at high pressure values. 6lthough measuring the amount of H trapped in defects remains difficult, fairly accurate predictions can be made regarding the molar fraction of incorporated H under equilibrium conditions using statistical thermodynamics. 21Davidson et al. 21conducted theoretical studies of a-Fe using DFT calculations and statistical mechanics, finding equilibrium H molar fractions of up to 10 À4 at partial pressures below 1 bar.Even higher molar fractions for fcc Fe have been reported by Nazarov et al. 22 However, Fe, Ni, and Pd have a much higher H solubility compared to Cu. 23 Furthermore, vacancies in these metals reach molar fractions of 10 À3 , 24 whereas in Cu, vacancy molar fractions are estimated to be one order of magnitude lower, even at high temperatures near the melting point. 25xperimental studies for ECD Cu have revealed H molar fractions of around 10 À4 near the melting temperature of Cu under extreme H partial pressures exceeding 10 kbar (1 GPa). 6,26ecent theoretical studies have investigated the structure of H complexes and vacancies in Cu. 15,27 However, it remains unclear how much H can be incorporated into perfect Cu lattice at more common pressures of 0.1 to 100 bar and whether SAVs can actually form under such conditions or much higher pressures and temperatures are required.Here we address these questions using theoretical simulations and statistical thermodynamics.
The first step of incorporating atomic H into Cu from the gas phase is the adsorption of H 2 on Cu surfaces.As described in previous theoretical studies, 28,29 this process leads to the dissociation of the H 2 molecule into atomic H, which will occupy hollow sites on the Cu surface. 28This reaction was found to have a barrier of 0.50 to 0.56 eV, depending on the exposed Cu surface. 29Atomic H can then diffuse from the Cu surface into the bulk with incorporation energies ranging between 0.64 and 0.90 eV. 30As described in previous works on H in metals, 9,31 the presence of H is expected to increase the concentration of metal vacancies.The increase in vacancy concentration due to the presence of H can be explained by the creation of lattice strain accomplished both by physical 32 and chemical 33 mechanisms.As a result of the incorporation of H, the bond strength between adjacent Cu atoms and that between neighboring host Cu atoms is expected to reduce, 34 facilitating the formation of vacancies.H also interacts with thermally formed Cu vacancies and can become trapped by them. 10,22,24However, for this interaction to be significant, H has to be closer than at the second nearest neighbor (2NN) distance from a Cu vacancy. 15he formation of a H i -V Cu complex can then facilitate the diffusivity of either H 31 or the whole complex. 15n this work, we neglect these dynamic processes and provide a thermodynamic analysis of the formation of H interstitials (H i ), Cu vacancies (V Cu ), and H i -V Cu complexes in perfect fcc Cu.By employing a combination of atomistic simulations and grand canonical statistical mechanics, we calculate the probabilities of formation of different configurations of clusters of up to 6 Cu vacancies as well as the equilibrium amount of H incorporated into perfect Cu bulk as a function of H 2 partial pressure (from 0.1 to 50 kbar (5 GPa)) and temperature.The effect of H i on the concentration of V Cu at various partial pressures and temperatures was then assessed.The results show that in Cu, unlike in other metals where H is more soluble, the molar fraction of absorbed atomic H at equilibrium is expected to be relatively low.Even at H partial pressures of 100 bar and temperatures close to the melting point of Cu (H absorption increases with temperature in this system), the equilibrium molar fraction of H is orders of magnitude lower than the equilibrium molar fraction of Cu vacancies.However, at partial pressures above 10 kbar and temperatures close to the melting point of Cu, the H molar fraction approaches that of vacancies.Only at these conditions does the presence of vacancies affect the equilibrium molar fraction of H in Cu.In other words, under normal temperatures and pressures there are much more vacancies than hydrogens.Therefore hydrogen atoms are unlikely to affect the number of vacancies.At high external pressures of hydrogen, more atoms incorporate and the presence of vacancies facilitates this process further.
The paper is organized as follows: in Section 2, the technical details of the DFT and force field calculations are summarized, and the grand canonical ensemble statistical model is outlined.The computed equilibrium molar fractions of vacancies and H in bulk Cu based on the employed thermodynamic model are presented in Sections 3.1 and 3.2, respectively.In Section 3.3, we present the DFT results for vacancy-H complexes and discuss how the presence of one type of defect affects the equilibrium concentration of the other.In Section 4, a summary and conclusions are provided.

Methods of calculation
2.1 DFT calculations DFT calculations were carried out for 32 and 108 atoms periodic cells using the Vienna Ab Initio Simulation Package (VASP) [35][36][37] and the Perdew-Burke-Ernzerhof (PBE) GGA exchange-correlation functional. 38The larger supercell was used in cases where 3 or more vacancies were introduced.For the 32 and 108 atoms cells, in line with previous works, 10,39 converged 6 Â 6 Â 6 and 4 Â 4 Â 4 k-point grids were used, respectively, with the 450 eV energy cutoff.Cu was treated using 11 valence electrons (3d 10 4s 1 ).The atomic positions were relaxed to the energy tolerance of 10 À5 eV.In accordance with previous DFT studies on H trapping in transition metals, 21,27 no self-interaction (+U) corrections have been used since, unlike in Cu oxides 40,41 and other transition metal oxides, 42 H at the examined concentrations is not affecting the nature of Cu-Cu bonding.
The formation energies (per vacancy) of the vacancy clusters were calculated as: where n is the number of Cu vacancies; E n and E 0 are the energies of the crystals with n and no vacancies, respectively; and m Cu is the chemical potential of a Cu atom in the bulk (approximated here as simply the DFT energy per atom of the pristine metal crystal).

Force field calculations
4][45][46] In particular, the EAM has been widely used to investigate vacancies in metals. 47,48The EAM potential developed by Mishin et al. 49 was used for all force field simulations.This potential has been fitted semiempirically and is being widely used for pure Cu systems. 50,51To carry out the EAM simulations of vacancy clusters in Cu, the Large-scale Atomic/ Molecular Massively Parallel Simulator (LAMMPS) 52 was used.We used the conjugate gradient method with the force tolerance of 10 À10 eV Å À1 for geometry optimization.Further cell equilibration was carried out under an isothermal-isobaric NPT ensemble at constant temperature (0 K) for 0.02 ns using a timestep of 1.0 fs.

Enumeration of defect configurations
One of the goals of this work is to identify molar fractions of various vacancy clusters of different sizes and configurations in bulk Cu.Any configuration of vacancies located in a cubic volume of 7.23 Å Â 7.23 Å Â 7.23 Å comprised of 32 atoms was considered a cluster.This volume was chosen in order to limit the number of considered vacancy configurations.The Site-Occupation Disorder (SOD) code 53 was used to identify nonequivalent V n configurations by accounting for the system's symmetry.If an isometric transformation or symmetry operation that can change a cluster configuration into one that has already been discovered is found, these configurations are considered to be equivalent.Following this process, 14 V 3 , 72 V 4 , 223 V 5 , and 870 V 6 non-equivalent configurations were identified.Using the EAM potential, 49 the geometries of 500atom cells containing each vacancy cluster configuration were fully optimized.The lowest energy configurations were then reoptimized using DFT in 108-atom simulation cells.

Grand-canonical analysis
To calculate equilibrium concentrations of complex defects including inter-defect interactions, one must compare energies of non-equivalent configurations using a suitable supercell.In the typical formulation, the probability, P m , of a given configuration (m) of lattice site substitutions or defects depends on the configuration energy, E m , and its degeneracy, O m , in the configurational space as follows: where is the canonical partition function, k B is the Boltzmann's constant, and the index m runs over all the symmetrically non-equivalent configurations.While this approach works well for the simulation of solid solutions, [55][56][57][58] it is not adequate for the treatment of dilute defects, because impractically large supercells would have to be considered.A more suitable approach in that case, still based on periodic DFT calculations and Boltzmann's probabilities, is to employ a grand-canonical ensemble of configurations, 59 where configurations with different compositions are included via the introduction of the relevant chemical potentials.For example, the probability of the mth configuration with n Cu vacancies in the supercell may be calculated as: where is the grand canonical partition function, E nm and O nm are the energy and degeneracy, respectively, of the specific configuration, n is the number of vacancies or interstitials, and m = Àm Cu is the chemical potential of vacancies, which is simply the negative value of the chemical potential of Cu atoms in the Cu lattice.By mixing the configurations with no vacancies and configurations with vacancies or clusters of vacancies, any defect composition can be represented, and the probability of any vacancy cluster can be calculated.The equilibrium total molar fraction of vacancies (x in Cu 1Àx & x , where & represents a vacancy) at any temperature is then obtained as: where N is the number of Cu sites in the supercell.
The same formalism can be used for the calculation of the equilibrium molar fraction of interstitial H atoms in bulk Cu.In this case, the chemical potential is a function of both temperature and of the external H 2 partial pressure (p), if the crystal is considered to be in equilibrium with H 2 gas: where g H 2 is the free energy per molecule of the H 2 gas.At low pressures, the ideal gas law for the pressure dependence can be used: where E DFT is the energy computed using DFT/PBE for an H 2 molecule in vacuum, and E ZP = 0.27 eV is the zero-point energy of the H 2 molecule, calculated from the computed vibrational frequency n = 4308 cm À1 , which is in reasonable agreement with the experimental frequency (n = 4401 cm À1 ). 60Dg(T, p 0 ), which corresponds to the change of the free energy per H 2 molecule from the absolute zero to temperature T at reference pressure p 0 = 1 bar, was extracted from the thermodynamic tables compiled by Chase. 61his is a common approach when evaluating the chemical potential of H 2 or O 2 gases, 59,62,63 to avoid the calculation of DFT free energies in the gas phase.For the analysis at high partial pressures (kbar range), where the ideal-gas approximation breaks, the H 2 free energies (relative to the zero point) were obtained from the high-pressure thermodynamic tables by Fukai. 64inally, an ensemble of supercell configurations including both n Cu vacancies and k H interstitial atoms is considered.The probability of the mth configuration with a given (n, k) composition is given by: where E knm are O knm are the energy and degeneracy of the configuration, respectively.Such calculations allow us to evaluate the probabilities of formation of H i -V Cu complexes.
The equilibrium molar fraction of H (d) in the crystal is then given by: In the equations above, the configuration energies are taken directly from the simulations, ignoring any vibrational contributions, to avoid the huge computational cost of vibrational analysis for all configurations in the ensembles.However, vibrational contributions can sometimes be significant in the thermodynamics of defect formation in metals, 24 especially when light atoms like H are involved.As a compromise, vibrational free energy corrections are introduced to the energies of the configurations involving interstitial H, by using the harmonic approximation to define a contribution from each frequency mode n involving H: where the vibrational frequencies are obtained from constrained DFT calculations in which only the H atoms were allowed to move.This approach significantly reduces the computational cost involved in phonon frequency evaluation while capturing the main contributions.Due to the low molar fractions of H incorporated (as will be discussed below), the statistical analysis in our calculations will be truncated to exclude configurations with k 4 1 (i.e. with more than one H atom per supercell), so vibrational calculations only had to be done for a few configurations.For an isolated H i in Cu, the DFT frequency of the triply degenerate mode was ca. 25 THz.In the configuration where H i is nearest neighbour to a V Cu , the mode degeneracy breaks into two 26 THz modes and one 21 THz mode.

Vacancy clusters and their equilibrium concentrations
The formation energies per vacancy for clusters with up to 6 vacancies are included in Fig. 1.The lowest-energy configuration for each cluster is shown in the inset images.In this and all the following structure figures, the same color code is used, with blue corresponding to Cu atoms, red representing Cu vacancies, and H atoms shown in white.Vacancies were found to form close-packed clusters, with V 4 forming a stacking fault tetrahedron and V 6 forming a symmetric diamond.The same lowest energy configurations have been reported for Cu vacancies by Ganchenkova et al. 27 Also, in agreement with Ganchenkova et al., 27 the formation energies differ within less than 0.1 eV/vacancy for the lowest energy configurations of each cluster.As the number of vacancies increases, so does the binding energy of the clusters while the formation energy per vacancy constantly reduces.Additionally, more degenerate configurations were discovered for the V 5 and V 6 clusters because the number of initial non-equivalent SOD-identified configurations grows as the number of vacancies increases.To calculate the total molar fraction of vacancies at various temperatures, all the configurations shown in Fig. 1 were considered along with their degeneracy.
The calculated molar fractions of the examined clusters are shown in Fig. 2(a).Molar fractions are calculated for temperatures ranging from 300 to 1400 K but Fig. 2(a) only shows molar fractions for temperatures higher than 500 K.Even at temperatures close to the melting point of Cu, one can see that V 1 dominates, followed by V 2 , which has a molar fraction at least 3 orders of magnitude lower compared to monovacancies.Alhough a number of studies discuss the role of di-vacancies in the curvature of the Arrhenius plot in Cu, 65,66 based on these results, the molar fraction of di-vacancies under equilibrium conditions will be negligible.This result is in good agreement with the di-vacancy molar fractions computed by Nordlund and Averback 67 using EAM and supports their conclusion that divacancies, unlike other defects like self-interstitial atoms, are not expected to contribute to the curvature of the Arrhenius plots in Cu.These findings suggest that, under equilibrium conditions, bigger vacancy clusters are unlikely to form in perfect fcc Cu.
The molar fractions of monovacancies at varying temperatures calculated using EAM and DFT are compared in Fig. 2(b).Experimental vacancy molar fractions by Hehenkamp et al. 25 and Simmons and Balluffi 54 are also provided for reference.The DFT values are reasonably close to the experimental ones, especially the ones by Simmons and Balluffi.Hehenkamp et al. attributed the discrepancy between the two experimental methods to the length of the examined samples by Simmons and Balluffi, which resulted in an inhomogeneity in the temperatures of the samples.We also note that the molar fractions found using the potential 49 are approximately one order of magnitude lower than those obtained using DFT.This potential overestimates the formation energy of a monovacancy Fig. 1 Formation energies per vacancy for V 1 to V 6 clusters obtained using the EAM potential. 49Inset images illustrate the lowest energy configurations for each cluster.by around 0.2 eV in comparison to experimental results, causing this discrepancy.

H concentration
First, we estimate the molar fraction of H in perfect Cu lattice at various pressures and temperatures.Unlike previous results from Nazarov et al. for other fcc metals, 22 our calculations show that H in Cu occupies interstitial octahedral rather than tetrahedral sites.This is in good agreement with the results of Ganchenkova et al. 27 The calculated equilibrium molar fractions (d) of H interstitials as functions of temperature and pressure are summarized in Fig. 3.As can be seen, in accordance with Sieverts' law, 68 the molar fraction of H at low pressures is proportional to the square root of the partial pressure (see the inset in Fig. 3(a)).Deviations from the Sieverts' law are seen at 1400 K.In the kbar range (Fig. 3(b)) the Sieverts's law is no lorger obeyed, in agreement with Evans. 69Because the H molar fraction is below 10 À9 at temperatures lower than 600 K, only higher temperatures are considered.The molar fraction of H at the partial pressure range of 0.1 to 100 bar is 3 orders of magnitude lower than that of vacancies at the same temperatures.At temperatures close to the melting point of Cu and external pressure values of 1 bar, the molar fraction of H is about 10 À6 and is too low to play any significant role in the formation of Cu vacancies.It increases significantly between 10 and 50 kbar (Fig. 3(b)) and for temperatures close to the melting point reaches the 10 À5 range.

Influence of H on the concentration of vacancies
We now turn to examining the effect of H on the concentration of Cu vacancies.To achieve this, we used SOD to identify all possible 2H-V 1 and 1H-V 2 complex configurations.The calculated energies for the 2H-V 1 and 1H-V 2 complexes are shown in Fig. 4(a) and (b), respectively.Once identified, all configurations were optimized using DFT.The highest energy configuration was chosen as a reference for all energies in both plots, thus negative energies correspond to more favorable configurations.In addition, the plots include inset images of the 2 lowest energy 2H-V 1 and 1H-V 2 configurations.In the lowest energy configuration of 2H-V 1 , two H atoms are located at the two octahedral sites closest to the vacancy and are separated by  2.6 Å.In the second lowest energy configuration, the two neighboring H-occupied octahedral sites, but are now separated by 3.6 Å.The 0.1 eV energy difference between the two configurations is attributed to the weak interaction between the two interstitial H atoms.
In 1H-V 2 complexes, only the configurations in which the two vacancies are separated by up to a 4NN distance were examined.As seen in Fig. 4(b), the configuration with two vacancies at a distance of 2.6 Å, where the H atom occupies an octahedral site at a V-H-V coordination, thus forming a 901 angle, was found to be the most favorable.The second lowest energy configuration, in which the two vacancies were at a 3.6 Å distance, was found to be 0.05 eV higher (Fig. 4(b)).As also noted by Nazarov et al., 10 the H-V Cu interaction is very short ranged.In cases where H occupied a second nearest available octahedral site from the vacancy, the gain in energy due to the H-V interaction was found to be negligible.
Under equilibrium conditions, monovacancies dominate over other vacancy clusters with molar fractions of 10 À4 at temperatures beyond 1300 K. On the other hand, the estimated molar fraction of H in non-defective Cu is as low as 10 À6 at the melting point and external partial pressure values of 1 bar.However, the interaction of H with vacancies discussed above will reduce formation energies compared to isolated H atoms and Cu vacancies.To assess the importance of this interaction, we calculated the molar fractions of H in the presence of monovacancies and di-vacancies and the molar fraction of vacancies in the presence of H.In Fig. 5(a), the molar fraction of vacancies without H is used as a reference.One can see that the difference between the plots for vacancies with and without H for partial pressures between 1 and 100 bar is indeed very small, within the 10 À11 to 10 À7 range.At high hydrogen pressures, between 10 and 50 kbar, the increase of the number of vacancies due to the presence of H reaches 10 À5 for temperatures close to the melting point (see Fig. 5(a)) and becomes more significant.
Since the number of thermally created vacancies is much larger than that of incorporated hydrogen at normal pressure and temperature, we can expect that at high pressure an extra amount of hydrogen can be incorporated into Cu via absorption by vacancies.Indeed, the differential molar fraction of H, d(p,T), with and without vacancies for partial pressures 1 and 100 bar in Fig. 5(b) demonstrates a moderate increase of 10 À5 close to the melting point.However, for a higher pressure of 10 and 50 kbar at temperatures close to the melting point the H molar fraction increased from 1.84 Â 10 À4 to 2.28 Â 10 À3 (see Fig. 5(c)).The latter value is in a very good agreement with the experimental 2.5 Â 10 À3 H molar fraction at 50 kbar reported by Fukai et al. 6 and reaches the same magnitude as the molar fraction of Cu vacancies.

Conclusions
We used a grand canonical statistical mechanics analysis to determine the molar fraction of vacancies and H in perfect bulk Cu at various temperatures and external pressure values.The study shows that V 1 monovacancies dominate over other vacancy clusters up to the melting temperature of Cu, followed by V 2 , with molar fractions at least three orders of magnitude  lower than that of the monovacancies.Without considering the presence of Cu vacancies, the molar fraction of H was found to be orders of magnitude lower than that of vacancies at the same temperatures and at hydrogen pressure values of around 1 bar.For hydrogen pressures within the kbar range and temperatures close to the melting point of Cu, H was found to reach and surpass the molar fraction of the equilibrium vacancies due to the interaction with vacancies.Further, the effect of H on the concentration of vacancies was investigated.Given the low molar fraction of H, it appeared to have little impact on the rise in vacancy molar fraction for external pressure lower than 100 bar.Likewise, vacancies were shown to have a negligible effect on the molar fraction of H.The impact of vacancies on the molar fraction of H started to become considerable at hydrogen pressure values higher than 10 kbar and temperatures close to the melting point.
These demonstrate that, even though the presence of H promotes the formation of vacancies, this effect is insignificant in Cu because of the small amount of incorporated H. Therefore, in perfect bulk Cu under equilibrium the presence of H will start to significantly increase the molar fraction of vacancies only at hydrogen pressures exceeding 10 kbar and temperatures close to the melting point.
It is important to stress that some of the commonly used deposition methods, such as the ECD process, take place in the presence of water-based electrolyte solutions.Thus, additional H is expected to be incorporated during the deposition and can become trapped in sinks such as vacancies, dislocations or grain boundaries increasing the H molar fraction and facilitating the formation of superabundant vacancies.Further analysis under non-equilibrium conditions, where defects such as dislocations and grain boundaries are taken into account, will be required in order to fully comprehend the role of H in the SAV formation in Cu.

Fig. 2
Fig.2(a) Molar fractions (x) of V 1 -V 6 vacancy clusters at varying temperatures obtained using EAM.(b) Comparison between the molar fraction obtained through DFT, EAM and experimental equilibrium molar fractions of vacancies from Hehenkamp et al.25 and Simmons and Balluffi.54

Fig. 3
Fig. 3 Molar fraction (d) of interstitial H atoms isotherms in pristine bulk Cu at varying temperatures with and without considering the vibrational contributions of H. Two pressure ranges, (a) 0.1 to 1 bar (inset: d against p 1/2 ) and (b) 10 to 50 kbar, were investigated.

Fig. 4
Fig. 4 Relative energies compared to the highest energy configuration for (a) 2H-V 1 and (b) 1H-V 2 complexes obtained using DFT and SOD.Inset images include the two lowest configurations for the examined complexes.The energies of the lowest energy configurations in both plots are highlighted in red.

Fig. 5
Fig. 5 (a) Additional molar fraction of vacancies (x) generated at the presence of H (at partial pressure p).(b) Differential plots for the molar fractions of H (d) with the presence of vacancies at varying temperatures of 1 and 100 bar partial pressure.(c) Molar fractions of H (d) with and without vacancies at 10 and 50 kbar.