Thermodynamics of impurity-enhanced vacancy formation in metals

Hydrogen induced vacancy formation in metals and metal alloys has been of great interest during the past couple of decades. The main reason for this phenomenon, often referred to as the superabundant vacancy formation, is the lowering of vacancy formation energy due to the trapping of hydrogen. By means of thermodynamics, we study the equilibrium vacancy formation in fcc metals (Pd, Ni, Co, and Fe) in correlation with the H amounts. The results of this study are compared and found to be in good agreement with experiments. For the accurate description of the total energy of the metal–hydrogen system, we take into account the binding energies of each trapped impurity, the vibrational entropy of defects, and the thermodynamics of divacancy formation. We demonstrate the effect of vacancy formation energy, the hydrogen binding, and the divacancy binding energy on the total equilibrium vacancy concentration. We show that the divacancy fraction gives the major contribution to the total vacancy fraction at high H fractions and cannot be neglected when studying superabundant vacancies. Our results lead to a novel conclusion that at high hydrogen fractions, superabundant vacancy formation takes place regardless of the binding energy between vacancies and hydrogen. We also propose the reason of superabundant vacancy formation mainly in the fcc phase. The equations obtained within this work can be used for any metal–impurity system, if the impurity occupies an interstitial site in the lattice. Published by AIP Publishing. [http://dx.doi.org/10.1063/1.4974530]


I. INTRODUCTION
The interaction of metals with hydrogen impurities has been broadly studied for technological and scientific purposes.Hydrogen (H) trapping to defects, such as vacancies, voids, dislocations, and grain boundaries, is one of the long standing problems in materials science.Whether the trapping is desired (H storage 1 ) or unwanted (tritium retention in the fusion reactor wall 2-4 and H embrittlement 5,6 ), the presence of H may alter material properties and introduce damage.Near the melting point of metals, the equilibrium vacancy fractions range from about 10 4 to 10 3 while due to the presence of large amounts of H this fraction can increase up to 0.1-0.3.This phenomenon, called superabundant vacancy (SAV) formation, has been of great interest since the experimental observation of H induced vacancy formation in nickel (Ni) and palladium (Pd) by Fukai et al. in 1993. 7,8][11][12][13][14][15][16][17] The origin of SAV formation lies in the lowering of the system's total energy.The main mechanism of impurity assisted vacancy formation is now commonly believed to be the vacancy formation energy decrease due to the impurity trapping.
Some admirable efforts to study SAV formation theoretically have been made previously.Kato et al. 18 studied SAV formation in tungsten (W) using a thermodynamics approach.The authors assumed a maximum of 6 trapped hydrogen (H) atoms in the monovacancy being bound with a constant binding energy.A similar study on H induced vacancy formation by Ohsawa et al. 19 showed a possibility for more than 6 H bound to the monovacancy and used a different binding energy for each trapped H into the monovacancy.The results of both studies were in good agreement, indicating a very low probability for W monovacancy to trap more than 6 H atoms at temperatures above 300 K.
Combining the first principles calculations with the thermodynamics approach, Nazarov et al. 20 have presented an extensive study of H interaction with vacancies in fcc iron (Fe).The authors found that a Fe monovacancy can accommodate up to 6 H atoms and showed that under rich H conditions the vacancy formation rate increases drastically.
None of the studies mentioned above include the divacancy formation thermodynamics.However, at high temperatures, where the monovacancy concentration is large, divacancies will be formed and give non-negligible contribution to the total vacancy concentration.Nevertheless, the commonly used assumption to calculate the divacancy concentration as a square of monovacancies 21 is largely overestimating the total vacancy concentration when interstitial H is present and, therefore, should not be used when studying SAV formation.In this article, we present a more accurate way to study SAVs.By the means of thermodynamics combined with already established material characteristics, we calculate the formation of each vacancy/vacancy complex taking into account the binding energies for multiple H to vacancies.In our model, we consider the defect vibrational entropies.In the previous studies, [18][19][20] its omission can be justified as the similarity of the compared structures (for example, defect free and containing only monovacancies), giving a very small entropy change.However, in some cases, a) laura.bukonte@helsinki.fi0021-8979/2017/121(4)/045102/11/$30.00 Published by AIP Publishing.121, 045102-1 the vibrational entropy effect has to be taken into account.This is crucial at high temperatures T, as the entropy term (TS vib ) becomes increasingly important.
We also discuss the atomic structure impact on SAV formation and the reason why SAVs are experimentally observed mainly in the fcc phase.
The results from the thermodynamics calculations are compared with the experimental study of Fukai et al., 22,23 where SAVs are found to be created in Co-H, Ni-H, Pd-H, and Fe-H systems.

II. THEORETICAL BACKGROUND: DEFECT EQUILIBRIUM
The vacancy equilibrium fraction ½V in a pure solid can be expressed as where N vac is the number of vacancies, N 0 -the number of host atoms, DH-the change in enthalpy, DS vib -the change in vibrational entropy, k B -the Boltzmann constant, and Ttemperature.
According to thermodynamics, [24][25][26] the change of the Gibbs energy DG at temperature T is associated with the change in enthalpy DH, and the change in the configurational DS conf and the vibrational DS vib entropy of the crystal as The equilibrium of a system is found when the Gibbs free energy reaches its minimum, i.e., the derivative of DG with respect to the changing variable x is zero To continue, we will assess each parameter presented in Equation (2) separately with regard to monovacancies and divacancies.

A. Defect enthalpies
To a first approximation, the total change in enthalpy DH is proportional to the number of monovacancies N v1 and divacancies N v 2 formed where E v1 f and E v2 f are the monovacancy and divacancy formation energies, respectively.
where N 0 is the number of atoms and E 0 -the total energy of a perfect crystal, and E v2 f can be obtained as E v 2 b being the divacancy binding energy.
The impurity atoms N imp present in a crystal can be located at interstitial sites (IS), also called solution sites, whereas others are trapped in the vacancies, i.e., trapping sites (TS).If there is an attraction between an impurity atom and a vacancy, vacancies can usually trap several impurity atoms.The binding energy E b for the nth trapped impurity is the energy needed to remove it from the vacancy and place it on an IS and is expressed as where Eðn 1Þ is the energy of the system, where (n 1) impurities are in the vacancy, and E 1 is the energy with one impurity in IS, E(n) is the energy of the system, where n impurities are in the vacancy, and E 0 is the energy of the reference system.
If the binding energy is positive, energy is gained when an additional impurity is trapped in the vacancy.The binding energy for each additional impurity usually depends on the number of impurities already occupying the vacancy.The cumulative binding energy for n impurities to a vacancy can be very generally be expressed by a polynomial function of n E cum b ðnÞ ¼ where a, b, and c are fitting constants (shown later in chapter IV).
To find the enthalpy change when impurity atoms are trapped in the vacancies DH imp , we define the number of trapped impurities as where n 1 is the number of trapped impurities in the monovacancy and n 2 -the number of trapped impurities in the divacancy.The enthalpy change of the system when N imp tr impurity atoms are trapped in N v 1 and N v 2 where A, B, C and E, F, G are the fitting parameters of H cumulative binding energy to a monovacancy and a divacancy, respectively.Then, the total enthalpy change of the system with N v1 þ N v2 vacancies and trapped impurity atoms N imp tr becomes The enthalpy of the vacancy formation is decreased due to the energy gained when impurity atoms from the IS are trapped in the vacancies.

B. Defect entropies
The total entropy term DS consists of the vibrational entropy DS vib and the configurational entropy DS conf .The vibrational entropy describes random vibrational motion in a defected crystal.The configurational entropy, on the other hand, has a probabilistic nature and is associated with the defect distribution in the lattice.

Configurational entropy
The configurational entropy is the main reason for defect presence in a crystal lattice.Having a crystal with no intrinsic defects at equilibrium is highly improbable, whereas, for example, vacancies give rise to many lattice configurations with different energy states.The increasing amount of defects results in higher configurational entropy and of the lowering of the Gibbs free energy of the system.The configurational entropy of the system depends on the thermodynamic probability, X j , that describes the permutation of defect distribution in the crystal lattice: (1) the number of possible ways to distribute divacancies among N latt lattice sites ( where z is the configurational factor (for fcc z ¼ 12, and for bcc z ¼ 8).
(2) The number of possible ways isolated monovacancies can reside in residual N latt -2N v2 sites (3) The number of ways to distribute n 1 N v1 trapped impurities among total number of trapping sites n max 1 N v1 , where n max 1 is the maximum number of impurities one monovacancy can accommodate (4) The number of ways to distribute n 2 N v2 trapped impurities among total number of trapping sites n max 2 N v2 , where n max 2 is the maximum number of impurities that can be trapped by one divacancy (5) The number of ways to distribute residual impurity atoms among all interstitial sites (m N 0 ) in the crystal lattice, where m is the number of IS per host atom As a result, the total thermodynamics probability to distribute defects in a crystal lattice is presented as Further, according to the Boltzmann's definition of entropy, we obtain for the total configurational entropy of the system as follows: For convenience, we assign new variables TS being the number of impurity atoms trapped in vacancies (TS) and N imp IS -residual impurities located at IS.Using Stirling's formula: logðA!Þ A logðAÞ A, and substituting we get 045102-3 Bukonte, Ahlgren, and Heinola J. Appl.Phys.121, 045102 (2017)

Vibrational entropy
The vibrational (or thermal) entropy S vib is associated with lattice vibrations.Since the lattice atoms around the vacancy are less bound, every vacancy slightly contributes to the total vibrational entropy.Additionally, interstitial or trapping site occupancy of impurity atoms can considerably change the vibrational entropy value.The vibrational entropy term is influential, especially, at high temperatures when the entropy term TDS vib is enhanced, as in Eq. ( 2), and increasingly contributes to the free energy of defect formation.The change of the total vibrational entropy DS vib when N v 1 þ N v 2 vacancies are formed in the system is where DS v 1 vib and DS v 2 vib stand for the change in the vibrational entropy of the crystal due to the formation of a monovacancy and divacancy, respectively, DS imp1 vib and DS imp2 vib are the differences in the vibrational entropy between IS and TS occupancy of impurity atoms, former being associated with monovacancy and latter -with divacancy.

III. GIBBS FREE ENERGY MINIMIZATION
The thermodynamic equilibrium of the system is achieved by minimizing the Gibbs free energy, i.e., by calculating the derivative of the Gibbs free energy with respect to the changing variables in the system, i.e., N v 1 ; N v 2 , n 1 and n 2 .
(1) Derivative with respect to the number of divacancies (2) Derivative with respect to the number of monovacancies (3) Derivative with respect to the number of trapped impurity atoms per monovacancy n 1 (4) Derivative with respect to the number of trapped impurity atoms per divacancy n 2 For convenience, we express the system parameters as fractions Fraction of impurities : Fraction of impurities in IS : where N 0 is the number of host atoms.Following the equilibrium condition and minimizing the Gibbs free energy of the system, we finally obtain the following system of equations: where , and b ¼ ln . This system of equations is numerically solved by Nelder-Mead algorithm 29 to obtain the four unknowns: the monovacancy fraction ½V 1 , the divacancy fraction ½V 2 , the number of trapped impurity atoms per monovacancy n 1 , and the number of trapped impurity atoms per divacancy n 2 at any given temperature T, where ½I is the total impurity fraction, m -interstitial sites per host atom, n max

IV. RESULTS AND DISCUSSION
Employing obtained Eqs. ( 33)-( 36), we calculate the total vacancy fraction as a function of H fraction in four metals: Pd, Ni, Co, and Fe.
Theoretical thermodynamics calculation results are compared with the experimental data by Fukai and co-workers, where authors using the lattice contraction measurements for Ni, Co, and Pd and X-ray diffraction measurements for Fe observed SAV formation. 22,23Vacancy formation energies, vacancy formation entropies, and H binding energies to the monovacancy for each metal-hydrogen system are found from the literature and used as the calculation parameters in our theoretical thermodynamics model.The H fractions and temperature of the system are taken from the corresponding experiments.
The equilibrium vacancy fraction is determined mainly by the vacancy formation energy.However, the effective vacancy formation energy of the material can be lowered due to the H trapping, and it is commonly believed to be the main reason of SAV formation.1][32] We assume that FIG. 1. Cumulative H binding energies, E cum b , to a metal monovacancy and divacancy.Monovacancies accommodate a maximum of 6 and divacancies a maximum of 10 H atoms.The H binding energies to metal monovacancies are taken from the literature and the choice of binding parameters is justified in Sections IV A 1-IV A 4. The binding energy for low H occupancy in the divacancy is assumed to be 0.16 eV higher than for monovacancies.The lines are the fits to the cumulative binding energy obtained using Eq. ( 8).A, B, C and E, F, G for monovacancy and divacancy, respectively, in Eq. (10).a Me stands for a hypothetical metal with bcc and fcc structures, for which the formation of vacancies is compared at the bcc and fcc phase (see Figure 2).

045102-5
Bukonte, Ahlgren, and Heinola a divacancy has one available site less per vacancy due to the proximity of two monovacancies resulting in 10 octahedral sites to accommodate H.The binding energies of H to a single vacancy in many metals vary from 0.1 to 1.2 eV. 33Theoretical calculations state that the binding energies of the H atom to the monovacancy are largest for single and double occupancies, for higher occupancies reducing by 0.2 eV in metals such as Ni and Fe (as a result of the repulsion between H atoms trapped in the vacancy), but remaining about approximately the same for Pd and Nb. 34An extensive study of H interactions with fcc metals by Nazarow et al. shows that H is stronger bound to the divacancy than to the monovacancy with the difference of 0.05-0.2eV. 35In Figure 1, we present the cumulative binding energy of H in monovacancies and divacancies and show that Eq. ( 8) can accurately describe the cumulative H binding to monovacancies.The binding energy parameters A, B, C and E, F, G, in Eq. ( 10), for monovacancy and divacancy, respectively, are listed in Table I.
The vibrational properties of the material are nonnegligible and become increasingly more important at high temperatures and high defect concentrations.However, the data found from the literature on defect entropies are limited.Vacancy formation entropies are reported to range from 1.80-3.0k B in most metals. 36,37Vibrational entropies of H occupying interstitial and trapped states for chosen metals are not available in the literature to our knowledge.However, our DFT calculations of vibrational properties of H in tungsten 38 suggest that the H atom has a larger vibrational frequency in the solute state at the tetrahedral lattice position than in the trapped state -octahedral site, leading to the negative change of entropy for H.We predict a similar trend also for other metals.
The metal (fcc: Pd, Co, Ni, and Fe) properties found in the literature are reviewed thoroughly and used as calculation parameters in our theoretical thermodynamics model.The summary of the utilised values is presented in Table II.
It is important to point out that SAV formation is experimentally usually observed in the fcc phase.The H configurational entropy factor plays an important role in this phenomenon: depending on the material, in an fcc structure, H can occupy one octahedral interstitial site or two tetrahedral interstitial sites per host atom, Eq. ( 16), (m ¼ 1 or m ¼ 2), whereas in a bcc crystal, H can be distributed over six tetrahedral interstitial sites or three octahedral interstitial sites per host atom (m ¼ 6 or m ¼ 3).In order to maximize the entropy and minimize the total energy of the system, vacancy formed in an fcc material would give rise to more configurations to distribute H than in the bcc phase.This is one of the reasons why larger fractions of H are needed for SAVs to be measurable in most bcc metals.The strength of the H configurational entropy contribution depends profoundly on the material.To show the effect of the material structure on SAV formation, in Figure 2, we present a The parameters in the last row describe a hypothetical metal (Me), for which the formation of vacancies is compared at the bcc and the fcc phase (see Figure 2).hypothetical metal (the parameters used in the thermodynamics calculations listed in Tables I and II) with two different structures-bcc and fcc.The hypothetical metal is assumed to accommodate H in six tetrahedral interstitial positions (m ¼ 6) in the bcc phase and one octahedral interstitial position, when having the fcc phase (m ¼ 1).We can see that a metal having the fcc structure has about two decades larger vacancy fraction than a metal with the bcc structure, with the same energy and entropy parameters.We will further discuss the importance of the H configurational entropy in Section V.

Ni
The superabundant vacancy formation is studied in the Ni-H system, having an fcc structure, and compared with the experimental data from the lattice contraction measurements. 22The reported monovacancy formation energy for Ni varies from 1.37 to 1.81 eV, [42][43][44] and the divacancy formation energy ranges from 2.4 to 3.0 eV. 41,45The H binding energy to the Ni monovacancy (E H b ) is reported to be 0.54 eV. 46In our model, we use a value of 0.54 eV for low occupancies (1st and 2nd trapped H) and 0.34 eV for the higher occupancies.9][50] In the thermodynamics calculations, we test different energies and entropy parameters (see Table II) and demonstrate their impact on the total vacancy fraction in Figure 3.We show, Figure 3(a), that the main contribution to the total vacancy fraction comes from divacancies.This is expected due to the low formation energy of the divacancy in Ni (2.42 eV 39 ).In Figure 3(b), we present the thermodynamics calculation results in comparison with the experimental data using different formation energies of monovacancy and divacancy (E v 1 f and E v 2 f , corresponding to the divacancy binding energy of 0.32 eV according to Eq. ( 6).The effect of H binding energy to the Ni vacancy is shown in Figure 3(c).The first value of E H b , in the legend of Figure 3(c), represents the H binding energy for low occupancies (1-2 H in vacancy) and the second-for higher occupancies.The binding energy for FIG. 3. Experimentally obtained and theoretically calculated vacancy fraction as a function of H fraction for the Ni-H system, the comparison between various parameters and effects on the total vacancy fraction.(a) The effect of monovacancy and divacancy fractions on the total vacancy fraction.A good agreement between theory and experiments is obtained with parameters from Refs.42 and 39 (E v1 f ¼ 1.37 eV, E v2 f ¼ 2.42 eV); (b) The effect of E v1 f and E v2 f on the total vacancy fraction, using divacancy binding energy 0.32 eV; (c) The effect of H binding energy E H b to the vacancy.The first value of E H b represents the binding energy for low occupancies (1-2 H in vacancy) and the second for higher occupancies.The binding energy of H for low occupancies in the divacancy is 0.16 eV higher than for monovacancies; (d) The effect of E v2 b on the total vacancy fraction.The monovacancy formation energy used here is 1.37 eV and the divacancy formation energy is obtained from Eq. ( 6), with E v2 b ¼ 0.03, 0.32, and 0.44 eV.

045102-7
Bukonte, Ahlgren, and Heinola low occupancies in the divacancy is assumed to be 0.16 eV higher than for monovacancies.The effect of divacancy binding energies (E v 2 b ) is presented in Figure 3(d).We also assess the effect of H vibrational entropy change on the total vacancy fraction using DS Imp vib values from 0.5 k B to 0 k B ; however, the impact was insignificant in this case.To summarise-the best agreement between theory and experiments for the Ni-H system is found with vacancy formation energy 39 and H binding energy (E H b ) to the Ni monovacancy 0.54 eV 46 for first two H in monovacancy and 0.34 for higher occupancies.

Co
The results of superabundant vacancy formation for the Co-H system are presented in Figure 4, where the experimental values are obtained from the lattice contraction measurements. 22Reported vacancy formation energy values for Co are fairly different 1.25 eV, 52 1.34 eV, 51 and 1.91 eV. 53No information in the literature is found on Co divacancies; therefore, we use the divacancy formation energy E v 2 f ¼ 2 E v 1 f , corresponding to zero divacancy binding energy.The H binding energies to Co monovacancy (E H b ) have been deduced from the thermal desorption data to be 0.47 eV for the first two trapped H atoms and 0.32 eV for three to six H atoms in the monovacancy. 54he best agreement with experiments is obtained in the case of the Co-H system at the smallest vacancy formation energy of 1.25 eV.

Pd
The experimentally measured and theoretically calculated equilibrium vacancy fraction for Pd-H is presented in Figure 5, where we show a good agreement with the experimental data at high H fractions.7][58] The attraction between two vacancies is known to be weak about 0.01 eV. 59Note that even if the binding energy of divacancies is small for Pd, the divacancy fraction at H fractions above 0.5 dominates over the monovacancy fraction.

fcc Fe
The thermodynamic calculation results and the data attained from lattice contraction measurements for the Fe-H system 23 are presented in Figure 6.The vacancy formation energy in fcc iron obtained from the positron annihilation experiment is 1.7 eV; 60 the DFT values are reported to be about 1.8 eV, 20 and we adopt the value of 1.7 eV in the thermodynamics calculations.The H binding energies to monovacancy are found from the literature for the bcc iron: 0.61, 0.7, 0.47, 0.35, ad 0.48 eV (for 1-5 H in the monovacancy) 61 (note: the H binding energy utilised in this work is for the bcc iron.The first principles study of fcc Fe and H interactions by Nazarov et al. 20 gives much smaller H binding energies: 0.28, 0.31, 0.19, 0.23, 0.27, and 0.27 eV (for 1-6 H in the monovacancy) for non-magnetic fcc Fe and even smaller having an anti-ferromagnetic double layer structure: 0.16, 0.19, 0.16, 0.15, 0.12, and 0.14 eV.These binding energies combined with our thermodynamics model resulted in a poor agreement with experiments.We believe that the firstprinciples calculations do not include the zero-point energies that are important to introduce in thermodynamics calculations for the realistic description of the studied system.)In fcc iron, the divacancy binding energy is known to be about 0 eV. 62n Figure 6, we obtain a good agreement with experimental data for the Fe-H system with the vacancy formation energy 1.7 eV, 60 divacancy formation energy 3.4 eV, and the H binding energies from Ref. 61.

V. COMPARISON BETWEEN DIFFERENT THERMODYNAMICS APPROACHES
The simplest way to calculate the equilibrium vacancy fraction is considering only monovacancies.Following the procedure described in Section II and neglecting divacancies, we obtain a much simpler system of equations: where M ¼ ln ½I IS m½IIS and a ¼ ln When the monovacancy fraction is large, e.g., at high temperatures, vacancy clustering becomes very important.The equilibrium divacancy fraction is often calculated as the square of the monovacancies 21 where z is the configurational factor that accounts for possible orientations of the divacancy.DS v 2 is the divacancy binding entropy for simplicity we assume e ðDS v 2 =kBÞ to be unity.The divacancy binding energy can be calculated from Eq. ( 6).The total vacancy fraction is then obtained by adding the divacancy contribution However, this approach is not sufficient to describe the SAV formation, especially, at high H fractions.In such case, both the monovacancy fraction and the divacancy fraction should be calculated taking into account H energy and entropy parameters, as shown in Section II.In Figure 7, we compare both approaches and results show that at high H fractions, vacancy fraction is overestimated if the divacancy fraction is assumed to be proportional to the square of monovacancy fraction as in Eq. (39).Moreover, the total vacancy fraction becomes increasingly overestimated with larger divacancy binding energies.
In Figure 8, we present the comparison of equilibrium vacancy fraction in the Ni-H system calculated using the FIG. 7. Comparison between theoretically calculated vacancy fraction as a function of H fraction for the Ni-H system using the thermodynamics model (see Section II) and a method described in Section V, where the divacancy fraction is calculated as a square of monovacancies.FIG. 8. Comparison between theoretically calculated vacancy fraction as a function of H fraction for the Ni-H system using the thermodynamics model of monovacancies and divacancies (see Section II, Eqs. ( 33)-( 36)) and the thermodynamics of only monovacancies Eqs. ( 37) and (38).In the inset, we show that the thermodynamics for only monovacancies is not sufficient enough to accurately describe the total vacancy fraction.Note the logarithmic scale.thermodynamics model of monovacancies and divacancies (see Section II, Eqs. ( 33)-( 36)) and the thermodynamics of only monovacancies Eqs. ( 37) and (38).A very good agreement is found for low H fractions below 0.001, however at larger H fractions, where the divacancy fraction becomes increasingly dominant over the monovacancy fraction, a deviation is observed between both approaches.This stems from the higher H binding energy to the divacancy.To prove this hypothesis, we demonstrate the case where H does not bind to vacancies, i.e., the fitting constants A, B, C and E, F, G to the cumulative binding curves of H are set to zero.The temperature of the system is chosen to be the melting point of pure Ni.In Figure 9, as anticipated, the monovacancy contribution to the total equilibrium vacancy fraction dominates over the divacancy contribution.Interestingly, even if the H binding energy to vacancies is zero, SAV formation is observed at H fractions above 10%.This is due to the H configurational entropy: when a vacancy is formed, more possible sites to accommodate H become available, increasing the entropy, and therefore minimizing the Gibbs free energy of the system.Vacancy formation at supersaturation of H is especially favourable in the fcc phase (as discussed previously in Section IV), due to more occupational sites available in the vacancy (six trapping sites in the monovacancy and 10 trapping sites in the divacancy) than in the interstitial position (one interstitial site per host atom in the case of Ni).
To conclude our findings: the thermodynamics of only monovacancies is capable to accurately describe the total equilibrium vacancy fraction below the H fraction 10 3 .Above this H fraction (if H is stronger bound to the divacancy than to the monovacancy), the divacancy contribution becomes dominant and the thermodynamics of both monovacancies and divacancies is a more sufficient way to calculate the total equilibrium vacancy fraction.At last, the approach where monovacancy fraction is obtained using the monovacancy thermodynamics and the divacancy fraction calculated from the square of monovacancies is not suitable for studying SAVs at large H fractions.The SAV formation occurs at large H fractions regardless of the binding energy between the H and a vacancy.The ultimate reason for SAV formation at supersaturation of H is the H configurational entropy, which makes SAV formation more favourable in the fcc phase compared to the bcc phase.

VI. CONCLUSIONS
Our general thermodynamics model quantitatively demonstrates the vacancy formation as a function of temperature and the impurity concentration.For the correct description of the free energy of the system, we take into account the binding energies of each trapped impurity, the vibrational entropy of defects, and the thermodynamics of divacancy formation.Our thermodynamics model shows that vacancies are formed in metals due to the presence of hydrogen impurities.The main reason for superabundant vacancy formation is the lowering of the total energy of the system due to hydrogen trapping.We demonstrate the effect of vacancy formation energy, the hydrogen binding, and the divacancy binding energy on the total equilibrium vacancy concentration, the hydrogen binding energy to the vacancies, being the most influential energy factor.We show that the divacancy fraction gives the major contribution to the total vacancy fraction at high H fractions and cannot be neglected when studying superabundant vacancies.Our study leads to an important conclusion that H configurational entropy is one of the reasons responsible for SAV formation observed most often in an fcc phase.Depending on the host material, in an fcc structure, H can occupy one octahedral or two tetrahedral interstitial sites per host atom, whereas in a bcc crystal H can be distributed over six tetrahedral or three octahedral interstitial positions positions per host atom.As a consequence, the formation of a vacancy in a fcc metal gives rise to more configurations to distribute H than in the bcc metal.Due to the increase of the configurational entropy when H is distributed among the available sites in vacancies, the vacancy formation is enhanced even if the binding energy between H and vacancy is zero.Hence, superabundant vacancies are always present at high H fractions.
For the first time, the theoretical thermodynamics calculations are directly compared to the experiments of superabundant vacancy formation and are found to be in good agreement for the fcc metals: Pd, Ni, Fe, and Co.The presented equations are general and can be used for any metal-impurity system, where the impurity occupies an interstitial site in the lattice.33)-( 36)) at temperature 1728 K. Here, the binding parameters A, B, C and E, F, G for monovacancies and divacancies, respectively, are zero.Note that even if the H binding energy to vacancies is zero, the vacancy fraction increases due to the presence of hydrogen.

1 and n max 1 -
the maximum number of impurities per monovacancy and divacancy, A; B; C; E v 1 fthe energy parameters for a monovacancy, E; F; G; E v2 f -the energy parameters associated with a divacancy, and DS v 1 vib ; DS imp1 vib ; DS v 2 vib ; DS imp2 vib -the vibrational entropy parameters for monovacancy and divacancy, respectively.

FIG. 2 .
FIG.2.The equilibrium vacancy fraction as a function of H fraction for a hypothetical metal with fcc and bcc structure, note the logarithmic y scale.The difference between both structures arises due to the H configurational entropy factor, see Section II B.

FIG. 6 .FIG. 4 .FIG. 5 .
FIG.6.Experimentally obtained and theoretically calculated vacancy fraction as a function of H concentration for the Fe-H system.E v1 f ¼ 1.7 eV and E v2 f ¼ 3.4 eV.

FIG. 9 .
FIG.9.The vacancy fraction as a function of H fraction for the Ni-H system using the thermodynamics model (see Section II, Eqs.(33)-(36)) at temperature 1728 K. Here, the binding parameters A, B, C and E, F, G for monovacancies and divacancies, respectively, are zero.Note that even if the H binding energy to vacancies is zero, the vacancy fraction increases due to the presence of hydrogen.

TABLE II .
The parameters utilised in the thermodynamics model: monovacancy formation energy E v1 f , divacancy formation energy E v2 f , monovacancy formation entropy DS v1 f , the change of entropy for impurity (H) occupying solute and trapped state DS Imp vib , and the temperature at which a thermodynamics calculations are performed.