Breakdown of the Arrhenius Law in Describing Vacancy Formation Energies : The Importance of Local Anharmonicity Revealed by Ab initio Thermodynamics

We study the temperature dependence of the Gibbs energy of vacancy formation in Al and Cu from T 1⁄4 0 K up to the melting temperature, fully taking into account anharmonic contributions. Our results show that the formation entropy of vacancies is not constant as often assumed but increases almost linearly with temperature. The resulting highly nonlinear temperature dependence in the Gibbs formation energy naturally explains the differences between positron annihilation spectroscopy and differential dilatometry data and shows that nonlinear thermal corrections are crucial to extrapolate high-temperature experimental data to T 1⁄4 0 K. Employing these corrections—rather than the linear Arrhenius extrapolation that is commonly assumed in analyzing experimental data—revised formation enthalpies are obtained that differ up to 20% from the previously accepted ones. Using the revised experimental formation enthalpies, we show that a large part of the discrepancies between DFT-GGA and unrevised experimental vacancy formation energies disappears. The substantial shift between previously accepted and the newly revised T 1⁄4 0 K formation enthalpies also has severe consequences in benchmarking ab initio methods against experiments, e.g., in deriving corrections that go beyond commonly used LDA and GGA exchangecorrelation functionals such as the AM05 functional.


I. MOTIVATION
Vacancies in a crystal are known to have a strong impact on mechanical strength and ductility, e.g., by enabling material transport, acting as pinning centers for dislocations, or by enabling a dislocation climb.A key quantity to characterize vacancies is their temperature-dependent Gibbs energy of formation, G f ðTÞ, since it provides direct information regarding thermodynamic stability, equilibrium concentration, and solubility.In the dilute limit, the concentration c is related to G f by c ¼ g expð−G f =k B TÞ; (1) with g a geometry factor (e.g., g ¼ 1 for monovacancies and g ¼ 6 for divacancies in fcc) and k B the Boltzmann constant.The experimental determination of c faces severe difficulties: (i) The vacancies must occur in concentrations well above the experimental detection limit, (ii) their detection should not be shadowed by other defects or impurities, and (iii) their concentration must have reached equilibrium.Particularly, conditions (i) and (iii) force experimentalists to go to high temperatures, where concentrations are high and defect kinetics is fast.In practice, measurements are restricted to a temperature range between ≈60% and 100% of the melting point, as indicated in Fig. 1 for Al and Cu by the gray shaded area.
Ab initio calculations, in particular, density-functional theory (DFT), have become the work horse to compute energies of defect formation.The majority of these studies, however, have been performed at T ¼ 0 K, and only recently, the quasiharmonic approximation or empirical potentials were used to estimate finite temperature effects in defect systems [13][14][15].The inclusion of explicit anharmonicity due to phonon-phonon interactions-relevant at high temperatures-was computationally too expensive to be evaluated on a DFT level.
Taking into account the theoretical restriction to low temperatures and the experimental restriction to high temperatures, it becomes evident that a direct and conclusive comparison of experiment and theory has so far been hampered by a large temperature gap.To bridge this gap and to provide a temperature dependence of G f , a common approach is the assumption of an Arrhenius-like behavior, with temperature-independent enthalpy and entropy of formation, H f and S f , that are obtained by a fit to experimental data.As shown in Fig. 1, the limited and scattered experimental data (black symbols) do not allow one to check the accuracy of this assumption.Indeed, Fig. 1 indicates deviations from the simple linear Arrhenius behavior [difference in slopesbetween positronannihilation spectroscopy (PAS) and differential dilatometry (DD) data], which have been extensively and controversially discussed over the last decades [1,[15][16][17][18][19].
We show here by means of highly accurate finite temperature DFT calculations that G f for the two prototype elements Al and Cu has a strong temperature dependence and that the common assumption of a linear Arrhenius extrapolation [Eq.(2)] may give rise to deviations of a few tenths of an eV in the formation enthalpies and an order of magnitude in the entropies.In particular, our results reveal that anharmonic phonon-phonon interactions-efficiently captured by the recently developed upsampled thermodynamic integration using Langevin dynamics (UP-TILD) method-explain the observed deviations of a few tenths of an eV compared to the quasiharmonic approximation.Because G f ðTÞ is found to be strongly nonlinear in temperature, we show that the almost universally accepted linear Arrhenius assumption needs to be replaced by a local Grüneisen theory (LGT) with a formation entropy linear in the temperature.Only the LGT accurately captures the ab initio computed temperature dependence.As will be discussed, these results have important implications for the interpretation of experimental data, and we provide revised T ¼ 0 K extrapolated vacancy formation enthalpies for Al and Cu to guide future studies that rely on the availability of highly accurate experimentally derived T ¼ 0 K data.Using these newly derived T ¼ 0 K data, we show that a large part of the previously reported discrepancies between DFT-GGA and experimental vacancy formation energies disappears.Consequently, previously introduced concepts such as surface corrections [16] or the AM05 functional that aim at correcting DFT errors and that have been justified by benchmarking against experimental T ¼ 0 K extrapolated vacancy formation enthalpies must be revisited.

A. General approach to compute the Gibbs energy of vacancy formation
Key to computing the temperature dependence of G f is the calculation of the bulk and vacancy supercell free energies containing the relevant excitation mechanisms, all being computed by DFT: the T ¼ 0 K energy E 0K and the electronic, quasiharmonic, and anharmonic free energy, F el , F qh , and F ah , respectively.The computation of the first three contributions is standard and described in detail, e.g., in Ref. [20].The computationally most challenging contribution is the anharmonic one, which only recently became accessible on a DFT level [21][22][23].Converging the anharmonic contribution of the defect formation energy G f to a precision similar to the contribution to the bulk free energy F is a completely new challenge since the relevant energy is scaled to the defect and not to the atom, as is the Experiments (PAS ¼ positron annihilation spectroscopy [1], DD ¼ differential dilatometry [1,2]) are limited to a region (gray shaded) close to the melting point, T melt Al=Cu .Extrapolations of available PAS [1,[3][4][5][6][7] and DD data [2,1,3,8-12] to T ¼ 0 K using the common Arrhenius ansatz, G f ðTÞ ¼ H f − TS f , introduce scatter in the reported values (filled/empty black bars mark corresponding intervals).Formation energies computed by common ab initio approximations such as the T ¼ 0 K (dotted line) and the electronic-plus-quasiharmonic (el þ qh; dashed line) approach are shown.The full curve (el þ qh þ ah) includes all free-energy contributions in particular anharmonicity.The error resulting when assuming the Arrhenius extrapolation, Δ Arr , is marked by the orange arrow at T ¼ 0 K. case for bulk properties.The targeted precision therefore needs to be much higher.To estimate the necessary accuracy, let us consider a precision of 1 meV (see Ref. [22] for a discussion on the required precision) to sample the anharmonic contribution of a supercell containing 100 atoms.Our molecular-dynamics runs are performed for the bulk and defect systems on a quasiharmonic reference, as described in Ref. [22].Since both reference systems are of similar quality, the standard deviation σ of the sampled energy differences is very similar in the bulk and defect cases, σ bulk ≈ σ defect .In a bulk calculation with 100 atoms, extensivity holds and the standard error σ n can be calculated as a per-atom quantity.In the case of defect formation energy, we are interested in energy differences of the two supercells (defect minus bulk), and the standard error will therefore be 100 times larger compared to the bulk free-energy calculation.Using the definition of the standard error σ n ≔ σ ffiffi n p , with n being the number of time steps in the corresponding moleculardynamics run, we see a dramatic effect on n defect =n bulk : To obtain the same precision in the defect formation energy as in the bulk free-energy calculation, we need a sampling time that is 10 4 times longer.Since converging the anharmonic bulk free energy to 1 meV/atom is already a difficult task on an ab initio level, it becomes clear that to converge defect formation energies to a similar quality is a considerable challenge.Here, we employ the UP-TILD [22] approach, which is based on a successive coarsening of configuration space, while simultaneously increasing accuracy in the energy calculation.At the final level, only a few hundred fully converged DFT configurations are needed to achieve statistical error bars of 0.1 meV/atom in FðV; TÞ and 10 meV/defect in a corresponding 108 (107) atom cell for G f .Using this approach, we calculate all contributions entering the free energy of the vacancy cell and perfect bulk cell, F vac and F bulk , as a function of volume and temperature.The temperature-and pressure-dependent Gibbs energy of formation is then given by The volume of the defect supercell Ω with N atoms and the volume per atom V of the perfect bulk are self-consistently determined to correspond to a given pressure P (standard atmospheric pressure).The volume of vacancy formation is given by v f ¼ Ω − NV.

B. Technical details
For our calculations, we employ the projector augmented wave method [24] as implemented in VASP [25,26].The exchange-correlation functional is described by the localdensity and generalized-gradient approximations (LDA/ GGA) within the scheme of Ceperley-Alder [27] as parametrized by Perdew and Zunger [28] for LDA and Perdew-Burke-Ernzerhof (PBE) [29] for GGA.Additional calculations were performed at T ¼ 0 K with the Perdew Wang (PW91) [30] parametrization to GGA and with the AM05 functional [31], which is assumed to largely overcome the deficiencies of the LDA and GGA functionals in describing vacancy formation energies.
Extensive convergence tests were conducted for all freeenergy contributions entering Eq. ( 3): E 0K , F el , F qh , and F ah .In general, convergence parameters are optimized such as to guarantee a precision of better than 10 meV/defect in the formation energies G f .This corresponds to converging the bulk free-energy contributions to ≈0.1 meV=atom in a 108-atom cell (3 × 3 × 3 ¼ 3 3 fcc supercell), as described in the previous section.For clarity, we consistently specify meV/atom when referring to free energies while using meV/defect for defect formation energies G f .The number of atoms per supercell given below refers to perfect bulk supercells, while corresponding monovacancy or divacancy cells have one or two atoms less.The supercells are given in units of the conventional fcc unit cell.To minimize errors, all calculations are performed using equal convergence parameters for the bulk and defect calculations.
The T ¼ 0 K contribution E 0K was investigated for Al up to a 4 3 supercell (256 atoms).Differences between a 3 3 (108 atoms) and a 4 3 supercell are below 2 meV/defect.In Cu, we find a difference of 6 meV/defect between a 2 3 (32 atoms) and a 3 3 supercell.We carefully checked the k-point convergence up to 1 million k-points times atom.Both Al and Cu converge at roughly 100,000 k-point times atom to our chosen convergence criterion.We tested planewave cutoffs up to 500 eV, and we find that 300 eV for Al and 400 eV for Cu are sufficient.
For the electronic contribution F el , a well-converged formation energy is achieved by a parametrization of the ðT; VÞ dependence on a grid including seven T steps and four volumes employing 32,000 k-points times atom.The exact details of the parametrization follow the scheme introduced in Ref. [32].By separating out the T ¼ 0 K contribution, an accuracy of better than 2 meV/defect is easily obtained.For Al and Cu, convergence tests were performed in a 2 3 supercell, and we find that the electronic contribution to G f is negligible in both elements.
Phonon calculations for obtaining the quasiharmonic free energyF qh weredonein 2 3 and 3 3 supercells, and ninevolume points were found for both elements to ensure converged vibrational free-energy contributions.We also carefully checked the k-point convergence up to 130,000 (90,000) k-points times atom for Al (Cu) and find converged results already at 55,000 (23,000) k-points times atom.For the Al phonon part, we tested plane-wave cutoffs of 300 eV and 400 eV and found the lower value to be well converged.For Cu, we tested cutoffs up to 500 eVand found a difference of 2 meV/defect compared to a cutoff of 400 eV.
The anharmonic free energy F ah was investigated in a 2 3 (32 atoms) and a 3 3 (108 atoms) supercell and treated with the UP-TILD method [22].The corresponding moleculardynamics simulations used a time step of 10 fs and a friction parameter of 0.01 for Al and 0.03 for Cu for the Langevin dynamics.We defined a convergence criterion for the statistical error to 0.1 meV/atom.This criterion resulted in molecular-dynamics trajectories (after equilibration) of 3,000 steps (30 ps) at the lowest considered temperatures of 250 K and of 15,000 steps (150 ps) at the melting temperature for both elements.For Al, we used a 3 3 kmesh and a plane-wave cutoff of 250 eV as low converged parameters for both cells but checked k-meshes up to 4 3  and cutoffs up to 300 eV explicitly at several selected hightemperature points.For our high converged runs of Al, we used a cutoff of 400 eV and a 6 3 k-mesh in the 2 3 supercell and a 4 3 k-mesh in the 3 3 supercell.For the low converged Cu molecular dynamics, we used a 3 3 k-mesh in the 2 3  supercell and a 2 3 k-mesh in the 3 3 cell.The high converged runs for Cu were performed at a 400-eV cutoff and a 6 3 k-mesh in the 32-atom cell (4 3 k-mesh for the 108 cell).For each temperature and volume, we sampled a mesh of five coupling parameters (0, 0.15, 0.5, 0.85, 1) and used the proposed cotangent fit [32] for parametrization.We used a dense mesh of > 25 volume-temperature points for both elements (five volumes and more than five temperatures) as fitting input for the parametrization based on a renormalized phonon frequency [22], , with a 0 …a 2 fitting coefficients.The resulting renormalized shift due to anharmonicity is used to fit the anharmonic free-energy surface according to the analytical model introduced in Ref. [22].Within this model, ω ah enters a nonlinear function for the free energy, which provides the capability to capture the strong nonlinearity observed in Fig. 1.

III. RESULTS
The resulting Gibbs energy of vacancy formation for Al and Cu is shown for the two commonly employed exchange-correlation functionals LDA and GGA-PBE in Fig. 1 (thick solid blue and orange lines).We consistently observe strong and clearly non-negligible deviations from linearity, particularly when considering the full and experimentally inaccessible temperature window.To quantify the error of the linear extrapolation which is commonly assumed in analyzing experimental data, we perform an Arrhenius fit through the GGA-PBE high-temperature data ("Arrhenius" lines in Fig. 1).The T ¼ 0 K extrapolated value of the energy of vacancy formation is 0.15 eV (≈23%) for Al and 0.22 eV (≈20%) for Cu above the original T ¼ 0 K value.Thus, strong non-Arrhenius behavior clearly results in non-negligible corrections to the T ¼ 0 K energy of vacancy formation.The large nonlinearity also explains the intensively debated difference in the slopes of the PAS and DD data [1,[15][16][17][18][19].The two methods are restricted to distinct temperature ranges, and they therefore probe the slope of the non-Arrhenius curve in different regions.
The possibility of non-Arrhenius behavior has been proposed in previous studies [33], but clarifying the microscopic origin remained impossible and has led to severe controversies between two schools of thought, suggesting either (a) strong finite-temperature contributions [16,17] or (b) significant concentrations of not one, but at least two types of point defects [18,8,34].Most experimental studies assumed, in their analyses, situation (b) with divacancies as a second type of point defects to explain the nonlinearities in G f .Fitting such a model of monovacancies and divacancies to experiment, about 40% of vacancytype defects at T melt have to be assumed to be divacancies for Al and 20% for Cu [34,8].A verification of these fits has been lacking so far since a separation of the total concentration into monovacancies and divacancies, which would allow one to prove or disprove the underlying assumptions, turned out to be experimentally infeasible.
Having our formalism at hand, we are able to clarify whether divacancies may give rise to the non-Arrhenius contributions observed in Fig. 1.For that purpose, we also apply the methodology described in Sec.II for computing the complete thermodynamics of divacancies in both Al and Cu.The divacancy-to-monovacancy concentration ratio at the melting point resulting from our calculations is 0.3% and 0.4% for Al and Cu, respectively, and not 40% and 20% as had to be assumed in the empirical model (b); i.e., they are two orders of magnitude smaller.A critical assumption that had to be made in the empirical model (b) is to invoke for divacancies a much larger entropy of formation than for monovacancies.Figure 2 shows a representative example of the ab initio computed and empirically proposed formation entropies for Cu.As can be seen, the actual entropy of divacancy formation (orange dashed line) is only slightly higher than that of the monovacancy (orange solid line).This result is in gross contrast to a factor of up to 5 higher entropy for the The thinner lines indicate the entropy of formation considering only the electronic and quasiharmonic free-energy contribution (el þ qh).The thick curves include all contributions in particular anharmonicity (el þ qh þ ah).Green lines represent numbers suggested in Ref. [8] for explaining non-Arrhenius behavior within a monovacancy þ divacancy model.Experimental PAS [37] and DD entropies (filled/empty black bars) are derived from the experimental data shown in Fig. 1.
divacancy, which had been postulated to fit the experimentally observed non-Arrhenius behavior (green dashed vs green solid line) [8,35,36].The consequence of these observations is that divacancies [model (b)] can be clearly ruled out as a source of the non-Arrhenius behavior for both elements [16].The strong temperature dependence of the formation energy of the monovacancy therefore remains an exclusive source.Thus, a crucial question to be answered is what type of physical excitation mechanism can lead to such a strong temperature dependence.
For answering this question, the fully ab initio approach described here is ideally suited.Its advantage is that all free-energy contributions are available, allowing one to identify the responsible one(s) for the strong non-Arrhenius behavior.The corresponding analysis reveals the explicit anharmonic contribution, i.e., the one describing excitations beyond the quasiharmonic approximation, as the main source: When it is excluded, we find an almost linear dependence in the Gibbs energy of formation (orange dashed lines labeled "el þ qh" in Fig. 1) and negligible contributions to the entropy of formation (thin orange dashed and solid lines in Fig. 2).The strong impact of explicit anharmonicity is surprising since for fcc crystals it is generally assumed to be small.To understand origin, we have analyzed in detail how temperature affects the distribution function ρ V;T ðx; yÞ of the metal atoms closest to the vacancy: Here, the sum runs over all time steps of a moleculardynamics run at a fixed volume V and temperature T; δðxÞ is equal to 1 for x ¼ 0 and otherwise 0. Further, X NN V;T;i and Y NN V;T;i are the coordinates of all the first nearest neighbors of the vacancy at the ith molecular-dynamics step transformed into the first quadrant of the xy plane by proper point group symmetry operations.Figure 3a shows an example of ρ V;T ðx; yÞ for Cu at the melting temperature, and Fig. 3b shows the Gauss broadened projection onto the [110] direction, i.e., along the line through the vacancy center and the neighboring atom.Using the ab initio computed distribution function, the temperature-dependent effective potential is constructed.Both ρ V;T and v eff V;T show an anisotropy and softening towards the vacancy [Fig.3b], which can be intuitively understood by the fact that bond compression is absent in this direction.As a consequence, the effective potential resembles a Morse potential and leads to a displacement of the time-averaged position towards the vacancy center with increasing temperature [orange diamonds in Fig. 3b].This net inward relaxation with increasing temperature leads to a local expansion of the host matrix at the expense of the vacancy volume.
The large anharmonicity in the effective potential is a direct consequence of destroying the inversion symmetry that an atom has in a perfect fcc crystal: While in the ideal bulk the effective potential will be symmetric and thus effectively cancel third-and higher-odd-order anharmonic contributions, the loss of inversion symmetry of an atom near the vacancy center gives rise to sizable odd-order contributions, as shown in the effective potential in Fig. 3b.The presence of odd, in particular, third-order asymmetric contributions near the vacancy (or any defect destroying inversion symmetry locally) naturally explains the surprisingly large anharmonic effects.
Based on this discussion, the largest anharmonic contributions should be along directions where inversion symmetry is destroyed locally.Indeed, this behavior is found in the distribution shown in Fig. 3a: Odd-order anharmonicity is large towards the vacancy center (along [110]) but absent for directions perpendicular to the [110] direction because of the presence of a mirror symmetry.FIG. 3 (color online).(a) Harmonic (black) and anharmonic (orange) distribution ρ V;T ðx; yÞ according to Eq. ( 5) for Cu at T melt ¼ 1360 K.The vacancy center is placed at (0, 0), and the equilibrium position of the first nearest neighbor is marked by a green cross.The points (black/orange) show the moleculardynamics trajectory of the atom at discrete time steps of 10 fs.The region close to the equilibrium position is densely populated, and thus the individual points are not resolvable on this scale.The harmonic data are obtained from thermodynamic integration runs at zero coupling constant.(b) Distribution function ρ V;T ðdÞ (dashed lines) according to Eq. ( 6), i.e., projection of ρ V;T ðx; yÞ onto the [110] direction indicated in (a), and corresponding effective potential according to Eq. ( 7) (solid lines).The zero line of the distribution function is shifted upwards by the energy k B T melt according to the temperature at which ρ V;T ðdÞ was calculated.The orange diamonds mark the shift of the center of mass of the anharmonic ρ V;T ðdÞ shown and additional anharmonic distributions (not explicitly shown) towards the vacancy at the following temperatures: 250, 450, 800, 1100, 1250 K (related to the energy axis by k B T).For clarity, the latter curve is scaled by a factor of 10 on the d axis.
Since we further find that asymmetries in the distributions of second-and further-neighbor atoms around the vacancy are small, the only relevant collective degree of freedom is the displacement of the nearest-neighbor atoms along a line through the vacancy center.
Being able to relate the anharmonicity to a single collective degree of freedom allows us to make a formal connection to the conventional Grüneisen theory that describes (quasi-)anharmonicity of an ideal bulk system with respect to a collective variable such as the lattice constant or the volume.Thus, anharmonicity of a vacancy (or, more generally, of any defect locally destroying the symmetry) may be regarded as local quasi-anharmonicity, where all physical quantities related to the bulk collective variable are replaced by an analogous variable related to the vacancy neighbor displacement; i.e., thermal expansion relates to the (inward) shift of the neighbors with increasing temperature, the bulk modulus to the second derivative of the effective potential shown in Fig. 3b.In the following, we call this new mechanism local Grüneisen theory (LGT).
The unveiled connection to the Grüneisen theory allows us to understand a further important feature observed in Fig. 2: The entropy of formation is not a temperatureindependent constant-the cornerstone of the Arrhenius approximation and a common assumption in defect models-but the actual dependence is approximately linear for both monovacancies and divacancies, resulting in an approximate T 2 scaling in the Gibbs energy of formation (Fig. 1).Our results also reveal for Al a linear temperature dependence of the monovacancy and divacancy entropies (not included in Fig. 2), suggesting that the origin of this behavior might be generic.Indeed, a linear temperature dependence of the entropy is also observed for the conventional (i.e., related to the volume expansion of bulk systems) Grüneisen model if a harmonic reference entropy (one at a fixed volume) is subtracted from the full entropy (containing the volume dependence).The situation for the LGT of a vacancy, discovered here, is similar: While the entropy of the defect supercell is equivalent to the quasiharmonic entropy in the conventional Grüneisen model, the entropy of the perfect bulk supercell is equivalent to the harmonic entropy in the conventional model.Thus, taking the difference between the entropy of the defect supercell and the entropy of the perfect bulk to obtain the entropy of formation within the LGT corresponds to the difference of the quasiharmonic entropy and a harmonic entropy within the conventional Grüneisen approachwhich is linear.
Empirical models based on a linear temperature dependence of the entropy of formation have in fact been considered earlier [33,38] but failed to receive sufficient attention since they were based on experimental diffusion rather than on concentration data.All reported experimental vacancy formation energies have instead used the conventional Arrhenius dependence with constant entropies of formation.Our present calculations reveal that, at least for the investigated elements, the Arrhenius model needs to be replaced by a model in which the formation entropy has to be expanded up to the first order in temperature S f ðTÞ ≈ S 0 þ S 0 T. For the dominant contribution of monovacancies, the S 0 term is negligible with respect to the S 0 T term (cf.Fig. 2), and we can further approximate S f ðTÞ ≈ S 0 T: This ab initio derived relation corresponds to fitting experimental point-defect formation energies to The new model based on Eq. ( 9) contains the same number of fitting coefficients as a linear Arrhenius model, requiring only H f 0K and the slope S 0 ¼ ∂S f =∂T to be fitted.Table I gives a concise comparison of the conventional linear Arrhenius approximation with the LGT proposed here, which has been inspired by the ab initio computed temperature dependence of the defect formation entropy.
Applying the LGT to experimental data, revised enthalpies of formation at T ¼ 0 K (H f 0K ) are given for Al and Cu in Table II and compared to values obtained by the linear Arrhenius extrapolation of PAS or DD data.For the latter, errors of up to 0.24 eV are observed that are consistent with the magnitude found for our explicitly computed data (Δ Arr in Fig. 1).For Al, the generally accepted value for the vacancy formation enthalpy is the one extracted from PAS data.As shown in Table II, the difference between our revised value and the PAS one is only 0.01 eV and thus unlikely to change any conclusions or benchmarks where

Fitting coefficients
Linear Arrhenius approximation Local Grüneisen theory (LGT) this value has been used.For Cu, the difference between the revised one and the one recommended by Landoldt-Börnstein is 0.22 eV, i.e., an order of magnitude larger, making it mandatory to use the revised value.Future studies also need to address the question of whether the validity of S 0 ≈ 0 that has been assumed in Eqs. ( 8) and ( 9) is specific to these two metals or whether it is of a more general behavior.Finally, we discuss the performance of the exchangecorrelation (xc) functionals employed in this study.Figure 1 clearly shows that the GGA-PBE-based curves lie almost on top of all experimental data sets in absolute value and curvature, while the LDA data predict systematically too-high formation energies (albeit with a similar curvature).This finding is in clear contrast to what has been recently concluded based on calculations without DFTbased anharmonicity and using Arrhenius-extrapolated experimental data as a reference [16,40,41].In these studies, LDA appeared to be the functional that predicts T ¼ 0 K energies of vacancy formation in best agreement to Arrhenius-extrapolated experimental data for Al and Cu, while the two commonly employed GGA functionals (PBE and PW91) predict the energies systematically lower.This behavior is illustrated in Fig. 1 and can also be seen in Table II by comparing the theoretical H f 0K values (LDA, GGA-PBE, GGA-PW91) to Arrhenius extrapolated PAS and DD data.In fact, the apparent underestimation of GGA appears to be a general trend for fcc metals when comparing DFT formation energies to experimentally derived T ¼ 0 K values reported in the literature [40].
The systematic underestimation of vacancy formation energies by GGA was explained by LDA describing the energetics of metal surfaces better than GGA [40,42].Since vacancies may be regarded as an inner surface, it has recently been proposed [16] that for vacancy calculations novel xc functionals should be used that partly correct this deficiency, such as the AM05 xc functional [31].
The present results clearly show that the good agreement of the LDA energies of vacancy formation with experimental data extrapolated to T ¼ 0 K is accidental and related to the neglect of the non-Arrhenius behavior in the temperature window not accessible by experiment.To further verify this conclusion, we performed additional T ¼ 0 K calculations including surface corrections and using the AM05 xc functional.The results are included in Table II and in the insets of Fig. 1.A clear trend is observed for the formation enthalpies (Table II): Both the AM05 xc functional and the surface-corrected PBE calculation give values that are significantly too large when comparing with the revised experimental reference, i.e., the one obtained by using the LGT rather than the Arrhenius-extrapolated data.It is interesting to note that for the unrevised "experimental" T ¼ 0 K enthalpies (see the Landoldt-Börnstein value in Table II), the surface-corrected GGA-PW91 functional (the one originally used in the work of Carling et al. [16]) shows an apparently better performance than the uncorrected LDA or GGA explaining why these xc corrections have been suggested, but it also highlights why it is so important to accurately extrapolate experimental data to T ¼ 0 K. Table II shows that the overall best performance is for the non-surface-corrected GGA-PBE functional.

IV. CONCLUSIONS
Based on an efficient technique to sample anharmonic free-energy contributions for point defects on a DFT level (UP-TILD), we were able to go beyond previous approaches that were limited to a harmonic or quasiharmonic description of the vibrational part of the Gibbs formation energy.Applying this technique to compute formation energies of prototype point defects, the monovacancy and divacancy in fcc-bulk Al and Cu, allowed to cover the entire temperature interval from T ¼ 0 K up to the melting temperature of the bulk systems.A detailed analysis of the computed temperature dependence revealed features in gross contrast to present textbook knowledge, which is based solely on experimental high-temperature concentration data.First, the formation entropy of all studied defects is not independent of temperature as is commonly assumed and as is the basis for the practically universally accepted Arrhenius model, but it shows in leading order a linear temperature dependence.This strong TABLE II.Comparison of experimental formation enthalpies from the compilation in the Landoldt-Börnstein series [39] and of Arrhenius-extrapolated PAS and DD formation enthalpies (averaged values of black filled/empty bars in Fig. 1) with experimental values obtained from the LGT proposed here (Table I) and with ab initio computed formation enthalpies at T ¼ 0 K.For computing the ab initio values, various flavors of exchange-correlation functionals have been employed: LDA, GGA-PBE, GGA-PW91, and AM05.Additionally, surface corrections based on the scheme from Ref. [40] have been used for the PBE and PW91 functionals.The columns labeled Δ exp give the difference from the LGT DD þ PAS value.
-----and unexpected linear dependence was traced back to explicit anharmonicity, i.e., the part of anharmonicity that is not included in the quasiharmonic treatment.Introducing a LGT with a collective variable that represents the distance between the vacancy center and its first nearest neighbor, we explained the origin of the large linear temperature dependence of the entropy of formation.Replacing the Arrhenius relation by a dependence following from the LGT, we showed that the highly controversially discussed discrepancies between PAS and DD data practically disappear and that the accepted T ¼ 0 K defect formation enthalpies-as compiled, e.g., in the Landoldt-Börnstein series [39]-must be revised.While the difference between revised and unrevised values is negligible for the Al vacancy (0.01 eV), it is sizable for the Cu vacancy (0.22 eV) and should be considered in future compilations and studies.Using the revised experimental T ¼ 0 K values and comparing them with DFT calculations employing different xc functionals, we showed that the conventionally accepted picture according to which LDA outperforms the various flavors of GGA in describing vacancy formation energies is, to a large part, a consequence of comparing with the original Arrheniusextrapolated experimental data and changes qualitatively when comparing with the temperature corrected (revised) values derived here.An important consequence is that the apparent improvement when employing surface corrections or using surface-corrected xc functionals such as AM05 vanishes.For example, using the revised values, GGA-PBE performs better than all studied xc functionals, while the previously described surface corrections that were needed to get an improved agreement with the unrevised values perform significantly worse.Finally, we note that the mechanisms behind the strong local anharmonicity are not restricted to vacancy defects but apply to any defect configuration that induces locally strong asymmetric displacement potentials, e.g., bulk systems with point defects such as doped semiconductors, alloys, or larger-scale defects like surfaces or interfaces and grain boundaries.Including anharmonicity in future studies of such systems is therefore expected to be crucial.

FIG. 2 (
FIG. 2 (color online).Entropy of formation, S f ¼ −dG f =dT, for the Cu monovacancy (solid line) and divacancy (dashed line).The thinner lines indicate the entropy of formation considering only the electronic and quasiharmonic free-energy contribution (el þ qh).The thick curves include all contributions in particular anharmonicity (el þ qh þ ah).Green lines represent numbers suggested in Ref.[8] for explaining non-Arrhenius behavior within a monovacancy þ divacancy model.Experimental PAS[37] and DD entropies (filled/empty black bars) are derived from the experimental data shown in Fig.1.

TABLE I .
Comparison of the previously applied linear Arrhenius approximation and the newly proposed and ab initio based LGT.