Asymptotic Behavior of Exact Exchange for Slabs:

Asymptotic Behavior of Exact Exchange for Slabs:

Beyond the Leading Order

Eberhard Engel ^ID

Center for Scientific Computing, J. W. Goethe-Universität Frankfurt, Max-von-Laue-Str. 1, D-60438 Frankfurt/Main, Germany; engel@th.physik.uni-frankfurt.de

Received: 26 February 2018; Accepted: 21 April 2018; Published: 29 April 2018

Abstract:Far outside the surface of slabs, the exact exchange (EXX) potentialvxfalls off as−1/z, ifz denotes the direction perpendicular to the surface and the slab is localized aroundz=0. Similarly, the EXX energy densityexbehaves as−n/(2z), wherenis 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 asymptoticvxandex. It turns out that in both cases, the corrections are proportional to 1/z²in general.

Keywords:density functional theory; exact exchange; slabs; asymptotic behavior PACS:71.15.Mb; 71.15.-m

1. 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–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.

However, slabs are also of interest in themselves. The most prominent example is single- and bi-layer graphene [7–10], but also other 2D materials have attracted considerable attention recently, such as silicene [11–13], hexagonal boron nitride nanosheets [14,15] and mono- and multi-layer transition metal dichalcogenides [16–18] (as well as structures obtained by a combination of these materials with graphene).

The most important quantity for the density functional description of surfaces and slabs is the exchange-correlation (xc) potentialvxc(for an overview of DFT, see [19]). Starting with Lang and Kohn [1–3], the xc-potential at surfaces has been studied extensively [20–43]. 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) potentialvxof slabs. Bothvxand the EXX energy densityexof 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–26] on semi-infinite matter), relying on the

Computation2018,6, 35; doi:10.3390/computation6020035 www.mdpi.com/journal/computation

jellium model (i.e., on full translational symmetry in the directions parallel to the surface). If the slab is parallel to thexy-plane and has its center atz= 0, the exactvxbehaves as−1/zfar outside the surface (atomic units are used throughout this work). Similarly, the EXX energy density of jellium slabs falls off as −n/(2z), wheren is the electron density [28]. This latter result is of particular interest, sinceex defines the Slater potential,v_Slater = 2ex/n, which is the core contribution to the often used Krieger–Li–Iafrate approximation [44] for the exactvxand 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 thexy-directions,

ex(r) −−→^zL −ⁿ(r)

2z (1)

vx(r) −−→^zL −¹

z . (2)

Here,Lcharacterizes the position of the surface, i.e., the slab extends from−LtoL.

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 bothexandvx. 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 asymptoticnas a function ofr.

The derivation is prepared by a brief review of the asymptotic form of the Kohn–Sham (KS) states (in Section2), 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 ofex(Section3),v_Slater(Section4) andvx(Section5) 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 inexhas the form−p(r)/(2z²)where pis a generalized density with weights determined by the first moments of the transition matrix.

pcan be explicitly expressed in terms ofn, if the asymptotic density is dominated by the states in the vicinity of a singlek-pointqin the interior of the first Brillouin zone (which is the standard situation for non-metallic slabs [43]). In this case, one findsp(r)−−−→^z→∞ mqβn(r), wheremqβdenotes the first moment of the most weakly-decaying stateβatq. In the same fashion, the next-to-leading order contribution tovx(r)is obtained as−m_qβ/z², and this correction also applies tov_Slater. It seems worthwhile to remark that all results can be generalized to the situation that there are several degenerate most weakly-decaying states atq. Moreover, the derivation of higher order corrections can basically follow the same route.

2. Asymptotic Behavior of States

In the case of slabs, the total KS potentialvs is invariant under translation by an arbitrary two-dimensional (2D) Bravais vector in thexy-directions, but confines the electrons in thez-direction,

vs(r) =

∑

G

e^iG·rkv(G,z). (3)

Here, r_k = (x,y), and G is a vector of the 2D reciprocal lattice in the xy-directions. The corresponding KS states are given by:

φ_kα(r) = ^e

ik·r_k

√A

∑

G

e^iG·rkckα(G,z), (4)

wherekis the 2D crystal momentum. The normalization is chosen so thatc_kαintegrates up to one for a single 2D unit cell with areaA,

δ_α,α⁰ = Z _∞

−∞dz

∑

G

c^∗_kα(G,z)_ckα0(G,z)_{. (5)} The coefficientsc_kαsatisfy the KS equations on the reciprocal lattice,

(G+k)²−∂²_z

2 c_kα(G,z) +

∑

G⁰

v(G−G⁰,z)c_kα(G⁰,z) = e_kαc_kα(G,z). (6) In [35], it has been shown that the exactvxbehaves as−1/zforz→∞. In the following,vs(r)is therefore assumed to have the asymptotic form:

vs(r) −−→^zL −^u

z +. . . , (7)

withubeing allowed to vanish. For the Fourier component withG=0, one then has:

v(0,z) −−→^zL −^u

z+. . . , (8)

while all otherv(G,z)decay at least as fast as 1/z²(typically even faster). It seems worthwhile to emphasize that the next-to-leading order contributions tov(0,z)not included in Equations (7) and (8) are irrelevant for the derivation of the asymptotic forms ofex,v_Slaterandvxin Sections3–5.

The asymptotic behavior of the solutions of the differential equation (6) with the potential (8) has been analyzed in detail in [47]. For largez, the solutions can be expressed as:

c(G,z) = c_h(G,z) +

∑

G⁰6=G

ci(G,G⁰,z) (9)

(the quantum numberskαare omitted for brevity in the rest of this section). Here, chdenotes the general solution of the homogeneous differential equation:

∂²_zch(G,z) = ^h(G+k)²−2e+2v(0,z)ⁱch(G,z). (10) For the asymptotic potential (8), this solution is to leading order given by:

c_h(G,z) −−→^zL f0(G)z^u/γ(G)e^−γ(G)z, (11) with:

γ(G) = ^h(G+k)²−2ei1/2

. (12)

On the other hand,ciaccounts for the coupling of all amplitudes in (6), ci(G,G⁰,z) = − ¹

γ(G) Z z

z₁dz⁰ z z⁰

u/γ(G)

e^{γ(G)(z0−z)}v(G−G⁰,z⁰)c(G⁰,z⁰)

− ¹ γ(G)

Z _∞ z dz⁰

z⁰ z

u/γ(G)

e^{γ(G)(z−z0)}v(G−G⁰,z⁰)c(G⁰,z⁰) (13) (z₁ < 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)² > k² for all G 6= 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 amplitudec_h(G,^¯ z)with(G^¯+k)² =k², which shows the same asymptotic decay asch(_0,z). In fact, for largez, these two amplitudes only differ by the constant f0in Equation (11), sinceγ(G^¯) =γ(0).

For all states and allG⁰6=0, the functionci(0,G⁰,z)falls off faster thanc_h(0,z). The contribution ofc_i(0,G⁰,z)is largest for states on the boundary of the first BZ, for which a second amplitudec_h(G,^¯ z) decays as slowly asch(0,z). For these states, one finds thatci(0, ¯G,z)is suppressed by:

Z _∞

z dz⁰v(−G,^¯ z⁰)

compared to c_h(0,z) (and an analogous statement holds for the relation between c_i(G,^¯ 0,z) and ch(G,^¯ z)). However,R_∞

z dz⁰v(−G,^¯ 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:

c(0,z) −−→^zL c_h(0,z) (14)

(and an analogous relation forc(G,^¯ z)for states on the boundary of the first BZ).

In the case of all otherc(G,z), the contribution ofci(G,G⁰,z)cannot be neglected, since, at least in general,c_h(G,z)decays faster than the productv(G,z)c_h(0,z). However, allci(G,G⁰,z)decay faster than the most weakly-decaying amplitudec_h(0,z), so that Equation (9) can be solved iteratively for largez. A single iteration of this equation, i.e., replacement of the full solutionc(G⁰,z)on the right-hand side of (13) bych(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 thev(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 6= 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 potentialv(G1,z)with non-zero G1 or byc_h(G,z)/c_h(0,z), if this ratio decays more slowly thanv(G₁,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 singlek-point well inside the first BZ; (ii)γ(0)is Taylor expandable at thisk-point; and (iii) all states in a complete neighborhood of thisk-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.

3. Asymptotic Behavior of Exchange Energy Density

In the case of spin-saturated slabs,ex(r)andn(r)are given by:

ex(r) = −A² Z

1BZ

d²k (2π)²

d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘk⁰α⁰ Z

d³r⁰ φ_kα^† (r)φ_k⁰_α⁰(r)φ^†_k0α⁰(r⁰)φ_kα(r⁰)

|r−r⁰| ⁽¹⁵⁾ n(r) = 2A

Z

1BZ

d²k (2π)²

∑

α

Θkα|φ_kα(r)|², (16)

whereΘkα denotes the occupation of the statekα and thek-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:

ex(r) =

∑

G

e^iG·rkex(G,z) (17)

ex(G,z) = − Z

1BZ

d²k (2π)²

d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘk⁰α⁰

∑

G⁰ Z _∞

−∞dz⁰ 2πe^{−|k−k0+G0||z−z0|}

|k−k⁰+G⁰|

×

∑

G⁰⁰

c^∗_k0α⁰(G⁰⁰−G⁰,z)c_kα(G⁰⁰+G,z)

∑

G⁰⁰⁰

c^∗_kα(G⁰⁰⁰,z⁰)c_k⁰_α⁰(G⁰⁰⁰−G⁰,z⁰). (18) In order to extend thek⁰-integration to the complete 2D k-space, one uses the fact that the coefficientsc_kα(G)only depend onk+G,

c_kα(G,z) = c_k+G,α(0,z). (19)

If one extends the occupation function accordingly,

Θk+G,α := Θkα, (20)

the exchange energy density can be expressed as:

ex(G,z) = − Z

1BZ

d²k (2π)²

Z d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘk⁰α⁰ Z _∞

−∞dz⁰2πe^{−|k0−k||z−z0|}

|k⁰−k|

×

∑

G⁰

c^∗_k0α⁰(G⁰,z)c_kα(G⁰+G,z)

∑

G⁰⁰

c^∗_kα(G⁰⁰,z⁰)c_k⁰_α⁰(G⁰⁰,z⁰). (21) At this point, one can choose the origin of thek⁰-integration to be atk,

ex(G,z) = −

Z d²k⁰ (2π)²

Z _∞

−∞dz⁰ 2πe^{−|k0||z−z0|}

|k⁰| ^f(k⁰,G,z,z⁰)_{, (22)} with:

f(k⁰,G,z,z⁰) := Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘ_k+k⁰_α⁰

∑

G⁰

c^∗_k+k0α⁰(G⁰,z)ckα(G⁰+G,z)

×

∑

G⁰⁰

c^∗_kα(G⁰⁰,z⁰)c_k+k⁰_α⁰(G⁰⁰,z⁰). (23) In the following, the behavior ofexin the regimezLis analyzed (withLcharacterizing 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,

ex(G,z) = e^L_x(G,z) +e^S_x(G,z) (24)

e^L_x(G,z) = − Z d²k⁰

2π Z _∞

−∞dz⁰e^{−|k0|(|z−z0|+µ)}

|k⁰| ^f(0,G,z,z⁰) (25) e^S_x(G,z) = −

Z d²k⁰ 2π

Z _∞

−∞dz⁰e^{−|k0||z−z0|}

|k⁰| h

f(k⁰,G,z,z⁰)−f(0,G,z,z⁰)e^−µ|k0|i

. (26)

Thek⁰-integral in (25) can be performed directly, e^L_x(G,z) = −

Z _z

−∞dz⁰ f(0,G,z,z⁰) z−z⁰+µ −

Z _∞

z dz⁰ f(0,G,z,z⁰)

z⁰−z+µ . (27)

For largez, the overlap function f(k⁰,G,z,z⁰)vanishes exponentially withz. This is immediately clear for the integrand in (23), since, in view of Equation (14),ckα(0,z)has the form (11) for largez, and allc_kα(G,z)decay even faster thanc_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,

e^L_x(G,z)−−→ −^zL Z z

−∞dz⁰ f(0,G,z,z⁰)

z−z⁰+µ . (28)

The function f(k⁰,G,z,z⁰)also vanishes exponentially forz⁰ →∞, the arguments being the same as in the case ofz. In fact, the BZ integration scheme of [43] is particularly appropriate if bothzandz⁰ are large. The exponential decay of f(k⁰,G,z,z⁰)_{with large}z⁰allows one to expand the denominator in powers ofz⁰/(z+µ). Subsequently, one can extend thez⁰-integration to infinity, without introducing additional power-law corrections,

e^L_x(G,z) −−→^zL − Z _∞

−∞dz⁰ f(0,G,z,z⁰) z+µ

1+ ^z

0

z+µ+. . .

. (29)

The leading term of this expansion can be evaluated by use of (5), Z _∞

−∞dz⁰ f(0,G,z,z⁰) = Z

1BZ

d²k (2π)²

∑

α

Θkα

∑

G⁰

c^∗_kα(G⁰,z)c_kα(G⁰+G,z) . This expression is easily identified as the Fourier component of the density (16),

n(r) = 2

∑

G

e^iG·rk Z

1BZ

d²k (2π)²

∑

α

Θkα

∑

G⁰

c^∗_kα(G⁰,z)ckα(G⁰+G,z) . (30) The long-range component of the exchange energy density therefore asymptotically behaves as:

e^L_x(G,z) −−→^zL − ⁿ(G,z)

2(z+µ)− ^p(G,z)

2(z+µ)^{2. (31)}

with:

p(G,z) = ₂ Z _∞

−∞dz⁰z⁰ f(_0,G,z,z⁰)_{. (32)}

Sinceµcan be chosen reasonably small, one obtains for the leading two orders:

e^L_x(G,z) −−→^zL −ⁿ(G,z) 2z

1−^µ

z

− ^p(G,z)

2z² +. . . . (33)

Now, consider the short-wavelength component ofex, Equation (26). In order to analyze its behavior for large z, thek⁰-integration is decomposed into integrals over an elementary (simply connected) cellV, which is periodically repeated to fill the complete 2D momentum space without leaving any void.Vcould, for instance, have the shape of the first BZ, but a reduced spatial extension.

In this case, the vectorsQ, which connect the centers of the cells, are correspondingly scaled reciprocal

lattice vectors. However, for the present purpose,Vcould also be a small square, so that the vectorsQ are the reciprocal lattice vectors of the associated 2D simple cubic lattice,

e^S_x(G,z) = −

∑

Q Z

V

d²k⁰ 2π

Z _∞

−∞dz⁰ e^{−|k0−Q||z−z0|}

|k⁰−Q| h

f(k⁰−Q,G,z,z⁰)− f(0,G,z,z⁰)e^{−µ|k0−Q|}i .

As soon asQ6=0, there is no Coulomb singularity in thek⁰-integral: in this case, thek⁰-integral is well defined for arbitraryz,z⁰. Forz→∞, these contributions decay exponentially faster than the term withQ=0, so that only the latter term remains to be analyzed,

e^S_x(G,z) −−→^zL − Z

V

d²k⁰ 2π

Z _∞

−∞dz⁰e^{−|k0||z−z0|}

|k⁰|

hf(k⁰,G,z,z⁰)−f(0,G,z,z⁰)e^−µ|k0|i

. (34)

For largez, an expansion of the expression in square brackets aboutk⁰ =₀is legitimate, since:

(i) thek⁰-integral is dominated by the vicinity of the Coulomb singularity; (ii) the integrand falls off rapidly for largez; and (iii)Vcan be chosen much smaller than the first BZ,

e^S_x(G,z) −−→^zL − Z

V

d²k⁰ 2π

Z _∞

−∞dz⁰e^{−|k0||z−z0|}

|k⁰|

×

∂f(k⁰,G,z,z⁰)

∂k⁰ _k0=0

·k⁰+ f(0,G,z,z⁰)µ|k⁰|+. . .

. (35)

Due to the inversion symmetry ofV, one has:

Z

V

d²k 2π

k

|k|^e

−|k||z−z⁰| = 0. (36)

Thek⁰-integral in the second term on the right-hand side of (35) can be evaluated further in polar coordinates,

Z

V

d²k

2π e^{−|k||z−z0|} = Z _2π

0

dϕ 2π

Z _S(ϕ)

0 kdk e^−|z−z0|k

= Z _2π

0

dϕ 2π

1

|z−z⁰|²−^e

−|z−z⁰|S

|z−z⁰|^{2 1}+|z−z⁰|S

, (37)

whereS(ϕ)defines the boundary of the integration areaV. If, for instance,Vis chosen to be a square with its corners at(±d,±d),S(ϕ)is given by:

S(ϕ) =



















d

cos(ϕ) 0≤ ϕ≤ ^π₄

d sin(ϕ) π

4 ≤ϕ≤ ^3π₄

−d cos(ϕ)

3π

4 ≤ ϕ≤ ^5π₄

−d sin(ϕ)

5π

4 ≤ ϕ≤ ^7π₄

d cos(ϕ)

7π

4 ≤ ϕ≤2π .

The leading contribution toe^S_xis thus obtained as e^S_x(G,z) −−→^zL −µ

Z _∞

−∞dz⁰ f(0,G,z,z⁰)

|z−z⁰|² Z _2π

0

dϕ 2π h

1−e^−|z−z0|S 1+|z−z⁰|Si

. (38)

It seems worthwhile to emphasize that thez⁰-integral is well defined, since for z⁰ close toz, one has:

1−e^−|z−z0|S 1+|z−z⁰|S

= 1−

1− |z−z⁰|S+(|z−z⁰|S)² 2 +. . .

1+|z−z⁰|S

= (|z−z⁰|S)² 2 +. . . .

Taking this into account, the leading term for largezis given by:

e^S_x(G,z) −−→^zL −µn(G,z)

2z² , (39)

sinceS >0 for allϕand f(0,G,z,z⁰)decays exponentially withz⁰. Upon addition of (33) and (39), one finds for the leading two orders of the asymptotic energy density,

ex(G,z) −−→ −^{zL n}(G,z)

2z − ^p(G,z)

2z² . (40)

The result is independent ofµ, as it should be. In real space, one thus obtains:

ex(r) −−→ −^{zL n}(r) 2z − ^p(r)

2z² . (41)

4. Asymptotic Behavior of Slater Potential

Equation (41) shows that the Slater potential asymptotically behaves as:

v_Slater(r) = 2ex(r) n(r)

−−→ −zL ¹

z− ^p(r)

n(r)z² . (42)

In order to evaluate the next-to-leading order term further, thek-integrations involved in p(r) andn(r)need to be analyzed. Asymptotically, the density (30) is dominated by the coefficientsckα

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 coefficientsc_kα(G=0,z)decay most slowly, as discussed in Section2,

n(r) −−→^zL 2 Z

1BZ

d²k (2π)²

∑

α

Θkα|ckα(0,z)|², (43) withc_kα(0,z)given by (11), due to (14). The same argument can be applied to the asymptoticp(r),

p(r) =

∑

G

e^iG·rkp(G,z)

−−→zL 2 Z _∞

−∞dz⁰z⁰ Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘkα⁰c^∗_kα0(0,z)c_kα(0,z)

∑

G⁰⁰

c^∗_kα(G⁰⁰,z⁰)c_kα⁰(G⁰⁰,z⁰). (44) In general, one particular bandαfalls off most slowly for eachk. On the other hand, if there are two degenerate highest occupied statesαandα⁰for somek, the asymptotically-leading terms in the corresponding coefficientsckα(0,z)andc_kα⁰(0,z)only differ by some constant, sinceγ_kα(0) =γ_kα0(0). It is thus sufficient to restrict the expression (44) to the terms withα⁰=α,

p(r) −−→^zL 2 Z

1BZ

d²k (_2π)²

∑

α

Θkαm_kα|c_kα(0,z)|². (45)

Here,m_kαdenotes the first moment of|c_kα(G,z)|², m_kα =

Z _∞

−∞dz z

∑

G

|c_kα(G,z)|², (46) 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) ofc_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 integrandAkαin Equation (18) of [43] corresponds to 2|f0,kα|², for the first momentAkαcorresponds to 2m_kα|f_0,kα|². Both (43) and (45) are asymptotically dominated by the state(s) with the minimum exponentγ_kα(G=₀), as discussed in [43]. If there is only a singlek-pointqand bandβfor which γ_kα(G=0)reaches its minimum value inside the integration region, one finds:

p(r) −−→^{zL mqβ}

2π |f_0,qβ(G=₀)|²z^2u/γqβ−1e^−2γqβz[det(Γ)]^−1/2 (47) n(r) −−→^{zL 1}

2π|f_0,qβ(G=₀)|²z^2u/γqβ−1e^−2γqβz[det(_Γ)]^−1/2, (48) whereΓis defined by:

Γ = ^∂

2γ_kβ(G=0)

∂k∂k (q). (49)

p(r)therefore has the same asymptotic behavior asn(r), so thatp/n −−−→^z→∞ mqβ. The ratiop/nalso approaches a constant for largez, iff0,qβ(G=0) =0, since this affects thez-dependence ofpandnin the same way. Moreover, the same statement applies, if there is more than onek-point and/or more than one band for whichγ_kα(G=₀)assumes its minimum value (see [43]). One thus finally arrives at:

v_Slater(r) = 2ex(r) n(r)

−−→ −zL ¹

z−^mqβ

z² +. . . , (50)

withm_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 atk = qunder the transformationz→ −z,mqβcan vanish. This is the case in particular, if one has reflection symmetry with respect to thexy-plane, so that|cqβ(G,−z)|=|cqβ(G,z)|.

5. Asymptotic Behavior of Exact Exchange Potential

For a discussion of the exactvx, one starts with the total exchange energy per unit cell, expressed in terms of the amplitudesckα(G,z),

Ex = −A Z

1BZ

d²k (2π)²

d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘk⁰α⁰

∑

G Z _∞

−∞dz Z _∞

−∞dz⁰2πe^{−|k−k0+G||z−z0|}

|k−k⁰+G|

×

∑

G⁰

c^∗_k0α⁰(G⁰,z)c_kα(G⁰+G,z)

∑

G⁰⁰

c^∗_kα(G⁰⁰,z⁰)c_k⁰_α⁰(G⁰⁰−G,z⁰) (51) (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 ofc_kα(G,z), since there exists a one-to-one correspondence between the KS states and the density (30),

δn(r) = 2

∑

G

e^iG·rk Z

1BZ

d²k (2π)²

∑

α

Θkα

∑

G⁰

δc^∗_kα(G⁰,z)ckα(G⁰+G,z) +c.c. . (52)

In addition, since only the asymptotic density is to be varied (but not the completen(r)),δc_kα(G,z) is non-zero only for largez[δc_kα(G,z)can be assumed to have the standard shape of a test function].

Now, consider the variation ofExunder the variation ofckα(G,z)_, δEx = −2A

Z

1BZ

d²k (2π)²

d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘk⁰α⁰

∑

G Z _∞

−∞dz Z _∞

−∞dz⁰2πe^{−|k−k0+G||z−z0|}

|k−k⁰+G|

×

∑

G⁰⁰

c^∗_k0α⁰(G⁰⁰,z⁰)ckα(G⁰⁰+G,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)ck⁰α⁰(G⁰−G,z) +c.c. . (53) Using (19) and (20), thek⁰-integration can be combined with the sum overG. If the origin of the resultingk⁰-integration over the completek-space is chosen to be atk, one obtains:

δEx = −2A Z

1BZ

d²k (_2π)²

Z d²k⁰ (_2π)²

∑

α,α⁰

ΘkαΘk+k⁰,α⁰ Z _∞

−∞dz Z _∞

−∞dz⁰ 2πe^{−|k0||z−z0|}

|k⁰|

×

∑

G⁰⁰

c^∗_k+k0,α⁰(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_k+k⁰_,α⁰(G⁰,z) +c.c. . (54) Next,δExis split into a long- and a short-wavelength component, in analogy to (24),

δE^L_x = −2A Z

1BZ

d²k (_2π)²

Z d²k⁰ (_2π)²

∑

α,α⁰

ΘkαΘkα⁰ Z _∞

−∞dz Z _∞

−∞dz⁰2πe^{−|k0|(|z−z0|+µ)}

|k⁰|

×

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z) +c.c. (55)

δE^S_x = δEx−δE_x^L. (56)

The evaluation ofδE^L_xproceeds in straightforward fashion. One first performs thek⁰-integral, δE^L_x = −2A

Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘkα⁰ Z _∞

−∞dz Z _∞

−∞dz⁰ 1

|z−z⁰|+µ

×

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)ckα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)ckα⁰(G⁰,z) +c.c. . (57) One can then follow the argument leading from (27) to (29), sincec_kα⁰(G⁰,z)falls off exponentially withz, irrespective ofδc^∗_kα(G⁰,z),

δE^L_x −−→^zL −2A Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘ_kα^{0 Z ∞}

−∞dz Z _∞

−∞dz⁰ 1 z+µ

1+ ^z

0

z+µ+. . .

×

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)ckα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z) +c.c. . (58) First, focus on the leading order term of the expansion (L0). This term can be evaluated further by use of the orthonormality relation (5),

δE^L0_x −−→^zL −2A Z

1BZ

d²k (2π)²

∑

α

Θkα Z _∞

−∞

dz z+µ

∑

G⁰

δc^∗_kα(G⁰,z)c_kα(G⁰,z) +c.c. . (59)

One can now re-introduce thex-y-integrations and identifyδn(r), δE^L0_x −−→^zL −

Z

Ad²rk Z _∞

−∞

dz z+µ

×2

∑

G

e^iG·rk Z

1BZ

d²k (2π)²

∑

α

Θkα

∑

G⁰

δc^∗_kα(G⁰,z)c_kα(G⁰+G,z) +c.c.

= −

Z

Ad²r_k Z _∞

−∞dzδn(r)

z+µ . (60)

Finally, an expansion in powers ofµ/zis legitimate, sinceµcan be chosen to be much smaller than thez-values for whichδc_kα(G,z)is non-zero,

δE^L0_x −−→ −^zL Z

Ad²r_k Z _∞

−∞dz δn(r) z

h1−^µ z +. . .i

. (61)

Next, consider the short-wavelength component (56), δE^S_x = −2A

Z

1BZ

d²k (2π)²

Z d²k⁰ (2π)²

∑

α,α⁰

Θkα Z _∞

−∞dz Z _∞

−∞dz⁰ 2πe^{−|k0||z−z0|}

|k⁰|

×

Θk+k⁰,α⁰

∑

G⁰⁰

c^∗_k+k0,α⁰(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_k+k⁰_,α⁰(G⁰,z)

−Θ_kα⁰

∑

c^∗_kα0(G⁰⁰,z⁰)ckα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z)

e^−|k0|µ

+c.c. . (62) The arguments that lead from (26) to (35) also apply to (62). Using (36), one arrives at:

δE_x^S −−→^zL −2Aµ Z

1BZ

d²k (2π)²

Z

V

d²k⁰ (2π)²

∑

α,α⁰

ΘkαΘkα⁰ Z _∞

−∞dz Z _∞

−∞dz⁰2πe^{−|k0||z−z0|}

×

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰)

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z) +c.c. . (63) One can then employ (37) for thek⁰-integration and expand in powers of 1/z, since the remaining z⁰-integral is free of any singularities,

δE^S_x −−→^zL −2Aµ Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘkα⁰ Z _∞

−∞dz⁰

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰)

× Z _∞

−∞

dz z²

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z) +c.c. . (64) 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), δE_x^L1 −−→^zL −2A

Z

1BZ

d²k (2π)²

∑

α,α⁰

ΘkαΘ_kα^{0 Z ∞}

−∞

dz z²

∑

G⁰

δc^∗_kα(G⁰,z)c_kα⁰(G⁰,z)

× Z _∞

−∞dz⁰z⁰

∑

G⁰⁰

c^∗_kα0(G⁰⁰,z⁰)c_kα(G⁰⁰,z⁰) +c.c. . (65)

A completely general representation ofδE^L1_x in terms ofδn, as achieved forδE_x^L0, is not obvious.

In the following, therefore, only the most important class of variations will be considered. For very largez, the BZ integral in the density (30) is dominated by the vicinity of thek-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 byqβ. If only the amplitudesckβin the vicinity ofqare varied, while all otherc_kαare kept unchanged,δnis given by:

δn(r) = 2

∑

G

e^iG·rk Z

1BZ

d²k (2π)²

∑

G⁰

δc^∗_kβ(G⁰,z)c_kβ(G⁰+G,z) +c.c. . (66)

Now, considerδE_x^L1for this variation. For any givenδc^∗_kα(G⁰,z), i.e., any given statekα, the sum overα⁰in (65) is dominated by the most weakly-decaying occupied amplitude among thec_kα⁰(G⁰,z). In the case of the variation (66),δE^L1_x thus reduces to:

δE^L1_x −−→^zL −2A Z

1BZ

d²k (_2π)^2mkβ

Z _∞

−∞

dz z²

∑

G⁰

δc^∗_kβ(G⁰,z)c_kβ(G⁰,z) +c.c. , (67) withmkβdefined by Equation (46). Sinceδckβ(G⁰,z)is non-zero only for largezandckβ(G⁰,z)decays exponentially for thesez-values, thek-integration in (67) can be simplified by a Taylor expansion of the smooth parts of the integrand aboutk=q,

δE^L1_x −−→^zL −2A m_qβ Z

1BZ

d²k (2π)²

Z _∞

−∞

dz z²

∑

G⁰

δc^∗_kβ(G⁰,z)c_kβ(G⁰,z) +c.c. . (68)

At this point, one can finally re-introduce thexy-integrations and utilize (66), δE^L1_x −−→^zL −mqβ

Z

Ad²r_k Z _∞

−∞dzδn(r)

z² . (69)

Thus, for the relevant class of variations δE_x^L1 can be explicitly expressed in terms of δn.

A completely general variation ofnincludes this particular class.

One thus concludes that, in general, the next-to-leading order contribution to the asymptotic exact exchange potential is proportional to 1/z²,

vx(r) = ^δEx δn(r)

−−→ −zL ¹

z−^mqβ

z² +. . . , (70)

where (61), (64) and (69) have finally been combined.

6. Concluding Remarks

In this work, the next-to-leading order contributions to the EXX energy density, the associated Slater potential and the exact EXX potential have been calculated for arbitrary slabs, including semi-conducting and insulating materials. If the asymptotic density is dominated by the states of the bandβin the vicinity of a singlek-pointqin the interior of the first BZ, one finds:

ex(r) −−→^{zL n}(r) 2

−¹ z−^mqβ

z² +Oz⁻³

(71) v_Slater(r) −−→^zL −¹

z−^mqβ

z² +Oz⁻³

(72) vx(r) −−→^zL −¹

z−^mqβ

z² +Oz⁻³

, (73)

wheremqβ denotes the first moment (46) of the most weakly-decaying occupied state. If there are several degenerate most weakly-decaying states at thek-pointqor states at severalk-points show