Semilocal density functional theory with correct surface asymptotics

Semilocal density functional theory is the most used computational method for electronic structure calculations in theoretical solid-state physics and quantum chemistry of large systems, providing good accuracy with a very attractive computational cost. Nevertheless, because of the non-locality of the exchange-correlation hole outside a metal surface, it was always considered inappropriate to describe the correct surface asymptotics. Here, we derive, within the semilocal density functional theory formalism, an exact condition for the image-like surface asymptotics of both the exchange-correlation energy per particle and potential. We show that this condition can be easily incorporated into a practical computational tool, at the simple meta-generalized-gradient approximation level of theory. Using this tool, we also show that the Airy-gas model exhibits asymptotic properties that are closely related to the ones at metal surfaces. This result highlights the relevance of the linear effective potential model to the metal surface asymptotics.


I. INTRODUCTION
The exact form of the potential felt by an electron leaving from or approaching a metal surface is of great importance for a variety of physical phenomena, including the interpretation of image states 1 , modeling quantum-transport 2 , low-energy electron diffraction (LEED) 3 , scanning tunneling microscopy 4,5 , and inverse or two-photon photoemission spectroscopy 6,7 . The asymptotic form of this image potential is −1/(4(z −z 0 )), with z being the distance from the surface, and z 0 representing the position of the so-called image plane 8 , and should be reproduced by any computational method aiming at an accurate description of the surface physics.
Within Kohn-Sham (KS) density-functional theory (DFT) 9,10 , which is the most used computational method for electronic structure calculations in theoretical solidstate physics, the shape of the image potential is dictated by the properties of the effective Kohn-Sham (KS) potential. This depends on the employed approximation for the exchange-correlation (XC) functional E xc [ρ], which gives the XC potential via the relation v xc (r) = δE xc [ρ] δρ(r) , where ρ is the electron density. It has been shown that the exact v xc asymptotically approaches the image potential 8,[11][12][13][14] , despite a different result has been obtained within the plasmon-pole approximation 15 . The popular local density approximation (LDA) 10 and the generalized-gradient approximation (GGA), however, fail in this task 16 showing either a too fast decay (e.g. exponential), or an inaccurate description of the surface energetics (as for the Becke exchange 17 ). Ad-hoc modifications of the LDA XC potential 18,19 have then been proposed to improve the asymptotic behavior of the XC potential, but such methods are not functional derivatives of any energy functional. Alternatively, non-local methods outside the KS framework 11,13,14,[20][21][22] are employed.
An accurate KS-DFT method with the correct surface asymptotics would be desirable for many reasons, including the local nature and the computational efficiency. However, a good functional shall yield not only the correct asymptotic XC potential, but also accurate energies. Thus, it is necessary to be defined by a realistic energy per particle ε xc (r). The latter is not a uniquely defined physical quantity, but an exact reference for it is the the conventional energy per particle, which is associated with the interaction of an electron at r with the coupling-constant-averaged charge of its XC hole [23][24][25] . The exact (conventional) ε xc (r) at metal surfaces decays as −1/(4(z − z 0 )), i.e. as the image potential 11,26 .
The simultaneous description of surface asymptotic and energy properties is anyway an ambitious objective, which is in fact not achievable at the GGA level [35][36][37][38] . In this article, we show that the issue can be instead solved at the meta-GGA level of theory, employing an exact condition which yields the correct image-like asymptotic behavior of both ε xc and v xc at metal surfaces. This condition can be easily implemented in any meta-GGA functional, keeping its original accuracy for ground-state properties not related to surface asymptotics. Hence, an accurate KS-DFT method with correct metal-surface asymptotic can be obtained for application in many surface science problems.

II. EXACT CONDITION FOR ASYMPTOTIC PROPERTIES
To start, we consider the simplest (and most used) model for a metal surface: the semi-infinite jellium surface. This model system is very important in surface science and solid-state physics, containing the physics of simple metal surfaces 8,39,40 .
The KS single-particle orbitals have the form where k || and r || are the two-dimensional wave-vector and position vector in the plane xy of the surface, k z and z are the corresponding components in the direction perpendicular to the surface, S and L are the normalization area and length, and φ kz (z) are the eigenfunctions of a one-dimensional KS Hamiltonian (for details see Appendix A).
In the vacuum region, far away from the surface (z → ∞), the single-particle orbitals φ kz (z) behave as 27 : where α KS > 0. Now, we consider a meta-GGA XC energy per particle of the form where η is a parameter to be fixed later, α(r) = τ (r) − τ W (r))/τ T F (r) is the well-known meta-GGA ingredient that measures the non-locality of the kineticenergy density 38 , with τ , τ W , and τ T F being the positivedefined exact KS, von Weizsäcker, and Thomas-Fermi kinetic-energy densities, respectively. We recall that 1/ 1 + α(r) 2 is the electron localization function, often used in the characterization of chemical bonds 41 . Eq. (5) yields the following asymptotic behaviors (see the Appendix B for details): Here, the KS potential has been obtained in the generalized KS framework using the formula 42 Equations. (6) and (7) show that, in contrast to previously developed XC functionals, both v xc and ǫ xc are proportional to the exact ones: if η = η 1 = 4 (η = η 2 = 6) then the exact energy-density (potential) is obtained. Unfortunately, η 1 = η 2 . Nevertheless, for both values Eq. (5) yields an asymptotic behavior qualitatively and quantitatively significantly beyond the current state-ofthe-art. We also remark that Eq. (5) is solely based on the properties of the reduced kinetic ingredient α. However, at the meta-GGA level of theory, other ingredients are also available (e.g. the gradient and the Laplacian of the density) so that the exact asymptotic description of both ǫ xc and v xc might be achieved.

III. PRACTICAL COMPUTATIONAL TOOL
As a first practical example, we consider the case η = η 1 and incorporate the condition of Eq. (5) into the popular TPSS meta-GGA functional 43 , using an approach similar to that of Ref. 44. The resulting XC functional will be termed surface-asymptotics (SA) TPSS. This functional is obtained by simply changing, in the TPSS exchange formula, the parameter κ (which determines the asymptotic behavior of the functional) from its original value of 0.804 to The correlation is left unchanged. (Note that the TPSS correlation decays exponentially, and thus our XC condition is incorporated in the TPSS exchange functional. This is a common procedure for semilocal functionals, which are based on a strong error cancellation between the exchange and correlation parts.) In Eq. (9), the parameters a = 2.413 and b = 0.348 have been fixed by imposing the constraints: κ = 0.804 for α = 1 and α = 0; whereas, when degenerate orbitals contribute to the tail of the density (α → ∞), κ → 2π (5)).  In Fig. 1, we show the TPSS and SA-TPSS exchange enhancement factors versus α for several values of the reduced gradient s = |∇ρ|/(2(3π 2 ) 1/3 ρ 4/3 )). When s is small, TPSS and SA-TPSS coincide for all values of α. As s increases, TPSS and SA-TPSS start to differ, especially at large values of α, as expected.
In Fig. 2, fully self-consistent KS-LDA orbitals were used to obtain zǫ xc (z) (upper panel) and zv xc (z) (lower panel) as a function of the scaled distance z for jellium surfaces with the electron-density parameters r s = 2 and r s = 4. The TPSS functional yields a wrong (exponential) asymptotic behavior for both quantities. Instead, the SA-TPSS functional gives the following asymptotic behaviors: ǫ xc (z) → −0.25/z and v xc (z) → −0.375/z. Thus, Fig. 2 provides a numerical proof of the validity of Eq. (5), as well as a validation of the simple construction used to obtain the SA-TPSS functional.
In Tables I and II we report the TPSS and SA-TPSS results for surface energies and various molecular tests, respectively. By construction, whenever the surface asymptotics plays a negligible role (e.g., covalent interactions in molecules: W4, OMRE, IP13, MGBL19) both functionals yield very similar results. In the case of jellium surface energies and non-covalent interactions (DI6, HB6), however, SA-TPSS improves over the standard TPSS. Finally, we have to stress that the SA-TPSS meta-GGA gives only the correct asymptotic decay of the XC energy per particle and potential at metal surfaces, but it can not provide the exact behavior for the exchange and correlation components, separately. This is a difficult task, that in our opinion, a simple semilocal functional can not obey.

IV. AIRY GAS ASYMPTOTIC PROPERTIES
As an additional example of the use of Eq. (5) and of the SA-TPSS functional, we consider the Airy gas 51,52 , which is the simplest possible model for an edge electron gas (for details see Appendix A). This model system plays an important role in DFT 52-54 , as it incorporates the correct physics of a semi-infinite metal surface and is simple enough to allow for analytical calculations. To our knowledge, the exact asymptotic behavior of the XC energy per particle and of KS XC potential of the Airy gas are unknown. Nevertheless, we can use the SA-TPSS functional to obtain some information about them.
The Airy-gas electron density and positive-defined kinetic-energy density are [51][52][53]55 q, and α), as a function of the scaled distance z. The Airy bulk is at z ≤ 0. The exponential decay of the electron density occurs at z ≥ 0.
where z is the scaled distance and Ai(z) is the Airy function. In Fig. 3, we show the Airy-gas semilocal ingredient s, the reduced Laplacian q = ∇ 2 ρ/[4(3π 2 ) 2/3 ρ 5/3 ], and α. In the bulk (z → −∞), both s and q are small, while α → 1 (showing that the Thomas-Fermi theory becomes exact). In the vacuum, all semilocal ingredients diverge (as in the case of the jellium surface).
The SA-TPSS XC functional gives the following analytical expressions: This result is interesting, since it suggests that, for the Airy gas, both ǫ xc and v xc decay as −1/z, as in the case of the jellium surface. Furthermore, because the coefficients (0.217 and 0.325) are close to, but smaller than, the ones for the jellium surface (0.25 and 0.375), we can conclude that at a metal surface the main contribution to the asymptotics comes from the region near the surface, where the effective potential is linear and well described by the Airy-gas model. We note that such a result cannot be obtained without the proper inclusion of exact surface conditions into the functional. This is shown in Fig.  4 where we plot, for several functionals, zǫ xc (z) (upper panel) and zv xc (z) (lower panel), versus the scaled distance z, in the vacuum region of the Airy gas. The TPSS functional shows a rather unphysical exponential decay for both ε xc and v xc . On the other hand, the AM05 functional 52 , whose energy density is fitted to the Airy gas, is close to our SA-TPSS functional for ǫ xc but displays a decay that is too fast for v xc . This latter feature represents a manifestation of the impossibility, at the GGA level, to describe correctly both ǫ xc and v xc , which can only be overcome at the meta-GGA level 38 .
The exchange-only asymptotic behaviors (ǫ x and v x ) at metal surfaces, are depending on the bulk and surface parameters (k F and W ) [27][28][29] , and thus they are created by bulk-and surface-electrons. When correlation is included, screening effects dump the bulk contribution, such that the asymptotic properties of ǫ xc and v xc are mainly created by a density-independent XC effect of the surface region, as proved by the SA-TPSS result for the Airy gas. Thus, the Airy gas model system can be efficiently used in modelling various phenomena outside metal surfaces, even being an alternative to the ad-hoc LDA XC potential modifications 18,19 .

V. CONCLUSIONS
In conclusion, we have derived an exact meta-GGA condition for the correct image-like surface asymptotics of the XC energy per particle ε xc and the KS XC potential v xc . Our formula [Eq. (5)] depends only on the semilocal ingredient α and takes advantage of the nonlocality of the kinetic energy density beyond the von Weizsäcker term 38 . The existence of this exact condition represents an important contribution in the framework of DFT, as it shows that surface asymptotics can be described by semilocal meta-GGA functionals. On the contrary, no GGA can be constructed that is able to describe correctly the asymptotics of both ε xc and v xc . In fact, there is, to our knowledge, no GGA functional that yields a realistic KS XC potential at metal surfaces.
We have demonstrated that our exact condition can be easily implemented in any meta-GGA functional, keeping its original accuracy for standard ground-state properties and providing, at the same time, a correct description of the surface asymptotics. We have constructed the SA-TPSS functional, which we have shown to perform as the TPSS for covalent chemistry and to improve over it for non-covalent interactions and surface-related problems. This new functional can thus be a promising tool for the investigation of surface-sensitive electronic-structure calculations, such as molecule/molecular complex-surface, cluster-surface, and surface-surface interactions.
We thank TURBOMOLE GmbH for providing the TURBOMOLE program package.

APPENDIX A
In case of semi-infinite jellium surfaces, the oneparticle eigenfunctions φ kz (z) have a continuous energy spectrum ǫ kz = V KS (∞) + k 2 z /2, and are solutions of the one-dimensional KS equation Here V KS (z) = V H (z) + V xc (z) is the sum of the total classical electrostatic potential (which incorporates the positive background), and the XC potential.
In case of the Airy gas, the effective potential has the linear form with F > 0 being the slope of the effective potential 51 . Thus, the one-particle normalized eigenfunctions φ j (Z) satisfy the one-dimensional equation being proportional to the Airy function. Here Z is the distance perpendicular to the surface. It is convenient to consider the scaled distance z = Z(2F ) 1/351 , as used in the Section IV.

APPENDIX B
The electron density and kinetic energy density of a jellim surface are with k F being the magnitude of the bulk Fermi wavevector. Using Eqs. (3) and (4) we obtain, when z → ∞, the following expressions ρ → 1/4 N 2 W k F π 2 z 2 − 1/8 where N is a normalization constant, W is the work function, and k F is the bulk Fermi wave-vector. Here F xc = ǫ MGGA xc /ǫ LDA x is the enhancement factor corre-sponding to the energy density of Eq. (5).