The density gradient expansion of correlation functions

We present a general scheme based on nonlinear response theory to calculate the expansion of correlation functions such as the pair-correlation function or the exchange-correlation hole of an inhomogeneous many-particle system in terms of density derivatives of arbitrary order. We further derive a consistency condition that is necessary for the existence of the gradient expansion. This condition is used to carry out an infinite summation of terms involving response functions up to infinite order from which it follows that the coefficient functions of the gradient expansion can be expressed in terms the local density profile rather than the background density around which the expansion is carried out. We apply the method to the calculation of the gradient expansion of the one-particle density matrix to second order in the density gradients and recover in an alternative manner the result of Gross and Dreizler (Z. Phys. A 302, 103 (1981)) which was derived using the Kirzhnits method. The nonlinear response method is more general and avoids the turning point problem of the Kirzhnits expansion. We further give a description of the exchange hole in momentum space and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B21, 5469 (1980)) for this case. This is used to derive that the second order gradient expansion of the system averaged exchange hole satisfies the hole sum rule and to calculate the gradient coefficient of the exchange energy without the need to regularize divergent integrals.


I. INTRODUCTION
Since the groundbreaking work of Hohenberg and Kohn [1] we know that the external potential of any inhomogeneous quantum many-body system is a functional of its ground state density n. This implies that the manybody ground state |Ψ[n] and hence any ground state expectation value

A[n] = Ψ[n]|Â|Ψ[n]
(1) for any operatorÂ is a functional of the density. This idea has inspired an enormous amount of work in a research field that is now known as density-functional theory. Density-functional theory [2][3][4] is currently one of the main approaches used in electronic structure theory. Over the years many extensions beyond its original formulation have been developed and currently it is widely applied in solid state physics and quantum chemistry, both for ground state and time-dependent systems [5]. Especially a large activity has gone into finding explicit expressions for the ground state energy E[n] as a functional of the density for many-electron systems. The ground state energy is obtained by minimizing this functional on an appropriate set of electronic densities, which is of great importance in determining structures and geometries of molecules and solids. By use of the Kohn-Sham method [6] this minimization problem is reduced to solving effective one-particle equations. The construction of the effective potential in these equations requires the knowledge of the so-called exchange-correlation energy functional E xc [n]. This functional can be expressed as [7,8] E xc [n] = 1 2 drdr ′ n(r)ρ xc (r ′ |r) |r − r ′ | where ρ xc (r ′ |r) is the coupling constant integrated exchange-correlation hole. The exchange-correlation hole has a physical interpretation as the difference between the conditional and the unconditional probability (which is simply the density n(r ′ )) to find a electron at r ′ , given that we know that there is an electron at r. Equation (2) is an important relation since it expresses the exchangecorrelation energy in terms of a quantity that has a local physical interpretation and can be studied by accurate wave function and many-body methods or modeled based on physical intuition. It has therefore played an important role in the development of approximate density functionals [8][9][10][11][12][13][14]. One of the simplest and already very successful approximations has been the local density approximation (LDA) in which the exchange-correlation hole is taken from the homogeneous electron gas and applied locally [8] to systems with arbitrary density profiles. The accuracy of the LDA has been considerably improved by means of the so-called generalized gradient approximations (GGAs) [15]. A large class of commonly used GGAs is based on the so-called real-space cutoff of the straightforward gradient expansion for the exchange and correlation hole. The first such functional for exchange was proposed by Perdew [9]. He noted that the gradient expansion of the exchange hole to second order in the density gradients violates the negativity condition and the hole sum rule by which the exchange hole integrates to a deficit of one electron. By enforcing these physical constraints a density functional was obtained that greatly improved on the LDA for binding energies of molecules. This procedure was later simplified [10] using the fact that the exchange correlation energy (2) can be written as where N is the number of electrons in the system and ρ xc (y) = 1 N dr n(r)ρ xc (r + y|r), is the so-called system-averaged exchange-correlation hole. This averaged hole still obeys the negativity condition as well as the sum rule and can therefore be subjected to the real-space cut-off procedure. One advantage of using the system averaged holes was to reduce the order of the derivatives in the gradient expansion. Furthermore it also allowed for the real-space cut-off procedure to be applied to the correlation hole since there is a known gradient expansion for the system averaged correlation hole in the random phase approximation (RPA) but no known expansion for the correlation hole itself. This was the basis for a GGA for the correlation energy [3,[16][17][18]. These GGAs relied heavily on the pioneering work by Langreth, Perdew and Mehl [19][20][21] who performed a wave-vector decomposition of the system averaged hole of Eq. (3) and calculated the Fourier transform ρ xc (k) = dy ρ xc (y) e −ik·y , from which the exchange correlation energy can be calculated as We are not aware of any work that goes beyond references [19][20][21] and improves on the straightforward gradient expansion of the system averaged correlation hole. The knowledge on the gradient expansion of the correlation hole is very limited indeed. As mentioned before, in contrast to the work on the exchange hole [2,9], there are no known expressions for the gradient expansion of the non-system averaged correlation hole. Part of the problem is that there has been no clear derivation on how to expand the expectation value of a general operator as in Eq.(1) in terms of density gradients. A well-known expansion method for a two-point function is the Kirzhnits expansion [2] which is, however, specifically adapted to expanding the noninteracting oneparticle density matrix in powers of density gradients and can not be used for more general correlation functions. The first goal of this paper is, therefore, to present a scheme based on nonlinear response theory that can be used to expand general correlation functions (such as the correlation hole) in terms of density gradients. The second goal of the paper is to clarify a point that is often overlooked in carrying out gradient expansions. We will illustrate the problem by considering the standard gradient expansion for a global quantity, namely the exchange-correlation energy E xc . The starting point of any gradient expansion is the consideration of a density profile n(r) = n 0 + δn(r) which consists of a small density variation δn(r) around a homogeneous density n 0 . By considering the limit of slow density variations one then finds that the lowest gradient correction to the exchange-correlation energy is an integral over a function of the form B xc (n 0 )(∇δn(r)) 2 where the coefficient B xc (n 0 ) is calculated from the static exchange-correlation kernel of the homogeneous electron gas of density n 0 [4]. Since ∇δn(r) = ∇n(r) we can replace the density variation under the derivative operator by the full density profile n(r). However, the coefficient B xc (n 0 ) still depends on the background density n 0 . This is a problem in application of the formula to general inhomogeneous systems such as molecules or surfaces in which a background density cannot be unambiguously defined (assuming that a low order gradient expansion makes sense in such systems [22,23]). The usual argument to get rid of the background dependence, is that the replacement B xc (n 0 ) → B xc (n 0 + δn(r)) can be made since the difference between these two quantities is or order δn(r). However, it is not clear that this is consistent with gradient coefficients that arise from terms of order (δn(r)) 3 . The point was discussed clearly by Svendsen and von Barth [22] who checked that the replacement of n 0 by the full density profile was consistent to third order in the density variation. This derivation was based on a specific relation between response functions that describe the change in the exchange energy to second and third order in the density variations. In reference [24] this derivation was extended by showing that the replacement of n 0 by the full density profile is consistent to all orders in the density variation. In this paper we will show that this is the case as well for the gradient expansion of correlation functions rather than scalar functions. This relies on certain relations between the response functions that must be satisfied for the gradient expansion to exist.
The paper is divided as follows. In section II we derive the basic equations and consistency conditions for the density gradient expansion of correlation functions. In section III we derive the gradient expansion of the one-particle density matrix for a noninteracting system with density n(r) (which is therefore equal to the density matrix of the Kohn-Sham system) and discuss its symmetry properties. We further calculate the gradient expansion of the exchange-hole to second order in the density gradients, both in real and in momentum space. The momentum space description is used to demonstrate that the system averaged exchange hole satisfies the sum rule (but not the negativity condition), and to calculate the gradient expansion of the exchange energy without the need to regularize divergent integrals. Finally in section IV we present our conclusions and outlook.

A. Expansion in density variations
Let us for a system with ground state density n(r) consider an arbitrary correlation function F [n](r ′ , r). Since the correlation function can be calculated from the manybody ground state, by virtue of the Hohenberg-Kohn theorem this function is a functional of the density. At this point it will not be important what the specific form of the correlation function is, as the structure of their gradient expansions will be independent of its specific form. Correlation functions of common interest are, for instance, the pair-correlation function, the exchangecorrelation hole or the one-particle density matrix.
When the density is constant in space, i.e. when n(r) = n 0 , we are describing the homogeneous electron gas and due to translational and rotational invariance we have that F [n](r ′ , r) = F 0 (n 0 , |r − r ′ |) where F 0 is a function of the homogeneous density and the distance between its spatial arguments. When the density n(r) = n 0 + δn(r) deviates from the constant density n 0 by a small amount δn(r) we can expand F in powers of this density variation. To derive compact expressions we first introduce the notation r m = (r 1 , . . . , r m ) dr m = dr 1 . . . dr m for the collection of m position vectors and the corresponding integration volume elements. We further define δn(r m ) = m j=1 δn(r j ).
With these conventions the expansion of F in density variations is given by where we defined The n 0 dependence of the functions F m will be important for the gradient expansion, but to shorten notation we suppress this dependence in the notation for the time being and will reintroduce it later. Here we assume that the derivatives F m exist, which means that we assume that F is a smooth functional of the density. When the derivatives do exist we can expect the Taylor series (4) to converge whenever |δn(r)|/n 0 ≪ 1. The density variations are therefore assumed to be small but we do not need to put any constraint on how rapid they can vary in space and therefore their spatial derivatives can be large. Let us now look at the symmetry properties of the functions F m . As follows directly from the definition (4) of the functions F m and the assumption of differentiability, we have the permutational symmetry F m (r ′ , r, r 1 . . . r m ) = F m (r ′ , r, r π(1) . . . r π(m) ) for all permutations π of the integers 1 . . . m. Since the functions F m are evaluated at the homogeneous density n 0 they also have the spatial symmetry properties of the homogeneous electron gas. These symmetries are the translational symmetry over an arbitrary vector a, rotational symmetry under arbitrary rotations by a rotation matrix R, and inversion symmetry. If we denote by a = (a, . . . , a) Rr m = (Rr 1 , . . . , Rr m ) the m-tuples of translation vectors and rotated position vectors we can write Since in Eq.(5) the vector a is arbitrary we can in particular choose a = −r and define a new function N m depending on m + 1 independent vectors by the relation The difference vector r ′ − r will appear several times in our equations and it will therefore be convenient to define the short notation y = r ′ − r. We will further use y = |y| to denote the length of this vector. Our interest will be in the functions N m and in particular their Fourier transforms (7) with Fourier inverse With these definitions we obtain the following expansion for F in powers of the density variations This equation forms the basis of the gradient expansion. To proceed we further derive a consistency condition that is necessary for the existence of the gradient expansion. It follows directly from the definition (4) of the functions F m that δF m (r ′ , r, r m ) = dr ′′ F m+1 (r ′ , r, r ′′ , r m )δn(r ′′ ).
Taking δn(r ′′ ) = δn 0 then to be a uniform shift in the density (which means that we look at a compression or decompression of the electron gas) yields ∂F m ∂n 0 (r ′ , r, r m ) = dr ′′ F m+1 (r ′ , r, r ′′ , r m ).
We can translate this condition to a condition on N m using its definition (6) and the definition (7) of its Fourier transform. This gives the condition ∂N m ∂n 0 (y, q 1 . . . q m ) = N m+1 (y, 0, q 1 . . . q m ). (9) This relation is of key importance for the existence of the gradient expansion. Without it we would not be able to eliminate the dependence of the gradient coefficients on the homogeneous density n 0 in favor of a dependence on the actual density profile. As we will see, this requires a summation over response functions F m to infinite order. There is another important advantage of this resummation. It will allow us to relax the constraint that the density variations be small (but varying arbitrarily fast) and replace it with the constraint that the density variations vary slowly but with an arbitrarily amplitude. This is discussed in the next section.

B. The gradient expansion
In order to expand F in density gradients we have to assume that the density gradients are small. This constraint is most easily phrased in momentum space by requiring that the Fourier coefficients δn(q) of the density variation have their main contribution from small wave vectors q. If this is the case then the main contribution to the integrals in Eq.(8) comes from this region of small vectors and we can expand the functions N m in powers of the wave vectors. Subsequently carrying out the integrals in Eq.(8) then leads to an expansion in density gradients, since powers in wave vectors correspond to orders of derivatives in real space. Since the response functions N m typically vary very rapidly for wave vectors close to the Fermi surface (or multiples thereof) we require that the Fourier coefficients δn(q) of the density variation have their main contributions from wave vectors that satisfy |q| < k F where k F is the Fermi wave vector.
The next task is then to expand the functions N m in powers of wave vectors. This task is simplified when we exploit the symmetry of the these functions. As follows directly from Eqs. (6) and (7) the functions N m in Fourier space inherit the symmetry properties of F m . They are symmetric functions of their wave vectors, i.e. N m (y, q 1 . . . q m ) = N m (y, q π(1) . . . q π(m) ) for any permutation π of the integer labels, and are invariant under rotations and inversion From these equations we therefore see that any power expansion of the functions N m must lead to an expansion in terms of the spatial vector y and the wave vectors q j which is invariant under rotations and inversions of these vectors. The mathematical question is then which polynomials have these properties. This was answered in the classic book by Weyl [25]; every such polynomial must be a function of the variables y · y, y · q i and q i · q j . We see that these inner products are indeed invariant under rotations and inversions. Since any polynomial in powers of y 2 = y · y can be re-summed to a function of y we find that the general expression for N m to second order in the wave vectors q j is given by where q 2 i = q i · q i and where we took into account that this expansion must be invariant under permutations of the wave vectors q j . This expansion is readily extended to higher order powers in the wave vectors. For a given choice of correlation function F the practical task will be to determine the explicit form of the coefficient functions N m j (y) as a function of y and the background density n 0 . If we insert Eq.(10) into Eq. (8) and Fourier transform back to real space we obtain the expansion Let us analyze the explicit form of the first five coefficients A m j for j = 0, .., 4 of this expansion, since these are exactly the terms that we need for a gradient expansion to second order in the density derivatives. The first term for j = 0 is simply given by The terms with j = 1, 2, 3 in Eq.(11) involve according to Eq.(10) the Fourier transforms of m symmetry equivalent terms and acquire the following forms in real space In the derivation of the last term we used that for an arbitrary function f where with r = (x 1 , x 2 , x 3 ) we used the notation ∂ i = ∂/∂x i as well as y i = x ′ i − x i . This expression is readily checked using Finally the terms with j = 4, 5 in Eq.(11) involve according to Eq.(10) the Fourier transforms of m(m − 1)/2 symmetry equivalent terms which yields We see that a general coefficient function A m j depends on the density variation and gradients thereof. The functions δn(r) that appear under the gradient operators can be replaced by the actual density profile n(r) = n 0 + δn(r). However, we see that we are still left with powers of δn(r) which are unprotected by derivative operators and which therefore can not be replaced by the full density profile n(r). It is precisely at this point that the consistency conditions (9) play a crucial role. If we use these conditions in Eq.(10) then we see that This equation relates certain coefficients coming from higher order response functions N m to those of lower order ones. In particular it tells us that by iteration If we insert these expressions together with the explicit form of the coefficient functions A m j back into Eq.(11) in Eq. (11) and shift some indices we obtain the expansion In this expression we recognize several Taylor series of the coefficient functions N m j (n 0 + δn(r), y) in powers of δn(r). These Taylor series can be re-summed to finally give F (r ′ , r) = F 0 (n(r), y) + iN 1 1 (n(r), y) y · ∇n(r) − N 1 2 (n(r), y) ∇ 2 n(r) − N 1 3 (n(r), y) y 2 ( where the implicit dependence on n 0 of the coefficient functions is now replaced by a dependence on the full density profile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstated the explicit density dependence of the coefficients in Eq. (14) to indicate that they are local functions of the density n(r) and therefore have a nontrivial spatial dependence on both r and y. We stress again that the derivative condition (9) was essential in eliminating the dependence of the gradient coefficient on the reference density n 0 in favor of a dependence on the full density profile.
Having obtained the general gradient expansion Eq. (14) it remains to find explicit expressions for the coefficient functions N m j . We will describe how to do this in detail for the coefficients needed for a gradient expansion to second order in the density derivatives. To calculate the coefficients N 1 1 , N 1 2 and N 1 3 we need to calculate the function N 1 (n 0 , y, q) for the homogeneous electron gas and expand it in powers of q: The determination of the coefficients N 2 4 and N 2 5 requires knowledge of the second order response function N 2 (n 0 , y, q 1 , q 2 ) = N 2 0 (n 0 , y) + N 2 1 (n 0 , y)y · (q 1 + q 2 ) + N 2 2 (n 0 , y)(q 2 1 + q 2 2 ) + N 2 3 (n 0 , y)((y · q 1 ) 2 + (y · q 2 ) 2 ) + N 2 4 (n 0 , y)(q 1 · q 2 ) + N 2 5 (n 0 , y)(y · q 1 )(y · q 2 ) + . . . (16) The expression (14) together with Eqs. (15) and (16) determine completely the gradient expansion of an arbitrary correlation function F to second order in the density gradients. These equations form the main result of the present work. If expansions to higher order gradients are required they can be derived without difficulty along the lines described above. To do this one needs to construct higher order symmetric polynomials in the wave vectors and carry out the required Fourier transforms. The consistency conditions Eq.(9) or equivalently Eq.(12) then allow for a complete re-summation and leads to gradient coefficients depending on n(r). What is needed to determine the explicit form of the gradient coefficients in practice is the determination of the functions N m . How to do this for N 1 and N 2 is described in the next section.
C. Determination of N 1 and N 2 A practical calculation of the coefficients of the gradient expansion of Eq. (14) requires the knowledge of the the functions N 1 and N 2 . According to Eqs.(4) and (6) these are defined as the first and second density derivatives of the correlation function F with respect to the density, i.e.
To derive useful expressions for these functions it is convenient to transform derivatives with respect to the density to derivatives with respect to the external potential as such derivatives appear naturally in perturbation theory. We can then use the well-developed tools of perturbation theory to calculate these functions. Let us there-fore define the new functions Since we evaluate these functions for the homogeneous electron gas they only depend on differences between the spatial coordinates. It will therefore be convenient to also define their Fourier transforms by as well as their partial Fourier transforms The chain rule of differentiation gives a connection between the derivatives with respect to the density and the derivatives with respect to the potential of the correlation function F . We have δ 2 F (r 1 , r 2 ) δv(r 3 )δv(r 4 ) = dr 5 δF (r 1 , r 2 ) δn(r 5 ) δ 2 n(r 5 ) δv(r 3 )δv(r 4 ) + dr 5 dr 6 δ 2 F (r 1 , r 2 ) δn(r 5 )δn(r 6 ) We see that in these expressions the first and second order derivative of the density with respect to the potential appear. We therefore define the linear density response function by χ as well as the second order density response function χ 2 by Using these definitions we find by Fourier transforming Eqs. (19) and (20) that N 2 (y, p, q) = F 2 (y, p, q) − N 1 (y, p + q)χ 2 (p, q) These equations give the desired relation between the density derivatives N 1 and N 2 and the potential derivatives F 1 and F 2 of the correlation function F . Relations between the higher order derivatives N m and F m can be derived in a completely analogous way. From Eqs. (21) and (22) we see that to calculate the functions N 1 and N 2 , and hence the gradient expansion coefficients of F to second order in the gradients, we need to calculate the density response and the change in F to second order in a perturbing potential δv. The problem is thereby reduced to doing perturbation theory on the electron gas. This is a well-developed field of research. One option is to use diagrammatic perturbation theory [26]. In the following section we will use standard perturbation theory to obtain these response function for the case that F is the one-particle density matrix and use that to calculate the gradient expansion of the exchange hole.

III. GRADIENT EXPANSION OF THE ONE-PARTICLE DENSITY MATRIX AND THE EXCHANGE HOLE
A. The one-particle density matrix As an application of our formalism we carry out the gradient expansion of the one-particle density matrix for a noninteracting system with density n(r). This problem has received large interest in the past since both the exchange energy E x [n] as well the Kohn-Sham kinetic energy T s [n] are explicit functionals of such a noninteracting density matrix [2,3]. Therefore a gradient expansion of this density matrix directly leads to a gradient expansion of these two functionals. Such a gradient expansion was first carried out by Gross and Dreizler [27] using the Kirzhnits expansion [2,28]. This expansion is adjusted to the specific form of the noninteracting oneparticle density matrix and can not be generalized to arbitrary correlation functions. The method presented in this work can, on the other hand, be applied to arbitrary correlation functions, but to demonstrate the soundness of our formalism we will use it to provide an alternative derivation of the result obtained earlier by Gross and Dreizler using the Kirzhnits expansion. An advantage of our derivation based on nonlinear response theory is that it avoids the turning point problem encountered in the Kirzhnits approach [2].
Let us start by defining the one-particle density matrix in second quantization as where σ is a spin coordinate. We will consider the case of spin-compensated systems. Then, in a noninteracting system with 2N electrons, the density matrix is given in terms of one-particle wave functions ϕ j (r) by where the pre-factor 2 is a spin factor. The one-particle states are eigenstates of a one-particle Schrödinger equation where v(r) is the external potential. To determine the gradient expansion coefficients of γ we need to determine how the density matrix changes if we make small changes v(r) → v(r) + δv(r) in the external potential. More precisely, we need to calculate the functional derivatives which are the equivalents of the functions F 1 and F 2 of Eqs. (17) and (18) for the specific choice of correlation function F = γ (note that they become functions of two and three independent vectors respectively when evaluated for a homogeneous system). To do this it is sufficient to know the functional derivatives of the orbitals ϕ j and eigenvalues ǫ j with respect to the potential v. These are simply given by doing first order perturbation theory on the one-particle Schrödinger Eq. (24). We find With these relations we can differentiate Eq.(23) twice with respect to the potential. Afterwards we then need to insert the plane wave states ϕ k (r) and their eigenvalues ǫ k appropriate for the homogeneous electron gas into the final expressions. In order not to interrupt the presentation these calculations are given in Appendix A and we simply present the results of these calculations here.
If we define the Fermi factors by where θ is the Heaviside function and where k F = (3π 2 n 0 ) 1/3 is the Fermi wave vector, then the resulting expressions are given by In this expression ǫ k = |k| 2 /2 are the one-particle energies. It remains to calculate first and second order density response functions χ(q) and χ 2 (p, q) to evaluate the functions N 1 and N 2 of Eqs. (21) and (22). Fortunately, for the specific choice F = γ that we made these density response functions are simply given by They are therefore obtained directly from Eqs. (27) and (28) so that we do not need to do any extra work to calculate them.
To calculate the gradient coefficient to second order in the gradients we now need to expand γ 1 and γ 2 to second order in the wave vectors p and q. Since these calculations are somewhat lengthy we do not present them here and refer to Appendix B for a description of the main steps. The resulting expansions are given by where we defined z = k F y and q = |q|. The expansions for the first and second order response functions χ and χ 2 of Eqs. (30) and (31) are obtained by taking y = 0 in these expressions. Together which Eqs. (32) and (33) these expression then immediately determine the functions N 1 and N 2 from Eqs. (21) and (22). The final result of this calculation to second order in the wave vectors is given by By comparison of these two expression to the expansions (15) and (16) we can directly read off the gradient coefficient functions of Eq. (14) to be where now we replaced n 0 → n(r) and hence the Fermi wave vector k F = (3π 2 n(r)) 1/3 that appears in these expressions is a spatially varying function. We further introduced the spherical Bessel functions By inserting the gradient coefficient function into Eq. (14) we find the following explicit gradient expansion for the one-particle density matrix This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent form as This expression becomes identical in form to that derived by Gross and Dreizler [27] using the Kirzhnits method after eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(19) of reference [27] ). We close this section with a final remark on our result. We note that the gradient expansion to finite order breaks the symmetry γ(r, r ′ ) = γ(r ′ , r) * . However, the full gradient expansion as described by Eq.(11) has this symmetry, which means that the symmetry is restored by taking all higher order gradients into account whenever the series converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such as the exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literature we will in the next section adopt the definition of the exchange hole that was used by Perdew [9].

B. The exchange hole and energy
We finally stress a few points on the properties of the exchange hole as derived from the gradient expansion and present an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler and Gross which has the advantage of avoiding the regularization of divergent integrals [27]. In doing so we confirm a result obtained by Langreth and Perdew [19]. The exchange hole can be calculated directly from the one-particle density matrix as Inserting the expansion (34) of for the one-particle density matrix and retaining terms to second order in the gradients yields the expression As was pointed out by Perdew [9] the second order gradient expansion of the exchange hole does not satisfy the negativity and sum rule constraints The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediately obvious from Eq. (36). The sum rule can, however, be conveniently analyzed in momentum space. To do this we define the Fourier transformed exchange hole to be This function has the form where we defined the unit vectork = k/k and k = |k|. The explicit form of the functions ρ 0 x and α j follow directly by Fourier transforming Eq. (36) and are given in Appendix C. The coefficient functions α j are tempered distributions [29] which have mathematically well-defined Fourier transforms. As follows directly from Eq.(37) the sum rule condition in momentum space translates to the requirement that ρ x (k = 0|r) = −1. Since ρ 0 x (k = 0|r) = −1 (see Eq.(C1) ) this implies that the sum rule would be fullfilled whenever α j (k = 0, n) = 0. However, as is clear from Eqs.(C2)-(C6) in the Appendix this is not satisfied for α 1 , whereas α 3 and α 5 even diverge for k → 0 when interpreted locally as functions rather than distributions. The gradient expansion of the exchange hole does therefore indeed not satisfy the sum rule. These divergencies, however, cancel when we calculate the system averaged exchange hole in momentum space which is given by the expression dr n(r)β ij (k, n(r))∂ i n(r)∂ j n(r).
where N is the number of electrons. In this expression ρ 0 x (k) is the system average of ρ 0 x (k|r) and in the last term under the integral sign we defined the tensor This tensor has a longitudinal part with coefficient α L and a transverse part with coefficient α T which describe the contributions to the system-averaged hole of density gradients ∇n parallel and perpendicular to the momentum vector k. These coefficients are calculated from the functions α j as as a short calculation on the basis of Eq.(38) will show. From these equations and the explicit expressions for α j given in Appendix C we can readily calculate the explicit expressions for α L,T . If we define the dimensionless vari-ablek = k/(2k F ) then we have where the functions Z L,T have the explicit form These expressions are in accordance with the results of Langreth and Perdew (see Eq.(3.55) of reference [19]) and this independent result therefore confirms the correctness of our alternative derivation. The interesting observation is now that and hence β ij (k, n) → 0 when k → 0. This implies that the last term in Eq. (39) does not contribute to the sum rule of the system averaged hole. The same conclusion can be derived for the second term in Eq. (39). Since ) is a local function of the density the integrand of the first integral in Eq.(39) is a total derivative and vanishes (assuming either vanishing densities at infinity or periodic boundary conditions). We therefore find that This implies that system averaged exchange hole obtained from the second order gradient expansion does fulfill the sum rule, although the exchange hole itself does not. We note, however, that this is only achieved by integrating over both positive and negative contributions. When the positive contributions to the exchange hole are removed the sum rule is only recovered for a finite hole radius [9,10]. We finally note that with expression Eq. (39) it is now a simple matter to calculate the exchange energy from We insert Eq. (39) and do the angular integrations first. It is therefore convenient to define the spherical and system averaged hole in momentum space by such that From Eq.(3) we then find that where we defined The exchange energy is obtained by inserting the expression for the averaged hole of Eq.(41) into Eq. (40). This then finally gives the expression where and where to calculate the first term in Eq.(43) we used Eq.(C1). The coefficient B x is the same gradient coefficient as obtained by Gross and Dreizler [27] and earlier by Sham [30]. However, the correct analytic exchange gradient coefficient is known [31][32][33][34][35] to be a factor 10/7 larger. The reason for this discrepancy is clearly described by Svendsen and von Barth [34] who showed that the Coulomb interaction is too singular to allow for the interchange of the operations of integration and the expansion in wave vectors. The problem does, for instance, not appear when Yukawa screened Coulomb interactions are used [34]. The conclusion is that there is nothing wrong with the Kirzhnits method, or the nonlinear response theory derivation of the gradient expansion of the exchange hole presented here, but one has be aware that for Coulomb interactions the procedures of expanding correlation functions in terms of density gradients followed by integrations may not yield the same result as directly expanding the integrated quantity in terms of density gradients. For this reason the original GGA of Perdew and Wang [10] based on the gradient expansion of the exchange hole was later reparametrized [16,36] to accomodate the correct gradient coefficient for the exchange energy.

IV. CONCLUSIONS AND OUTLOOK
We derived a general scheme based on nonlinear response theory to calculate the density gradient expansion of general correlation functions and showed that in order to express the gradient coefficients in terms of the full density profile summations to infinite order must be carried out over response functions of arbitrarily large order. A consistency condition was derived to do this. The formalism was used to derive the gradient expansion of the one-particle density matrix and the exchange hole to second order in the gradients. We confirm the derivation of Dreizler and Gross that used the Kirzhnits expansion method. We further analyzed the exchange hole in momentum space to derive that the system averaged hole to second order in the gradients satisfies the sum rule and to derive the gradient expansion of the exchange energy without the need to regularize divergent integrals.
The scheme that we presented is very general and can be applied to more general correlation functions. A natural application of the scheme would be to calculate the gradient expansion of the correlation hole ρ c . Regarding the gradient terms of the correlation hole little is known. We essentially only have some detailed information on the long-range properties of the spherically and system averaged hole. This information comes from the work of Langreth and Perdew that calculated ρ c (k) within the RPA. In our notation their results (see also Appendix C of [19]) can be summarized as where the function z c has been parametrized by Langreth and Mehl [20,21] as where k s = (4k F /π) 1/2 is the Thomas-Fermi or screening wave vector. We can transform to real space to obtain the following expression for the spherically and system averaged correlation hole We see from Eq.(44) that z c (n, k = 0) = 0 and as a consequence the sum rule ρ c (0) = 0 for the correlation hole is not satisfied. We know, however, that the RPA becomes accurate in the high density limit and since from Eq.(45) the gradient coefficient of the correlation hole depends on k s y we see that high densities are equivalent to large separations y. Similarly, low density corresponds to short-distance behavior. However, RPA can not be trusted in this region and we have no precise knowledge on the small distance behavior of the gradient coefficient of ρ c (y) . A model for the general short-range behavior was proposed [16][17][18] on physical arguments (it should not affect the Coulomb cusp of and the on-top value of the LDA hole) after which the real-space cutoff procedure was applied (see e.g. Fig. 4 in [17] and Fig. 5 in [18]) to obtain a GGA for correlation. It would, however, be desirable to have a first principles approach to the calculation of the short range properties of the gradient coefficient of the correlation hole. As follows from our derivation this requires the knowledge of the functions (17) and (18), or at least their expansion to second order in wave vectors, for the case that F represents the pair correlation function or the correlation hole of the electron gas. It may therefore be worthwhile to use our current scheme to explore these response functions beyond the RPA. For short range correlation an approach based on diagrammatic summation of ladder diagrams suggests itself. Of course, for the development of density functionals for general systems beyond the weakly inhomogeneous regime it is not sufficient to use the straightforward gradient expansion. However, general short-range features, such as the way the exchange-correlation hole deforms close to the reference electron in the presence of a density gradient can be transferred to more general systems than the weakly inhomogeneous ones. Perhaps in the combination with truly nonlocal functionals for the exchangecorrelation hole [37] this can lead to the development of more accurate density functionals. This will be a topic of our future research. Acknowledgments I would like to thank the students in the course "Density Functional Theory" at the University of Jyväskylä for inspiration and Esa Räsänen for useful discussions and Klaas Giesbertz for checking a large part of the equations .
Appendix A: Calculation of γ 1 and γ 2 In this Appendix we will derive in the expression for the first and second derivative of the one-particle density matrix γ with respect to the potential v. For a clearer interpretation and compactness of notation we will only put a comma between the variables representing the original coordinates of the density matrix and the variables that appear as argument of the potential variations. Direct differentiation of Eq.(23) using Eq.(26) then gives where we used the short notation j = r j and introduced the occupation numbers f j = 1 for an occupied state and f j = 0 for an unoccupied state. Inserting plane wave states appropriate for the electron gas we find that γ 1 (r 2 r 3 , r 1 ) = 2 dqdp (2π) 6 n p+q − n p ǫ p+q − ǫ p e ip·(r2−r3)+iq·(r1−r3) .
To obtain the explicit form of this function we have to differentiate both the orbitals and the eigenvalues using Eqs. (25) and (26). Let us denote the term that arises from the differentiating the orbitals by X and the term that arises from differentiating the eigenvalues by Y . Then we find that and where we further defined The function Φ has the useful properties Φ(ijj) = 0 and Φ(ijk) = Φ(ikj). The function γ 2 (23, 14) = γ 2 (23, 41) is symmetric in the indices 4 and 1 due to the fact that the differentiations with respect to the potential commute. This is, however, not obvious from Eqs.(A2) and (A3). We therefore want to rewrite the form of the function γ 2 in such a way that this symmetry is obvious. To do this we first note that since Φ(iji) = −Φ(jij) the terms with k = i and k = j in Eq.(A2) sum up to a function Z of a similar form as Y in Eq.(A3), namely We can therefore write γ 2 (23,14) as It is easily seen that the sum of Y and Z is symmetric under interchange of 1 and 4. However, this is not obvious in the first term of the equation above since at first sight Φ(jik) does not appear to be equal to Φ(ijk) since which seems to be different from expression (A4). However, this is just appearance. The reason is that the occupations can only attain the values zero and one. For the case f i = f j = 1 or f i = f j = 0 we see directly that Eq.(A4) and (A7) attain the same value. A little inspection shows that this is also true for cases f i = 0, f j = 1 and f i = 1, f j = 0. We therefore find for k = i, j that Φ(ijk) = Φ(jik). Therefore Eq.(A6) can be simplified to This expression is now explicitly symmetric in the indices 1 and 4. Let us now evaluate this expression for the homogeneous electron gas. In the electron gas the oneparticle states are plane waves where V is the volume of the system. Since |ϕ k | 2 = 1/V it follows from Eqs.(A3) and (A5) that the terms Y and Z are identically zero and the function γ(23, 14) of Eq.(A8) therefore attains the form where we defined In this expression n p = θ(ǫ F −ǫ p ) is the zero-temperature Fermi function and ǫ F = k 2 F /2 is the Fermi energy. When we replace the sum by an integral the restriction q = k, p is a region of measure zero and we can therefore freely integrate over all wave vectors. After a few substitutions we obtain where we defined This gives the integrand of Eq.(28).

Expansion of γ 1
To determine γ 1 to second order in the wave vectors we have to expand the function to second order in powers of q. This can be done by expanding the integrand in powers of q. If we denote ∆ = ǫ p+q − ǫ p = p · q + q 2 /2 with q = |q|, then we can expand the Fermi function n p+q in a distributional Taylor series as Since n p = θ(ǫ F − ǫ p ) is the Heaviside function we have Inserting this into Eq.(B1) then gives We see that in these integrals several derivatives of the delta function appear. We now use the general mathematical expression where the sum runs over all zeros of the function y(x), i.e. y(x i ) = 0. Using this equation in spherical coordinates with y(p) = (k 2 F − p 2 )/2 we find where in Eq.(B6) we also evaluated the third derivative of the delta function as we will need it later in the expansion of γ 2 . With Eqs.(B3)-(B5) we can readily evaluate the three integrals in Eq.(B2). Let us denote these integrals by A, B and C. Then we find for the first integral where in the second step we used Eq.(B4). It therefore remains to calculate the last integral However, up to order q 2 it is suffices to calculate where in the second step we used Eq.(B5). Adding the results γ 1 = A + B + C gives the expression (32).

Expansion of γ 2
Analogously to γ 1 we can expand γ 2 (y, p, q) of Eq.(28) in powers of p and q. To reduce our effort we note that we can write γ 2 (y, p, q) = A(y, p, q) + A(y, q, p) where A(y, p, q) = 2 dk (2π) 3 e ik·y Φ(k, k + p + q, k + p) We therefore only need to expand A(y, p, q) in powers of the wave vectors and symmetrize with respect to p and q afterwards. To do this we first need to expand the function If we denote we find by expanding the Fermi functions in a distributional Taylor series that With this expression we find that the function A(y, p, q) can be calculated as The three integrals appearing on the left hand side of this expression are sufficient to extract the powers to second order in p and q. Let us call these three integrals A 1 , A 2 and A 3 in order of appearance. The first integral gives The second integral is given by where we used Eq.(B5) in the second step. Finally the third integral has the explicit form However, we only need the terms to second order in p and q. Therefore, it is sufficient to calculate (p · ∇ y ) 2 + (p · ∇ y )((p + q) · ∇ y ) + ((p + q) · ∇ y ) 2 3 cos(z) − z 2 cos(z) + 3z sin(z) 2π 2 k 5 F where to evaluate the integral over the delta function we used Eq.(B6). To perform the derivatives in the last term we use that for an arbitrary function f (z) With this expression we find that Collecting our results A = A 1 + A 2 + A 3 we therefore obtain the following expression for γ 2 (y, p, q) γ 2 (y, p, q) = A(y, p, q) + A(y, q, p) This is equal to expression (33).

Appendix C: Exchange hole in momentum space
Here we present the explicit expressions for the coefficient functions α j of Eq. (38) which are calculated by Fourier transforming Eq.(36). If we define the dimensionless quantityk = |k|/(2k F (r)) then ρ 0 x is given by which is simply the Fourier transform of the LDA hole. We see that ρ x (0|r) = −1 and therefore the LDA hole satisfies the sum rule. The functions α j are given by We present a general scheme based on nonlinear response theory to calculate the expansion of correlation functions such as the pair-correlation function or the exchange-correlation hole of an inhomogeneous many-particle system in terms of density derivatives of arbitrary order. We further derive a consistency condition that is necessary for the existence of the gradient expansion. This condition is used to carry out an infinite summation of terms involving response functions up to infinite order from which it follows that the coefficient functions of the gradient expansion can be expressed in terms the local density profile rather than the background density around which the expansion is carried out. We apply the method to the calculation of the gradient expansion of the one-particle density matrix to second order in the density gradients and recover in an alternative manner the result of Gross and Dreizler (Z. Phys. A 302, 103 (1981)) which was derived using the Kirzhnits method. The nonlinear response method is more general and avoids the turning point problem of the Kirzhnits expansion. We further give a description of the exchange hole in momentum space and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B21, 5469 (1980)) for this case. This is used to derive that the second order gradient expansion of the system averaged exchange hole satisfies the hole sum rule and to calculate the gradient coefficient of the exchange energy without the need to regularize divergent integrals.

I. INTRODUCTION
Since the groundbreaking work of Hohenberg and Kohn [1] we know that the external potential of any inhomogeneous quantum many-body system is a functional of its ground state density n. This implies that the manybody ground state |Ψ[n] and hence any ground state expectation value

A[n] = Ψ[n]|Â|Ψ[n]
(1) for any operatorÂ is a functional of the density. This idea has inspired an enormous amount of work in a research field that is now known as density-functional theory. Density-functional theory [2][3][4] is currently one of the main approaches used in electronic structure theory. Over the years many extensions beyond its original formulation have been developed and currently it is widely applied in solid state physics and quantum chemistry, both for ground state and time-dependent systems [5]. Especially a large activity has gone into finding explicit expressions for the ground state energy E[n] as a functional of the density for many-electron systems. The ground state energy is obtained by minimizing this functional on an appropriate set of electronic densities, which is of great importance in determining structures and geometries of molecules and solids. By use of the Kohn-Sham method [6] this minimization problem is reduced to solving effective one-particle equations. The construction of the effective potential in these equations requires the knowledge of the so-called exchange-correlation energy functional E xc [n]. This functional can be expressed as [7,8] where ρ xc (r ′ |r) is the coupling constant integrated exchange-correlation hole. The exchange-correlation hole has a physical interpretation as the difference between the conditional and the unconditional probability (which is simply the density n(r ′ )) to find a electron at r ′ , given that we know that there is an electron at r. Equation (2) is an important relation since it expresses the exchangecorrelation energy in terms of a quantity that has a local physical interpretation and can be studied by accurate wave function and many-body methods or modeled based on physical intuition. It has therefore played an important role in the development of approximate density functionals [8][9][10][11][12][13][14]. One of the simplest and already very successful approximations has been the local density approximation (LDA) in which the exchange-correlation hole is taken from the homogeneous electron gas and applied locally [8] to systems with arbitrary density profiles. The accuracy of the LDA has been considerably improved by means of the so-called generalized gradient approximations (GGAs) [15]. A large class of commonly used GGAs is based on the so-called real-space cutoff of the straightforward gradient expansion for the exchange and correlation hole. The first such functional for exchange was proposed by Perdew [9]. He noted that the gradient expansion of the exchange hole to second order in the density gradients violates the negativity condition and the hole sum rule by which the exchange hole integrates to a deficit of one electron. By enforcing these physical constraints a density functional was obtained that greatly improved on the LDA for binding energies of molecules. This procedure was later simplified [10] using the fact that the exchange correlation energy (2) can be written as where N is the number of electrons in the system and ρ xc (y) = 1 N dr n(r)ρ xc (r + y|r), is the so-called system-averaged exchange-correlation hole. This averaged hole still obeys the negativity condition as well as the sum rule and can therefore be subjected to the real-space cut-off procedure. One advantage of using the system averaged holes was to reduce the order of the derivatives in the gradient expansion. Furthermore it also allowed for the real-space cut-off procedure to be applied to the correlation hole since there is a known gradient expansion for the system averaged correlation hole in the random phase approximation (RPA) but no known expansion for the correlation hole itself. This was the basis for a GGA for the correlation energy [3,[16][17][18]. These GGAs relied heavily on the pioneering work by Langreth, Perdew and Mehl [19][20][21] who performed a wave-vector decomposition of the system averaged hole of Eq. (3) and calculated the Fourier transform ρ xc (k) = dy ρ xc (y) e −ik·y , from which the exchange correlation energy can be calculated as We are not aware of any work that goes beyond references [19][20][21] and improves on the straightforward gradient expansion of the system averaged correlation hole. The knowledge on the gradient expansion of the correlation hole is very limited indeed. As mentioned before, in contrast to the work on the exchange hole [2,9], there are no known expressions for the gradient expansion of the non-system averaged correlation hole. Part of the problem is that there has been no clear derivation on how to expand the expectation value of a general operator as in Eq.(1) in terms of density gradients. A well-known expansion method for a two-point function is the Kirzhnits expansion [2] which is, however, specifically adapted to expanding the noninteracting oneparticle density matrix in powers of density gradients and can not be used for more general correlation functions. The first goal of this paper is, therefore, to present a scheme based on nonlinear response theory that can be used to expand general correlation functions (such as the correlation hole) in terms of density gradients. The second goal of the paper is to clarify a point that is often overlooked in carrying out gradient expansions. We will illustrate the problem by considering the standard gradient expansion for a global quantity, namely the exchange-correlation energy E xc . The starting point of any gradient expansion is the consideration of a density profile n(r) = n 0 + δn(r) which consists of a small density variation δn(r) around a homogeneous density n 0 . By considering the limit of slow density variations one then finds that the lowest gradient correction to the exchange-correlation energy is an integral over a function of the form B xc (n 0 )(∇δn(r)) 2 where the coefficient B xc (n 0 ) is calculated from the static exchange-correlation kernel of the homogeneous electron gas of density n 0 [4]. Since ∇δn(r) = ∇n(r) we can replace the density variation under the derivative operator by the full density profile n(r). However, the coefficient B xc (n 0 ) still depends on the background density n 0 . This is a problem in application of the formula to general inhomogeneous systems such as molecules or surfaces in which a background density cannot be unambiguously defined (assuming that a low order gradient expansion makes sense in such systems [22,23]). The usual argument to get rid of the background dependence, is that the replacement B xc (n 0 ) → B xc (n 0 + δn(r)) can be made since the difference between these two quantities is or order δn(r). However, it is not clear that this is consistent with gradient coefficients that arise from terms of order (δn(r)) 3 . The point was discussed clearly by Svendsen and von Barth [22] who checked that the replacement of n 0 by the full density profile was consistent to third order in the density variation. This derivation was based on a specific relation between response functions that describe the change in the exchange energy to second and third order in the density variations. In reference [24] this derivation was extended by showing that the replacement of n 0 by the full density profile is consistent to all orders in the density variation. In this paper we will show that this is the case as well for the gradient expansion of correlation functions rather than scalar functions. This relies on certain relations between the response functions that must be satisfied for the gradient expansion to exist.
The paper is divided as follows. In section II we derive the basic equations and consistency conditions for the density gradient expansion of correlation functions. In section III we derive the gradient expansion of the one-particle density matrix for a noninteracting system with density n(r) (which is therefore equal to the density matrix of the Kohn-Sham system) and discuss its symmetry properties. We further calculate the gradient expansion of the exchange-hole to second order in the density gradients, both in real and in momentum space. The momentum space description is used to demonstrate that the system averaged exchange hole satisfies the sum rule (but not the negativity condition), and to calculate the gradient expansion of the exchange energy without the need to regularize divergent integrals. Finally in section IV we present our conclusions and outlook.

A. Expansion in density variations
Let us for a system with ground state density n(r) consider an arbitrary correlation function F [n](r ′ , r). Since the correlation function can be calculated from the manybody ground state, by virtue of the Hohenberg-Kohn theorem this function is a functional of the density. At this point it will not be important what the specific form of the correlation function is, as the structure of its gradient expansion will be independent of its specific form. Correlation functions of common interest are, for instance, the pair-correlation function, the exchange-correlation hole or the one-particle density matrix.
When the density is constant in space, i.e. when n(r) = n 0 , we are describing the homogeneous electron gas and due to translational and rotational invariance we have that F [n](r ′ , r) = F 0 (n 0 , |r − r ′ |) where F 0 is a function of the homogeneous density and the distance between its spatial arguments. When the density n(r) = n 0 + δn(r) deviates from the constant density n 0 by a small amount δn(r) we can expand F in powers of this density variation. To derive compact expressions we first introduce the notation With these conventions the expansion of F in density variations is given by where we defined The n 0 dependence of the functions F m will be important for the gradient expansion, but to shorten notation we suppress this dependence in the notation for the time being and will reintroduce it later. Here we assume that the derivatives F m exist, which means that we assume that F is a smooth functional of the density. In the electron gas two Taylor expansions around two different values of n 0 have very likely the same mathematical convergence properties since the two systems have the same physical properties (unless we are close to a very low density corresponding to a Wigner crystal transition). The point n 0 is therefore unlikely to be a non-analytic point. When the derivatives F m do exist we can therefore expect the Taylor series (4) to converge whenever |δn(r)|/n 0 ≪ 1 or when they are integrable with a small integral norm. The density variations are therefore assumed to be small but we do not need to put any constraint on how rapid they can vary in space and therefore their spatial derivatives can be large. We note, however, that the notion of the existence of the functional derivative is closely related to the allowed domain of densities and the norm defined on this function space. The smallness or integrability of |δn(r)| implies the use of a supremum norm or possibly an L p -integral norm (i.e. requiring |δn(r)| p to be integrable). Other norms may induce weaker constraints on the derivative functions F m but stronger constraints on the density variations δn(r). A rigorous discussion on the issue of functional differentiability for the case of the energy functional is given in Refs. [25][26][27][28].
Let us now look at the symmetry properties of the functions F m . As follows directly from the definition (4) of the functions F m and the assumption of differentiability, we have the permutational symmetry The difference vector r ′ − r will appear several times in our equations and it will therefore be convenient to define the short notation y = r ′ − r. We will further use y = |y| to denote the length of this vector. Our interest will be in the functions N m and in particular their Fourier transforms N m (y, q m ) = dr m N m (y, r m ) e −iq1·r1−...−iqm·rm (7) with Fourier inverse N m (n 0 ; y, r m ) = dq m (2π) 3m N m (y, q m )e iq1·r1+...+iqm·rm .
With these definitions we obtain the following expansion for F in powers of the density variations This equation forms the basis of the gradient expansion. To proceed we further derive a consistency condition that is necessary for the existence of the gradient expansion. It follows directly from the definition (4) of the functions F m that δF m (r ′ , r, r m ) = dr ′′ F m+1 (r ′ , r, r ′′ , r m )δn(r ′′ ).
Taking δn(r ′′ ) = δn 0 then to be a uniform shift in the density (which means that we look at a compression or decompression of the electron gas) yields We can translate this condition to a condition on N m using its definition (6) and the definition (7) of its Fourier transform. This gives the condition ∂N m ∂n 0 (y, q 1 . . . q m ) = N m+1 (y, 0, q 1 . . . q m ). (9) This relation is of key importance for the existence of the gradient expansion. Without it we would not be able to eliminate the dependence of the gradient coefficients on the homogeneous density n 0 in favor of a dependence on the actual density profile. As we will see, this requires a summation over response functions F m to infinite order. There is another important advantage of this resummation. It will allow us to relax the constraint that the density variations be small (but varying arbitrarily fast) and replace it with the constraint that the density variations vary slowly but with an arbitrary amplitude. This is discussed in the next section.

B. The gradient expansion
In order to expand F in density gradients we have to assume that the density gradients are small. This constraint is most easily phrased in momentum space by requiring that the Fourier coefficients δn(q) of the density variation have their main contribution from small wave vectors q. If this is the case then the main contribution to the integrals in Eq.(8) comes from this region of small vectors and we can expand the functions N m in powers of the wave vectors. Subsequently carrying out the integrals in Eq.(8) then leads to an expansion in density gradients, since powers in wave vectors correspond to orders of derivatives in real space. Since the response functions N m typically vary very rapidly for wave vectors close to the Fermi surface (or multiples thereof) we require that the Fourier coefficients δn(q) of the density variation have their main contributions from wave vectors that satisfy |q| < k F where k F is the Fermi wave vector. Note that the procedure of interchanging integration and summation requires absolute convergence of the wave vector expansion. It is hard to prove this property for the general response functions that we consider here and we therefore have to assume its validity.
The next task is then to expand the functions N m in powers of wave vectors. This task is simplified when we exploit the symmetry of the these functions. As follows directly from Eqs. (6) and (7)  From these equations we therefore see that any power expansion of the functions N m must lead to an expansion in terms of the spatial vector y and the wave vectors q j which is invariant under rotations and inversions of these vectors. The mathematical question is then which polynomials have these properties. This was answered in the classic book by Weyl [29]; every such polynomial must be a function of the variables y · y, y · q i and q i · q j . We see that these inner products are indeed invariant under rotations and inversions. Since any polynomial in powers of y 2 = y · y can be re-summed to a function of y we find that the general expression for N m to second order in the wave vectors q j is given by where q 2 i = q i · q i and where we took into account that this expansion must be invariant under permutations of the wave vectors q j . This expansion is readily extended to higher order powers in the wave vectors. For a given choice of correlation function F the practical task will be to determine the explicit form of the coefficient functions N m j (y) as a function of y and the background density n 0 . If we insert Eq.(10) into Eq.(8) and Fourier transform back to real space we obtain the expansion Let us analyze the explicit form of the first six coefficients A m j for j = 0, .., 5 of this expansion, since these are exactly the terms that we need for a gradient expansion to second order in the density derivatives. The first term for j = 0 is simply given by The terms with j = 1, 2, 3 in Eq.(11) involve according to Eq.(10) the Fourier transforms of m symmetry equivalent terms and acquire the following forms in real space In the derivation of the last term we used that for an arbitrary function f where with r = (x 1 , x 2 , x 3 ) we used the notation ∂ i = ∂/∂x i as well as y i = x ′ i − x i . This expression is readily checked using Finally the terms with j = 4, 5 in Eq.(11) involve according to Eq.(10) the Fourier transforms of m(m − 1)/2 symmetry equivalent terms which yields We see that a general coefficient function A m j depends on the density variation and gradients thereof. The functions δn(r) that appear under the gradient operators can be replaced by the actual density profile n(r) = n 0 + δn(r). However, we see that we are still left with powers of δn(r) which are unprotected by derivative operators and which therefore can not be replaced by the full density profile n(r). It is precisely at this point that the consistency conditions (9) play a crucial role. If we use these conditions in Eq.(10) then we see that This equation relates certain coefficients coming from higher order response functions N m to those of lower order ones. In particular it tells us that by iteration If we insert these expressions together with the explicit form of the coefficient functions A m j back into Eq.(11) in Eq. (11) and shift some indices we obtain the expansion ∂n m 0 (y · ∇n(r)) 2 + . . .
In this expression we recognize several Taylor series of the coefficient functions N m j (n 0 + δn(r), y) in powers of δn(r). These Taylor series can be re-summed to finally give F (r ′ , r) = F 0 (n(r), y) + iN 1 1 (n(r), y) y · ∇n(r) − N 1 2 (n(r), y) ∇ 2 n(r) − N 1 3 (n(r), y) y 2 ( y y · ∇) 2 n(r) where the implicit dependence on n 0 of the coefficient functions is now replaced by a dependence on the full density profile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstated the explicit density dependence of the coefficients in Eq. (14) to indicate that they are local functions of the density n(r) and therefore have a nontrivial spatial dependence on both r and y. We stress again that the derivative condition (9) was essential in eliminating the dependence of the gradient coefficient on the reference density n 0 in favor of a dependence on the full density profile. Having obtained the general form of gradient expansion we can wonder about its convergence. In order for the gradient expansion to be useful it at least needs to be an asymptotic series. This is known to be the case for the gradient expansion of some commonly studied functionals. For instance, for small densities variations very accurate results for the exchange energy are obtained by summing all terms up to fourth order in the density derivatives [4,23].
Having obtained the general gradient expansion Eq. (14) it remains to find explicit expressions for the coefficient functions N m j . We will describe how to do this in detail for the coefficients needed for a gradient expansion to second order in the density derivatives. To calculate the coefficients N 1 1 , N 1 2 and N 1 3 we need to calculate the function N 1 (n 0 , y, q) for the homogeneous electron gas and expand it in powers of q: N 1 (n 0 , y, q) = N 1 0 (n 0 , y) + N 1 1 (n 0 , y)y · q + N 1 2 (n 0 , y)q 2 + N 1 3 (n 0 , y)(y · q) 2 + . . .
The determination of the coefficients N 2 4 and N 2 5 requires knowledge of the second order response function N 2 (n 0 , y, q 1 , q 2 ) = N 2 0 (n 0 , y) + N 2 1 (n 0 , y)y · (q 1 + q 2 ) + N 2 2 (n 0 , y)(q 2 1 + q 2 2 ) + N 2 3 (n 0 , y)((y · q 1 ) 2 + (y · q 2 ) 2 ) + N 2 4 (n 0 , y)(q 1 · q 2 ) + N 2 5 (n 0 , y)(y · q 1 )(y · q 2 ) + . . . (16) The expression (14) together with Eqs. (15) and (16) determine completely the gradient expansion of an arbitrary correlation function F to second order in the density gradients. These equations form the main result of the present work. If expansions to higher order gradients are required they can be derived without difficulty along the lines described above. To do this one needs to construct higher order symmetric polynomials in the wave vectors and carry out the required Fourier transforms. The consistency conditions Eq.(9) or equivalently Eq.(12) then allow for a complete re-summation and leads to gradient coefficients depending on n(r). What is needed to determine the explicit form of the gradient coefficients in practice is the determination of the functions N m . How to do this for N 1 and N 2 is described in the next section.

C. Determination of N 1 and N 2
A practical calculation of the coefficients of the gradient expansion of Eq. (14) requires the knowledge of the the functions N 1 and N 2 . According to Eqs.(4) and (6) these are defined as the first and second density derivatives of the correlation function F with respect to the density, i.e.
To derive useful expressions for these functions it is convenient to transform derivatives with respect to the den-sity to derivatives with respect to the external potential as such derivatives appear naturally in perturbation theory. We can then use the well-developed tools of perturbation theory to calculate these functions. Let us therefore define the new functions Since we evaluate these functions for the homogeneous electron gas they only depend on differences between the spatial coordinates. It will therefore be convenient to also define their Fourier transforms by as well as their partial Fourier transforms The chain rule of differentiation gives a connection between the derivatives with respect to the density and the derivatives with respect to the potential of the correlation function F . We have δ 2 F (r 1 , r 2 ) δv(r 3 )δv(r 4 ) = dr 5 δF (r 1 , r 2 ) δn(r 5 ) δ 2 n(r 5 ) δv(r 3 )δv(r 4 ) + dr 5 dr 6 δ 2 F (r 1 , r 2 ) δn(r 5 )δn(r 6 ) We see that in these expressions the first and second order derivative of the density with respect to the potential appear. We therefore define the linear density response function by χ as well as the second order density response function χ 2 by δ 2 n(r 1 ) δv(r 2 )δv(r 3 ) | n0 = χ 2 (r 2 − r 1 , r 3 − r 1 ) = dpdq (2π) 6 χ 2 (p, q)e ip·(r2−r1)+iq·(r3−r1) .
Using these definitions we find by Fourier transforming Eqs. (19) and (20) that These equations give the desired relation between the density derivatives N 1 and N 2 and the potential derivatives F 1 and F 2 of the correlation function F . Relations between the higher order derivatives N m and F m can be derived in a completely analogous way. From Eqs. (21) and (22) we see that to calculate the functions N 1 and N 2 , and hence the gradient expansion coefficients of F to second order in the gradients, we need to calculate the density response and the change in F to second order in a perturbing potential δv. The problem is thereby reduced to doing perturbation theory on the electron gas. This is a well-developed field of research. One option is to use diagrammatic perturbation theory [30]. In the following section we will use standard perturbation theory to obtain these response function for the case that F is the one-particle density matrix and use that to calculate the gradient expansion of the exchange hole.

III. GRADIENT EXPANSION OF THE ONE-PARTICLE DENSITY MATRIX AND THE EXCHANGE HOLE
A. The one-particle density matrix As an application of our formalism we carry out the gradient expansion of the one-particle density matrix for a noninteracting system with density n(r). This problem has received large interest in the past since both the exchange energy E x [n] as well the Kohn-Sham kinetic energy T s [n] are explicit functionals of such a noninteracting density matrix [2,3]. Therefore a gradient expansion of this density matrix directly leads to a gradient expansion of these two functionals. Such a gradient expansion was first carried out by Gross and Dreizler [31] using the Kirzhnits expansion [2,32]. This expansion is adjusted to the specific form of the noninteracting oneparticle density matrix and can not be generalized to arbitrary correlation functions. The method presented in this work can, on the other hand, be applied to arbitrary correlation functions, but to demonstrate the soundness of our formalism we will use it to provide an alternative derivation of the result obtained earlier by Gross and Dreizler using the Kirzhnits expansion. An advantage of our derivation based on nonlinear response theory is that it avoids the turning point problem encountered in the Kirzhnits approach [2].
Let us start by defining the one-particle density matrix in second quantization as where σ is a spin coordinate. We will consider the case of spin-compensated systems. Then, in a noninteracting system with 2N electrons, the density matrix is given in terms of one-particle wave functions ϕ j (r) by where the pre-factor 2 is a spin factor. The one-particle states are eigenstates of a one-particle Schrödinger equation where v(r) is the external potential. To determine the gradient expansion coefficients of γ we need to determine how the density matrix changes if we make small changes v(r) → v(r) + δv(r) in the external potential. More precisely, we need to calculate the functional derivatives which are the equivalents of the functions F 1 and F 2 of Eqs. (17) and (18) for the specific choice of correlation function F = γ (note that they become functions of two and three independent vectors respectively when evaluated for a homogeneous system). To do this it is sufficient to know the functional derivatives of the orbitals ϕ j and eigenvalues ǫ j with respect to the potential v. These are simply given by doing first order perturbation theory on the one-particle Schrödinger Eq. (24). We find With these relations we can differentiate Eq.(23) twice with respect to the potential. Afterwards we then need to insert the plane wave states ϕ k (r) and their eigenvalues ǫ k appropriate for the homogeneous electron gas into the final expressions. In order not to interrupt the presentation these calculations are given in Appendix A and we simply present the results of these calculations here.
If we define the Fermi factors by , where θ is the Heaviside function and where k F = (3π 2 n 0 ) 1/3 is the Fermi wave vector, then the resulting expressions are given by where In this expression ǫ k = |k| 2 /2 are the one-particle energies. It remains to calculate the first and second order density response functions χ(q) and χ 2 (p, q) to evaluate the functions N 1 and N 2 of Eqs. (21) and (22). Fortunately, for the specific choice F = γ that we made these density response functions are simply given by χ 2 (p, q) = γ 1 (y = 0, p, q).
They are therefore obtained directly from Eqs. (27) and (28) so that we do not need to do any extra work to calculate them.
To calculate the gradient coefficient to second order in the gradients we now need to expand γ 1 and γ 2 to second order in the wave vectors p and q. Since these calculations are somewhat lengthy we do not present them here and refer to Appendix B for a description of the main steps. The resulting expansions are given by where we defined z = k F y. The expansions for the first and second order response functions χ and χ 2 of Eqs. (30) and (31) are obtained by taking y = 0 in these expressions. Together which Eqs. (32) and (33) these expression then immediately determine the functions N 1 and N 2 from Eqs. (21) and (22). The final result of this calculation to second order in the wave vectors is given by By comparison of these two expression to the expansions (15) and (16) we can directly read off the gradient coefficient functions of Eq. (14) to be where now we replaced n 0 → n(r) and hence the Fermi wave vector k F = (3π 2 n(r)) 1/3 that appears in these expressions is a spatially varying function. We further introduced the spherical Bessel functions j 0 (z) = sin z z j 1 (z) = sin z − z cos z z 2 .
By inserting the gradient coefficient functions into Eq.(14) we find the following explicit gradient expansion for the one-particle density matrix This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent form as This expression becomes identical in form to that derived by Gross and Dreizler [31] using the Kirzhnits method after eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(19) of reference [31] ). We close this section with a final remark on our result. We note that the gradient expansion to finite order breaks the symmetry γ(r, r ′ ) = γ(r ′ , r) * . However, the full gradient expansion as described by Eq.(11) has this symmetry, which means that the symmetry is restored by taking all higher order gradients into account whenever the series converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such as the exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literature we will in the next section adopt the definition of the exchange hole that was used by Perdew [9].

B. The exchange hole and energy
We finally stress a few points on the properties of the exchange hole as derived from the gradient expansion and present an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler and Gross which has the advantage of avoiding the regularization of divergent integrals [31]. In doing so we confirm a result obtained by Langreth and Perdew [19]. The exchange hole can be calculated directly from the one-particle density matrix as ρ x (r + y|r) = − |γ(r + y, r)| 2 2n(r) .
Inserting the expansion (34) of for the one-particle density matrix and retaining terms to second order in the gradients yields the expression As was pointed out by Perdew [9] the second order gradient expansion of the exchange hole does not satisfy the negativity and sum rule constraints The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediately obvious from Eq. (36). The sum rule can, however, be conveniently analyzed in momentum space. To do this we define the Fourier transformed exchange hole to be This function has the form ρ x (k|r) =ρ 0 x (k|r) + α 1 (k, n)k · ∇n(r) + α 2 (k, n) ∇ 2 n(r) + α 3 (k, n) (k · ∇n(r)) 2 + α 4 (k, n)(∇n(r)) 2 + α 5 (k, n) (k · ∇) 2 n(r) + . . .
where we defined the unit vectork = k/k and k = |k|. The explicit form of the functions ρ 0 x and α j follow directly by Fourier transforming Eq. (36) and are given in Appendix C. The coefficient functions α j are tempered distributions [33] which have mathematically well-defined Fourier transforms. As follows directly from Eq.(37) the sum rule condition in momentum space translates to the requirement that ρ x (k = 0|r) = −1. Since ρ 0 x (k = 0|r) = −1 (see Eq.(C1) ) this implies that the sum rule would be fullfilled whenever α j (k = 0, n) = 0. However, as is clear from Eqs.(C2)-(C6) in the Appendix this is not satisfied for α 1 , whereas α 3 and α 5 even diverge for k → 0 when interpreted locally as functions rather than distributions. The gradient expansion of the exchange hole does therefore indeed not satisfy the sum rule. These divergencies, however, cancel when we calculate the system averaged exchange hole in momentum space which is given by the expression N ρ x (k) = dr n(r)ρ x (k|r) = N ρ 0 x (k) +k · dr n(r)α 1 (k, n(r))∇n(r) dr n(r)β ij (k, n(r))∂ i n(r)∂ j n(r).
where N is the number of electrons. In this expression ρ 0 x (k) is the system average of ρ 0 x (k|r) and in the last term under the integral sign we defined the tensor This tensor has a longitudinal part with coefficient α L and a transverse part with coefficient α T which describe the contributions to the system-averaged hole of density gradients ∇n parallel and perpendicular to the momentum vector k. These coefficients are calculated from the functions α j as as a short calculation on the basis of Eq.(38) will show. From these equations and the explicit expressions for α j given in Appendix C we can readily calculate the explicit expressions for α L,T . If we define the dimensionless variablek = k/(2k F ) then we have where the functions Z L,T have the explicit form These expressions are in accordance with the results of Langreth and Perdew (see Eq.(3.55) of reference [19]) and this independent result therefore confirms the correctness of our alternative derivation. The interesting observation is now that lim k→0 Z L,T (k) = 0 and hence β ij (k, n) → 0 when k → 0. This implies that the last term in Eq. (39) does not contribute to the sum rule of the system averaged hole. The same conclusion can be derived for the second term in Eq.(39). Since α 1 (k = 0, n) = i/(4nk F ) (see Eq.(C2)) is a local function of the density the integrand of the first integral in Eq.(39) is a total derivative and vanishes (assuming either vanishing densities at infinity or periodic boundary conditions). We therefore find that This implies that system averaged exchange hole obtained from the second order gradient expansion does fulfill the sum rule, although the exchange hole itself does not. We note, however, that this is only achieved by integrating over both positive and negative contributions. When the positive contributions to the exchange hole are removed the sum rule is only recovered for a finite hole radius [9,10]. We finally note that with expression Eq. (39) it is now a simple matter to calculate the exchange energy from We insert Eq. (39) and do the angular integrations first. It is therefore convenient to define the spherical and system averaged hole in momentum space by such that From Eq.(3) we then find that where we defined The exchange energy is obtained by inserting the expression for the averaged hole of Eq. where and where to calculate the first term in Eq.(43) we used Eq.(C1). The coefficient B x is the same gradient coefficient as obtained by Gross and Dreizler [31] and earlier by Sham [34]. However, the correct analytic exchange gradient coefficient is known [35][36][37][38][39] to be a factor 10/7 larger. The reason for this discrepancy is clearly described by Svendsen and von Barth [38] who showed that the Coulomb interaction is too singular to allow for the interchange of the operations of integration and the expansion in wave vectors. The problem does, for instance, not appear when Yukawa-screened Coulomb interactions are used [38]. The conclusion is that there is nothing wrong with the Kirzhnits method, or the nonlinear response theory derivation of the gradient expansion of the exchange hole presented here, but one has be aware that for Coulomb interactions the procedures of expanding correlation functions in terms of density gradients followed by integrations involving Coulomb potentials may not yield the same result as directly expanding the integrated quantity in terms of density gradients. For this reason the original GGA of Perdew and Wang [10] based on the gradient expansion of the exchange hole was later reparametrized [16,40] to accomodate the correct gradient coefficient for the exchange energy.

IV. CONCLUSIONS AND OUTLOOK
We derived a general scheme based on nonlinear response theory to calculate the density gradient expansion of general correlation functions and showed that in order to express the gradient coefficients in terms of the full density profile summations to infinite order must be carried out over response functions of arbitrarily large order. A consistency condition was derived to do this. The formalism was used to derive the gradient expansion of the one-particle density matrix and the exchange hole to second order in the gradients. We confirm the derivation of Dreizler and Gross that used the Kirzhnits expansion method. We further analyzed the exchange hole in momentum space to derive that the system averaged hole to second order in the gradients satisfies the sum rule and to derive the gradient expansion of the exchange energy without the need to regularize divergent integrals.
The scheme that we presented is very general and can be applied to more general correlation functions. A natural application of the scheme would be to calculate the gradient expansion of the correlation hole ρ c . Regarding the gradient terms of the correlation hole little is known. We essentially only have some detailed information on the long-range properties of the spherically and system averaged hole. This information comes from the work of Langreth and Perdew who calculated ρ c (k) within the RPA. In our notation their results (see also Appendix C of [19]) can be summarized as N ρ c (k) = N ρ 0 c (k) + dr π 2 16k 4 F z c (n, k)(∇n) 2 where the function z c has been parametrized by Langreth and Mehl [20,21] as where k s = (4k F /π) 1/2 is the Thomas-Fermi or screening wave vector. We can transform to real space to obtain the following expression for the spherically and system averaged correlation hole N ρ c (y) = N ρ 0 c (y) + dr 1 96k 4 F k 2 s (1 + (k s y) 2 /12) 2 (∇n) 2 (45) We see from Eq.(44) that z c (n, k = 0) = 0 and as a consequence the sum rule ρ c (0) = 0 for the correlation hole is not satisfied. We know, however, that the RPA becomes accurate in the high density limit and since from Eq.(45) the gradient coefficient of the correlation hole depends on k s y we see that high densities are equivalent to large separations y. Similarly, low density corresponds to short-distance behavior. However, RPA can not be trusted in this region and we have no precise knowledge on the small distance behavior of the gradient coefficient of ρ c (y) . A model for the general short-range behavior was proposed [16][17][18] on physical arguments (it should not affect the Coulomb cusp of and the on-top value of the LDA hole) after which the real-space cutoff procedure was applied (see e.g. Fig. 4 in [17] and Fig. 5 in [18]) to obtain a GGA for correlation. It would, however, be desirable to have a first principles approach to the calculation of the short range properties of the gradient coefficient of the correlation hole. As follows from our derivation this requires the knowledge of the functions (17) and (18), or at least their expansion to second order in wave vectors, for the case that F represents the pair correlation function or the correlation hole of the electron gas. It may therefore be worthwhile to use our current scheme to explore these response functions beyond the RPA. For short range correlation an approach based on diagrammatic summation of ladder diagrams suggests itself. Of course, for the development of density functionals for general systems beyond the weakly inhomogeneous regime it is not sufficient to use the straightforward gradient expansion [23]. However, general shortrange features, such as the way the exchange-correlation hole deforms close to the reference electron in the presence of a density gradient can be transferred to more general systems than the weakly inhomogeneous ones. Perhaps in combination with truly nonlocal functionals for the exchange-correlation hole [41] this can lead to the development of more accurate density functionals. This will be a topic of our future research.
Finally, we would like to make some general remarks on the extension of our derivations to temperature-or time-dependent systems. In the case of temperaturedependent systems the expectation values of observables are traces over a grand canonical ensemble. By Mermin's theorem [42] such expectation values are still functionals of the density and therefore the derivations in this work can be carried out unchanged. In fact, the dependence on temperature is probably going to improve the convergence properties of the gradient series by smoothening of the sharp Fermi functions which, for instance, were the cause of the oscillatory nature of the gradient coefficients of the one-particle density matrix. The situation of time-dependent systems is more delicate. Although the existence of a density-potential mapping has been well-established [43][44][45][46] and hence a density functional theory can be formulated, the time-variable induces severe complications. As was shown by Vignale and Kohn [47,48] the temporal and spatial non-locality of time-dependent density functionals are intimately correlated. Any temporal non-locality (or frequency dependence in the language of equilibrium response functions) induces long-range spatial properties that prevents the gradient expansion from existing. A way out of this problem is to use new variables such as the current-density [47,48] or the local density deformation density in a Lagrangian frame [49][50][51] for which a local density approximation and a corresponding gradient expansion does exist. This has been exploited in Ref. [52] to construct a GGA functional within time-dependent current density functional theory. The study of correlation functions in the same fashion remains a challenge for the future. which seems to be different from expression (A4). However, this is just appearance. The reason is that the occupations can only attain the values zero and one. For the case f i = f j = 1 or f i = f j = 0 we see directly that Eq.(A4) and (A7) attain the same value. A little inspection shows that this is also true for cases f i = 0, f j = 1 and f i = 1, f j = 0. We therefore find for k = i, j that Φ(ijk) = Φ(jik). Therefore Eq.(A6) can be simplified to γ 2 (23, 14) = 2 ∞ ijk,k =(i,j) Φ(ijk) ϕ k (1)ϕ * i (1)ϕ * k (4)ϕ j (4) + ϕ j (1)ϕ * k (1)ϕ * i (4)ϕ k (4) + Y (23, 14) + Z (23,14) (A8) This expression is now explicitly symmetric in the indices 1 and 4. Let us now evaluate this expression for the homogeneous electron gas. In the electron gas the oneparticle states are plane waves where V is the volume of the system. Since |ϕ k | 2 = 1/V it follows from Eqs.(A3) and (A5) that the terms Y and Z are identically zero and the function γ(23, 14) of Eq.(A8) therefore attains the form γ 2 (23, 14) = 2 V 3 k,p,q,q =(k,p) Φ(k, p, q)e i(k·r2−p·r3) e i(q−k)·r1+i(p−q)·r4 + e i(p−q)·r1+i(q−k)·r4 where we defined Φ(k, p, q) = 1 ǫ q − ǫ p n q − n k ǫ q − ǫ k − n p − n k ǫ p − ǫ k .
which is simply the Fourier transform of the LDA exchange hole. We see that ρ x (0|r) = −1 and therefore the LDA hole satisfies the sum rule. The functions α j are given by