Asymptotic Behavior of Exact Exchange for Slabs : Beyond the Leading Order

Far outside the surface of slabs, the exact exchange (EXX) potential vx falls off as −1/z, if z denotes the direction perpendicular to the surface and the slab is localized around z = 0. Similarly, the EXX energy density ex behaves as −n/(2z), where n is the electron density. Here, an alternative proof of these relations is given, in which the Coulomb singularity in the EXX energy is treated in a particularly careful fashion. This new approach allows the derivation of the next-to-leading order contributions to the asymptotic vx and ex. It turns out that in both cases, the corrections are proportional to 1/z2 in general.


Introduction and Summary of Results
The application of Kohn-Sham density functional theory (DFT) to surfaces was initiated by the seminal work of Lang and Kohn [1][2][3] (see also [4]).While the initial focus was on metals and the jellium model, soon, also semi-conducting materials were investigated [5].In practice, DFT calculations for surfaces often rely on slabs, in which a limited number of atomic layers simulates the bulk material [5,6].The number of layers required for an accurate description of the surface primarily depends on the electronic structure, in particular on the localization of the states.Obviously, a decoupling of the surface states on the two sides of the slab is desirable.If relaxation of the atomic positions at the surface or surface reconstruction is of interest, one has to make sure that some bulk-like layers remain in the middle of the slab.Nevertheless, often a rather small number of layers is sufficient in the case of non-metallic solids.Particularly accurate results for quantities such as the work function and the surface energy can be obtained by extrapolating data from slabs with different thicknesses to the limit of infinite thickness.
The most important quantity for the density functional description of surfaces and slabs is the exchange-correlation (xc) potential v xc (for an overview of DFT, see [19]).Starting with Lang and Kohn [1][2][3], the xc-potential at surfaces has been studied extensively .While the discussion of the xc-potential turned out to be difficult for semi-infinite matter (see [31,33]), definitive information is available for the exact exchange (EXX) potential v x of slabs.Both v x and the EXX energy density e x of metallic slabs have been investigated by Horowitz and collaborators [28,29,32,36] and by Qian [33] (based on extensive work by Sahni and collaborators [23][24][25][26] on semi-infinite matter), relying on the jellium model (i.e., on full translational symmetry in the directions parallel to the surface).If the slab is parallel to the xy-plane and has its center at z = 0, the exact v x behaves as −1/z far outside the surface (atomic units are used throughout this work).Similarly, the EXX energy density of jellium slabs falls off as −n/(2z), where n is the electron density [28].This latter result is of particular interest, since e x defines the Slater potential, v Slater = 2e x /n, which is the core contribution to the often used Krieger-Li-Iafrate approximation [44] for the exact v x and its refinement, the localized Hartree-Fock/common energy denominator approximation [45,46].Both results have subsequently been generalized [34,35] to non-jellium slabs, for which one has a Bravais lattice in the xy-directions, Here, L characterizes the position of the surface, i.e., the slab extends from −L to L.
In this contribution, a particularly careful derivation of the relations (1) and (2) is given, in which the crucial Coulomb singularity in the EXX energy is handled in a mathematically unambiguous fashion.In contrast to the approaches used in [34,35], this derivation allows one to extract the next-to-leading order contributions to both e x and v x .Moreover, unlike the argument in [35], the present proof of (2) does not employ the ultimate asymptotic form of the density from the very outset, but rather is based on a completely general variation of the asymptotic n as a function of r.
The derivation is prepared by a brief review of the asymptotic form of the Kohn-Sham (KS) states (in Section 2), relying on the analysis of [47].In particular, the coupling of the Fourier amplitudes for the states is discussed, in order to show which amplitudes are asymptotically dominant.On this basis, the asymptotic behavior of e x (Section 3), v Slater (Section 4) and v x (Section 5) is analyzed.In this analysis, the Coulomb singularity is kept under control by addition and subtraction of a suitably-screened exchange kernel.One finds that the next-to-leading order term in e x has the form −p(r)/(2z 2 ) where p is a generalized density with weights determined by the first moments of the transition matrix.p can be explicitly expressed in terms of n, if the asymptotic density is dominated by the states in the vicinity of a single k-point q in the interior of the first Brillouin zone (which is the standard situation for non-metallic slabs [43]).In this case, one finds p(r) z→∞ − −− → m qβ n(r), where m qβ denotes the first moment of the most weakly-decaying state β at q.In the same fashion, the next-to-leading order contribution to v x (r) is obtained as −m qβ /z 2 , and this correction also applies to v Slater .It seems worthwhile to remark that all results can be generalized to the situation that there are several degenerate most weakly-decaying states at q.Moreover, the derivation of higher order corrections can basically follow the same route.

Asymptotic Behavior of States
In the case of slabs, the total KS potential v s is invariant under translation by an arbitrary two-dimensional (2D) Bravais vector in the xy-directions, but confines the electrons in the z-direction, Here, r = (x, y), and G is a vector of the 2D reciprocal lattice in the xy-directions.The corresponding KS states are given by: where k is the 2D crystal momentum.The normalization is chosen so that c kα integrates up to one for a single 2D unit cell with area A, The coefficients c kα satisfy the KS equations on the reciprocal lattice, In [35], it has been shown that the exact v x behaves as −1/z for z → ∞.In the following, v s (r) is therefore assumed to have the asymptotic form: with u being allowed to vanish.For the Fourier component with G = 0, one then has: while all other v(G, z) decay at least as fast as 1/z 2 (typically even faster).It seems worthwhile to emphasize that the next-to-leading order contributions to v(0, z) not included in Equations ( 7) and ( 8) are irrelevant for the derivation of the asymptotic forms of e x , v Slater and v x in Sections 3-5.
The asymptotic behavior of the solutions of the differential equation ( 6) with the potential (8) has been analyzed in detail in [47].For large z, the solutions can be expressed as: (the quantum numbers kα are omitted for brevity in the rest of this section).Here, c h denotes the general solution of the homogeneous differential equation: For the asymptotic potential (8), this solution is to leading order given by: with: On the other hand, c i accounts for the coupling of all amplitudes in (6), (z 1 < z must be chosen sufficiently large, so that the z -integration only extends over the asymptotic region).
The asymptotic behavior of the solutions is primarily controlled by the exponent γ(G), Equation (12).As long as (G + k) 2 > k 2 for all G = 0, i.e., for states inside the first Brillouin zone (BZ), there is a unique most slowly-decaying amplitude, c h (0, z).For states on the boundary of the first BZ, on the other hand, one has a second amplitude c h ( Ḡ, z) with ( Ḡ + k) 2 = k 2 , which shows the same asymptotic decay as c h (0, z).In fact, for large z, these two amplitudes only differ by the constant f 0 in Equation (11), since γ( Ḡ) = γ(0).
For all states and all G = 0, the function c i (0, G , z) falls off faster than c h (0, z).The contribution of c i (0, G , z) is largest for states on the boundary of the first BZ, for which a second amplitude c h ( Ḡ, z) decays as slowly as c h (0, z).For these states, one finds that c i (0, Ḡ, z) is suppressed by: ∞ z dz v(− Ḡ, z ) compared to c h (0, z) (and an analogous statement holds for the relation between c i ( Ḡ, 0, z) and c h ( Ḡ, z)).However, ∞ z dz v(− Ḡ, z ) vanishes at least as fast as 1/z, so that the second term on the right-hand side of ( 9) is irrelevant for the asymptotically-leading amplitude(s) in the case of all states.One thus has: (and an analogous relation for c( Ḡ, z) for states on the boundary of the first BZ).
In the case of all other c(G, z), the contribution of c i (G, G , z) cannot be neglected, since, at least in general, c h (G, z) decays faster than the product v(G, z)c h (0, z).However, all c i (G, G , z) decay faster than the most weakly-decaying amplitude c h (0, z), so that Equation ( 9) can be solved iteratively for large z.A single iteration of this equation, i.e., replacement of the full solution c(G , z) on the right-hand side of ( 13) by c h (G , z), is sufficient to obtain the correct asymptotic behavior of the next-to-most weakly-decaying amplitudes.For all yet faster decaying amplitudes, further iteration can (but need not) be required, depending on the asymptotic behavior of the v(G, z) involved.
In summary, the asymptotic behavior the states well inside the first BZ is determined by the amplitude with G = 0, which, in turn, approaches the homogeneous solution (11) for large z.All amplitudes with G = 0 are suppressed relative to c h (0, z), with the suppression factor c(G, z)/c h (0, z) depending on the structure of the potential.In the case of the next-to-most weakly-decaying amplitudes, the suppression factor is either given by the most weakly-decaying Fourier component of the potential v(G 1 , z) with non-zero G 1 or by c h (G, z)/c h (0, z), if this ratio decays more slowly than v(G 1 , z).
In the following, it is assumed for simplicity that: (i) the most weakly-decaying occupied state in the first BZ, i.e., the minimum exponent γ(0), is found for a single k-point well inside the first BZ; (ii) γ(0) is Taylor expandable at this k-point; and (iii) all states in a complete neighborhood of this k-point are also occupied.In this situation, the statements of the previous paragraph apply to all asymptotically-relevant states, and the simplest version of the BZ integration scheme of [43] can be used.

Asymptotic Behavior of Exchange Energy Density
In the case of spin-saturated slabs, e x (r) and n(r) are given by: where Θ kα denotes the occupation of the state kα and the k-integrations extend over the first BZ.Insertion of (4) into the exchange energy density (15) and use of the 2D Fourier transform of the Coulomb interaction yields: In order to extend the k -integration to the complete 2D k-space, one uses the fact that the coefficients c kα (G) only depend on k + G, If one extends the occupation function accordingly, the exchange energy density can be expressed as: At this point, one can choose the origin of the k -integration to be at k, with: In the following, the behavior of e x in the regime z L is analyzed (with L characterizing the extension of the slab whose center is at the origin).In order to isolate the asymptotically-dominating term without any mathematical ambiguity, the integrand is first split into a short-and a long-wavelength component, The k -integral in ( 25) can be performed directly, For large z, the overlap function f (k , G, z, z ) vanishes exponentially with z.This is immediately clear for the integrand in (23), since, in view of Equation ( 14), c kα (0, z) has the form (11) for large z, and all c kα (G, z) decay even faster than c kα (0, z).Moreover, this exponential decay is preserved by the k-integration, which may be performed using the approach of [43].Consequently, the second term in (27) is exponentially smaller than the first term in the asymptotic regime, The function f (k , G, z, z ) also vanishes exponentially for z → ∞, the arguments being the same as in the case of z.In fact, the BZ integration scheme of [43] is particularly appropriate if both z and z are large.The exponential decay of f (k , G, z, z ) with large z allows one to expand the denominator in powers of z /(z + µ).Subsequently, one can extend the z -integration to infinity, without introducing additional power-law corrections, The leading term of this expansion can be evaluated by use of ( 5), This expression is easily identified as the Fourier component of the density (16), The long-range component of the exchange energy density therefore asymptotically behaves as: with: Since µ can be chosen reasonably small, one obtains for the leading two orders: Now, consider the short-wavelength component of e x , Equation (26).In order to analyze its behavior for large z, the k -integration is decomposed into integrals over an elementary (simply connected) cell V, which is periodically repeated to fill the complete 2D momentum space without leaving any void.V could, for instance, have the shape of the first BZ, but a reduced spatial extension.In this case, the vectors Q, which connect the centers of the cells, are correspondingly scaled reciprocal lattice vectors.However, for the present purpose, V could also be a small square, so that the vectors Q are the reciprocal lattice vectors of the associated 2D simple cubic lattice, As soon as Q = 0, there is no Coulomb singularity in the k -integral: in this case, the k -integral is well defined for arbitrary z, z .For z → ∞, these contributions decay exponentially faster than the term with Q = 0, so that only the latter term remains to be analyzed, For large z, an expansion of the expression in square brackets about k = 0 is legitimate, since: (i) the k -integral is dominated by the vicinity of the Coulomb singularity; (ii) the integrand falls off rapidly for large z; and (iii) V can be chosen much smaller than the first BZ, Due to the inversion symmetry of V, one has: The k -integral in the second term on the right-hand side of ( 35) can be evaluated further in polar coordinates, where S(ϕ) defines the boundary of the integration area V. If, for instance, V is chosen to be a square with its corners at (±d, ±d), S(ϕ) is given by:

The leading contribution to e S
x is thus obtained as It seems worthwhile to emphasize that the z -integral is well defined, since for z close to z, one has: Taking this into account, the leading term for large z is given by: since S > 0 for all ϕ and f (0, G, z, z ) decays exponentially with z .Upon addition of ( 33) and ( 39), one finds for the leading two orders of the asymptotic energy density, The result is independent of µ, as it should be.In real space, one thus obtains:

Asymptotic Behavior of Slater Potential
Equation (41) shows that the Slater potential asymptotically behaves as: In order to evaluate the next-to-leading order term further, the k-integrations involved in p(r) and n(r) need to be analyzed.Asymptotically, the density (30) is dominated by the coefficients c kα with the lowest exponents γ kα (G) for which Θ kα > 0. If the lowest value of γ kα (G) is found well inside the first BZ (which is usually the case [43]), the coefficients c kα (G = 0, z) decay most slowly, as discussed in Section 2, with c kα (0, z) given by (11), due to (14).The same argument can be applied to the asymptotic p(r), In general, one particular band α falls off most slowly for each k.On the other hand, if there are two degenerate highest occupied states α and α for some k, the asymptotically-leading terms in the corresponding coefficients c kα (0, z) and c kα (0, z) only differ by some constant, since γ kα (0) = γ kα (0).It is thus sufficient to restrict the expression (44) to the terms with α = α, Here, m kα denotes the first moment of |c kα (G, z)| 2 , in the non-degenerate case and a combination of the moments of the degenerate states otherwise.The combination of Equations ( 43) and ( 45) with the asymptotic form (11) of c kα shows that the asymptotic behavior of both the density (30) and the first moment (32) is controlled by a BZ integral of the form discussed in [43], i.e., Equation (18) of this reference.In the case of the density, the generic integrand A kα in Equation ( 18) of [43] corresponds to 2| f 0,kα | 2 , for the first moment A kα corresponds to 2m kα | f 0,kα | 2 .Both ( 43) and ( 45) are asymptotically dominated by the state(s) with the minimum exponent γ kα (G = 0), as discussed in [43].If there is only a single k-point q and band β for which γ kα (G = 0) reaches its minimum value inside the integration region, one finds: where Γ is defined by: p(r) therefore has the same asymptotic behavior as n(r), so that p/n z→∞ − −− → m qβ .The ratio p/n also approaches a constant for large z, if f 0,qβ (G = 0) = 0, since this affects the z-dependence of p and n in the same way.Moreover, the same statement applies, if there is more than one k-point and/or more than one band for which γ kα (G = 0) assumes its minimum value (see [43]).One thus finally arrives at: with m qβ being replaced by a superposition of the moments of all asymptotically-relevant states in the more general situation.Depending on the symmetry of the relevant states at k = q under the transformation z → −z, m qβ can vanish.This is the case in particular, if one has reflection symmetry with respect to the xy-plane, so that |c qβ (G, −z)| = |c qβ (G, z)|.

Asymptotic Behavior of Exact Exchange Potential
For a discussion of the exact v x , one starts with the total exchange energy per unit cell, expressed in terms of the amplitudes c kα (G, z), (which is obtained by integration of (17) over the complete unit cell).In the following, the variation of this expression under a variation δn(r) of the asymptotic density is studied.However, a variation of the asymptotic density implies a variation of c kα (G, z), since there exists a one-to-one correspondence between the KS states and the density (30), In addition, since only the asymptotic density is to be varied (but not the complete n(r)), δc kα (G, z) is non-zero only for large z [δc kα (G, z) can be assumed to have the standard shape of a test function].Now, consider the variation of E x under the variation of c kα (G, z), Using ( 19) and ( 20), the k -integration can be combined with the sum over G.If the origin of the resulting k -integration over the complete k-space is chosen to be at k, one obtains: Next, δE x is split into a long-and a short-wavelength component, in analogy to (24), The evaluation of δE L x proceeds in straightforward fashion.One first performs the k -integral, One can then follow the argument leading from ( 27) to (29), since c kα (G , z) falls off exponentially with z, irrespective of δc * kα (G , z), First, focus on the leading order term of the expansion (L0).This term can be evaluated further by use of the orthonormality relation (5), One can now re-introduce the x-y-integrations and identify δn(r), Finally, an expansion in powers of µ/z is legitimate, since µ can be chosen to be much smaller than the z-values for which Next, consider the short-wavelength component (56), The arguments that lead from ( 26) to ( 35) also apply to (62).Using (36), one arrives at: One can then employ (37) for the k -integration and expand in powers of 1/z, since the remaining z -integral is free of any singularities, Use of the orthonormality relation (5) and the density variation (52) then shows that this expression cancels exactly with the first order contribution to (61).
Finally, consider the first order contribution (L1) to the expansion in (58), A completely general representation of δE L1 x in terms of δn, as achieved for δE L0 x , is not obvious.In the following, therefore, only the most important class of variations will be considered.For very large z, the BZ integral in the density (30) is dominated by the vicinity of the k-point(s) for which the exponent γ kα (G), Equation (12), assumes its lowest value in the integration region [43].If one is only interested in a variation of the asymptotic density, it is therefore sufficient to vary only the amplitudes of the states in the vicinity of the states with minimal γ kα (G).Assume, for simplicity, that the lowest value of the exponent is obtained only for a single state, denoted by qβ.If only the amplitudes c kβ in the vicinity of q are varied, while all other c kα are kept unchanged, δn is given by: Now, consider δE L1 x for this variation.For any given δc * kα (G , z), i.e., any given state kα, the sum over α in (65) is dominated by the most weakly-decaying occupied amplitude among the c kα (G , z).In the case of the variation (66), δE L1 x thus reduces to: the same asymptotic behavior (due to symmetry), this coefficient has to be suitably generalized (compare [43]).
The results (71)-(73) show exactly the required behavior under a translation of the coordinate system.Consider two coordinate systems whose origin differs by a shift ∆ in the z-direction, z = z + ∆.At the same physical point far from the surface, Equation (73) then yields in the two coordinate systems, , z ) c kα (G , z ) + c.c. .