Exchange-Correlation Hole of a Generalized Gradient Approximation for Solids and Surfaces

We propose a generalized gradient approximation (GGA) for the angle- and system-averaged exchange-correlation hole of a many-electron system. This hole, which satisfies known exact constraints, recovers the PBEsol (Perdew-Burke-Ernzerhof for solids) exchange-correlation energy functional, a GGA that accurately describes the equilibrium properties of densely packed solids and their surfaces. We find that our PBEsol exchange-correlation hole describes the wavevector analysis of the jellium exchange-correlation surface energy in agreement with a sophisticated time-dependent density-functional calculation (whose three-dimensional wavevector analysis we report here).


I. INTRODUCTION
In the Kohn-Sham (KS) density functional theory 1 for the ground-state energy of a many-electron system, only the exchange-correlation (xc) energy E xc [n] has to be approximated. The exact xc energy of an arbitrary inhomogeneous system of density n(r), which incorporates all the quantum many-body effects beyond the Hartree approximation, can be obtained from the spherical averagē n xc (r, u) of the coupling-constant averaged xc hole densityn xc (r, r ′ ) at r ′ around an electron at r as follows 2,3 E xc [n] = dr n(r) ε xc [n](r), (1) where ε xc [n](r) is the xc energy per particle at point r: withn xc (r, u) = 1 4π dΩn xc (r, r ′ ), dΩ being a differential solid angle around the direction of u = r ′ − r, and k representing the magnitude of the wavevector. (Unless otherwise stated, atomic units are used throughout, i.e., e 2 = = m e = 1.) The "Jacob's ladder" classification of the widelyused ground-state density-functional approximations for E xc [n] andn xc (r, u) has three complete non-empirical rungs: the local-spin-density approximation (LSDA), 1 the generalized-gradient approximation (GGA), 4,5,6 and the meta-GGA. 7,8 Due to its simplicity and accuracy, one of the most commonly used xc density-functional approximation in solid-state physics and quantum chemistry calculations is nowadays the semilocal PBE-GGA. 4 Recent work 9 has shown, however, that the exchange density-functional approximations should recover, in the limit of slowly-varying densities, the universal secondorder gradient-expansion (GE) approximation of the exchange energy, where ǫ unif x is the exchange energy per particle of the uniform electron gas, µ GE x = 10/81 is the GE exchange coefficient, 10 and s = |∇n|/(2k F n) is the reduced density gradient which measures the variation of the electron density over a Fermi wavelength λ F = 2π/k F , with k F = (3π 2 n) 1/3 representing the magnitude of the local Fermi wavevector. Recovery of the correct secondorder gradient expansion for correlation 11 in the slowlyvarying limit is much less important for the construction of density-functional approximations. (See Table 1 of  Ref. 9). A GGA, which has as ingredients only the spin densities n ↑ and n ↓ and their gradients ∇n ↑ and ∇n ↓ , cannot recover, in the limit of slowly varying densities, the GE approximation of the exchange energy and at the same time be accurate for atoms. 9,12 The semilocal PBE has the correct correlation GE coefficient in the high-density limit, and is accurate for atoms, but its exchange GE coefficient is almost twice as large as the exact coefficient, i.e., µ P BE x ≈ 2µ GE x . Because of this, 12 PBE overestimates the equilibrium lattice constants of solids and yields surface energies that are too low.
Following the ideas of Ref. 9, PBEsol (PBE for solids) was constructed 12 . PBEsol is a GGA that has the same form as PBE but restores the density-gradient expansion for exchange by replacing µ P BE x = 0.2195 with µ P BEsol x . By fitting the jellium xc surface energies (as had been done previously in Ref. 13 in the construction of a GGA which relies on the Airy-gas approximation 14 ), the PBEsol correlation GE coefficient was set to µ P BEsol c = 0.046. (For PBE, µ P BE c = 0.0667). Thus, PBEsol can easily be applied in solid-state calculations (just by changing the coefficients in a PBE code) and yields good equilibrium lattice constants and jellium surface energies. 12 Several other applications of PBEsol have already proved the accuracy of this GGA for solids. In particular, PBEsol considerably improves the structure of gold clusters 15 and works better than PBE for isomerization energies and isodesmic stabilization energies of hydrocarbon molecules. 16 PBEsol also describes ferro-and anti-ferro-electric ABO 3 crystals 17 much better than LSDA or PBE-GGA.
In this paper, we first construct the PBEsol angleaveraged xc hole densityn P BEsol xc (r, u). A nonempirical derivation of the PBE xc hole was reported in Ref. 5, starting from the second-order density-gradient expansion of the xc hole and cutting off the spurious largeu contributions to satisfy exact constraints according to which (i) the exchange-hole density must be negative, (ii) the exchange hole must integrate to -1, and (iii) the correlation hole must integrate to zero. Later on, a fully smoothed analytic model was constructed for the PBE exchange hole. 6 Our construction of the PBEsol xc hole begins with and appropriately modifies the sharp cutoff correlation hole of Ref. 5 and the smooth exchange hole of Ref. 6. It should be recalled that, because of an integration by parts that occurs in the underlying gradient expansion, a GGA hole is meaningful only after averaging over the electron density n(r) (as in our tests and applications), and this system-averaging itself smooths sharp cutoffs.
Finally, we use our PBEsol xc hole to carry out a three-dimensional (3D) wavevector analysis of the jellium xc surface energy. This wavevector analysis was carried out in Ref. 3 in the random-phase approximation (RPA). Here, we go beyond the RPA in the framework of timedependent density-functional theory (TDDFT), and we compare these calculations with the results we obtain from our PBEsol xc hole density.
The paper is organized as follows. In Sec. II, we present the PBEsol angle-averaged xc-hole densitȳ n P BEsol xc (r, u). In section III, we perform the wavevector analysis of the jellium xc surface energy. In Sec. IV, we summarize our conclusions.

II. PBESOL-GGA ANGLE-AVERAGED EXCHANGE-CORRELATION HOLE
We assume here that the PBEsol correlation energy can be constructed from a gradient expansion for the correlation hole in much the same way that the PBE correlation energy was so constructed 5 . The GGA angleaveraged correlation hole is 5 where r s = (9π/4) 1/3 /k F is a local density parameter, ζ = (n ↑ − n ↓ )/(n ↑ + n ↓ ) is the relative spin polariza-tion, φ = [(1 + ζ) 2/3 + (1 − ζ) 2/3 ]/2 is a spin-scaling factor, v = φk s u with k s = (4k F /π) 1/2 is the reduced electron-electron separation on the scale of the screening length, and t = |∇n|/(2φk s n) is the reduced density gradient measuring the variation of the electron density over the screening length. The sharp cutoff v c is found such that Eq. (5) satisfies the correlation hole sum rule dr n c (r, is the LSDA correlation hole 18 given by Eq. (45) of Ref. 5, and the gradient correction to the correlation hole is given by the following expression: 5 is constructed so that the second-order gradient expansion of the PBEsol correlation energy is recovered. Here E 1 (y) = ye y ∞ y dt e −t /t is between 0 and 1, and µ GGA c is the GGA gradient coefficient in the slowly-varying limit (µ P BE c for PBE and µ P BEsol c for PBE-sol). In  (6) can be negative at small v and small r s ; however, because of the energy sum rule both ∞ 0 du uB c (r s , ζ, v) and ∞ 0 du u 2 B c (r s , ζ, v) are positive, which ensures that the cutoff procedure is correct and for every value of r s , ζ, and t there is a v c such that Eq. (5) satisfies the correlation-hole sum rule.
The exchange energy and exchange-hole density for a spin-polarized system may be evaluated from their spin-unpolarized counterparts by using the spin-scaling and thus, we only need to consider the spin-unpolarized system. As in the case of the analytical PBE exchange hole of Ref. 6, we choose the following ansatz for the nonoscillatory dimensionless exchange-hole shape: and also the small-u behavior of the exchange hole is recovered by 6 : Here F P BEsol x (s) is the PBEsol enhancement factor. 12 The integrals of Eqs. (11) and (12) can be solved analytically, 6 and Eqs. (11)-(13) Fig. 2) can be fitted to the following analytic expression: where  In Fig. 3, we plot the dimensionless exchange hole shape J P BEsol (s, y) [using the analytical fit of Eq. (14)] versus y = k F u for several values of the reduced gradient s. Our J P BEsol (s, y) looks similar to the J P BE (s, y) of Ref. 6, but J P BE (s, y) is deeper because µ P BE x = 0.2195 > µ P BEsol x = 0.1235. Finally, we look at the xc enhancement factor, which displays the nonlocality: 22 ǫ unif x (n) being the exchange energy per particle of a spinunpolarized uniform electron gas. For a spin-unpolarized system in the high-density limit (r s → 0) the exchange energy is dominant and Eq. (15) defines the exchange

III. WAVEVECTOR ANALYSIS OF THE JELLIUM XC SURFACE ENERGY
The xc surface energy, σ xc , can be defined as the xc energy cost per per unit area to create a planar surface by cutting the bulk. In a jellium model, in which the electron system is translationally invariant in the plane of the surface, and assuming the surface to be normal to the z-axis, the surface energy can be written as follows 3  where 24 is the wavevector-resolved xc surface energy, and b xc (k, z) = 4π ∞ 0 du u 2 sin ku ku n xc (z, u) −n unif xc (u) . (18) Equations (16)- (18) comprise an angle-averaged threedimensional wavevector analysis of the surface xc energy into contributions from density fluctuations of various wavevectors k, following from the Fourier transform of the Coulomb interaction in Eq. (2). In these and subsequent equations, k F = (3π 2n ) 1/3 is the bulk (not the local) Fermi wavevector, and r s is the bulk (not the local) density parameter.
The exact low-wavevector limit of σ xc is known to be 2 where ω p = (4πn) 1/2 and ω s = ω p / √ 2 are the bulkand surface-plasmon energies, andn is the bulk density. Eq. (19) was used in the wavevector-interpolation approach reported in Refs. 2, 7, and 25, and it was naturally recovered by the RPA approach reported in Ref. 3.
In the calculations presented below, we have considered, as in Ref. 3, a jellium slab of background thickness a = 2.23λ F (where λ F = 2π/k F ) and bulk parameter r s = 2.07. This slab corresponds to about four atomic layers of Al(100).
For the GGA calculations of γ xc (k), the function b xc (k, z) entering Eq. (17)  The wavevector-resolved exact-exchange surface energy γ x (k) is shown in Fig. 6. Figures 7-8 and Fig. 9 show, respectively, the wavevector-resolved xc and correlationonly surface energies γ xc (k) and γ c (k). Figure 6 shows that γ P BEsol In Fig. 7, we compare the wavevector-resolved exact-RPA surface energy (as reported in Ref. 3) with its TDDFT counterpart (which we have not reported elsewhere, and which required three months of computation). In the long-wavelength limit (k → 0), both RPA and TDDFT calculations approach the exact lowwavevector limit of Eq. (19). In the large-k limit, where RPA is wrong, our TDDFT approach (which reproduces accurately the xc energy of the uniform electron gas) is expected to be accurate. Furthermore, the uniform-gas-based isotropic xc kernel that we use in our TDDFT calculation has been shown recently to yield essentially the same two-dimensional (2D) wavevector analysis as a more sophisticated high-level correlated approach (the inhomogeneous Singwi-Tosi-Land-Sjölander method) which does not use an isotropic kernel derived from the uniform gas. 34 Hence, we take the TDDFT wavevector-resolved surface energy represented in Fig. 7 by a solid line as the benchmark curve against which we compare various GGA's. Figure 8 shows our wavevector-resolved TDDFT surface energy together with its PBE and PBEsol counter- parts. γ P BEsol xc is nearly exact at small wavevectors, where it matches the exact initial slope of Eq. (19). At intermediate wavevectors, γ P BEsol xc is very close to our benchmark TDDFT calculation. As in the case of wavevector-resolved exchange-only surface energies, PBE and PBEsol calculations correctly recover the nonoscillatory LSDA and differ from the more accurate TDDFT calculation due to the inaccuracy of the nonoscillatory model of the exchange-hole shape [see, e.g., Eq. (10)]. Figure 8 shows that γ P BEsol xc nicely matches our benchmark TDDFT calculation at low and intermediate wavevectors, so we recommend using the model PBEsol xc hole, in all solid-state calculations where the hole is needed but full nonlocality is not important, instead of the more expensive TDDFT. We recall that the system-averaged hole, unlike the energy density, is uniquely defined, and is an observable at full coupling strength 5 . Figure 9 exhibits the wavevector-resolved correlationonly PBE, PBEsol, exact-RPA, and TDDFT surface energies. We observe that at long wavelengths (k → 0), where the LSDA is known to fail badly, and at short wavelengths (large k), where RPA is wrong, the GGA's under consideration are considerably close to our benchmark TDDFT calculations. At intermediate wavevectors, however, GGA's cannot describe γ c accurately, although PBEsol has been shown in Fig. 8 to give a very good description of γ xc . This is due to a cancellation of the errors introduced at these wavevectors within the exchange and correlation contributions to γ xc , which almost cancel each other 35 .

IV. CONCLUSIONS
We have constructed a PBEsol angle-averaged xc holē n xc (r, u) that satisfies known exact constraints and recovers the recently reported PBEsol xc energy functional. Our construction of the PBEsol xc hole begins from and appropriately modifies the sharp cutoff correlation hole of Ref. 5 and the smooth exchange hole of Ref. 6. We also generalize [see Eq. (7)] the sharp cutoff procedure for the correlation hole to any GGA which has a positive gradient expansion coefficient.
We have found that our PBEsol xc hole describes accurately the wavevector-resolved xc jellium surface energy for all values of the wavevector, thus providing support for the PBEsol GGA for solids and surfaces.