The density gradient expansion of correlation functions
aa r X i v : . [ c ond - m a t . m t r l - s c i ] A p r The density gradient expansion of correlation functions
Robert van Leeuwen
1, 2 Department of Physics, Nanoscience Center, University of Jyv¨askyl¨a, 40014 Jyv¨askyl¨a, Finland European Theoretical Spectroscopy Facility (ETSF) (Dated: October 30, 2018)We present a general scheme based on nonlinear response theory to calculate the expansion ofcorrelation functions such as the pair-correlation function or the exchange-correlation hole of aninhomogeneous many-particle system in terms of density derivatives of arbitrary order. We furtherderive a consistency condition that is necessary for the existence of the gradient expansion. Thiscondition is used to carry out an infinite summation of terms involving response functions up toinfinite order from which it follows that the coefficient functions of the gradient expansion can beexpressed in terms the local density profile rather than the background density around which theexpansion is carried out. We apply the method to the calculation of the gradient expansion of theone-particle density matrix to second order in the density gradients and recover in an alternativemanner the result of Gross and Dreizler (Z. Phys. A , 103 (1981)) which was derived using theKirzhnits method. The nonlinear response method is more general and avoids the turning pointproblem of the Kirzhnits expansion. We further give a description of the exchange hole in momentumspace and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B , 5469 (1980))for this case. This is used derive that the second order gradient expansion of the system averagedexchange hole satisfies the hole sum rule and to calculate the gradient coefficient of the exchangeenergy without the need to regularize divergent integrals. PACS numbers: 31.15.ee, 71.15.Mb, 31.10.+z
I. INTRODUCTION
Since the groundbreaking work of Hohenberg andKohn [1] we know that the external potential of any inho-mogeneous quantum many-body system is a functional ofits ground state density n . This implies that the many-body ground state | Ψ[ n ] i and hence any ground stateexpectation value A [ n ] = h Ψ[ n ] | ˆ A | Ψ[ n ] i (1)for any operator ˆ A is a functional of the density. Thisidea has inspired an enormous amount of work in a re-search field that is now known as density-functional the-ory. Density-functional theory [2–4] is currently one ofthe main approaches used in electronic structure the-ory. Over the years many extensions beyond its origi-nal formulation have been developed and currently it iswidely applied in solid state physics and quantum chem-istry, both for ground state and time-dependent systems[5]. Especially a large activity has gone into finding ex-plicit expressions for the ground state energy E [ n ] as afunctional of the density for many-electron systems. Theground state energy is obtained by minimizing this func-tional on an appropriate set of electronic densities, whichis of great importance in determining structures and ge-ometries of molecules and solids. By use of the Kohn-Sham method [6] this minimization problem is reducedto solving effective one-particle equations. The construc-tion of the effective potential in these equations requiresthe knowledge of the so-called exchange-correlation en-ergy functional E xc [ n ]. This functional can be expressed as [7, 8] E xc [ n ] = 12 Z d r d r ′ n ( r ) ρ xc ( r ′ | r ) | r − r ′ | (2)where ρ xc ( r ′ | r ) is the coupling constant integratedexchange-correlation hole. The exchange-correlation holehas a physical interpretation as the difference betweenthe conditional and the unconditional probability (whichis simply the density n ( r ′ )) to find a electron at r ′ , giventhat we know that there is an electron at r . Equation (2)is an important relation since it expresses the exchange-correlation energy in terms of a quantity that has a lo-cal physical interpretation and can be studied by accu-rate wave function and many-body methods or modeledbased on physical intuition. It has therefore played animportant role in the development of approximate den-sity functionals [8–14]. One of the simplest and alreadyvery successful approximations has been the local densityapproximation (LDA) in which the exchange-correlationhole is taken from the homogeneous electron gas and ap-plied locally [8] to systems with arbitrary density profiles.The accuracy of the LDA has been considerably improvedby means of the so-called generalized gradient approxi-mations (GGAs) [15]. A large class of commonly usedGGAs is based on the so-called real-space cutoff of thestraightforward gradient expansion for the exchange andcorrelation hole. The first such functional for exchangewas proposed by Perdew [9]. He noted that the gradientexpansion of the exchange hole to second order in thedensity gradients violates the negativity condition andthe hole sum rule by which the exchange hole integratesto a deficit of one electron. By enforcing these physicalconstraints a density functional was obtained that greatlyimproved on the LDA for binding energies of molecules.This procedure was later simplified [10] using the factthat the exchange correlation energy (2) can be writtenas E xc [ n ] = N Z d y h ρ xc ( y ) i| y | where N is the number of electrons in the system and h ρ xc ( y ) i = 1 N Z d r n ( r ) ρ xc ( r + y | r ) , (3)is the so-called system-averaged exchange-correlationhole. This averaged hole still obeys the negativity con-dition as well as the sum rule and can therefore be sub-jected to the real-space cut-off procedure. One advan-tage of using the system averaged holes was to reducethe order of the derivatives in the gradient expansion.Furthermore it also allowed for the real-space cut-off pro-cedure to be applied to the correlation hole since there isa known gradient expansion for the system averaged cor-relation 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–18]. These GGAs relied heavily on the pioneeringwork by Langreth, Perdew and Mehl [19–21] who per-formed a wave-vector decomposition of the system aver-aged hole of Eq. (3) and calculated the Fourier transform h ρ xc ( k ) i = Z d y h ρ xc ( y ) i e − i k · y , from which the exchange correlation energy can be cal-culated as E xc [ n ] = N Z d k (2 π ) h ρ xc ( k ) i π | k | . We are not aware of any work that goes beyond refer-ences [19–21] and improves on the straightforward gra-dient expansion of the system averaged correlation hole.The knowledge on the gradient expansion of the correla-tion hole is very limited indeed. As mentioned before, incontrast to the work on the exchange hole [2, 9], thereare no known expressions for the gradient expansion ofthe non-system averaged correlation hole.Part of the problem is that there has been no clear deriva-tion on how to expand the expectation value of a gen-eral operator as in Eq.(1) in terms of density gradients.A well-known expansion method for a two-point func-tion is the Kirzhnits expansion [2] which is, however,specifically adapted to expanding the noninteracting one-particle density matrix in powers of density gradients andcan not be used for more general correlation functions.The first goal of this paper is, therefore, to present ascheme based on nonlinear response theory that can beused to expand general correlation functions (such as thecorrelation 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 stan-dard gradient expansion for a global quantity, namelythe exchange-correlation energy E xc . The starting pointof any gradient expansion is the consideration of a den-sity profile n ( r ) = n + δn ( r ) which consists of a smalldensity variation δn ( r ) around a homogeneous density n . By considering the limit of slow density variationsone then finds that the lowest gradient correction to theexchange-correlation energy is an integral over a func-tion of the form B xc ( n )( ∇ δn ( r )) where the coefficient B xc ( n ) is calculated from the static exchange-correlationkernel of the homogeneous electron gas of density n [4].Since ∇ δn ( r ) = ∇ n ( r ) we can replace the density vari-ation under the derivative operator by the full densityprofile n ( r ). However, the coefficient B xc ( n ) still de-pends on the background density n . This is a problemin application of the formula to general inhomogeneoussystems such as molecules or surfaces in which a back-ground density cannot be unambiguously defined (assum-ing that a low order gradient expansion makes sense insuch systems [22, 23]). The usual argument to get ridof the background dependence, is that the replacement B xc ( n ) → B xc ( n + δn ( r )) can be made since the differ-ence between these two quantities is or order δn ( r ). How-ever, it is not clear that this is consistent with gradientcoefficients that arise from terms of order ( δn ( r )) . Thepoint was discussed clearly by Svendsen and von Barth[22] who checked that the replacement of n by the fulldensity profile was consistent to third order in the densityvariation. This derivation was based on a specific rela-tion between response functions that describe the changein the exchange energy to second and third order in thedensity variations. In reference [24] this derivation wasextended by showing that the replacement of n by thefull density profile is consistent to all orders in the den-sity variation. In this paper we will show that this isthe case as well for the gradient expansion of correlationfunctions rather than scalar functions. This relies on cer-tain relations between the response functions that mustbe satisfied for the gradient expansion to exist.The paper is divided as follows. In section II we de-rive the basic equations and consistency conditions forthe density gradient expansion of correlation functions.In section III we derive the gradient expansion of theone-particle density matrix for a noninteracting systemwith density n ( r ) (which is therefore equal to the den-sity matrix of the Kohn-Sham system) and discuss itssymmetry properties. We further calculate the gradientexpansion of the exchange-hole to second order in thedensity gradients, both in real and in momentum space.The momentum space description is used to demonstratethat the system averaged exchange hole satisfies the sumrule (but not the negativity condition), and to calculatethe gradient expansion of the exchange energy withoutthe need to regularize divergent integrals. Finally in sec-tion IV we present our conclusions and outlook. II. GRADIENT EXPANSION OFCORRELATION FUNCTIONSA. Expansion in density variations
Let us for a system with ground state density n ( r ) con-sider an arbitrary correlation function F [ n ]( r ′ , r ). Sincethe correlation function can be calculated from the many-body ground state, by virtue of the Hohenberg-Kohn the-orem this function is a functional of the density. Atthis point it will not be important what the specificform of the correlation function is, as the structure oftheir gradient expansions will be independent of its spe-cific 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 , we are describing the homogeneous electrongas and due to translational and rotational invariancewe have that F [ n ]( r ′ , r ) = F ( n , | r − r ′ | ) where F is a function of the homogeneous density and the dis-tance between its spatial arguments. When the density n ( r ) = n + δn ( r ) deviates from the constant density n by a small amount δn ( r ) we can expand F in powers ofthis density variation. To derive compact expressions wefirst introduce the notation r m = ( r , . . . , r m ) d r m = d r . . . d r m for the collection of m position vectors and the corre-sponding integration volume elements. We further define δn ( r m ) = m Y j =1 δn ( r j ) . With these conventions the expansion of F in densityvariations is given by F [ n ]( r ′ , r ) = F ( n , | r ′ − r | )+ ∞ X m =1 m ! Z d r m F m ( n ; r ′ , r , r m ) δn ( r m )where we defined F m ( n ; r ′ , r , r m ) = δ m F ( r ′ , r ) δn ( r ) . . . δn ( r m ) | n . (4)The n dependence of the functions F m will be impor-tant for the gradient expansion, but to shorten notationwe suppress this dependence in the notation for the timebeing and will reintroduce it later. Here we assume thatthe derivatives F m exist, which means that we assumethat 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 ≪
1. The density vari-ations are therefore assumed to be small but we do notneed to put any constraint on how rapid they can vary inspace and therefore their spatial derivatives can be large.Let us now look at the symmetry properties of the func-tions F m . As follows directly from the definition (4) ofthe functions F m and the assumption of differentiability,we have the permutational symmetry F m ( r ′ , r , r . . . r m ) = F m ( r ′ , r , r π (1) . . . r π ( m ) )for all permutations π of the integers 1 . . . m . Since thefunctions F m are evaluated at the homogeneous density n they also have the spatial symmetry properties of thehomogeneous electron gas. These symmetries are thetranslational symmetry over an arbitrary vector a , rota-tional symmetry under arbitrary rotations by a rotationmatrix R , and inversion symmetry. If we denote by a = ( a , . . . , a ) R r m = ( R r , . . . , R r m )the m -tuples of translation vectors and rotated positionvectors we can write F m ( r ′ , r , r m ) = F m ( r ′ + a , r + a , r m + a ) (5) F m ( r ′ , r , r m ) = F m ( R r ′ , R r , R r m ) F m ( r ′ , r , r m ) = F m ( − r ′ , − r , − r m ) . Since in Eq.(5) the vector a is arbitrary we can in par-ticular choose a = − r and define a new function N m depending on m + 1 independent vectors by the relation F m ( r ′ , r , r . . . r m )= F m ( r ′ − r , , r − r , . . . , r m − r )= N m ( r ′ − r , r − r , . . . , r m − r ) . (6)The difference vector r ′ − r will appear several times inour equations and it will therefore be convenient to definethe short notation y = r ′ − r . We will further use y = | y | to denote the length of this vector. Our interest willbe in the functions N m and in particular their Fouriertransforms N m ( y , q m ) = Z d r m N m ( y , r m ) e − i q · r − ... − i q m · r m (7)with Fourier inverse N m ( n ; y , r m ) = Z d q m (2 π ) m N m ( y , q m ) e i q · r + ... + i q m · r m . With these definitions we obtain the following expansionfor F in powers of the density variations F [ n ]( r ′ , r ) = F ( n , y ) + ∞ X m =1 m ! Z d r m Z d q m (2 π ) m N m ( y , q m ) e i q · ( r − r )+ ... + i q m · ( r m − r ) δn ( r m )= F ( n , y ) + ∞ X m =1 m ! Z d q m (2 π ) m N m ( y , q m ) e − i ( q + ... + q m ) · r δn ( − q m ) (8)This equation forms the basis of the gradient expansion.To proceed we further derive a consistency condition thatis 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 ) = Z d r ′′ F m +1 ( r ′ , r , r ′′ , r m ) δn ( r ′′ ) . Taking δn ( r ′′ ) = δn then to be a uniform shift in thedensity (which means that we look at a compression ordecompression of the electron gas) yields ∂F m ∂n ( r ′ , r , r m ) = Z d r ′′ 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 Fouriertransform. This gives the condition ∂N m ∂n ( y , q . . . q m ) = N m +1 ( y , , q . . . q m ) . (9)This relation is of key importance for the existence of thegradient expansion. Without it we would not be able toeliminate the dependence of the gradient coefficients onthe homogeneous density n in favor of a dependence onthe actual density profile. As we will see, this requiresa summation over response functions F m to infinite or-der. There is another important advantage of this re-summation. It will allow us to relax the constraint thatthe density variations be small (but varying arbitrarilyfast) and replace it with the constraint that the densityvariations 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 toassume that the density gradients are small. This con-straint is most easily phrased in momentum space by re-quiring that the Fourier coefficients δn ( q ) of the density variation have their main contribution from small wavevectors q . If this is the case then the main contribu-tion to the integrals in Eq.(8) comes from this region ofsmall vectors and we can expand the functions N m inpowers of the wave vectors. Subsequently carrying outthe integrals in Eq.(8) then leads to an expansion in den-sity gradients, since powers in wave vectors correspondto orders of derivatives in real space. Since the responsefunctions N m typically vary very rapidly for wave vec-tors close to the Fermi surface (or multiples thereof) werequire that the Fourier coefficients δn ( q ) of the densityvariation have their main contributions from wave vec-tors that satisfy | q | < k F where k F is the Fermi wavevector.The next task is then to expand the functions N m inpowers of wave vectors. This task is simplified when weexploit the symmetry of the these functions. As followsdirectly from Eqs.(6) and (7) the functions N m in Fourierspace inherit the symmetry properties of F m . They aresymmetric functions of their wave vectors, i.e. N m ( y , q . . . q m ) = N m ( y , q π (1) . . . q π ( m ) )for any permutation π of the integer labels, and are in-variant under rotations and inversion N m ( y , q . . . q m ) = N m ( R y , R q . . . R q m ) N m ( y , q . . . q m ) = N m ( − y , − q . . . − q m ) . From these equations we therefore see that any powerexpansion of the functions N m must lead to an expansionin terms of the spatial vector y and the wave vectors q j which is invariant under rotations and inversions ofthese vectors. The mathematical question is then whichpolynomials have these properties. This was answeredin the classic book by Weyl [25]; every such polynomialmust be a function of the variables y · y , y · q i and q i · q j .We see that these inner products are indeed invariantunder rotations and inversions. Since any polynomial inpowers of y = y · y can be re-summed to a function of y we find that the general expression for N m to secondorder in the wave vectors q j is given by N m ( y , q . . . q m ) = N m ( y ) + N m ( y ) y · m X i =1 q i + N m ( y ) m X i =1 q i + N m ( y ) m X i =1 ( y · q i ) + N m ( y ) m X i>j ( q i · q j ) + N m ( y ) m X i>j ( y · q i )( y · q j ) + . . . (10)where q i = q i · q i and where we took into account thatthis expansion must be invariant under permutations ofthe wave vectors q j . This expansion is readily extendedto higher order powers in the wave vectors. For a givenchoice of correlation function F the practical task will beto determine the explicit form of the coefficient functions N mj ( y ) as a function of y and the background density n .If we insert Eq.(10) into Eq.(8) and Fourier transformback to real space we obtain the expansion F ( r ′ , r ) = F ( n , y ) + ∞ X m =1 m ! ∞ X j =0 N mj ( y ) A mj ( y , r ) . (11)Let us analyze the explicit form of the first five coeffi-cients A mj for j = 0 , .., j = 0 is simply given by A m = δn ( r ) m . The terms with j = 1 , , m symmetry equivalentterms and acquire the following forms in real space A m = i m ( δn ( r )) m − y · ∇ δn ( r ) A m = − m ( δn ( r )) m − ∇ δn ( r ) A m = − m ( δn ( r )) m − y ( y y · ∇ ) δn ( r )In the derivation of the last term we used that for anarbitrary function f ( y y · ∇ ) f = X ij y i y ∂ i ( y j y ∂ j ) f ( r ) = X ij y i y j y ∂ i ∂ j f where with r = ( x , x , x ) we used the notation ∂ i = ∂/∂x i as well as y i = x ′ i − x i . This expression is readilychecked using ∂ i ( y j y ) = − δ ij y + y i y j y . Finally the terms with j = 4 , m ( m − / A m = − m ( m − δn ( r )) m − ( ∇ δn ( r )) A m = − m ( m − δn ( r )) m − ( y · ∇ δn ( r )) . We see that a general coefficient function A mj dependson the density variation and gradients thereof. Thefunctions δn ( r ) that appear under the gradient opera-tors can be replaced by the actual density profile n ( r ) = n + δn ( r ). However, we see that we are still left withpowers of δn ( r ) which are unprotected by derivative op-erators and which therefore can not be replaced by thefull density profile n ( r ). It is precisely at this point thatthe consistency conditions (9) play a crucial role. If weuse these conditions in Eq.(10) then we see that N m +1 j ( y ) = ∂N mj ∂n ( y ) . (12)This equation relates certain coefficients coming fromhigher order response functions N m to those of lowerorder ones. In particular it tells us that by iteration N m ( y ) = ∂ m F ∂n m ( y ) N mj ( y ) = ∂ m − N j ∂n m − ( y ) ( j = 1 , , N mj ( y ) = ∂ m − N j ∂n m − ( y ) ( j = 4 , . If we insert these expressions together with the explicitform of the coefficient functions A mj back into Eq.(11) inEq.(11) and shift some indices we obtain the expansion F ( r ′ , r ) = F ( y ) + ∞ X m =1 m ! ∂ m F ∂n m ( y ) δn ( r ) m + ∞ X m =0 m ! δn ( r ) m (cid:2) i ∂ m N ∂n m ( y · ∇ n ( r )) − ∂ m N ∂n m ∇ n ( r ) − ∂ m N ∂n m y ( y y · ∇ ) n ( r ) (cid:3) + − ∞ X m =0 m ! δn ( r ) m (cid:2) ∂ m N ∂n m ( ∇ n ( r )) + ∂ m N ∂n m ( y · ∇ n ( r )) (cid:3) + . . . (13)In this expression we recognize several Taylor series of the coefficient functions N mj ( n + δn ( r ) , y ) in powers of δn ( r ).These Taylor series can be re-summed to finally give F ( r ′ , r ) = F ( n ( r ) , y ) + iN ( n ( r ) , y ) y · ∇ n ( r ) − N ( n ( r ) , y ) ∇ n ( r ) − N ( n ( r ) , y ) y ( y y · ∇ ) n ( r ) − N ( n ( r ) , y ) ( ∇ n ( r )) − N ( n ( r ) , y ) ( y · ∇ n ( r )) + . . . , (14)where the implicit dependence on n of the coefficient functions is now replaced by a dependence on the full densityprofile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstatedthe 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 in favor of adependence on the full density profile.Having obtained the general gradient expansion Eq.(14) it remains to find explicit expressions for the coefficientfunctions N mj . We will describe how to do this in detail for the coefficients needed for a gradient expansion tosecond order in the density derivatives. To calculate the coefficients N , N and N we need to calculate the function N ( n , y , q ) for the homogeneous electron gas and expand it in powers of q : N ( n , y , q ) = N ( n , y ) + N ( n , y ) y · q + N ( n , y ) q + N ( n , y )( y · q ) + . . . (15)The determination of the coefficients N and N requires knowledge of the second order response function N ( n , y , q , q ) = N ( n , y ) + N ( n , y ) y · ( q + q ) + N ( n , y )( q + q ) + N ( n , y )(( y · q ) + ( y · q ) ) + N ( n , y )( q · q ) + N ( n , y )( y · q )( y · q ) + . . . (16)The expression (14) together with Eqs.(15) and (16) de-termine completely the gradient expansion of an arbi-trary correlation function F to second order in the den-sity gradients. These equations form the main result ofthe present work. If expansions to higher order gradientsare required they can be derived without difficulty alongthe lines described above. To do this one needs to con-struct higher order symmetric polynomials in the wavevectors 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 gra-dient coefficients depending on n ( r ). What is needed todetermine the explicit form of the gradient coefficients inpractice is the determination of the functions N m . Howto do this for N and N is described in the next section. C. Determination of N and N A practical calculation of the coefficients of the gradi-ent expansion of Eq.(14) requires the knowledge of thethe functions N and N . According to Eqs.(4) and (6)these are defined as the first and second density deriva-tives of the correlation function F with respect to thedensity, i.e. δF ( r ′ , r ) δn ( r ′′ ) | n = N ( y , r ′′ − r ) δF ( r ′ , r ) δn ( r ′′ ) δn ( r ′′′ ) | n = N ( y , r ′′ − r , r ′′′ − r ) . To derive useful expressions for these functions it is con-venient to transform derivatives with respect to the den-sity to derivatives with respect to the external potentialas such derivatives appear naturally in perturbation the-ory. We can then use the well-developed tools of pertur-bation theory to calculate these functions. Let us there-fore define the new functions F ( y , r ′′ − r ) = δF ( r ′ , r ) δv ( r ′′ ) | n (17) F ( y , r ′′ − r , r ′′′ − r ) = δ F ( r ′ , r ) δv ( r ′′ ) δv ( r ′′′ ) | n . (18)Since we evaluate these functions for the homogeneouselectron gas they only depend on differences between thespatial coordinates. It will therefore be convenient to alsodefine their Fourier transforms by F ( r , r ) = Z d p d q (2 π ) F ( p , q ) e i p · r + i q · r F ( r , r , r ) = Z d k d p d q (2 π ) F ( k , p , q ) × e i k · r + i p · r + i q · r as well as their partial Fourier transforms F ( y , q ) = Z d p (2 π ) F ( p , q ) e i p · y F ( y , p , q ) = Z d k (2 π ) F ( k , p , q ) e i k · y . The chain rule of differentiation gives a connection be-tween the derivatives with respect to the density and thederivatives with respect to the potential of the correlationfunction F . We have δF ( r , r ) δv ( r ) = Z d r δF ( r , r ) δn ( r ) δn ( r ) δv ( r ) (19) δ F ( r , r ) δv ( r ) δv ( r ) = Z d r δF ( r , r ) δn ( r ) δ n ( r ) δv ( r ) δv ( r )+ Z d r d r δ F ( r , r ) δn ( r ) δn ( r ) δn ( r ) δv ( r ) δn ( r ) δv ( r ) . (20)We see that in these expressions the first and second orderderivative of the density with respect to the potentialappear. We therefore define the linear density responsefunction by χδn ( r ) δv ( r ) | n = χ ( r − r ) = Z d q (2 π ) χ ( q ) e i q · ( r − r ) as well as the second order density response function χ by δ n ( r ) δv ( r ) δv ( r ) | n = χ ( r − r , r − r )= Z d p d q (2 π ) χ ( p , q ) e i p · ( r − r )+ i q · ( r − r ) . Using these definitions we find by Fourier transformingEqs.(19) and (20) that N ( y , q ) = F ( y , q ) χ ( q ) (21) N ( y , p , q ) = (cid:2) F ( y , p , q ) − N ( y , p + q ) χ ( p , q ) (cid:3) χ ( p ) χ ( q ) (22) These equations give the desired relation between thedensity derivatives N and N and the potential deriva-tives F and F of the correlation function F . Relationsbetween the higher order derivatives N m and F m can bederived in a completely analogous way. From Eqs.(21)and (22) we see that to calculate the functions N and N , and hence the gradient expansion coefficients of F to second order in the gradients, we need to calculatethe density response and the change in F to secondorder in a perturbing potential δv . The problem isthereby reduced to doing perturbation theory on theelectron 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 standardperturbation theory to obtain these response functionfor the case that F is the one-particle density matrixand use that to calculate the gradient expansion of theexchange hole. III. GRADIENT EXPANSION OF THEONE-PARTICLE DENSITY MATRIX AND THEEXCHANGE HOLEA. The one-particle density matrix
As an application of our formalism we carry out thegradient expansion of the one-particle density matrix fora noninteracting system with density n ( r ). This prob-lem has received large interest in the past since both theexchange energy E x [ n ] as well the Kohn-Sham kineticenergy T s [ n ] are explicit functionals of such a noninter-acting density matrix [2, 3]. Therefore a gradient expan-sion of this density matrix directly leads to a gradientexpansion of these two functionals. Such a gradient ex-pansion was first carried out by Gross and Dreizler [27]using the Kirzhnits expansion [2, 28]. This expansion isadjusted to the specific form of the noninteracting one-particle density matrix and can not be generalized to ar-bitrary correlation functions. The method presented inthis work can, on the other hand, be applied to arbitrarycorrelation functions, but to demonstrate the soundnessof our formalism we will use it to provide an alternativederivation of the result obtained earlier by Gross andDreizler using the Kirzhnits expansion. An advantage ofour derivation based on nonlinear response theory is thatit avoids the turning point problem encountered in theKirzhnits approach [2].Let us start by defining the one-particle density matrixin second quantization as γ ( r ′ , r ) = X σ h Ψ | ˆ ψ † ( r σ ) ˆ ψ ( r ′ σ ) | Ψ i where σ is a spin coordinate. We will consider the caseof spin-compensated systems. Then, in a noninteractingsystem with 2 N electrons, the density matrix is given interms of one-particle wave functions ϕ j ( r ) by γ ( r ′ , r ) = N X j =1 ϕ j ( r ′ ) ϕ ∗ j ( r ) (23)where the pre-factor 2 is a spin factor. The one-particlestates are eigenstates of a one-particle Schr¨odinger equa-tion (cid:18) − ∇ + v ( r ) (cid:19) ϕ j ( r ) = ǫ j ϕ j ( r ) (24)where v ( r ) is the external potential. To determine thegradient expansion coefficients of γ we need to determinehow the density matrix changes if we make small changes v ( r ) → v ( r ) + δv ( r ) in the external potential. More pre-cisely, we need to calculate the functional derivatives γ ( r , r , r ) = δγ ( r , r ) δv ( r ) γ ( r , r , r , r ) = δ γ ( r , r ) δv ( r ) δv ( r )which are the equivalents of the functions F and F ofEqs.(17) and (18) for the specific choice of correlationfunction F = γ (note that they become functions of twoand three independent vectors respectively when evalu-ated for a homogeneous system). To do this it is sufficientto know the functional derivatives of the orbitals ϕ j andeigenvalues ǫ j with respect to the potential v . These aresimply given by doing first order perturbation theory onthe one-particle Schr¨odinger Eq.(24). We find δǫ j δv ( r ) = | ϕ j ( r ) | (25) δϕ j ( r ) δv ( r ′ ) = ∞ X k = j ϕ k ( r ) ϕ ∗ k ( r ′ ) ϕ j ( r ′ ) ǫ j − ǫ k . (26)With these relations we can differentiate Eq.(23) twicewith respect to the potential. Afterwards we then needto insert the plane wave states ϕ k ( r ) and their eigenval-ues ǫ k appropriate for the homogeneous electron gas into the final expressions. In order not to interrupt the pre-sentation these calculations are given in Appendix A andwe simply present the results of these calculations here.If we define the Fermi factors by n k = θ ( k F − | k | ) , where θ is the Heaviside function and where k F =(3 π n ) / is the Fermi wave vector, then the resultingexpressions are given by γ ( y , q ) = 2 Z d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · y (27) γ ( y , p , q ) = 2 Z d k (2 π ) e i k · y [Φ( k , k + p + q , k + p )+ Φ( k , k + p + q , k + q )] (28)whereΦ( k , p , q ) = 1 ǫ q − ǫ p (cid:18) n q − n k ǫ q − ǫ k − n p − n k ǫ p − ǫ k (cid:19) . (29)In this expression ǫ k = | k | / χ ( q ) and χ ( p , q ) to evaluate thefunctions N and N of Eqs. (21) and (22). Fortunately,for the specific choice F = γ that we made these densityresponse functions are simply given by χ ( q ) = γ ( y = 0 , q ) (30) χ ( p , q ) = γ ( y = 0 , p , q ) . (31)They are therefore obtained directly from Eqs. (27) and(28) so that we do not need to do any extra work tocalculate them.To calculate the gradient coefficient to second order inthe gradients we now need to expand γ and γ to secondorder in the wave vectors p and q . Since these calcula-tions are somewhat lengthy we do not present them hereand refer to Appendix B for a description of the mainsteps. The resulting expansions are given by γ ( y , q ) = − k F π (cid:18) sin zz (cid:20) − i q · y − ( q · y ) (cid:21) − q k cos( z ) (cid:19) (32) γ ( y , p , q ) = 1 π k F (cid:18) cos( z ) − i p + q ) · y cos( z ) + 112 k ( p + q + p · q )(cos( z ) + z sin( z )) −
124 [( p · y ) + ( q · y ) + 3(( p + q ) · y ) ] cos( z ) (cid:19) (33)where we defined z = k F y and q = | q | . The expansions for the first and second order response functions χ and χ of Eqs.(30) and (31) are obtained by taking y = 0 in these expressions. Together which Eqs.(32) and (33) theseexpression then immediately determine the functions N and N from Eqs.(21) and (22). The final result of thiscalculation to second order in the wave vectors is given by N ( y , q ) = sin zz (cid:20) − i q · y − ( q · y ) (cid:21) + q k sin( z ) − z cos( z ) zN ( y , p , q ) = π k (cid:20)(cid:18) cos( z ) − sin( z ) z (cid:19) (cid:18) − i p + q ) · y −
16 [( p · y ) + ( q · y ) ] (cid:19) + 112 k ( p + q + p · q )(3 cos( z ) + z sin( z ) − z ) z ) + 112 ( p · y )( q · y )[4 sin( z ) z − z )] (cid:21) By comparison of these two expression to the expansions (15) and (16) we can directly read off the gradient coefficientfunctions of Eq.(14) to be N ( y ) = − i j ( z ) N ( y ) = 112 k z j ( z ) N ( y ) = − j ( z ) N ( y ) = π k [ z j ( z ) − zj ( z )] N ( y ) = π k [ j ( z ) + 3 zj ( z )]where now we replaced n → n ( r ) and hence the Fermi wave vector k F = (3 π n ( r )) / that appears in these expressionsis a spatially varying function. We further introduced the spherical Bessel functions j ( z ) = sin zz j ( z ) = sin z − z cos zz . By inserting the gradient coefficient function into Eq.(14) we find the following explicit gradient expansion for theone-particle density matrix γ ( r ′ , r ) = k π j ( z ) z + 12 j ( z ) y · ∇ n ( r ) − k z j ( z ) ∇ n ( r ) + 16 k z j ( z )( y y · ∇ ) n ( r ) − π k ( z j ( z ) − zj ( z )) ( ∇ n ( r )) − π k ( j ( z ) + 3 zj ( z )) ( y · ∇ n ( r )) + . . . (34)This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent formas γ ( r ′ , r ) = 1 π h k j ( z ) z − ∇ k k F zj ( z ) + 112 1 k F (cid:16) ∇ · ( ∇ k · y y ) (cid:17) · y y z j ( z ) + 14 ∇ k · y y zj ( z ) −
196 ( ∇ k ) k ( j ( z ) z − zj ( z )) + 132 1 k (cid:16) ∇ k · y y (cid:17) ( z j ( z ) − z j ( z )) i (35)This expression becomes identical in form to that derived by Gross and Dreizler [27] using the Kirzhnits methodafter eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(19) ofreference [27] ). We close this section with a final remark on our result. We note that the gradient expansion to finiteorder breaks the symmetry γ ( r , r ′ ) = γ ( r ′ , r ) ∗ . However, the full gradient expansion as described by Eq.(11) has thissymmetry, which means that the symmetry is restored by taking all higher order gradients into account whenever theseries converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such asthe exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literaturewe 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 andpresent an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler andGross which has the advantage of avoiding the regularization of divergent integrals [27]. In doing so we confirm aresult obtained by Langreth and Perdew [19]. The exchange hole can be calculated directly from the one-particledensity matrix as ρ x ( r + y | r ) = − | γ ( r + y , r ) | n ( r ) . ρ x ( r + y | r ) = − n ( r ) j ( z ) z − j ( z ) j ( z ) z ( y · ∇ n ( r )) + 14 k j ( z ) ∇ n ( r ) − k zj ( z ) j ( z )( y y · ∇ ) n ( r )+ π k j ( z )( zj ( z ) − j ( z ))( ∇ n ( r )) + π k ( j ( z ) j ( z ) z + 3 j ( z ) − j ( z ) )( y · ∇ n ( r )) + . . . (36)As was pointed out by Perdew [9] the second order gradient expansion of the exchange hole does not satisfy thenegativity and sum rule constraints ρ x ( r + y | r ) ≤ Z d y ρ x ( r + y | r ) = − . The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediatelyobvious from Eq.(36). The sum rule can, however, be conveniently analyzed in momentum space. To do this we definethe Fourier transformed exchange hole to be ρ x ( k | r ) = Z d y ρ x ( r + y | r ) e − i k · y (37)This function has the form ρ x ( k | r ) = ρ ( k | r ) + α ( k, n ) ˆ k · ∇ n ( r ) + α ( k, n ) ∇ n ( r ) + α ( k, n ) (ˆ k · ∇ n ( r )) + α ( k, n )( ∇ n ( r )) + α ( k, n ) (ˆ k · ∇ ) n ( r ) + . . . (38)where we defined the unit vector ˆ k = k /k and k = | k | . The explicit form of the functions ρ and α j follow directly byFourier 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 conditionin momentum space translates to the requirement that ρ x ( k = 0 | r ) = −
1. Since ρ ( k = 0 | r ) = − α j ( k = 0 , n ) = 0. However, as is clear from Eqs.(C2)-(C6)in the Appendix this is not satisfied for α , whereas α and α even diverge for k → N h ρ x ( k ) i = Z d r n ( r ) ρ x ( k | r )= N h ρ ( k ) i + ˆ k · Z d r n ( r ) α ( k, n ( r )) ∇ n ( r ) + X i,j =1 Z d r n ( r ) β ij ( k, n ( r )) ∂ i n ( r ) ∂ j n ( r ) . (39)where N is the number of electrons. In this expression h ρ ( k ) i is the system average of ρ ( k | r ) and in the lastterm under the integral sign we defined the tensor β ij ( n, k ) = α L ( n, k ) k i k j k + α T ( n, k ) (cid:16) δ ij − k i k j k (cid:17) . This tensor has a longitudinal part with coefficient α L and a transverse part with coefficient α T which describethe contributions to the system-averaged hole of densitygradients ∇ n parallel and perpendicular to the momen-tum vector k . These coefficients are calculated from the functions α j as α T = α − n ∂ ( nα ) ∂nα L = α T + α − n ∂ ( nα ) ∂n , 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 explicitexpressions for α L , T . If we define the dimensionless vari-1able ¯ k = k/ (2 k F ) then we have α T = π nk Z T (¯ k ) α L = π nk Z L (¯ k ) , where the functions Z L , T have the explicit form Z T ( x ) = − x θ (1 − x ) + 43 δ ( x − Z L ( x ) = − x θ (1 − x ) + δ ( x −
1) + 13 δ ′ ( x − k → Z L , T (¯ k ) = 0and hence β ij ( k, n ) → k →
0. This implies thatthe last term in Eq.(39) does not contribute to the sumrule of the system averaged hole. The same conclusioncan be derived for the second term in Eq.(39). Since α ( k = 0 , n ) = i/ (4 nk F ) (see Eq.(C2)) is a local func-tion of the density the integrand of the first integral inEq.(39) is a total derivative and vanishes (assuming ei-ther vanishing densities at infinity or periodic boundaryconditions). We therefore find thatlim k → h ρ x ( k ) i = lim k → h ρ ( k ) i = − . This implies that system averaged exchange hole ob-tained from the second order gradient expansion doesfulfill the sum rule, although the exchange hole itself doesnot. We note, however, that this is only achieved by in-tegrating over both positive and negative contributions.When the positive contributions to the exchange hole areremoved the sum rule is only recovered for a finite hole ra-dius [9, 10]. We finally note that with expression Eq.(39)it is now a simple matter to calculate the exchange energyfrom E x [ n ] = N Z d k (2 π ) h ρ x ( k ) i π | k | . We insert Eq. (39) and do the angular integrations first.It is therefore convenient to define the spherical and sys-tem averaged hole in momentum space by hh ρ x ( k ) ii = 14 π Z d Ω k h ρ x ( k ) i such that E x [ n ] = Nπ Z ∞ dk hh ρ x ( k ) ii . (40) From Eq.(3) we then find that N hh ρ x ( k ) ii = N hh ρ ( k ) ii + Z d r n ( r )[ 13 α L ( k, n ) + 23 α T ( k, n )]( ∇ n ) = N hh ρ ( k ) ii + Z d r π k Z x (¯ k )( ∇ n ) , (41)where we defined Z x (¯ k ) = 13 Z L (¯ k ) + 23 Z T (¯ k ) . (42)The exchange energy is obtained by inserting the expres-sion for the averaged hole of Eq.(41) into Eq.(40). Thisthen finally gives the expression E x [ n ] = − (cid:16) π (cid:17) Z d r n ( r )+ Z d r B x ( n ( r ))( ∇ n ( r )) (43)where B x ( n ) = π k Z ∞ d ¯ k Z x (¯ k )= − π (3 π ) / n / and where to calculate the first term in Eq.(43) we usedEq.(C1). The coefficient B x is the same gradient coeffi-cient as obtained by Gross and Dreizler [27] and earlierby Sham [30]. However, the correct analytic exchangegradient coefficient is known [31–35] to be a factor 10/7larger. The reason for this discrepancy is clearly de-scribed by Svendsen and von Barth [34] who showed thatthe Coulomb interaction is too singular to allow for theinterchange of the operations of integration and the ex-pansion in wave vectors. The problem does, for instance,not appear when Yukawa screened Coulomb interactionsare used [34]. The conclusion is that there is nothingwrong with the Kirzhnits method, or the nonlinear re-sponse theory derivation of the gradient expansion of theexchange hole presented here, but one has be aware thatfor Coulomb interactions the procedures of expandingcorrelation functions in terms of density gradients fol-lowed by integrations may not yield the same result asdirectly expanding the integrated quantity in terms ofdensity gradients. For this reason the original GGA ofPerdew and Wang [10] based on the gradient expansionof the exchange hole was later reparametrized [16, 36]to accomodate the correct gradient coefficient for the ex-change energy. IV. CONCLUSIONS AND OUTLOOK
We derived a general scheme based on nonlinear re-sponse theory to calculate the density gradient expan-sion of general correlation functions and showed that in2order to express the gradient coefficients in terms of thefull density profile summations to infinite order must becarried out over response functions of arbitrarily large or-der. A consistency condition was derived to do this. Theformalism was used to derive the gradient expansion ofthe one-particle density matrix and the exchange hole tosecond order in the gradients. We confirm the derivationof Dreizler and Gross that used the Kirzhnits expansionmethod. We further analyzed the exchange hole in mo-mentum space to derive that the system averaged hole tosecond order in the gradients satisfies the sum rule andto derive the gradient expansion of the exchange energywithout the need to regularize divergent integrals.The scheme that we presented is very general and canbe applied to more general correlation functions. A nat-ural application of the scheme would be to calculate thegradient expansion of the correlation hole ρ c . Regardingthe gradient terms of the correlation hole little is known.We essentially only have some detailed information onthe long-range properties of the spherically and systemaveraged hole. This information comes from the work ofLangreth and Perdew that calculated hh ρ c ( k ) ii within theRPA. In our notation their results (see also Appendix Cof [19]) can be summarized as N hh ρ c ( k ) ii = N hh ρ ( k ) ii + Z d r π k z c ( n, k )( ∇ n ) where the function z c has been parametrized by Langrethand Mehl [20, 21] as z c ( n, k ) = 4 √ k s e − √ kk s (44)where k s = (4 k F /π ) / is the Thomas-Fermi or screeningwave vector. We can transform to real space to obtainthe following expression for the spherically and systemaveraged correlation hole N hh ρ c ( y ) ii = N hh ρ ( y ) ii + Z d r k k (1 + ( k s y ) / ( ∇ n ) (45)We see from Eq.(44) that z c ( n, k = 0) = 0 and as a con-sequence the sum rule hh ρ c (0) ii = 0 for the correlationhole is not satisfied. We know, however, that the RPAbecomes accurate in the high density limit and since fromEq.(45) the gradient coefficient of the correlation hole de-pends on k s y we see that high densities are equivalent tolarge separations y . Similarly, low density correspondsto short-distance behavior. However, RPA can not betrusted in this region and we have no precise knowledgeon the small distance behavior of the gradient coefficientof hh ρ c ( y ) ii . A model for the general short-range behav-ior was proposed [16–18] on physical arguments (it shouldnot affect the Coulomb cusp of and the on-top value ofthe LDA hole) after which the real-space cutoff proce-dure 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 thecalculation of the short range properties of the gradientcoefficient of the correlation hole. As follows from ourderivation this requires the knowledge of the functions(17) and (18), or at least their expansion to second orderin wave vectors, for the case that F represents the paircorrelation function or the correlation hole of the electrongas. It may therefore be worthwhile to use our currentscheme to explore these response functions beyond theRPA. For short range correlation an approach based ondiagrammatic summation of ladder diagrams suggests it-self. Of course, for the development of density function-als for general systems beyond the weakly inhomogeneousregime it is not sufficient to use the straightforward gra-dient expansion. However, general short-range features,such as the way the exchange-correlation hole deformsclose to the reference electron in the presence of a densitygradient can be transferred to more general systems thanthe weakly inhomogeneous ones. Perhaps in the combi-nation with truly nonlocal functionals for the exchange-correlation hole [37] this can lead to the development ofmore accurate density functionals. This will be a topicof our future research. Acknowledgments
I would like to thank the students in the course ”Den-sity Functional Theory” at the University of Jyv¨askyl¨afor inspiration and Esa R¨as¨anen for useful discussionsand Klaas Giesbertz for checking a large part of the equa-tions .
Appendix A: Calculation of γ and γ In this Appendix we will derive in the expression forthe first and second derivative of the one-particle densitymatrix γ with respect to the potential v . For a clearerinterpretation and compactness of notation we will onlyput a comma between the variables representing the orig-inal coordinates of the density matrix and the variablesthat appear as argument of the potential variations. Di-rect differentiation of Eq.(23) using Eq.(26) then gives γ (23 ,
1) = δγ (23) δv (1)= 2 ∞ X i,j ( f j − f i ) ϕ j (1) ϕ ∗ i (1) ϕ i (2) ϕ ∗ j (3) ǫ j − ǫ i (A1)where we used the short notation j = r j and introducedthe occupation numbers f j = 1 for an occupied state and f j = 0 for an unoccupied state. Inserting plane wavestates appropriate for the electron gas we find that γ ( r r , r ) = 2 Z d q d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · ( r − r )+ i q · ( r − r ) . γ has the simpleform γ ( p , q ) = 2 n p + q − n p ǫ p + q − ǫ p in Fourier space. This is precisely the integrand ofEq.(27). The function γ can be calculated by straight-forward differentiation of Eq.(A1) with respect to the po-tential γ (23 ,
14) = δγ (23 , δv (4) . To obtain the explicit form of this function we have to dif-ferentiate both the orbitals and the eigenvalues using Eqs.(25) and (26). Let us denote the term that arises fromthe differentiating the orbitals by X and the term thatarises from differentiating the eigenvalues by Y . Then wefind that γ (23 ,
14) = X (23 ,
14) + Y (23 , , where X (23 ,
14) = 2 ∞ X i,j,k ϕ i (2) ϕ ∗ j (3) h Φ( ijk ) ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ Φ( jik ) ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i (A2) Y (23 ,
14) = − ∞ X i,j ϕ i (2) ϕ ∗ j (3)Φ( iji )( | ϕ j (4) | − | ϕ i (4) | ) ϕ j (1) ϕ ∗ i (1) (A3)and where we further definedΦ( ijk ) = 1 ǫ k − ǫ j (cid:18) f k − f i ǫ k − ǫ i − f j − f i ǫ j − ǫ i (cid:19) . (A4)The function Φ has the useful properties Φ( ijj ) = 0 andΦ( ijk ) = Φ( ikj ). The function γ (23 ,
14) = γ (23 ,
41) issymmetric in the indices 4 and 1 due to the fact that thedifferentiations 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 γ in such a way that this symmetry is obvious. To do thiswe first note that since Φ( iji ) = − Φ( jij ) the terms with k = i and k = j in Eq.(A2) sum up to a function Z of asimilar form as Y in Eq.(A3), namely Z (23 ,
14) = − ∞ X i,j ϕ i (2) ϕ ∗ j (3)Φ( iji )( | ϕ j (1) | − | ϕ i (1) | ) ϕ j (4) ϕ ∗ i (4) . (A5) We can therefore write γ (23 ,
14) as γ (23 ,
14) = 2 ∞ X ijk,k =( i,j ) ϕ i (2) ϕ ∗ j (3) h Φ( ijk ) ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ Φ( jik ) ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i + Y (23 ,
14) + Z (23 ,
14) (A6)It is easily seen that the sum of Y and Z is symmetricunder interchange of 1 and 4. However, this is not obvi-ous in the first term of the equation above since at firstsight Φ( jik ) does not appear to be equal to Φ( ijk ) sinceΦ( jik ) = 1 ǫ k − ǫ i (cid:18) f k − f j ǫ k − ǫ j − f i − f j ǫ i − ǫ j (cid:19) (A7)which seems to be different from expression (A4). How-ever, this is just appearance. The reason is that the oc-cupations can only attain the values zero and one. Forthe case f i = f j = 1 or f i = f j = 0 we see directly thatEq.(A4) and (A7) attain the same value. A little inspec-tion shows that this is also true for cases f i = 0 , f j = 1and f i = 1 , f j = 0. We therefore find for k = i, j thatΦ( ijk ) = Φ( jik ). Therefore Eq.(A6) can be simplified to γ (23 ,
14) = 2 ∞ X ijk,k =( i,j ) ϕ i (2) ϕ ∗ j (3)Φ( ijk ) h ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i + Y (23 ,
14) + Z (23 ,
14) (A8)This expression is now explicitly symmetric in the indices1 and 4. Let us now evaluate this expression for thehomogeneous electron gas. In the electron gas the one-particle states are plane waves ϕ k ( r ) = e i k · r √ V where V is the volume of the system. Since | ϕ k | = 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 γ (23 ,
14) = 2 V X k , p , q , q =( k , p ) Φ( k , p , q ) e i ( k · r − p · r ) h e i ( q − k ) · r + i ( p − q ) · r + e i ( p − q ) · r + i ( q − k ) · r i where we definedΦ( k , p , q ) = 1 ǫ q − ǫ p (cid:18) n q − n k ǫ q − ǫ k − n p − n k ǫ p − ǫ k (cid:19) . In this expression n p = θ ( ǫ F − ǫ p ) is the zero-temperatureFermi function and ǫ F = k / q = k , p is a region of measure zero and we can therefore freelyintegrate over all wave vectors. After a few substitutionswe obtain γ (23 ,
14) = Z d k d p d q (2 π ) γ ( k , p , q ) e i k · ( r − r )+ i p · ( r − r )+ i q · ( r − r ) where we defined γ ( k , p , q ) = 2[Φ( k , k + p + q , k + p )+Φ( k , k + p + q , k + q )]This gives the integrand of Eq.(28). Appendix B: Expansion of γ and γ
1. Expansion of γ To determine γ to second order in the wave vectorswe have to expand the function γ ( y , q ) = 2 Z d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · y (B1)to second order in powers of q . This can be done byexpanding the integrand in powers of q . If we denote∆ = ǫ p + q − ǫ p = p · q + q / q = | q | , then wecan expand the Fermi function n p + q in a distributionalTaylor series as n p + q = n p + ∆ dndǫ | ǫ p + ∆ d ndǫ | ǫ p + ∆ d ndǫ | ǫ p + . . . Since n p = θ ( ǫ F − ǫ p ) is the Heaviside function we have n p + q = n p − ∆ δ ( ǫ F − ǫ p ) + ∆ δ ′ ( ǫ F − ǫ p ) − ∆ δ ′′ ( ǫ F − ǫ p ) + . . . Inserting this into Eq.(B1) then gives γ ( y , q ) = − Z d p (2 π ) δ ( ǫ F − ǫ p ) e i p · y + Z d p (2 π ) ∆ δ ′ ( ǫ F − ǫ p ) e i p · y − Z d p (2 π ) ∆ δ ′′ ( ǫ F − ǫ p ) e i p · y + . . . (B2)We see that in these integrals several derivatives of thedelta function appear. We now use the general mathe-matical expression δ ( n ) ( y ( x )) = dydx ddx ! n X i δ ( x − x i ) | dydx ( x i ) | where the sum runs over all zeros of the function y ( x ), i.e. y ( x i ) = 0. Using this equation in spherical coordinateswith y ( p ) = ( k − p ) / δ ( ǫ F − ǫ p ) = 1 k F δ ( p − k F ) (B3) δ ′ ( ǫ F − ǫ p ) = − pk F δ ′ ( p − k F ) (B4) δ ′′ ( ǫ F − ǫ p ) = − δ ′ ( p − k F ) k F p + δ ′′ ( p − k F ) k F p (B5) δ ′′′ ( ǫ F − ǫ p ) = − k F h p δ ′′′ ( p − k F ) − p δ ′′ ( p − k F )+ 3 p δ ′ ( p − k F ) i (B6)where in Eq.(B6) we also evaluated the third derivative ofthe delta function as we will need it later in the expansionof γ . With Eqs.(B3)-(B5) we can readily evaluate thethree integrals in Eq.(B2). Let us denote these integralsby A, B and C . Then we find for the first integral A = − k F Z d p (2 π ) δ ( p − k F ) e i p · y = − k F π sin zz where z = k F y . The second integral gives B = Z d p (2 π ) ( p · q + q / δ ′ ( ǫ F − ǫ p ) e i p · y = ( q / − i q · ∇ y ) Z d p (2 π ) δ ′ ( ǫ F − ǫ p ) e i p · y = ( q / − i q · ∇ y ) 12 π k F cos( k F y )= q π k F cos z + i q · y π k F sin zz where in the second step we used Eq.(B4). It thereforeremains to calculate the last integral C = − Z d p (2 π ) ( p · q + q / δ ′′ ( ǫ F − ǫ p ) e i p · y However, up to order q it is suffices to calculate C = − Z d p (2 π ) ( p · q ) δ ′′ ( ǫ F − ǫ p ) e i p · y = 13 ( q · ∇ y ) Z d p (2 π ) δ ′′ ( ǫ F − ǫ p ) e i p · y = −
13 ( q · ∇ y ) π k (cos( z ) + z sin( z ))= − q π k F cos( z ) + ( q · y ) π k F sin( z ) z where in the second step we used Eq.(B5). Adding theresults γ = A + B + C gives the expression (32).5
2. Expansion of γ Analogously to γ we can expand γ ( y , p , q ) of Eq.(28)in powers of p and q . To reduce our effort we note thatwe can write γ ( y , p , q ) = A ( y , p , q ) + A ( y , q , p )where A ( y , p , q ) = 2 Z dk (2 π ) e i k · y Φ( k , k + p + q , k + p )We therefore only need to expand A ( y , p , q ) in powersof the wave vectors and symmetrize with respect to p and q afterwards. To do this we first need to expand thefunctionΦ( k , k + p + q , k + p ) =1 ǫ k + p − ǫ k + p + q (cid:20) n k + p − n k ǫ k + p − ǫ k − n k + p + q − n k ǫ k + p + q − ǫ k (cid:21) . If we denote∆ = ǫ k + p − ǫ k = p p · k ∆ = ǫ k + p + q − ǫ k = ( p + q ) p + q ) · k , we find by expanding the Fermi functions in a distribu-tional Taylor series thatΦ( k , k + p + q , k + p ) = 12 δ ′ ( ǫ F − ǫ k ) −
16 (∆ + ∆ ) δ ′′ ( ǫ F − ǫ k )+ 124 (∆ + ∆ ∆ + ∆ ) δ ′′′ ( ǫ F − ǫ k ) + . . . With this expression we find that the function A ( y , p , q )can be calculated as A ( y , p , q )= Z d k (2 π ) δ ′ ( ǫ F − ǫ k ) e i k · y − Z d k (2 π ) δ ′′ ( ǫ F − ǫ k )(∆ + ∆ ) e i k · y + 112 Z d k (2 π ) δ ′′′ ( ǫ F − ǫ k )(∆ + ∆ ∆ + ∆ ) e i k · y + . . . The three integrals appearing on the left hand side of thisexpression are sufficient to extract the powers to secondorder in p and q . Let us call these three integrals A , A and A in order of appearance. The first integral gives A = Z d k (2 π ) δ ′ ( ǫ F − ǫ k ) e i k · y = cos z π k F . The second integral is given by A = − Z dk (2 π ) δ ′′ ( ǫ F − ǫ k ) (cid:20) p + ( p + q ) p + q ) · k (cid:21) e i k · y = (cid:20) − p + ( p + q ) i p + q ) · ∇ y (cid:21) Z dk (2 π ) δ ′′ ( ǫ F − ǫ k ) e i k · y = (cid:20) − p + ( p + q ) i p + q ) · ∇ y (cid:21) cos( z ) + z sin( z ) − π k = cos( z ) + z sin( z )12 π k ( p + ( p + q ) ) − i π k F (2 p + q ) · y cos( z )where we used Eq.(B5) in the second step. Finally the third integral has the explicit form A = 112 Z dk (2 π ) h ( p p · k ) + ( p p · k )( ( p + q ) p + q ) · k )+ ( ( p + q ) p + q ) · k ) i δ ′′′ ( ǫ F − ǫ k ) e i k · y p and q . Therefore, it is sufficient to calculate A = 112 Z dk (2 π ) (cid:2) ( p · k ) + ( p · k )(( p + q ) · k ) + (( p + q ) · k ) (cid:3) δ ′′′ ( ǫ F − ǫ k ) e i k · y = − (cid:2) ( p · ∇ y ) + ( p · ∇ y )(( p + q ) · ∇ y ) + (( p + q ) · ∇ y ) (cid:3) Z dk (2 π ) δ ′′′ ( ǫ F − ǫ k ) e i k · y = − (cid:2) ( p · ∇ y ) + ( p · ∇ y )(( p + q ) · ∇ y ) + (( p + q ) · ∇ y ) (cid:3) (cid:18) z ) − z cos( z ) + 3 z sin( z )2 π k (cid:19) where to evaluate the integral over the delta function we used Eq.(B6). To perform the derivatives in the last termwe use that for an arbitrary function f ( z )( p · ∇ y )( q · ∇ y ) f ( z ) = ( p · y )( q · y ) k y (cid:20) d fdz − z dfdz (cid:21) + ( p · q ) k F y dfdz With this expression we find that A = − π k ( p + p · ( p + q ) + ( p + q ) )(cos( z ) + z sin( z )) − π k F [( p · y ) + ( p · y )(( p + q ) · y ) + (( p + q ) · y ) ] cos( z )Collecting our results A = A + A + A we therefore obtain the following expression for γ ( y , p , q ) γ ( y , p , q ) = A ( y , p , q ) + A ( y , q , p )= 1 π k F (cid:18) cos( z ) − i p + q ) · y cos( z ) + 112 k ( p + q + p · q )(cos( z ) + z sin( z )) −
124 [( p · y ) + ( q · y ) + 3(( p + q ) · y ) ] cos( z ) (cid:21) 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 Fouriertransforming Eq.(36). If we define the dimensionless quantity ¯ k = | k | / (2 k F ( r )) then ρ x is given by ρ x ( k | r ) = ( − k −
12 ¯ k ) θ (1 − ¯ k ) (C1)which is simply the Fourier transform of the LDA hole. We see that ρ x (0 | r ) = − α j are given by α = i nk F θ (1 − ¯ k ) (C2) α = − ¯ k nk θ (1 − ¯ k ) (C3) α = − n k (cid:16) θ (1 − ¯ k )¯ k + 14 δ (¯ k −
1) + 34 δ ′ (¯ k − (cid:17) (C4) α = 172 n k (cid:16) k θ (1 − ¯ k ) − δ (¯ k − (cid:17) (C5) α = 124 nk (cid:16) θ (1 − ¯ k )¯ k + δ (¯ k − (cid:17) (C6) [1] P. Hohenberg and W. Kohn,Phys. Rev. , B864 (1964) [2] R. M. Dreizler and E.K.U. Gross, Density Functional Theory: An Approach to the Quantum Many-Body Prob-lem
Springer (Berlin) (1990)[3] E. Engel and R. M. Dreizler,
Density Functional Theory:An Advanced Course
Springer (Heidelberg) (2011)[4] U. von Barth,
Basic Density-Functional Theory- anOverview , Phys. Scr. T , 9 (2004)[5] C. A. Ullrich,
Time-Dependent Density-Functional The-ory: Concepts and Applications
Oxford University Press,Oxford (2012)[6] W. Kohn and L. J. Sham,Phys. Rev. , A1133 (1965)[7] C. - O. Almbladh, Technical Report, University of Lund(1972)[8] O. Gunnarsson and B. I. Lundqvist,Phys. Rev. B , 4274 (1976)[9] J. P. Perdew, Phys. Rev. Lett. , 1665 (1985), Erra-tum: Phys. Rev. Lett. , 2370 (1985)[10] J. P. Perdew and W. Yue, Phys. Rev. B , 8800 (1986),Erratum: Phys. Rev. B40, 3399 (1989)[11] A. D. Becke, J. Chem. Phys. , 7184 (1986).[12] A. D. Becke, J. Chem. Phys. , 1053 (1988).[13] A. D. Becke and M. R. Roussel,Phys. Rev. A , 3761 (1989)[14] V. Tschinke and T. Ziegler,Can. J. Chem. , 460 (1989).[15] J. P. Perdew and S. Kurth, ”Density Functionals forNon-relativistic Coulomb Systems in the New Century” inLecture Notes in Physics , pp 1-55, (2003), Editors:C. Fiolhais, F. Nogueira, M. A. L. Marques, Springer(Heidelberg).[16] J. P. Perdew, in Electronic Structure of Solids ’91 , pp.11-20, Eds. P. Ziesche and H. Eschrig, Akademie Verlag,Berlin (1991)[17] K. Burke, J. P. Perdew and Y. Wang, in
Electronic Den-sity Functional Theory: Recent Progress and New Direc-tions , p.81 , Eds. J.F. Dobson, G. Vignale and M. P.Das, Plenum, New York, (1997)[18] J. P. Perdew, K. Burke and Y. Wang,Phys. Rev. B , 16533 (1996), Erratum:Phys. Rev. B57, 14999 (1996) [19] D. C. Langreth and J. P. Perdew,Phys. Rev. B , 5469 (1980)[20] D. C. Langreth and M. J. Mehl,Phys. Rev. Lett. , 446 (1981)[21] D. C. Langreth and M. J. Mehl,Phys. Rev. B , 1809 (1983), ErratumPhys. Rev. B , 2310 (1984)[22] P. S. Svendsen and U. von Barth,Phys. Rev. B , 17402 (1996)[23] M. Springer, P. S. Svendsen and U. von Barth,Phys. Rev. B , 17392 (1996)[24] R. van Leeuwen, Adv. Quant. Chem. , 25 (2003)[25] H. Weyl, The classical groups: Their invariants and rep-resentations , Princeton University Press (1939)[26] G. Stefanucci and R. van Leeuwen,
NonequilibriumMany-Body Theory of Quantum Systems , CambridgeUniversity Press, Cambridge (2013)[27] E. K. U. Gross and R. M. Dreizler,Z. Phys. A , 103 (1981)[28] A. Putaja, E. R¨as¨anen, R. van Leeuwen, J. G. Vilhenaand M. A. L. Marques, Phys. Rev. B , 165101 (2012)[29] J. J. Duistermaat and J. A. C. Kolk, Distributions: The-ory and Applications , Springer, Dordrecht (2010)[30] L. J. Sham in
Computational methods in Band Theory ,eds. P. M. Marcus, J. F. Janak and A. R. Williams(Plenum, New York, 1971)[31] P. R. Antoniewicz and L. Kleinman,Phys. Rev. B , 6779 (1985)[32] L. Kleinman and S. Lee, Phys. Rev. B , 4634 (1988)[33] E. Engel and S. H. Vosko, Phys. Rev. B , 4940 (1990),Erratum: Phys. Rev. B , 1446 (1991)[34] P. S. Svendsen and U. von Barth,Int. J. Quant. Chem. , 351 (1995)[35] D. C. Langreth and S. H. Vosko,Adv. Quant. Chem. , 175 (1990)[36] J. P. Perdew, Physica B , 1 (1991)[37] K. J. H. Giesbertz, R. van Leeuwen and U. von Barth,Phys. Rev. A , 022514 (2013) r X i v : . [ c ond - m a t . m t r l - s c i ] A p r The density gradient expansion of correlation functions
Robert van Leeuwen
1, 2 Department of Physics, Nanoscience Center, University of Jyv¨askyl¨a, Survontie 9, 40014 Jyv¨askyl¨a, Finland European Theoretical Spectroscopy Facility (ETSF) (Dated: October 30, 2018)We present a general scheme based on nonlinear response theory to calculate the expansion ofcorrelation functions such as the pair-correlation function or the exchange-correlation hole of aninhomogeneous many-particle system in terms of density derivatives of arbitrary order. We furtherderive a consistency condition that is necessary for the existence of the gradient expansion. Thiscondition is used to carry out an infinite summation of terms involving response functions up toinfinite order from which it follows that the coefficient functions of the gradient expansion can beexpressed in terms the local density profile rather than the background density around which theexpansion is carried out. We apply the method to the calculation of the gradient expansion of theone-particle density matrix to second order in the density gradients and recover in an alternativemanner the result of Gross and Dreizler (Z. Phys. A , 103 (1981)) which was derived usingthe Kirzhnits method. The nonlinear response method is more general and avoids the turningpoint problem of the Kirzhnits expansion. We further give a description of the exchange hole inmomentum space and confirm the wave vector analysis of Langreth and Perdew (Phys. Rev. B ,5469 (1980)) for this case. This is used to derive that the second order gradient expansion of thesystem averaged exchange hole satisfies the hole sum rule and to calculate the gradient coefficientof the exchange energy without the need to regularize divergent integrals. PACS numbers: 31.15.E-, 71.15.Mb, 31.15.eg
I. INTRODUCTION
Since the groundbreaking work of Hohenberg andKohn [1] we know that the external potential of any inho-mogeneous quantum many-body system is a functional ofits ground state density n . This implies that the many-body ground state | Ψ[ n ] i and hence any ground stateexpectation value A [ n ] = h Ψ[ n ] | ˆ A | Ψ[ n ] i (1)for any operator ˆ A is a functional of the density. Thisidea has inspired an enormous amount of work in a re-search field that is now known as density-functional the-ory. Density-functional theory [2–4] is currently one ofthe main approaches used in electronic structure the-ory. Over the years many extensions beyond its origi-nal formulation have been developed and currently it iswidely applied in solid state physics and quantum chem-istry, both for ground state and time-dependent systems[5]. Especially a large activity has gone into finding ex-plicit expressions for the ground state energy E [ n ] as afunctional of the density for many-electron systems. Theground state energy is obtained by minimizing this func-tional on an appropriate set of electronic densities, whichis of great importance in determining structures and ge-ometries of molecules and solids. By use of the Kohn-Sham method [6] this minimization problem is reducedto solving effective one-particle equations. The construc-tion of the effective potential in these equations requiresthe knowledge of the so-called exchange-correlation en-ergy functional E xc [ n ]. This functional can be expressed as [7, 8] E xc [ n ] = 12 Z d r d r ′ n ( r ) ρ xc ( r ′ | r ) | r − r ′ | (2)where ρ xc ( r ′ | r ) is the coupling constant integratedexchange-correlation hole. The exchange-correlation holehas a physical interpretation as the difference betweenthe conditional and the unconditional probability (whichis simply the density n ( r ′ )) to find a electron at r ′ , giventhat we know that there is an electron at r . Equation (2)is an important relation since it expresses the exchange-correlation energy in terms of a quantity that has a lo-cal physical interpretation and can be studied by accu-rate wave function and many-body methods or modeledbased on physical intuition. It has therefore played animportant role in the development of approximate den-sity functionals [8–14]. One of the simplest and alreadyvery successful approximations has been the local densityapproximation (LDA) in which the exchange-correlationhole is taken from the homogeneous electron gas and ap-plied locally [8] to systems with arbitrary density profiles.The accuracy of the LDA has been considerably improvedby means of the so-called generalized gradient approxi-mations (GGAs) [15]. A large class of commonly usedGGAs is based on the so-called real-space cutoff of thestraightforward gradient expansion for the exchange andcorrelation hole. The first such functional for exchangewas proposed by Perdew [9]. He noted that the gradientexpansion of the exchange hole to second order in thedensity gradients violates the negativity condition andthe hole sum rule by which the exchange hole integratesto a deficit of one electron. By enforcing these physicalconstraints a density functional was obtained that greatlyimproved on the LDA for binding energies of molecules.This procedure was later simplified [10] using the factthat the exchange correlation energy (2) can be writtenas E xc [ n ] = N Z d y h ρ xc ( y ) i| y | where N is the number of electrons in the system and h ρ xc ( y ) i = 1 N Z d r n ( r ) ρ xc ( r + y | r ) , (3)is the so-called system-averaged exchange-correlationhole. This averaged hole still obeys the negativity con-dition as well as the sum rule and can therefore be sub-jected to the real-space cut-off procedure. One advan-tage of using the system averaged holes was to reducethe order of the derivatives in the gradient expansion.Furthermore it also allowed for the real-space cut-off pro-cedure to be applied to the correlation hole since there isa known gradient expansion for the system averaged cor-relation 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–18]. These GGAs relied heavily on the pioneeringwork by Langreth, Perdew and Mehl [19–21] who per-formed a wave-vector decomposition of the system aver-aged hole of Eq. (3) and calculated the Fourier transform h ρ xc ( k ) i = Z d y h ρ xc ( y ) i e − i k · y , from which the exchange correlation energy can be cal-culated as E xc [ n ] = N Z d k (2 π ) h ρ xc ( k ) i π | k | . We are not aware of any work that goes beyond refer-ences [19–21] and improves on the straightforward gra-dient expansion of the system averaged correlation hole.The knowledge on the gradient expansion of the correla-tion hole is very limited indeed. As mentioned before, incontrast to the work on the exchange hole [2, 9], thereare no known expressions for the gradient expansion ofthe non-system averaged correlation hole.Part of the problem is that there has been no clear deriva-tion on how to expand the expectation value of a gen-eral operator as in Eq.(1) in terms of density gradients.A well-known expansion method for a two-point func-tion is the Kirzhnits expansion [2] which is, however,specifically adapted to expanding the noninteracting one-particle density matrix in powers of density gradients andcan not be used for more general correlation functions.The first goal of this paper is, therefore, to present ascheme based on nonlinear response theory that can beused to expand general correlation functions (such as thecorrelation 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 stan-dard gradient expansion for a global quantity, namelythe exchange-correlation energy E xc . The starting pointof any gradient expansion is the consideration of a den-sity profile n ( r ) = n + δn ( r ) which consists of a smalldensity variation δn ( r ) around a homogeneous density n . By considering the limit of slow density variationsone then finds that the lowest gradient correction to theexchange-correlation energy is an integral over a func-tion of the form B xc ( n )( ∇ δn ( r )) where the coefficient B xc ( n ) is calculated from the static exchange-correlationkernel of the homogeneous electron gas of density n [4].Since ∇ δn ( r ) = ∇ n ( r ) we can replace the density vari-ation under the derivative operator by the full densityprofile n ( r ). However, the coefficient B xc ( n ) still de-pends on the background density n . This is a problemin application of the formula to general inhomogeneoussystems such as molecules or surfaces in which a back-ground density cannot be unambiguously defined (assum-ing that a low order gradient expansion makes sense insuch systems [22, 23]). The usual argument to get ridof the background dependence, is that the replacement B xc ( n ) → B xc ( n + δn ( r )) can be made since the differ-ence between these two quantities is or order δn ( r ). How-ever, it is not clear that this is consistent with gradientcoefficients that arise from terms of order ( δn ( r )) . Thepoint was discussed clearly by Svendsen and von Barth[22] who checked that the replacement of n by the fulldensity profile was consistent to third order in the densityvariation. This derivation was based on a specific rela-tion between response functions that describe the changein the exchange energy to second and third order in thedensity variations. In reference [24] this derivation wasextended by showing that the replacement of n by thefull density profile is consistent to all orders in the den-sity variation. In this paper we will show that this isthe case as well for the gradient expansion of correlationfunctions rather than scalar functions. This relies on cer-tain relations between the response functions that mustbe satisfied for the gradient expansion to exist.The paper is divided as follows. In section II we de-rive the basic equations and consistency conditions forthe density gradient expansion of correlation functions.In section III we derive the gradient expansion of theone-particle density matrix for a noninteracting systemwith density n ( r ) (which is therefore equal to the den-sity matrix of the Kohn-Sham system) and discuss itssymmetry properties. We further calculate the gradientexpansion of the exchange-hole to second order in thedensity gradients, both in real and in momentum space.The momentum space description is used to demonstratethat the system averaged exchange hole satisfies the sumrule (but not the negativity condition), and to calculatethe gradient expansion of the exchange energy withoutthe need to regularize divergent integrals. Finally in sec-tion IV we present our conclusions and outlook. II. GRADIENT EXPANSION OFCORRELATION FUNCTIONSA. Expansion in density variations
Let us for a system with ground state density n ( r ) con-sider an arbitrary correlation function F [ n ]( r ′ , r ). Sincethe correlation function can be calculated from the many-body ground state, by virtue of the Hohenberg-Kohn the-orem this function is a functional of the density. At thispoint it will not be important what the specific form ofthe correlation function is, as the structure of its gradientexpansion will be independent of its specific form. Corre-lation functions of common interest are, for instance, thepair-correlation function, the exchange-correlation holeor the one-particle density matrix.When the density is constant in space, i.e. when n ( r ) = n , we are describing the homogeneous electrongas and due to translational and rotational invariancewe have that F [ n ]( r ′ , r ) = F ( n , | r − r ′ | ) where F is a function of the homogeneous density and the dis-tance between its spatial arguments. When the density n ( r ) = n + δn ( r ) deviates from the constant density n by a small amount δn ( r ) we can expand F in powers ofthis density variation. To derive compact expressions wefirst introduce the notation r m = ( r , . . . , r m ) d r m = d r . . . d r m for the collection of m position vectors and the corre-sponding integration volume elements. We further define δn ( r m ) = m Y j =1 δn ( r j ) . With these conventions the expansion of F in densityvariations is given by F [ n ]( r ′ , r ) = F ( n , | r ′ − r | )+ ∞ X m =1 m ! Z d r m F m ( n ; r ′ , r , r m ) δn ( r m )where we defined F m ( n ; r ′ , r , r m ) = δ m F ( r ′ , r ) δn ( r ) . . . δn ( r m ) | n . (4)The n dependence of the functions F m will be importantfor the gradient expansion, but to shorten notation wesuppress this dependence in the notation for the time be-ing and will reintroduce it later. Here we assume that thederivatives F m exist, which means that we assume that F is a smooth functional of the density. In the electrongas two Taylor expansions around two different values of n have very likely the same mathematical convergenceproperties since the two systems have the same physicalproperties (unless we are close to a very low density cor-responding to a Wigner crystal transition). The point n is therefore unlikely to be a non-analytic point. Whenthe derivatives F m do exist we can therefore expect theTaylor series (4) to converge whenever | δn ( r ) | /n ≪ | δn ( r ) | implies the use of a supremum norm or possiblyan L p -integral norm (i.e. requiring | δn ( r ) | p to be inte-grable). Other norms may induce weaker constraints onthe derivative functions F m but stronger constraints onthe density variations δn ( r ). A rigorous discussion onthe issue of functional differentiability for the case of theenergy functional is given in Refs.[25–28].Let us now look at the symmetry properties of the func-tions F m . As follows directly from the definition (4) ofthe functions F m and the assumption of differentiability,we have the permutational symmetry F m ( r ′ , r , r . . . r m ) = F m ( r ′ , r , r π (1) . . . r π ( m ) )for all permutations π of the integers 1 . . . m . Since thefunctions F m are evaluated at the homogeneous density n they also have the spatial symmetry properties of thehomogeneous electron gas. These symmetries are thetranslational symmetry over an arbitrary vector a , rota-tional symmetry under arbitrary rotations by a rotationmatrix R , and inversion symmetry. If we denote by a = ( a , . . . , a ) R r m = ( R r , . . . , R r m )the m -tuples of translation vectors and rotated positionvectors we can write F m ( r ′ , r , r m ) = F m ( r ′ + a , r + a , r m + a ) (5) F m ( r ′ , r , r m ) = F m ( R r ′ , R r , R r m ) F m ( r ′ , r , r m ) = F m ( − r ′ , − r , − r m ) . Since in Eq.(5) the vector a is arbitrary we can in par-ticular choose a = − r and define a new function N m depending on m + 1 independent vectors by the relation F m ( r ′ , r , r . . . r m )= F m ( r ′ − r , , r − r , . . . , r m − r )= N m ( r ′ − r , r − r , . . . , r m − r ) . (6)The difference vector r ′ − r will appear several times inour equations and it will therefore be convenient to definethe short notation y = r ′ − r . We will further use y = | y | to denote the length of this vector. Our interest willbe in the functions N m and in particular their Fouriertransforms N m ( y , q m ) = Z d r m N m ( y , r m ) e − i q · r − ... − i q m · r m (7)with Fourier inverse N m ( n ; y , r m ) = Z d q m (2 π ) m N m ( y , q m ) e i q · r + ... + i q m · r m . With these definitions we obtain the following expansionfor F in powers of the density variations F [ n ]( r ′ , r ) = F ( n , y ) + ∞ X m =1 m ! Z d r m Z d q m (2 π ) m N m ( y , q m ) e i q · ( r − r )+ ... + i q m · ( r m − r ) δn ( r m )= F ( n , y ) + ∞ X m =1 m ! Z d q m (2 π ) m N m ( y , q m ) e − i ( q + ... + q m ) · r δn ( − q m ) (8)This equation forms the basis of the gradient expansion.To proceed we further derive a consistency condition thatis 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 ) = Z d r ′′ F m +1 ( r ′ , r , r ′′ , r m ) δn ( r ′′ ) . Taking δn ( r ′′ ) = δn then to be a uniform shift in thedensity (which means that we look at a compression ordecompression of the electron gas) yields ∂F m ∂n ( r ′ , r , r m ) = Z d r ′′ 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 Fouriertransform. This gives the condition ∂N m ∂n ( y , q . . . q m ) = N m +1 ( y , , q . . . q m ) . (9)This relation is of key importance for the existence of thegradient expansion. Without it we would not be able toeliminate the dependence of the gradient coefficients onthe homogeneous density n in favor of a dependence onthe actual density profile. As we will see, this requiresa summation over response functions F m to infinite or-der. There is another important advantage of this re-summation. It will allow us to relax the constraint thatthe density variations be small (but varying arbitrarilyfast) and replace it with the constraint that the densityvariations 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 toassume that the density gradients are small. This con-straint is most easily phrased in momentum space by re-quiring that the Fourier coefficients δn ( q ) of the densityvariation have their main contribution from small wave vectors q . If this is the case then the main contribu-tion to the integrals in Eq.(8) comes from this region ofsmall vectors and we can expand the functions N m inpowers of the wave vectors. Subsequently carrying outthe integrals in Eq.(8) then leads to an expansion in den-sity gradients, since powers in wave vectors correspondto orders of derivatives in real space. Since the responsefunctions N m typically vary very rapidly for wave vec-tors close to the Fermi surface (or multiples thereof) werequire that the Fourier coefficients δn ( q ) of the densityvariation have their main contributions from wave vec-tors that satisfy | q | < k F where k F is the Fermi wavevector. Note that the procedure of interchanging inte-gration and summation requires absolute convergence ofthe wave vector expansion. It is hard to prove this prop-erty for the general response functions that we considerhere and we therefore have to assume its validity.The next task is then to expand the functions N m inpowers of wave vectors. This task is simplified when weexploit the symmetry of the these functions. As followsdirectly from Eqs.(6) and (7) the functions N m in Fourierspace inherit the symmetry properties of F m . They aresymmetric functions of their wave vectors, i.e. N m ( y , q . . . q m ) = N m ( y , q π (1) . . . q π ( m ) )for any permutation π of the integer labels, and are in-variant under rotations and inversion N m ( y , q . . . q m ) = N m ( R y , R q . . . R q m ) N m ( y , q . . . q m ) = N m ( − y , − q . . . − q m ) . From these equations we therefore see that any powerexpansion of the functions N m must lead to an expansionin terms of the spatial vector y and the wave vectors q j which is invariant under rotations and inversions ofthese vectors. The mathematical question is then whichpolynomials have these properties. This was answeredin the classic book by Weyl [29]; every such polynomialmust be a function of the variables y · y , y · q i and q i · q j .We see that these inner products are indeed invariantunder rotations and inversions. Since any polynomial inpowers of y = 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 N m ( y , q . . . q m ) = N m ( y ) + N m ( y ) y · m X i =1 q i + N m ( y ) m X i =1 q i + N m ( y ) m X i =1 ( y · q i ) + N m ( y ) m X i>j ( q i · q j ) + N m ( y ) m X i>j ( y · q i )( y · q j ) + . . . (10)where q i = q i · q i and where we took into account thatthis expansion must be invariant under permutations ofthe wave vectors q j . This expansion is readily extendedto higher order powers in the wave vectors. For a givenchoice of correlation function F the practical task will beto determine the explicit form of the coefficient functions N mj ( y ) as a function of y and the background density n .If we insert Eq.(10) into Eq.(8) and Fourier transformback to real space we obtain the expansion F ( r ′ , r ) = F ( n , y ) + ∞ X m =1 m ! ∞ X j =0 N mj ( y ) A mj ( y , r ) . (11)Let us analyze the explicit form of the first six coefficients A mj for j = 0 , .., j = 0 is simply given by A m = δn ( r ) m . The terms with j = 1 , , m symmetry equivalentterms and acquire the following forms in real space A m = i m ( δn ( r )) m − y · ∇ δn ( r ) A m = − m ( δn ( r )) m − ∇ δn ( r ) A m = − m ( δn ( r )) m − y ( y y · ∇ ) δn ( r )In the derivation of the last term we used that for anarbitrary function f ( y y · ∇ ) f = X ij y i y ∂ i ( y j y ∂ j ) f ( r ) = X ij y i y j y ∂ i ∂ j f where with r = ( x , x , x ) we used the notation ∂ i = ∂/∂x i as well as y i = x ′ i − x i . This expression is readilychecked using ∂ i ( y j y ) = − δ ij y + y i y j y . Finally the terms with j = 4 , m ( m − / A m = − m ( m − δn ( r )) m − ( ∇ δn ( r )) A m = − m ( m − δn ( r )) m − ( y · ∇ δn ( r )) . We see that a general coefficient function A mj dependson the density variation and gradients thereof. Thefunctions δn ( r ) that appear under the gradient opera-tors can be replaced by the actual density profile n ( r ) = n + δn ( r ). However, we see that we are still left withpowers of δn ( r ) which are unprotected by derivative op-erators and which therefore can not be replaced by thefull density profile n ( r ). It is precisely at this point thatthe consistency conditions (9) play a crucial role. If weuse these conditions in Eq.(10) then we see that N m +1 j ( y ) = ∂N mj ∂n ( y ) . (12)This equation relates certain coefficients coming fromhigher order response functions N m to those of lowerorder ones. In particular it tells us that by iteration N m ( y ) = ∂ m F ∂n m ( y ) N mj ( y ) = ∂ m − N j ∂n m − ( y ) ( j = 1 , , N mj ( y ) = ∂ m − N j ∂n m − ( y ) ( j = 4 , . If we insert these expressions together with the explicitform of the coefficient functions A mj back into Eq.(11) inEq.(11) and shift some indices we obtain the expansion F ( r ′ , r ) = F ( y ) + ∞ X m =1 m ! ∂ m F ∂n m ( y ) δn ( r ) m + ∞ X m =0 m ! δn ( r ) m (cid:2) i ∂ m N ∂n m ( y · ∇ n ( r )) − ∂ m N ∂n m ∇ n ( r ) − ∂ m N ∂n m y ( y y · ∇ ) n ( r ) (cid:3) + − ∞ X m =0 m ! δn ( r ) m (cid:2) ∂ m N ∂n m ( ∇ n ( r )) + ∂ m N ∂n m ( y · ∇ n ( r )) (cid:3) + . . . (13)In this expression we recognize several Taylor series of the coefficient functions N mj ( n + δn ( r ) , y ) in powers of δn ( r ).These Taylor series can be re-summed to finally give F ( r ′ , r ) = F ( n ( r ) , y ) + iN ( n ( r ) , y ) y · ∇ n ( r ) − N ( n ( r ) , y ) ∇ n ( r ) − N ( n ( r ) , y ) y ( y y · ∇ ) n ( r ) − N ( n ( r ) , y ) ( ∇ n ( r )) − N ( n ( r ) , y ) ( y · ∇ n ( r )) + . . . , (14)where the implicit dependence on n of the coefficient functions is now replaced by a dependence on the full densityprofile. We see that this is achieved by summing over response functions N m to infinite order. For clarity we reinstatedthe 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 in favor of adependence on the full density profile. Having obtained the general form of gradient expansion we can wonder aboutits convergence. In order for the gradient expansion to be useful it at least needs to be an asymptotic series. Thisis known to be the case for the gradient expansion of some commonly studied functionals. For instance, for smalldensities variations very accurate results for the exchange energy are obtained by summing all terms up to fourthorder in the density derivatives [4, 23].Having obtained the general gradient expansion Eq.(14) it remains to find explicit expressions for the coefficientfunctions N mj . We will describe how to do this in detail for the coefficients needed for a gradient expansion tosecond order in the density derivatives. To calculate the coefficients N , N and N we need to calculate the function N ( n , y , q ) for the homogeneous electron gas and expand it in powers of q : N ( n , y , q ) = N ( n , y ) + N ( n , y ) y · q + N ( n , y ) q + N ( n , y )( y · q ) + . . . (15)The determination of the coefficients N and N requires knowledge of the second order response function N ( n , y , q , q ) = N ( n , y ) + N ( n , y ) y · ( q + q ) + N ( n , y )( q + q ) + N ( n , y )(( y · q ) + ( y · q ) ) + N ( n , y )( q · q ) + N ( n , y )( y · q )( y · q ) + . . . (16)The expression (14) together with Eqs.(15) and (16) de-termine completely the gradient expansion of an arbi-trary correlation function F to second order in the den-sity gradients. These equations form the main result ofthe present work. If expansions to higher order gradientsare required they can be derived without difficulty alongthe lines described above. To do this one needs to con-struct higher order symmetric polynomials in the wavevectors 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 gra-dient coefficients depending on n ( r ). What is needed todetermine the explicit form of the gradient coefficients inpractice is the determination of the functions N m . Howto do this for N and N is described in the next section. C. Determination of N and N A practical calculation of the coefficients of the gradi-ent expansion of Eq.(14) requires the knowledge of thethe functions N and N . According to Eqs.(4) and (6)these are defined as the first and second density deriva-tives of the correlation function F with respect to thedensity, i.e. δF ( r ′ , r ) δn ( r ′′ ) | n = N ( y , r ′′ − r ) δF ( r ′ , r ) δn ( r ′′ ) δn ( r ′′′ ) | n = N ( y , r ′′ − r , r ′′′ − r ) . To derive useful expressions for these functions it is con-venient to transform derivatives with respect to the den-sity to derivatives with respect to the external potentialas such derivatives appear naturally in perturbation the-ory. We can then use the well-developed tools of pertur-bation theory to calculate these functions. Let us there-fore define the new functions F ( y , r ′′ − r ) = δF ( r ′ , r ) δv ( r ′′ ) | n (17) F ( y , r ′′ − r , r ′′′ − r ) = δ F ( r ′ , r ) δv ( r ′′ ) δv ( r ′′′ ) | n . (18)Since we evaluate these functions for the homogeneouselectron gas they only depend on differences between thespatial coordinates. It will therefore be convenient to alsodefine their Fourier transforms by F ( r , r ) = Z d p d q (2 π ) F ( p , q ) e i p · r + i q · r F ( r , r , r ) = Z d k d p d q (2 π ) F ( k , p , q ) × e i k · r + i p · r + i q · r as well as their partial Fourier transforms F ( y , q ) = Z d p (2 π ) F ( p , q ) e i p · y F ( y , p , q ) = Z d k (2 π ) F ( k , p , q ) e i k · y . The chain rule of differentiation gives a connection be-tween the derivatives with respect to the density and thederivatives with respect to the potential of the correlationfunction F . We have δF ( r , r ) δv ( r ) = Z d r δF ( r , r ) δn ( r ) δn ( r ) δv ( r ) (19) δ F ( r , r ) δv ( r ) δv ( r ) = Z d r δF ( r , r ) δn ( r ) δ n ( r ) δv ( r ) δv ( r )+ Z d r d r δ F ( r , r ) δn ( r ) δn ( r ) δn ( r ) δv ( r ) δn ( r ) δv ( r ) . (20)We see that in these expressions the first and second orderderivative of the density with respect to the potentialappear. We therefore define the linear density responsefunction by χδn ( r ) δv ( r ) | n = χ ( r − r ) = Z d q (2 π ) χ ( q ) e i q · ( r − r ) as well as the second order density response function χ by δ n ( r ) δv ( r ) δv ( r ) | n = χ ( r − r , r − r )= Z d p d q (2 π ) χ ( p , q ) e i p · ( r − r )+ i q · ( r − r ) . Using these definitions we find by Fourier transformingEqs.(19) and (20) that N ( y , q ) = F ( y , q ) χ ( q ) (21) N ( y , p , q ) = (cid:2) F ( y , p , q ) − N ( y , p + q ) χ ( p , q ) (cid:3) χ ( p ) χ ( q ) (22)These equations give the desired relation between thedensity derivatives N and N and the potential deriva-tives F and F of the correlation function F . Relationsbetween the higher order derivatives N m and F m can bederived in a completely analogous way. From Eqs.(21)and (22) we see that to calculate the functions N and N , and hence the gradient expansion coefficients of F to second order in the gradients, we need to calculatethe density response and the change in F to secondorder in a perturbing potential δv . The problem isthereby reduced to doing perturbation theory on theelectron 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 standardperturbation theory to obtain these response functionfor the case that F is the one-particle density matrixand use that to calculate the gradient expansion of theexchange hole. III. GRADIENT EXPANSION OF THEONE-PARTICLE DENSITY MATRIX AND THEEXCHANGE HOLEA. The one-particle density matrix
As an application of our formalism we carry out thegradient expansion of the one-particle density matrix fora noninteracting system with density n ( r ). This prob-lem has received large interest in the past since both theexchange energy E x [ n ] as well the Kohn-Sham kineticenergy T s [ n ] are explicit functionals of such a noninter-acting density matrix [2, 3]. Therefore a gradient expan-sion of this density matrix directly leads to a gradientexpansion of these two functionals. Such a gradient ex-pansion was first carried out by Gross and Dreizler [31]using the Kirzhnits expansion [2, 32]. This expansion isadjusted to the specific form of the noninteracting one-particle density matrix and can not be generalized to ar-bitrary correlation functions. The method presented inthis work can, on the other hand, be applied to arbitrarycorrelation functions, but to demonstrate the soundnessof our formalism we will use it to provide an alternativederivation of the result obtained earlier by Gross andDreizler using the Kirzhnits expansion. An advantage ofour derivation based on nonlinear response theory is thatit avoids the turning point problem encountered in theKirzhnits approach [2].Let us start by defining the one-particle density matrixin second quantization as γ ( r ′ , r ) = X σ h Ψ | ˆ ψ † ( r σ ) ˆ ψ ( r ′ σ ) | Ψ i where σ is a spin coordinate. We will consider the caseof spin-compensated systems. Then, in a noninteractingsystem with 2 N electrons, the density matrix is given interms of one-particle wave functions ϕ j ( r ) by γ ( r ′ , r ) = N X j =1 ϕ j ( r ′ ) ϕ ∗ j ( r ) (23)where the pre-factor 2 is a spin factor. The one-particlestates are eigenstates of a one-particle Schr¨odinger equa-tion (cid:18) − ∇ + v ( r ) (cid:19) ϕ j ( r ) = ǫ j ϕ j ( r ) (24)where v ( r ) is the external potential. To determine thegradient expansion coefficients of γ we need to determinehow the density matrix changes if we make small changes v ( r ) → v ( r ) + δv ( r ) in the external potential. More pre-cisely, we need to calculate the functional derivatives γ ( r , r , r ) = δγ ( r , r ) δv ( r ) γ ( r , r , r , r ) = δ γ ( r , r ) δv ( r ) δv ( r )which are the equivalents of the functions F and F ofEqs.(17) and (18) for the specific choice of correlationfunction F = γ (note that they become functions of twoand three independent vectors respectively when evalu-ated for a homogeneous system). To do this it is sufficientto know the functional derivatives of the orbitals ϕ j andeigenvalues ǫ j with respect to the potential v . These aresimply given by doing first order perturbation theory onthe one-particle Schr¨odinger Eq.(24). We find δǫ j δv ( r ) = | ϕ j ( r ) | (25) δϕ j ( r ) δv ( r ′ ) = ∞ X k = j ϕ k ( r ) ϕ ∗ k ( r ′ ) ϕ j ( r ′ ) ǫ j − ǫ k . (26) With these relations we can differentiate Eq.(23) twicewith respect to the potential. Afterwards we then needto insert the plane wave states ϕ k ( r ) and their eigenval-ues ǫ k appropriate for the homogeneous electron gas intothe final expressions. In order not to interrupt the pre-sentation these calculations are given in Appendix A andwe simply present the results of these calculations here.If we define the Fermi factors by n k = θ ( k F − | k | ) , where θ is the Heaviside function and where k F =(3 π n ) / is the Fermi wave vector, then the resultingexpressions are given by γ ( y , q ) = 2 Z d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · y (27) γ ( y , p , q ) = 2 Z d k (2 π ) e i k · y [Φ( k , k + p + q , k + p )+ Φ( k , k + p + q , k + q )] (28)whereΦ( k , p , q ) = 1 ǫ q − ǫ p (cid:18) n q − n k ǫ q − ǫ k − n p − n k ǫ p − ǫ k (cid:19) . (29)In this expression ǫ k = | k | / χ ( q ) and χ ( p , q ) to evaluate thefunctions N and N of Eqs. (21) and (22). Fortunately,for the specific choice F = γ that we made these densityresponse functions are simply given by χ ( q ) = γ ( y = 0 , q ) (30) χ ( p , q ) = γ ( y = 0 , p , q ) . (31)They are therefore obtained directly from Eqs. (27) and(28) so that we do not need to do any extra work tocalculate them.To calculate the gradient coefficient to second order inthe gradients we now need to expand γ and γ to secondorder in the wave vectors p and q . Since these calcula-tions are somewhat lengthy we do not present them hereand refer to Appendix B for a description of the mainsteps. The resulting expansions are given by γ ( y , q ) = − k F π (cid:18) sin zz (cid:20) − i q · y − ( q · y ) (cid:21) − q k cos( z ) (cid:19) (32) γ ( y , p , q ) = 1 π k F (cid:18) cos( z ) − i p + q ) · y cos( z ) + 112 k ( p + q + p · q )(cos( z ) + z sin( z )) −
124 [( p · y ) + ( q · y ) + 3(( p + q ) · y ) ] cos( z ) (cid:19) (33)where we defined z = k F y . The expansions for the first and second order response functions χ and χ of Eqs.(30)and (31) are obtained by taking y = 0 in these expressions. Together which Eqs.(32) and (33) these expression thenimmediately determine the functions N and N from Eqs.(21) and (22). The final result of this calculation to secondorder in the wave vectors is given by N ( y , q ) = sin zz (cid:20) − i q · y − ( q · y ) (cid:21) + q k sin( z ) − z cos( z ) zN ( y , p , q ) = π k (cid:20)(cid:18) cos( z ) − sin( z ) z (cid:19) (cid:18) − i p + q ) · y −
16 [( p · y ) + ( q · y ) ] (cid:19) + 112 k ( p + q + p · q )(3 cos( z ) + z sin( z ) − z ) z ) + 112 ( p · y )( q · y )[4 sin( z ) z − z )] (cid:21) By comparison of these two expression to the expansions (15) and (16) we can directly read off the gradient coefficientfunctions of Eq.(14) to be N ( y ) = − i j ( z ) N ( y ) = 112 k z j ( z ) N ( y ) = − j ( z ) N ( y ) = π k [ z j ( z ) − zj ( z )] N ( y ) = π k [ j ( z ) + 3 zj ( z )]where now we replaced n → n ( r ) and hence the Fermi wave vector k F = (3 π n ( r )) / that appears in these expressionsis a spatially varying function. We further introduced the spherical Bessel functions j ( z ) = sin zz j ( z ) = sin z − z cos zz . By inserting the gradient coefficient functions into Eq.(14) we find the following explicit gradient expansion for theone-particle density matrix γ ( r ′ , r ) = k π j ( z ) z + 12 j ( z ) y · ∇ n ( r ) − k z j ( z ) ∇ n ( r ) + 16 k z j ( z )( y y · ∇ ) n ( r ) − π k ( z j ( z ) − zj ( z )) ( ∇ n ( r )) − π k ( j ( z ) + 3 zj ( z )) ( y · ∇ n ( r )) + . . . (34)This expression is the main result of this section. After some manipulations we can rewrite it in an equivalent formas γ ( r ′ , r ) = 1 π h k j ( z ) z − ∇ k k F zj ( z ) + 112 1 k F (cid:16) ∇ · ( ∇ k · y y ) (cid:17) · y y z j ( z ) + 14 ∇ k · y y zj ( z ) −
196 ( ∇ k ) k ( j ( z ) z − zj ( z )) + 132 1 k (cid:16) ∇ k · y y (cid:17) ( z j ( z ) − z j ( z )) i (35)This expression becomes identical in form to that derived by Gross and Dreizler [31] using the Kirzhnits methodafter eliminating the effective Fermi vector of their work in favor of the density (this requires inversion of Eq.(19) ofreference [31] ). We close this section with a final remark on our result. We note that the gradient expansion to finiteorder breaks the symmetry γ ( r , r ′ ) = γ ( r ′ , r ) ∗ . However, the full gradient expansion as described by Eq.(11) has thissymmetry, which means that the symmetry is restored by taking all higher order gradients into account whenever theseries converges. The symmetry violation has clearly consequences for the quantities that are derived from it, such asthe exchange hole, as it introduces ambiguity in defining such functions. To be in accordance with existing literaturewe will in the next section adopt the definition of the exchange hole that was used by Perdew [9].0 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 andpresent an alternative derivation of the gradient expansion for the exchange energy compared to that of Dreizler andGross which has the advantage of avoiding the regularization of divergent integrals [31]. In doing so we confirm aresult obtained by Langreth and Perdew [19]. The exchange hole can be calculated directly from the one-particledensity matrix as ρ x ( r + y | r ) = − | γ ( r + y , r ) | n ( r ) . Inserting the expansion (34) of for the one-particle density matrix and retaining terms to second order in the gradientsyields the expression ρ x ( r + y | r ) = − n ( r ) j ( z ) z − j ( z ) j ( z ) z ( y · ∇ n ( r )) + 14 k j ( z ) ∇ n ( r ) − k zj ( z ) j ( z )( y y · ∇ ) n ( r )+ π k j ( z )( zj ( z ) − j ( z ))( ∇ n ( r )) + π k ( j ( z ) j ( z ) z + 3 j ( z ) − j ( z ) )( y · ∇ n ( r )) + . . . (36)As was pointed out by Perdew [9] the second order gradient expansion of the exchange hole does not satisfy thenegativity and sum rule constraints ρ x ( r + y | r ) ≤ Z d y ρ x ( r + y | r ) = − . The violations of the negativity constraint is readily verified. However, the violation of the sum rule is not immediatelyobvious from Eq.(36). The sum rule can, however, be conveniently analyzed in momentum space. To do this we definethe Fourier transformed exchange hole to be ρ x ( k | r ) = Z d y ρ x ( r + y | r ) e − i k · y (37)This function has the form ρ x ( k | r ) = ρ ( k | r ) + α ( k, n ) ˆ k · ∇ n ( r ) + α ( k, n ) ∇ n ( r ) + α ( k, n ) (ˆ k · ∇ n ( r )) + α ( k, n )( ∇ n ( r )) + α ( k, n ) (ˆ k · ∇ ) n ( r ) + . . . (38)where we defined the unit vector ˆ k = k /k and k = | k | . The explicit form of the functions ρ and α j follow directly byFourier 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 conditionin momentum space translates to the requirement that ρ x ( k = 0 | r ) = −
1. Since ρ ( k = 0 | r ) = − α j ( k = 0 , n ) = 0. However, as is clear from Eqs.(C2)-(C6)in the Appendix this is not satisfied for α , whereas α and α even diverge for k → N h ρ x ( k ) i = Z d r n ( r ) ρ x ( k | r )= N h ρ ( k ) i + ˆ k · Z d r n ( r ) α ( k, n ( r )) ∇ n ( r ) + X i,j =1 Z d r n ( r ) β ij ( k, n ( r )) ∂ i n ( r ) ∂ j n ( r ) . (39)where N is the number of electrons. In this expression h ρ ( k ) i is the system average of ρ ( k | r ) and in the last term under the integral sign we defined the tensor β ij ( n, k ) = α L ( n, k ) k i k j k + α T ( n, k ) (cid:16) δ ij − k i k j k (cid:17) . α L and a transverse part with coefficient α T which describethe contributions to the system-averaged hole of densitygradients ∇ n parallel and perpendicular to the momen-tum vector k . These coefficients are calculated from thefunctions α j as α T = α − n ∂ ( nα ) ∂nα L = α T + α − n ∂ ( nα ) ∂n , 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 explicitexpressions for α L , T . If we define the dimensionless vari-able ¯ k = k/ (2 k F ) then we have α T = π nk Z T (¯ k ) α L = π nk Z L (¯ k ) , where the functions Z L , T have the explicit form Z T ( x ) = − x θ (1 − x ) + 43 δ ( x − Z L ( x ) = − x θ (1 − x ) + δ ( x −
1) + 13 δ ′ ( x − k → Z L , T (¯ k ) = 0and hence β ij ( k, n ) → k →
0. This implies thatthe last term in Eq.(39) does not contribute to the sumrule of the system averaged hole. The same conclusioncan be derived for the second term in Eq.(39). Since α ( k = 0 , n ) = i/ (4 nk F ) (see Eq.(C2)) is a local func-tion of the density the integrand of the first integral inEq.(39) is a total derivative and vanishes (assuming ei-ther vanishing densities at infinity or periodic boundaryconditions). We therefore find thatlim k → h ρ x ( k ) i = lim k → h ρ ( k ) i = − . This implies that system averaged exchange hole ob-tained from the second order gradient expansion doesfulfill the sum rule, although the exchange hole itself doesnot. We note, however, that this is only achieved by in-tegrating over both positive and negative contributions.When the positive contributions to the exchange hole areremoved the sum rule is only recovered for a finite hole ra-dius [9, 10]. We finally note that with expression Eq.(39) it is now a simple matter to calculate the exchange energyfrom E x [ n ] = N Z d k (2 π ) h ρ x ( k ) i π | k | . We insert Eq. (39) and do the angular integrations first.It is therefore convenient to define the spherical and sys-tem averaged hole in momentum space by hh ρ x ( k ) ii = 14 π Z d Ω k h ρ x ( k ) i such that E x [ n ] = Nπ Z ∞ dk hh ρ x ( k ) ii . (40)From Eq.(3) we then find that N hh ρ x ( k ) ii = N hh ρ ( k ) ii + Z d r n ( r )[ 13 α L ( k, n ) + 23 α T ( k, n )]( ∇ n ) = N hh ρ ( k ) ii + Z d r π k Z x (¯ k )( ∇ n ) , (41)where we defined Z x (¯ k ) = 13 Z L (¯ k ) + 23 Z T (¯ k ) . (42)The exchange energy is obtained by inserting the expres-sion for the averaged hole of Eq.(41) into Eq.(40). Thisthen finally gives the expression E x [ n ] = − (cid:16) π (cid:17) Z d r n ( r )+ Z d r B x ( n ( r ))( ∇ n ( r )) (43)where B x ( n ) = π k Z ∞ d ¯ k Z x (¯ k )= − π (3 π ) / n / and where to calculate the first term in Eq.(43) we usedEq.(C1). The coefficient B x is the same gradient coeffi-cient as obtained by Gross and Dreizler [31] and earlierby Sham [34]. However, the correct analytic exchangegradient coefficient is known [35–39] to be a factor 10/7larger. The reason for this discrepancy is clearly de-scribed by Svendsen and von Barth [38] who showed thatthe Coulomb interaction is too singular to allow for theinterchange of the operations of integration and the ex-pansion in wave vectors. The problem does, for instance,not appear when Yukawa-screened Coulomb interactionsare used [38]. The conclusion is that there is nothingwrong with the Kirzhnits method, or the nonlinear re-sponse theory derivation of the gradient expansion of the2exchange hole presented here, but one has be aware thatfor Coulomb interactions the procedures of expandingcorrelation functions in terms of density gradients fol-lowed by integrations involving Coulomb potentials maynot yield the same result as directly expanding the inte-grated quantity in terms of density gradients. For thisreason the original GGA of Perdew and Wang [10] basedon the gradient expansion of the exchange hole was laterreparametrized [16, 40] to accomodate the correct gradi-ent coefficient for the exchange energy. IV. CONCLUSIONS AND OUTLOOK
We derived a general scheme based on nonlinear re-sponse theory to calculate the density gradient expan-sion of general correlation functions and showed that inorder to express the gradient coefficients in terms of thefull density profile summations to infinite order must becarried out over response functions of arbitrarily large or-der. A consistency condition was derived to do this. Theformalism was used to derive the gradient expansion ofthe one-particle density matrix and the exchange hole tosecond order in the gradients. We confirm the derivationof Dreizler and Gross that used the Kirzhnits expansionmethod. We further analyzed the exchange hole in mo-mentum space to derive that the system averaged hole tosecond order in the gradients satisfies the sum rule andto derive the gradient expansion of the exchange energywithout the need to regularize divergent integrals.The scheme that we presented is very general and canbe applied to more general correlation functions. A nat-ural application of the scheme would be to calculate thegradient expansion of the correlation hole ρ c . Regardingthe gradient terms of the correlation hole little is known.We essentially only have some detailed information onthe long-range properties of the spherically and systemaveraged hole. This information comes from the work ofLangreth and Perdew who calculated hh ρ c ( k ) ii within theRPA. In our notation their results (see also Appendix Cof [19]) can be summarized as N hh ρ c ( k ) ii = N hh ρ ( k ) ii + Z d r π k z c ( n, k )( ∇ n ) where the function z c has been parametrized by Langrethand Mehl [20, 21] as z c ( n, k ) = 4 √ k s e − √ kk s (44)where k s = (4 k F /π ) / is the Thomas-Fermi or screeningwave vector. We can transform to real space to obtainthe following expression for the spherically and systemaveraged correlation hole N hh ρ c ( y ) ii = N hh ρ ( y ) ii + Z d r k k (1 + ( k s y ) / ( ∇ n ) (45) We see from Eq.(44) that z c ( n, k = 0) = 0 and as a con-sequence the sum rule hh ρ c (0) ii = 0 for the correlationhole is not satisfied. We know, however, that the RPAbecomes accurate in the high density limit and since fromEq.(45) the gradient coefficient of the correlation hole de-pends on k s y we see that high densities are equivalent tolarge separations y . Similarly, low density correspondsto short-distance behavior. However, RPA can not betrusted in this region and we have no precise knowledgeon the small distance behavior of the gradient coefficientof hh ρ c ( y ) ii . A model for the general short-range behav-ior was proposed [16–18] on physical arguments (it shouldnot affect the Coulomb cusp of and the on-top value ofthe LDA hole) after which the real-space cutoff proce-dure was applied (see e.g. Fig. 4 in [17] and Fig. 5 in[18]) to obtain a GGA for correlation. It would, how-ever, be desirable to have a first principles approach tothe calculation of the short range properties of the gra-dient coefficient of the correlation hole. As follows fromour derivation this requires the knowledge of the func-tions (17) and (18), or at least their expansion to secondorder in wave vectors, for the case that F represents thepair correlation function or the correlation hole of theelectron gas. It may therefore be worthwhile to use ourcurrent scheme to explore these response functions be-yond the RPA. For short range correlation an approachbased on diagrammatic summation of ladder diagramssuggests itself. Of course, for the development of densityfunctionals for general systems beyond the weakly inho-mogeneous regime it is not sufficient to use the straight-forward gradient expansion [23]. However, general short-range features, such as the way the exchange-correlationhole deforms close to the reference electron in the pres-ence of a density gradient can be transferred to moregeneral systems than the weakly inhomogeneous ones.Perhaps in combination with truly nonlocal functionalsfor the exchange-correlation hole [41] this can lead to thedevelopment of more accurate density functionals. Thiswill be a topic of our future research.Finally, we would like to make some general remarkson the extension of our derivations to temperature- ortime-dependent systems. In the case of temperature-dependent systems the expectation values of observablesare traces over a grand canonical ensemble. By Mermin’stheorem [42] such expectation values are still functionalsof the density and therefore the derivations in this workcan be carried out unchanged. In fact, the dependenceon temperature is probably going to improve the conver-gence properties of the gradient series by smootheningof the sharp Fermi functions which, for instance, werethe cause of the oscillatory nature of the gradient coeffi-cients of the one-particle density matrix. The situationof time-dependent systems is more delicate. Althoughthe existence of a density-potential mapping has beenwell-established [43–46] and hence a density functionaltheory can be formulated, the time-variable induces se-vere complications. As was shown by Vignale and Kohn[47, 48] the temporal and spatial non-locality of time-3dependent density functionals are intimately correlated.Any temporal non-locality (or frequency dependence inthe language of equilibrium response functions) induceslong-range spatial properties that prevents the gradientexpansion from existing. A way out of this problem isto use new variables such as the current-density [47, 48]or the local density deformation density in a Lagrangianframe [49–51] for which a local density approximationand a corresponding gradient expansion does exist. Thishas been exploited in Ref.[52] to construct a GGA func-tional within time-dependent current density functionaltheory. The study of correlation functions in the samefashion remains a challenge for the future. Acknowledgments
I would like to thank the students in the course ”Den-sity Functional Theory” at the University of Jyv¨askyl¨afor inspiration. I would further like to thank Esa R¨as¨anenfor useful discussions and Klaas Giesbertz for checking alarge part of the equations .
Appendix A: Calculation of γ and γ In this Appendix we will derive the expression for thefirst and second derivative of the one-particle density ma-trix γ with respect to the potential v . For a clearer inter-pretation and compactness of notation we will only puta comma between the variables representing the originalcoordinates of the density matrix and the variables thatappear as argument of the potential variations. Directdifferentiation of Eq.(23) using Eq.(26) then gives γ (23 ,
1) = δγ (23) δv (1)= 2 ∞ X i,j ( f j − f i ) ϕ j (1) ϕ ∗ i (1) ϕ i (2) ϕ ∗ j (3) ǫ j − ǫ i (A1)where we used the short notation j = r j and introducedthe occupation numbers f j = 1 for an occupied state and f j = 0 for an unoccupied state. Inserting plane wavestates appropriate for the electron gas we find that γ ( r r , r ) = 2 Z d q d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · ( r − r )+ i q · ( r − r ) . We therefore find that the function γ has the simpleform γ ( p , q ) = 2 n p + q − n p ǫ p + q − ǫ p in Fourier space. This is precisely the integrand ofEq.(27). The function γ can be calculated by straight-forward differentiation of Eq.(A1) with respect to the po-tential γ (23 ,
14) = δγ (23 , δv (4) . To obtain the explicit form of this function we have to dif-ferentiate both the orbitals and the eigenvalues using Eqs.(25) and (26). Let us denote the term that arises fromthe differentiating the orbitals by X and the term thatarises from differentiating the eigenvalues by Y . Then wefind that γ (23 ,
14) = X (23 ,
14) + Y (23 , , where X (23 ,
14) = 2 ∞ X i,j,k ϕ i (2) ϕ ∗ j (3) h Φ( ijk ) ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ Φ( jik ) ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i (A2) Y (23 ,
14) = − ∞ X i,j ϕ i (2) ϕ ∗ j (3)Φ( iji )( | ϕ j (4) | − | ϕ i (4) | ) ϕ j (1) ϕ ∗ i (1) (A3)and where we further definedΦ( ijk ) = 1 ǫ k − ǫ j (cid:18) f k − f i ǫ k − ǫ i − f j − f i ǫ j − ǫ i (cid:19) . (A4)The function Φ has the useful properties Φ( ijj ) = 0 andΦ( ijk ) = Φ( ikj ). The function γ (23 ,
14) = γ (23 ,
41) issymmetric in the indices 4 and 1 due to the fact that thedifferentiations 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 γ in such a way that this symmetry is obvious. To do thiswe first note that since Φ( iji ) = − Φ( jij ) the terms with k = i and k = j in Eq.(A2) sum up to a function Z of asimilar form as Y in Eq.(A3), namely Z (23 ,
14) = − ∞ X i,j ϕ i (2) ϕ ∗ j (3)Φ( iji )( | ϕ j (1) | − | ϕ i (1) | ) ϕ j (4) ϕ ∗ i (4) . (A5)We can therefore write γ (23 ,
14) as γ (23 ,
14) = 2 ∞ X ijk,k =( i,j ) ϕ i (2) ϕ ∗ j (3) h Φ( ijk ) ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ Φ( jik ) ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i + Y (23 ,
14) + Z (23 ,
14) (A6)It is easily seen that the sum of Y and Z is symmetricunder interchange of 1 and 4. However, this is not obvi-ous in the first term of the equation above since at firstsight Φ( jik ) does not appear to be equal to Φ( ijk ) sinceΦ( jik ) = 1 ǫ k − ǫ i (cid:18) f k − f j ǫ k − ǫ j − f i − f j ǫ i − ǫ j (cid:19) (A7)4which seems to be different from expression (A4). How-ever, this is just appearance. The reason is that the oc-cupations can only attain the values zero and one. Forthe case f i = f j = 1 or f i = f j = 0 we see directly thatEq.(A4) and (A7) attain the same value. A little inspec-tion shows that this is also true for cases f i = 0 , f j = 1and f i = 1 , f j = 0. We therefore find for k = i, j thatΦ( ijk ) = Φ( jik ). Therefore Eq.(A6) can be simplified to γ (23 ,
14) = 2 ∞ X ijk,k =( i,j ) ϕ i (2) ϕ ∗ j (3)Φ( ijk ) h ϕ k (1) ϕ ∗ i (1) ϕ ∗ k (4) ϕ j (4)+ ϕ j (1) ϕ ∗ k (1) ϕ ∗ i (4) ϕ k (4) i + Y (23 ,
14) + Z (23 ,
14) (A8)This expression is now explicitly symmetric in the indices1 and 4. Let us now evaluate this expression for thehomogeneous electron gas. In the electron gas the one-particle states are plane waves ϕ k ( r ) = e i k · r √ V where V is the volume of the system. Since | ϕ k | = 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 γ (23 ,
14) = 2 V X k , p , q , q =( k , p ) Φ( k , p , q ) e i ( k · r − p · r ) h e i ( q − k ) · r + i ( p − q ) · r + e i ( p − q ) · r + i ( q − k ) · r i where we definedΦ( k , p , q ) = 1 ǫ q − ǫ p (cid:18) n q − n k ǫ q − ǫ k − n p − n k ǫ p − ǫ k (cid:19) . In this expression n p = θ ( ǫ F − ǫ p ) is the zero-temperatureFermi function and ǫ F = k / q = k , p is a region of measure zero and we can therefore freelyintegrate over all wave vectors. After a few substitutionswe obtain γ (23 ,
14) = Z d k d p d q (2 π ) γ ( k , p , q ) e i k · ( r − r )+ i p · ( r − r )+ i q · ( r − r ) where we defined γ ( k , p , q ) = 2[Φ( k , k + p + q , k + p )+Φ( k , k + p + q , k + q )]This gives the integrand of Eq.(28). Appendix B: Expansion of γ and γ
1. Expansion of γ To determine γ to second order in the wave vectorswe have to expand the function γ ( y , q ) = 2 Z d p (2 π ) n p + q − n p ǫ p + q − ǫ p e i p · y (B1)to second order in powers of q . This can be done byexpanding the integrand in powers of q . If we denote∆ = ǫ p + q − ǫ p = p · q + q / q = | q | , then wecan expand the Fermi function n p + q in a distributionalTaylor series as n p + q = n p + ∆ dndǫ | ǫ p + ∆ d ndǫ | ǫ p + ∆ d ndǫ | ǫ p + . . . Since n p = θ ( ǫ F − ǫ p ) is the Heaviside function we have n p + q = n p − ∆ δ ( ǫ F − ǫ p ) + ∆ δ ′ ( ǫ F − ǫ p ) − ∆ δ ′′ ( ǫ F − ǫ p ) + . . . Inserting this into Eq.(B1) then gives γ ( y , q ) = − Z d p (2 π ) δ ( ǫ F − ǫ p ) e i p · y + Z d p (2 π ) ∆ δ ′ ( ǫ F − ǫ p ) e i p · y − Z d p (2 π ) ∆ δ ′′ ( ǫ F − ǫ p ) e i p · y + . . . (B2)We see that in these integrals several derivatives of thedelta function appear. We now use the general mathe-matical expression δ ( n ) ( y ( x )) = dydx ddx ! n X i δ ( x − x i ) | dydx ( x i ) | where the sum runs over all zeros of the function y ( x ), i.e. y ( x i ) = 0. Using this equation in spherical coordinateswith y ( p ) = ( k − p ) / δ ( ǫ F − ǫ p ) = 1 k F δ ( p − k F ) (B3) δ ′ ( ǫ F − ǫ p ) = − pk F δ ′ ( p − k F ) (B4) δ ′′ ( ǫ F − ǫ p ) = − δ ′ ( p − k F ) k F p + δ ′′ ( p − k F ) k F p (B5) δ ′′′ ( ǫ F − ǫ p ) = − k F h p δ ′′′ ( p − k F ) − p δ ′′ ( p − k F )+ 3 p δ ′ ( p − k F ) i (B6)5where in Eq.(B6) we also evaluated the third derivative ofthe delta function as we will need it later in the expansionof γ . With Eqs.(B3)-(B5) we can readily evaluate thethree integrals in Eq.(B2). Let us denote these integralsby A, B and C . Then we find for the first integral A = − k F Z d p (2 π ) δ ( p − k F ) e i p · y = − k F π sin zz where z = k F y . The second integral gives B = Z d p (2 π ) ( p · q + q / δ ′ ( ǫ F − ǫ p ) e i p · y = ( q / − i q · ∇ y ) Z d p (2 π ) δ ′ ( ǫ F − ǫ p ) e i p · y = ( q / − i q · ∇ y ) 12 π k F cos( k F y )= q π k F cos z + i q · y π k F sin zz where in the second step we used Eq.(B4). It thereforeremains to calculate the last integral C = − Z d p (2 π ) ( p · q + q / δ ′′ ( ǫ F − ǫ p ) e i p · y However, up to order q it is suffices to calculate C = − Z d p (2 π ) ( p · q ) δ ′′ ( ǫ F − ǫ p ) e i p · y = 13 ( q · ∇ y ) Z d p (2 π ) δ ′′ ( ǫ F − ǫ p ) e i p · y = −
13 ( q · ∇ y ) π k (cos( z ) + z sin( z ))= − q π k F cos( z ) + ( q · y ) π k F sin( z ) z where in the second step we used Eq.(B5). Adding theresults γ = A + B + C gives the expression (32).
2. Expansion of γ Analogously to γ we can expand γ ( y , p , q ) of Eq.(28)in powers of p and q . To reduce our effort we note thatwe can write γ ( y , p , q ) = A ( y , p , q ) + A ( y , q , p )where A ( y , p , q ) = 2 Z d k (2 π ) e i k · y Φ( k , k + p + q , k + p ) We therefore only need to expand A ( y , p , q ) in powersof the wave vectors and symmetrize with respect to p and q afterwards. To do this we first need to expand thefunctionΦ( k , k + p + q , k + p ) =1 ǫ k + p − ǫ k + p + q (cid:20) n k + p − n k ǫ k + p − ǫ k − n k + p + q − n k ǫ k + p + q − ǫ k (cid:21) . If we denote∆ = ǫ k + p − ǫ k = p p · k ∆ = ǫ k + p + q − ǫ k = ( p + q ) p + q ) · k , we find by expanding the Fermi functions in a distribu-tional Taylor series thatΦ( k , k + p + q , k + p ) = 12 δ ′ ( ǫ F − ǫ k ) −
16 (∆ + ∆ ) δ ′′ ( ǫ F − ǫ k )+ 124 (∆ + ∆ ∆ + ∆ ) δ ′′′ ( ǫ F − ǫ k ) + . . . With this expression we find that the function A ( y , p , q )can be calculated as A ( y , p , q )= Z d k (2 π ) δ ′ ( ǫ F − ǫ k ) e i k · y − Z d k (2 π ) δ ′′ ( ǫ F − ǫ k )(∆ + ∆ ) e i k · y + 112 Z d k (2 π ) δ ′′′ ( ǫ F − ǫ k )(∆ + ∆ ∆ + ∆ ) e i k · y + . . . The three integrals appearing on the left hand side of thisexpression are sufficient to extract the powers to secondorder in p and q . Let us call these three integrals A , A and A in order of appearance. The first integral gives A = Z d k (2 π ) δ ′ ( ǫ F − ǫ k ) e i k · y = cos z π k F . The second integral is given by6 A = − Z d k (2 π ) δ ′′ ( ǫ F − ǫ k ) (cid:20) p + ( p + q ) p + q ) · k (cid:21) e i k · y = (cid:20) − p + ( p + q ) i p + q ) · ∇ y (cid:21) Z d k (2 π ) δ ′′ ( ǫ F − ǫ k ) e i k · y = (cid:20) − p + ( p + q ) i p + q ) · ∇ y (cid:21) cos( z ) + z sin( z ) − π k = cos( z ) + z sin( z )12 π k ( p + ( p + q ) ) − i π k F (2 p + q ) · y cos( z )where we used Eq.(B5) in the second step. Finally the third integral has the explicit form A = 112 Z d k (2 π ) h ( p p · k ) + ( p p · k )( ( p + q ) p + q ) · k )+ ( ( p + q ) p + q ) · k ) i δ ′′′ ( ǫ F − ǫ k ) e i k · y However, we only need the terms to second order in p and q . Therefore, it is sufficient to calculate A = 112 Z d k (2 π ) (cid:2) ( p · k ) + ( p · k )(( p + q ) · k ) + (( p + q ) · k ) (cid:3) δ ′′′ ( ǫ F − ǫ k ) e i k · y = − (cid:2) ( p · ∇ y ) + ( p · ∇ y )(( p + q ) · ∇ y ) + (( p + q ) · ∇ y ) (cid:3) Z dk (2 π ) δ ′′′ ( ǫ F − ǫ k ) e i k · y = − (cid:2) ( p · ∇ y ) + ( p · ∇ y )(( p + q ) · ∇ y ) + (( p + q ) · ∇ y ) (cid:3) (cid:18) z ) − z cos( z ) + 3 z sin( z )2 π k (cid:19) where to evaluate the integral over the delta function we used Eq.(B6). To perform the derivatives in the last termwe use that for an arbitrary function f ( z )( p · ∇ y )( q · ∇ y ) f ( z ) = ( p · y )( q · y ) k y (cid:20) d fdz − z dfdz (cid:21) + ( p · q ) k F y dfdz With this expression we find that A = − π k ( p + p · ( p + q ) + ( p + q ) )(cos( z ) + z sin( z )) − π k F [( p · y ) + ( p · y )(( p + q ) · y ) + (( p + q ) · y ) ] cos( z )Collecting our results A = A + A + A we therefore obtain the following expression for γ ( y , p , q ) γ ( y , p , q ) = A ( y , p , q ) + A ( y , q , p )= 1 π k F (cid:18) cos( z ) − i p + q ) · y cos( z ) + 112 k ( p + q + p · q )(cos( z ) + z sin( z )) −
124 [( p · y ) + ( q · y ) + 3(( p + q ) · y ) ] cos( z ) (cid:21) 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 Fouriertransforming Eq.(36). If we define the dimensionless quantity ¯ k = | k | / (2 k F ( r )) then ρ x is given by ρ x ( k | r ) = ( − k −
12 ¯ k ) θ (1 − ¯ k ) (C1)7which is simply the Fourier transform of the LDA exchange hole. We see that ρ x (0 | r ) = − α j are given by α = i nk F θ (1 − ¯ k ) (C2) α = − ¯ k nk θ (1 − ¯ k ) (C3) α = − n k (cid:16) θ (1 − ¯ k )¯ k + 14 δ (¯ k −
1) + 34 δ ′ (¯ k − (cid:17) (C4) α = 172 n k (cid:16) k θ (1 − ¯ k ) − δ (¯ k − (cid:17) (C5) α = 124 nk (cid:16) θ (1 − ¯ k )¯ k + δ (¯ k − (cid:17) (C6) [1] P. Hohenberg and W. Kohn,Phys. Rev. , B864 (1964)[2] R. M. Dreizler and E.K.U. Gross, Density FunctionalTheory: An Approach to the Quantum Many-Body Prob-lem
Springer (Berlin) (1990)[3] E. Engel and R. M. Dreizler,
Density Functional Theory:An Advanced Course
Springer (Heidelberg) (2011)[4] U. von Barth,
Basic Density-Functional Theory- anOverview , Phys. Scr. T , 9 (2004)[5] C. A. Ullrich,
Time-Dependent Density-Functional The-ory: Concepts and Applications
Oxford University Press,Oxford (2012)[6] W. Kohn and L. J. Sham,Phys. Rev. , A1133 (1965)[7] C. - O. Almbladh, Technical Report, University of Lund(1972)[8] O. Gunnarsson and B. I. Lundqvist,Phys. Rev. B , 4274 (1976)[9] J. P. Perdew, Phys. Rev. Lett. , 1665 (1985), Erra-tum: Phys. Rev. Lett. , 2370 (1985)[10] J. P. Perdew and Y. Wang, Phys. Rev. B , 8800 (1986),Erratum: Phys. Rev. B40, 3399 (1989)[11] A. D. Becke, J. Chem. Phys. , 7184 (1986).[12] A. D. Becke, J. Chem. Phys. , 1053 (1988).[13] A. D. Becke and M. R. Roussel,Phys. Rev. A , 3761 (1989)[14] V. Tschinke and T. Ziegler,Can. J. Chem. , 460 (1989).[15] J. P. Perdew and S. Kurth, ”Density Functionals forNon-relativistic Coulomb Systems in the New Century” inLecture Notes in Physics , pp 1-55, (2003), Editors:C. Fiolhais, F. Nogueira, M. A. L. Marques, Springer(Heidelberg).[16] J. P. Perdew, in Electronic Structure of Solids ’91 , pp.11-20, Eds. P. Ziesche and H. Eschrig, Akademie Verlag,Berlin (1991)[17] K. Burke, J. P. Perdew and Y. Wang, in
Electronic Den-sity Functional Theory: Recent Progress and New Direc-tions , p.81 , Eds. J.F. Dobson, G. Vignale and M. P.Das, Plenum, New York, (1997)[18] J. P. Perdew, K. Burke and Y. Wang, Phys. Rev. B , 16533 (1996), Erratum:Phys. Rev. B57, 14999 (1996)[19] D. C. Langreth and J. P. Perdew,Phys. Rev. B , 5469 (1980)[20] D. C. Langreth and M. J. Mehl,Phys. Rev. Lett. , 446 (1981)[21] D. C. Langreth and M. J. Mehl,Phys. Rev. B , 1809 (1983), ErratumPhys. Rev. B , 2310 (1984)[22] P. S. Svendsen and U. von Barth,Phys. Rev. B , 17402 (1996)[23] M. Springer, P. S. Svendsen and U. von Barth,Phys. Rev. B , 17392 (1996)[24] R. van Leeuwen, Adv. Quant. Chem. , 25 (2003)[25] P. E. Lammert, Int. J. Quant. Chem. , 1943 (2007)[26] P. E. Lammert, Phys. Rev. A , 012109 (2010)[27] H. Eschrig, The Fundamentals of Density FunctionalTheory , Teubner, Stuttgart (1996)[28] H. Eschrig, Phys. Rev. B , 205120 (2010)[29] H. Weyl, The classical groups: Their invariants and rep-resentations , Princeton University Press (1939)[30] G. Stefanucci and R. van Leeuwen,
NonequilibriumMany-Body Theory of Quantum Systems , CambridgeUniversity Press, Cambridge (2013)[31] E. K. U. Gross and R. M. Dreizler,Z. Phys. A , 103 (1981)[32] A. Putaja, E. R¨as¨anen, R. van Leeuwen, J. G. Vilhenaand M. A. L. Marques, Phys. Rev. B , 165101 (2012)[33] J. J. Duistermaat and J. A. C. Kolk, Distributions: The-ory and Applications , Springer, Dordrecht (2010)[34] L. J. Sham in
Computational methods in Band Theory ,eds. P. M. Marcus, J. F. Janak and A. R. Williams(Plenum, New York, 1971)[35] P. R. Antoniewicz and L. Kleinman,Phys. Rev. B , 6779 (1985)[36] L. Kleinman and S. Lee, Phys. Rev. B , 4634 (1988)[37] E. Engel and S. H. Vosko, Phys. Rev. B , 4940 (1990),Erratum: Phys. Rev. B , 1446 (1991)[38] P. S. Svendsen and U. von Barth,Int. J. Quant. Chem. , 351 (1995)[39] D. C. Langreth and S. H. Vosko, Adv. Quant. Chem. , 175 (1990)[40] J. P. Perdew, Physica B , 1 (1991)[41] K. J. H. Giesbertz, R. van Leeuwen and U. von Barth,Phys. Rev. A , 022514 (2013)[42] N. D. Mermin, Phys. Rev. , A1441, (1965)[43] E. Runge and E. K. U. Gross,Phys. Rev. Lett. , 997 (1984)[44] R. van Leeuwen, Phys. Rev. Lett. , 3863 (1999)[45] M. Ruggenthaler and R. van Leeuwen,EPL , 13001, (2011)[46] M. Ruggenthaler, K. J. H. Giesbertz, M. Penz and R.van Leeuwen, Phys. Rev. A , 052504 (2012) [47] G. Vignale and W. Kohn, Phys. Rev. , 2037 (1996)[48] G. Vignale and W. Kohn, in Electronic Density Func-tional Theory: Recent Progress and New Directions ,edited by J. Dobson, M. P. Das, and G. Vignale (PlenumPress, New York, 1998).[49] I. V. Tokatly, Phys. Rev. B , 165104 (2005)[50] I. V. Tokatly, Phys. Rev. B , 165105 (2005)[51] I. V. Tokatly, Phys. Rev. B , 125105 (2007)[52] J. Tao and G. Vignale,Phys. Rev. Lett.97