Temperature dependence of the Gibbs energy of vacancy formation of fcc Ni

Quantum-mechanical calculations are used to determine the temperature dependence of the Gibbs energy of vacancy formation in nickel. Existing data reveal a discrepancy between the high-temperature estimates from experiments and low-temperature approximations from density functional theory. Our ﬁnite-temperature calculations—which include the effects of magnetism and fully interacting phonon vibrations—demonstrate that this discrepancy is mostly caused by the previously neglected explicit anharmonic contribution.


I. INTRODUCTION
Vacancies are of fundamental importance to the properties of metals and alloys. For example, they exacerbate the degradation of nickel-based superalloys at elevated temperatures, since vacancy-mediated diffusion can cause (i) vacancy absorption/emission (atomic migration) during climbing of dislocations (dislocation creep) [1], (ii) short-range atomic reordering during the lengthening of microtwins [2][3][4], (iii) imbalance of the diffusive flux of substitutional atoms which causes the formation of pores (Kirkendall porosity) [5], (iv) ionic diffusion within surface protective oxides during environmental attack by corrosive species [6], (v) degradation of thermal barrier and bond coat by interdiffusion with an alloy substrate [7], and (vi) nucleation and coarsening of precipitates such as topologically close-packed (TCP) phases [8].
The key quantity determining the vacancy concentration c is the Gibbs energy of vacancy formation G f according to [9] where g is a geometry factor (e.g., g = 1 for monovacanies and g = 6 for divacancies in fcc), k B the Boltzmann constant, and T the temperature. At present, the magnitude of G f and in particular its variation with temperature is the subject of significant uncertainty, even for pure Ni. For differential dilatometry (DD), which is able to obtain absolute vacancy concentrations and generally referred to as the most accurate method [10], the experimental detection limit for vacancies is 10 −5 . DD is therefore limited to temperatures close to the melting point where vacancy concentrations are high enough (see Fig. 1). Applying indirect techniques, such as residual electrical resistivity ( ρ) or calorimetry (DSC), the experimental detection limit can be decreased to about 10 −6 , widening the experimentally accessible temperature range down to 1000 K for Ni. However, the scatter in the data at lower temperatures is large and a temperature dependence impossible to determine. Enthalpies of vacancy formation extracted from positron annihilation spectroscopy (PAS) likewise show a large scatter of more than 0.3 eV (black solid bar in Fig. 1 and Table I).
On the theory side, numerous density functional theory (DFT) studies have been reported to estimate the monovacancy formation enthalpy in fcc Ni, but they considered only the T = 0 K limit (Table II). The temperature dependence of the ab initio Gibbs energy of vacancy formation was determined only within the low-temperature quasiharmonic approximation [11,12]. Corresponding results [11,12] suggest the entropy of formation to be negative (Table II), in clear contrast to the positive DD data (Table I).
We will show that explicit anharmonicity beyond quasiharmonicity, originating from the breakdown of the inversion symmetry close the vacancy, is mainly responsible for the missing entropy contribution. For that purpose, we will perform finite-temperature calculations-which include the effects of magnetism and phonon vibrations-using an adapted version of the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) method.

A. General free-energy approach
The temperature-and pressure-dependent Gibbs energy of vacancy formation is given by [42] (2) where F vac ( ,T ; N ) is the free-energy surface of the vacancy supercell with volume and N atoms, F bulk (v,T ) the freeenergy surface of the perfect bulk for a single atom with volume 0 300 600 900 1200 1500  [12,19,31,41] and DFT calculation for local density approximation (LDA) at T = 0 K [40]. The gray dotted line indicates the quality of the TU-TILD potential (cf. Sec. II C).
v per atom, P the pressure, and v f = ( − Nv) the volume of vacancy formation.
For both of the free-energy surfaces F vac (V ,T ) and F bulk (V ,T ) (where V denotes here a generic volume dependence), the set of relevant contributions for a magnetic material as Ni reads where E 0 K is the total energy at T = 0 K, F el the electronic free energy of a static lattice, F qh and F ah the quasiharmonic and explicitly anharmonic free-energy contribution from atomic vibrations, F el-vib the adiabatic coupling between electrons and atomic vibrations, and F mag-vib the magnetic contribution to the free energy containing the impact of atomic vibrations. The term E 0 K (V ) is meant to include the total electronic and magnetic energy at T = 0 K, which implies that F el and F mag-vib contain only the explicitly temperaturedependent electronic and magnetic contribution. We are mainly interested in G f (T ) in the experimentally accessible range, i.e., 60%-100% of the melting temperature of Ni (1728 K [43]). This temperature range is far above the Curie temperature of Ni (631 K [44]), meaning that experimentally the paramagnetic state is measured. In principle, this fact would require one to compute the terms entering Eq. (3) using paramagnetic DFT calculations. Such calculations are particularly challenging for Ni because (in contrast to, for example, Fe) paramagnetic spin configurations treated by conventional unconstrained spin DFT generally converge into the nonmagnetic state [45,46]. Inclusion of longitudinal spin fluctuations is required to stabilize finite local magnetic moments [45]. Modeling these fluctuations necessitates constrained spin DFT calculations [47], which increase drastically the computational expenses.

B. Total energy, electronic, and quasiharmonic contribution
We calculated E 0 K (V ) at 12 volumes with lattice parameters ranging from 3.42 to 3.64Å, fully relaxing the atomic coordinates. The volume dependence was parametrized using the Vinet equation of state [53]. Table III shows the convergence parameters and the resulting enthalpies of formation revealing only a small change of −0.02 eV upon increasing the supercell from 2 × 2 × 2 to 3 × 3 × 3, and no significant change when increasing further to a 4 × 4 × 4 supercell. The additional nonmagnetic calculations for the 3 × 3 × 3 supercell reveal only a small impact of magnetism on the enthalpy of vacancy formation (≈0.03 eV).
The electronic free energy F el for the ferromagnetic static lattice was calculated using the self-consistent field (SCF) approach [54] within finite-temperature DFT [55]. We used a mesh of 20 volume-temperature points and convergence parameters as given in Table IV. Based on the mesh, the difference between the vacancy and bulk electronic free energy was fitted using a polynomial in volume and tem- Table IV (last column) shows no difference in the Gibbs energy of formation when using a 2 × 2 × 2 or a 3 × 3 × 3 supercell for the electronic free-energy calculations.
The quasiharmonic free-energy contribution F qh to the vacancy formation energy was obtained by the small displacement method [56,57]. We used displacements of 0.03 Bohr radius (∼0.016Å) and treated all atoms spin polarized. Table V lists the tested convergence parameters as well as the considered volumes. The volume dependence was fitted with a second-order polynomial. The difference in the Gibbs energy of vacancy formation between the 2 × 2 × 2 and 3 × 3 × 3 supercells in F qh is below 0.01 eV at the melting point (Table V,  last column).   TABLE III. Convergence parameters for the T = 0 K calculations and the resulting vacancy formation enthalpy H f at atmospheric pressure, for both ferromagnetic (FM) and nonmagnetic (NM) cases. E cut denotes the plane-wave energy cutoff. The notation for the lattice constants (a start -a end ,δa; a 1 ,a 2 ) means that the mesh is going from a start to a end in steps of δa, with additional lattice constants a 1 and a 2 .  The term F ah was calculated using an adapted version of the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) method [59,60]. We had to adapt the method because of frequent jumps of the vacancy at temperatures above ≈1000 K, which made the use of the quasiharmonic reference in the first thermodynamic integration of the TU-TILD method impracticable. The thermodynamic integration to the harmonic reference was therefore performed only at a lower temperature (800 K) to obtain once the absolute free energy. The free energy at higher temperatures was obtained by integrating the internal energy over the temperature [61]. This was straightforward to do (i.e., long enough simulations, dense enough temperature mesh) because only the TU-TILD potential is required at this stage and no DFT. In particular, we used a large enough supercell of 6 × 6 × 6 (864 atoms). We used an embedded atom method (EAM) potential fitted to energies of atomic configurations resulting from DFT molecular dynamics (MD) simulations for the perfect bulk at 1728 K and several volumes.
The TU-TILD potential reproduces the DFT energies very well; the standard deviation of the energy difference between DFT and TU-TILD is only a few meV/atom even at the highest considered temperature of 1728 K (5 meV/atom in the 2 × 2 × 2 and 2 meV/atom in the 3 × 3 × 3 supercell). The quality of the TU-TILD potential is best reflected by the resulting temperature dependence of the Gibbs energy of vacancy formation (gray dotted line in Fig. 1). Note, however, that since the TU-TILD potential is only used as an optimized reference for the following thermodynamic integration to DFT (see below), its quality affects only the statistical convergence of the MD simulations but not the final result.
The thermodynamic integration from the TU-TILD potential to DFT was performed in a 2 × 2 × 2 and a 3 × 3 × 3 supercell with convergence parameters as given in Table VI. Note that vacancy jumps are unproblematic at this stage. MD simulations were performed with a time step of 5 fs using the Langevin thermostat with a friction parameter of 0.01 fs −1 for several thousands of steps such as to obtain a standard error below ≈0.2 meV/atom. We used a dense set of explicitly calculated volume and temperature points (Table VI) and parametrized the free-energy surface with a polynomial a 00 + a 10 T + a 01 V + a 11 T V + a 20 T 2 + a 02 V 2 with fitting coefficients a ij .
In some of the MD simulations for the vacancy supercell at high temperatures and large volumes, we found spin flips close to the vacancy, changing the energetics of the system. In order to stay within the ferromagnetic regime (see the discussion in Sec. II A), we sorted the corresponding MD parts out and averaged only over ferromagnetic solutions. Future studies are required to investigate this issue in more detail.

D. Coupling of electrons with vibrations
The free energy due to the coupling between electrons and lattice vibrations F el-vib (V ,T ) has been calculated by averaging the electronic free energy (referenced to the electronic internal energy at T = 0 K) of all upsampled atomic configurations in the TU-TILD scheme; this method is detailed by Zhang et al. [54]. TABLE V. Convergence parameters for the quasiharmonic calculations and the resulting Gibbs energy of vacancy formation G f at the melting point (1728 K) and at atmospheric pressure, including the converged T = 0 K and electronic contributions (i.e., only the quasiharmonic contribution varies in this table for a convenient comparison). E cut denotes the plane-wave energy cutoff and "Augmentation" refers to the grid used for the calculation of the augmentation charges (gp/atom = grid points; see Ref. [58] for details). The notation for the lattice constants (a start -a end ,δa; a 1 ,a 2 ) means that the mesh is going from a start to a end in steps of δa, with additional lattice constants a 1 and a 2 .

Supercell
Size

E. Magnetic excitations
As already discussed in Sec. II A, the inclusion of magnetic excitations from ab initio is challenging (see also, e.g., Ref. [62] for a recent overview). Since we are predominantly interested in the high-temperature regime, the magnetic free energy can be modeled by neglecting the magnetic short-range order (SRO) as [63,64] where S mag-vib denotes the magnetic entropy. The two main contributions to S mag-vib are due to transverse and longitudinal spin fluctuations. Transverse spin fluctuations, i.e., excitations of noncollinear spin states, are commonly taken into account via the hightemperature limit as [65][66][67] where M(V ,T ) denotes the local magnetic moment (at specific V and T ). Equation (5) can be derived from the hightemperature limit of the quantum-mechanical spin model, but note that it applies rigorously only for a localized magnetic system. For Ni, the application of localized spin models has previously provided good agreement with experiment for the magnetization shape and for thermodynamic properties such as, e.g., specific heat capacities [68]. Some recently proposed alternative entropy equations, S mag−vib = −3k B ln (M) and −k B ln (M) [69,70], have been so far evaluated only for very few selected systems and their applicability for itinerant magnetic systems has not been rigorously proven yet. Longitudinal spin fluctuations can cause variations in the magnitude of M(V ,T ). This can be due to explicit temperaturedependent spin excitations such as single-particle spin-flip transitions between bands of opposite spin (Stoner excitations) as well as indirectly via the coupling of M(V ,T ) to the overall magnetic state, the atomic motion, or the volume V . The impact of atomic motion on M(V ,T ) can be significant, as we showed previously for Fe [71]. It can be taken into account via ab initio MD simulations. The impact of the volume can be straightforwardly accounted for by coupling Eq. (5) to the derived thermal expansion. The impact of the magnetic state on M(V ,T ) is very challenging to capture because it requires constrained spin DFT (unconstrained calculations converge into a nonmagnetic solution; see, e.g., Ref. [46]) which need to be further coupled to MD simulations.
In the present work, we have employed Eq. (5) to capture transverse magnetic excitations. The magnetic moment M was obtained from our ab initio MD simulations by averaging over all N upsampled atomic configurations at specific V and T and in the ferromagnetic state (as discussed in Sec. II A), With this procedure the impact of longitudinal spin fluctuations due to the variation of V , the atomic motion, as well as the electronic temperature is taken into account. Single spin-flip transitions and the impact of the magnetic state on the local magnetic moments are neglected.

III. RESULTS AND DISCUSSION
A. Influence of excitation mechanisms on the vacancy formation free energy Figure 1 shows the calculated temperature-dependent monovacancy formation Gibbs energy and equilibrium vacancy concentration in fcc Ni, contributed by the various excitation mechanisms, in comparison with experimental measurements. The electronic contribution (pink line) gives a formation energy larger than the T = 0 K calculation (red dashed), while quasiharmonicity (blue line) shifts down the curve to slightly lower values. However, the shift is by far not enough to reproduce the experimental data and a discrepancy of almost 0.3 eV remains at the melting temperature. Anharmonicity (green line and shading) is the dominant contribution, shifting down the curve by about 0.2 eV and leading to a strong non-Arrhenius behavior. The contributions from the electron-phonon coupling (purple line) and from the magnetic excitations (red line) further decrease the formation energy at high temperatures.
The final theory curve including all excitation mechanisms (red line) is in reasonable agreement with experimental data. The remaining discrepancy may stem from the error associated with both experiments [30,31] and modeling; the strong interactions between vacancies and tracer elements such as oxygen [72,73], hydrogen [74,75], nitrogen [72], and carbon [76] have been rationalized recently. In terms of the modeling, it should be emphasized that the calculation of the formation energy of a point defect puts extreme demands on numerical accuracy and precision. Although we have tried to keep numerical errors well below 1 meV/atom in our calculations, even a tiny error of 0.1 meV/atom scales up to a 10 meV/defect error in the formation energy for a 3 × 3 × 3 supercell. Further, we stress again the challenge in capturing the full magnetic contribution discussed above.  Figure 2 shows the calculated vacancy formation entropy, i.e., the negative derivative of the Gibbs energy of formation, in comparison with differential dilatometry measurements. The quasiharmonic approximation (blue line) gives only a very small entropy of formation. Considering all excitation terms listed in Eq. (3), particularly the anharmonic contribution, gives a much larger entropy of formation in consistency with experimental measurements: Strong anharmonicity is suggested here. As discussed by Glensk et al. [42] for Al and Cu, this is largely due to the softening of the effective potential towards the vacancy, which locally destroys the inversion symmetry and introduces sizable odd-order anharmonic contributions, although no linear temperature dependency-which is suggested in the local Grüneisen theory (LGT) [42]-is seen in the present case.
To highlight the softening of the effective potential, Fig. 3(a) shows an example of the trajectories of the 12 nearest neighbors (1NNs) of a vacancy (mapped to the first quadrant of the xy plane with coordinates of X NN V ,T ;i ,Y NN V ,T ;i ) based on quasiharmonic (black) and fully anharmonic (red) simulations at  1600 K. In the anharmonic data, one can clearly see a high possibility for the neighbors to move towards the vacancy. Figure 3(b) shows the corresponding distribution function of the distance between 1NNs and the vacancy where δ(x) equals 1 for x = 0 and otherwise 0, and the resulting effective potentials according to One can clearly see an anisotropy and softening towards the vacancy which cannot-by definition-be reproduced within the quasiharmonic approximation. These results are in nice consistency with previous findings for nonmagnetic Al and Cu [42].
To rationalize the impact of the magnetic contribution at high temperatures (Fig. 2), we have studied the averaged magnetic moments in our MD simulations. Overall, we observe a similar behavior for the vacancy and perfect bulk supercell (Fig. 4): The magnetic moment decreases with increasing temperature, and increases with increasing volume. Both dependencies are small with changes of about 0.05μ B and tend to compensate each other if the effect of thermal expansion is considered. The difference between the vacancy and perfect bulk supercell is an order of magnitude smaller and becomes only apparent at the highest temperatures and volumes (red shading in Fig. 4). The slightly higher magnetic moment of the vacancy supercell is mainly caused by the 1NNs of the vacancy. Despite being small, this difference in magnetic moment between the vacancy and perfect bulk supercell is responsible for the steep increase in the entropy of formation at high temperatures.

IV. CONCLUSIONS
The Gibbs energy of vacancy formation in fcc Ni has been calculated using ab initio calculations and by accounting for different excitation mechanisms, including the electronic, quasiharmonic, anharmonic, and magnetic contribution as well as the coupling between electrons/magnetic moments and atomic vibrations. Anharmonicity and its coupling with electrons has been calculated by using an adapted version of the two-stage upsampled thermodynamic integration using Langevin dynamics (TU-TILD) scheme.
In contrast with purely quasiharmonic calculations which suggest a very small entropy of vacancy formation, a much larger entropy of vacancy formation due to locally asymmetric displacement potentials is revealed by the anharmonic calculations carried out here. The calculated results which account for all excitation mechanisms agree reasonably well with recent differential dilatometry (DD) measurements close to the melting temperature, and reconcile the experimental data with T = 0 K ab initio calculations due to a strong non-Arrhenius behavior.
We expect these results to be valuable for the prediction of high-temperature diffusion and defect evolution in nickelbased alloys and thus for the design of different grades of superalloys.