The Kohn-Sham system in one-matrix functional theory
TThe Kohn-Sham system in one-matrix functional theory
Ryan Requist ∗ and Oleg Pankratov Lehrstuhl f¨ur theoretische Festk¨orperphysikFriedrich-Alexander-Universit¨at Erlangen-N¨urnberg (Dated: November 21, 2018)A system of electrons in a local or nonlocal external potential can be studied with 1-matrixfunctional theory (1MFT), which is similar to density functional theory (DFT) but takes the one-particle reduced density matrix (1-matrix) instead of the density as its basic variable. Within 1MFT,Gilbert derived [PRB , 2111 (1975)] effective single-particle equations analogous to the Kohn-Sham (KS) equations in DFT. The self-consistent solution of these 1MFT-KS equations reproducesnot only the density of the original electron system but also its 1-matrix. While in DFT it isusually possible to reproduce the density using KS orbitals with integer (0 or 1) occupancy, in1MFT reproducing the 1-matrix requires in general fractional occupancies. The variational principleimplies that the KS eigenvalues of all fractionally occupied orbitals must collapse at self-consistencyto a single level, equal to the chemical potential. We show that as a consequence of the degeneracythe iteration of the KS equations is intrinsically divergent. Fortunately, the level shifting method,commonly introduced in Hartree-Fock calculations, is always able to force convergence. We introducean alternative derivation of the 1MFT-KS equations that allows control of the eigenvalue collapse byconstraining the occupancies. As an explicit example, we apply the 1MFT-KS scheme to calculatethe ground state 1-matrix of an exactly solvable two-site Hubbard model. PACS numbers: 71.15.Mb,31.15.xr
I. INTRODUCTION
Density functional theory (DFT) benefits from op-erating with the electron density, which as a functionof just three coordinates is much easier to work withthan the full many-body wavefunction. According to theHohenberg-Kohn (HK) theorem, the density of an elec-tron system in a local external potential v ( (cid:126)r ) may befound by minimizing a universal energy functional E v [ n ],whose basic variable is the density. Remarkably, the den-sity uniquely determines the ground state wavefunction(if it is nondegenerate), i.e., there can be only one groundstate wavefunction yielding a given density, no matterwhat v ( (cid:126)r ) is. However, if the external potential is non-local, then the density alone is generally not sufficient touniquely determine the ground state (see Appendix A fora simple example). Gilbert extended the HK theoremto systems with nonlocal and spin dependent externalpotential v ( x, x (cid:48) ), where x = ( (cid:126)r, σ ). It was proved thati) the ground state wavefunction is uniquely determinedby the ground state 1-matrix (one-particle reduced den-sity matrix) and ii) there is a universal energy functional E v [ γ ] of the 1-matrix, which attains its minimum at theground state 1-matrix. The 1-matrix is defined as γ ( x, x (cid:48) ) = N (cid:90) dx . . . dx N ρ ( x, x , . . . x N ; x (cid:48) , x , . . . x N ) , (1)where (cid:82) dx = (cid:80) σ (cid:82) d r and ˆ ρ = (cid:80) i w i | Ψ i (cid:105) (cid:104) Ψ i | is thefull N -electron density matrix with ensemble weights w i such that (cid:80) i w i = 1. An external potential may benonlocal with respect to the space coordinates and/orthe spin coordinates. For example, pseudopotentials arenonlocal in space, and Zeeman coupling − ( (cid:126) | e | /mc ) (cid:126)B · (cid:126)σ , where (cid:126)σ is the vector of Pauli matrices, is nonlocal inspin space. The coupling of electron motion and an ex-ternal vector potential, ( | e | / mc )( (cid:126)p · (cid:126)A + (cid:126)A · (cid:126)p ), may alsobe treated as a nonlocal potential because (cid:126)p is a differen-tial operator. It is rather intuitive that for such externalpotentials, which couple to the system in more complexways than the local potential v ( (cid:126)r ), it is necessary — inorder to permit statements analogous to the HK theorem— to refine the basic variable accordingly. Hence spin-DFT, whose basic variables are the density and themagnetization density, applies to systems with Zeemancoupling. Current-DFT, whose basic variables are thedensity and the paramagnetic current density, has thescope to treat systems in which the current is coupled toan external magnetic field. Generally, if one considers anexternal potential that is nonlocal in space and spin, thenecessary basic variable is the one-matrix, which con-tains all of the single-particle information of the system,including the density, magnetization density and param-agnetic current density.The DFT-type approach that takes the 1-matrix as ba-sic variable will be referred to here as 1-matrix functionaltheory (1MFT). As in DFT, an exact and explicit energyfunctional is generally unknown. An important differencebetween 1MFT and DFT is that the kinetic energy is asimple linear functional of the 1-matrix, while it is not aknown functional of the density. Thus, in 1MFT the onlypart of the energy not known explicitly is the electron-electron interaction energy W [ γ ]. Several approximate 1-matrix energy functionals have been proposed and testedrecently (see Refs. 7 and 8 and references therein.) No-tably, the so-called BBC n approximations, which aremodifications of the Buijse-Baerends functional, havegiven fairly accurate results for the potential energy a r X i v : . [ c ond - m a t . s t r- e l ] N ov curves of diatomic molecules and the momentum distri-bution and correlation energy of the homogeneous elec-tron gas. In Ref. 8, a density dependent fitting param-eter was introduced into the BBC1 functional such thatthe resulting functional yields the correct correlation en-ergy of the homogeneous electron gas at all values of den-sity. There is also the prospect of using 1MFT to obtainaccurate estimates for the band gaps of non-highly corre-lated insulators. Many of the approximate functionalsthat have been proposed are similar to an early approxi-mation by M¨uller. Actual calculations in 1MFT are more difficult thanin DFT. The energy functional E v [ γ ] must be minimizedin a space of higher dimension because the 1-matrix is amore complex quantity than the density. In the calcu-lations cited above, the energy has been minimized di-rectly by standard methods, e.g., the conjugate gradientmethod. In DFT the energy is generally not minimizedby such direct methods. Instead, the Kohn-Sham (KS)scheme provides an efficient way to find the groundstate density. In this scheme, one introduces an aux-iliary system of N noninteracting electrons, called theKS system, which experiences an effective local potential v s ( (cid:126)r ). This effective potential is a functional of the den-sity such that the self-consistent ground state of the KSsystem reproduces the ground state density of the inter-acting system. It is interesting to ask whether there isalso a KS scheme in 1MFT. The question may be statedas follows: does there exist a 1-matrix dependent effec-tive potential v s ( x, x (cid:48) ) such that, at self-consistency, asystem of noninteracting electrons experiencing this po-tential reproduces the exact ground state 1-matrix of theinteracting system? Although Gilbert derived such aneffective potential, the implications were thought to be“paradoxical” because the KS system was found to have ahigh (probably infinite) degree of degeneracy. Evidently,the KS eigenvalues in 1MFT do not have the meaning ofapproximate single-particle energy levels, in contrast toDFT and other self-consistent-field theories, where theeigenvalues may often be interpreted as the negative ofionization energies, owing to Koopmans’ theorem. Thestatus of the KS scheme in 1MFT appears to have re-mained unresolved, and recently it has been arguedthat the KS scheme does not exist in 1MFT. Gilbert derived the KS equations from the stationaryprinciple for the energy. The KS potential was foundto be v s ( x, x (cid:48) ) = v ( x, x (cid:48) ) + δW/δγ ( x (cid:48) , x ) . (2)In this article, we propose an alternative derivation of theKS equations, which, in our view, gives insight into thenature of the “paradoxical” degeneracy of the KS system.One-matrix energy functionals are often expressed interms of the so-called natural orbitals and occupationnumbers. This makes them similar to “orbital depen-dent” functionals in DFT. The natural orbitals are theeigenfunctions of the 1-matrix, and the occupation num-bers are the corresponding eigenvalues. These quanti- ties play a central role in 1MFT. Recently, it was shownthat when a given energy functional is expressed in termsof the natural orbitals and occupation numbers, the KSpotential can be found by using a chain rule to evaluatethe functional derivative in Eq. 2. Although the concept of the KS system can indeedbe extended to 1MFT, it has in this setting some veryunusual properties. In particular, the KS orbitals mustbe fractionally occupied, for otherwise the KS systemcould not reproduce the 1-matrix of the interacting sys-tem, which always has noninteger eigenvalues (occupa-tion numbers). This is different from the situation inDFT, where it is usually possible to reproduce the den-sity using only integer (0 or 1) occupation numbers, orin any case, only a finite number of fractionally occu-pied states. Due to the necessity of fractional occupationnumbers, the 1MFT-KS system cannot be described bya single Slater determinant. However, we find that it canbe described by an ensemble of Slater determinants, i.e.,a mixed state. In order that the variational principle isnot violated, all the states that comprise the ensemblemust be degenerate. This implies that the eigenvalues ofall fractionally occupied orbitals collapse to a single level,equal to the chemical potential. The degeneracy has im-portant consequences for the solution of the KS equationsby iteration. We prove that the iteration of the KS equa-tions is intrinsically divergent because the KS system hasa divergent response function χ s = δγ/δv s at the groundstate. Fortunately, convergence can always be obtainedwith the level shifting method. To illustrate explicitlythe unique properties of the 1MFT-KS system, we applyit to a simple Hubbard model with two sites. The modeldescribes approximately systems which have two local-ized orbitals with a strong on-site interaction, e.g., thehydrogen molecule with large internuclear separation. The Schr¨odinger equation for this model is exactly solv-able, and we find that the KS equations in 1MFT and inDFT can be derived analytically. It is interesting to com-pare 1MFT and DFT in this context. We demonstratethat divergent behavior will appear also in DFT whenthe operator 1 − χ s χ − , where χ and χ s are the densityresponse functions of the interacting and KS systems, re-spectively, has any eigenvalue with modulus greater than1. In this expression the null space of χ is assumed to beexcluded.This article is organized as follows. In Sec. II, we derivethe KS equations in 1MFT and discuss how to solve themself-consistently by iteration. In Sec. III, we comparethree approaches to ground state quantum mechanics —direct solution of the Schr¨odinger equation, 1MFT andDFT — by using them to solve the two-site Hubbardmodel. II. KOHN-SHAM SYSTEM IN 1MFT
It is not obvious that a KS-type scheme exists in 1MFTfor the following reason. Recall that in DFT the KS sys-tem consists of N noninteracting particles and reproducesthe density of the interacting system. The density of theKS system, if it is nondegenerate, is the sum of contri-butions of the N lowest energy occupied orbitals n ( (cid:126)r ) = occ (cid:88) i | φ i ( (cid:126)r ) | . (3)On the other hand, in 1MFT the KS system should repro-duce the 1-matrix of the interacting system. The eigen-functions of the 1-matrix are the so-called natural or-bitals, and the eigenvalues are the corresponding occupa-tion numbers. Occupying the N lowest energy orbitalsin analogy to (3), one obtains γ ( x, x (cid:48) ) = occ (cid:88) i φ i ( x ) φ ∗ i ( x (cid:48) ) . (4)Such an expression, in which the orbitals have only inte-ger (0 or 1) occupation, cannot reproduce the 1-matrix ofan interacting system because the orbitals of an interact-ing system have generally fractional occupation (see thediscussion in the following section.) The difference be-tween the 1-matrix in (4) and the 1-matrix of an interact-ing system is clearly demonstrated by the so-called idem-potency property. The 1-matrix in (4) is idempotent, i.e., (cid:82) dx (cid:48) γ ( x, x (cid:48) ) γ ( x (cid:48) , x (cid:48)(cid:48) ) = γ ( x, x (cid:48)(cid:48) ), while the 1-matrix ofan interacting system is never idempotent. However, ifthe KS system is degenerate and its ground state is anensemble state, the 1-matrix becomes γ ( x, x (cid:48) ) = (cid:88) i f i φ i ( x ) φ ∗ i ( x (cid:48) ) (5)with fractional occupation numbers f i . The N -particleground state density matrix of the KS system is ˆ ρ s = (cid:80) i w i | Φ i (cid:105) (cid:104) Φ i | , where the Φ i are Slater determinantseach formed from N degenerate KS orbitals. The oc-cupation numbers f i are related to the ensemble weights w i by f i = (cid:88) j w j Θ ji , (6)where Θ ji equals 1 if φ i is one of the orbitals in thedeterminant Φ j and 0 otherwise. A. Derivation of the 1MFT Kohn-Sham equations
In this section, we discuss Gilbert’s derivation of theKS equations in 1MFT and propose an alternative deriva-tion. We begin by reviewing the definition of the univer-sal 1-matrix energy functional E v [ γ ].One-matrix functional theory describes the groundstate of a system of N electrons with the Hamiltonianˆ H = (cid:80) Ni =1 (ˆ t i + ˆ v i ) + ˆ W , where ˆ t = −∇ r / v is the local or nonlocal external poten-tial operator, and ˆ W = (cid:80) i 1) orbitals, for otherwise the energy could belowered. When q i = q gsi for all i , γ min = γ gs , and (19) im-plies (cid:15) j ( { q gsi } ) = 0 for all fractionally occupied orbitals.Thus, we find that the KS eigenvalues of all orbitals thatare fractionally occupied in the ground state must col-lapse to a single level when the chosen set of occupa-tion numbers approach their ground state values, i.e., as q i → q gsi . In Gilbert’s derivation of the 1MFT-KS equa-tions, all fractionally occupied KS orbitals were found tohave the eigenvalue (cid:15) i = µ , where µ is the chemical po-tential. As we have not introduced the chemical potentialin our derivation (we consider a system with a fixed num-ber of electrons), the eigenvalues collapse to 0 instead of µ . The above arguments do not apply to orbitals withoccupation numbers exactly 0 or 1 because these values lie on the boundary of the allowed interval [0 , 1] specifiedby ENR condition (iii). All that can be concluded fromthe fact that E v has a minimum in the ENR space is (cid:15) j ≥ f j = 0 and (cid:15) j ≤ f j = 1. States with occupation numbers exactly 0 or1 have been called “pinned states.” Instances of suchstates in real systems have been reported, though theiroccurrence is generally considered to be exceptional. Due to the collapse of the eigenvalues at the groundstate, the KS Hamiltonian becomes the null operatorˆ h [ γ gs ] = ˆ0 (20)in the subspace of fractionally occupied orbitals. This isof course analogous to the familiar condition df /dx = 0for the extremum of a function f ( x ). Gilbert described a similar result (with µ ˆ I replacing ˆ0) as “paradoxical,” astatement that has been repeated. The problem with(20) is that while we expect the KS Hamiltonian to de-fine the natural orbitals, any state is an eigenstate of thenull operator. However, the KS Hamiltonian is a func-tional of the 1-matrix, and when the occupation num-bers are perturbed from their ground state values, thedegeneracy is lifted and the KS Hamiltonian does defineunique orbitals. In the KS scheme outlined above, thiscorresponds to the optimization of the orbitals with oc-cupation numbers fixed to values q i , perturbed from theground states values. In the limit that the occupationnumbers approach their ground state values, the opti-mal orbitals approach the ground state natural orbitals.The degenerate eigenvalues generally split linearly withrespect to perturbations away from the ground state. Inparticular, ∂(cid:15) i ∂q j = (cid:90) dydy (cid:48) (cid:104) φ i | δ ˆ hδγ ( y, y (cid:48) ) ∂γ ( y, y (cid:48) ) ∂q j | φ i (cid:105) = − (cid:90) dxdx (cid:48) dydy (cid:48) φ ∗ i ( x ) φ ∗ j ( y (cid:48) ) χ − ( xx (cid:48) , yy (cid:48) ) × φ i ( x (cid:48) ) φ j ( y ) . (21)Here, χ is the static response function defined as χ ( x, x (cid:48) ; y, y (cid:48) ) = δγ ( x, x (cid:48) ) δv ( y, y (cid:48) ) . (22)The relation δh/δγ = − χ − used in (21) is derived inthe following section. If χ has a null space, its inverseis defined only on a restricted space. For example, (21)does not apply to pinned states as there is a null vectorassociated with each pinned state (see below).Our derivation of the 1MFT-KS equations in fact as-sumes that the static response function χ of the interact-ing system has no null vectors except for those associatedwith pinned states and a constant shift of the potential.We now show that this guarantees G v is stationary. Ifthe interacting system has any other null vectors we canno longer be certain that G v is stationary and the deriva-tion should be modified as described below. We haveremarked already (see Ref. 34) that G v is stationaryin the VR space, i.e., it satisfies the stationary condi-tion δ G v = 0 with respect to an arbitrary variation ofthe 1-matrix in the VR space. However, our derivationof the KS equations requires G v to be stationary in theENR space. As the VR space is a subspace of the ENRspace, this is a stronger condition. The assumption that χ has no null vectors (apart from those associated withpinned states) is equivalent to assuming that any ENRvariation (apart from variations of the pinned occupa-tion numbers) is also a VR variation. For if χ has nonull vectors, then it is invertible and any ENR variation δ ˆ γ can be induced by the perturbation δ ˆ v = χ − δ ˆ γ . Hence, with the above assumption, G v is guaranteed tobe stationary with respect to any ENR variation. Theabove arguments do not apply to variations of pinnedoccupation numbers because there are null vectors asso-ciated with such variations; nevertheless, G v is station-ary with respect to such variations as this is maintainedby the Lagrange multipliers (cid:15) i . It is now clear how tomodify the derivation of the KS equations when χ hasadditional null vectors. By introducing new Lagrangemultipliers, the additional null vectors can be treated inanalogy with the pinned states. For example, suppose χ has one additional null vector ˆ u = (cid:80) ij u ij | φ i (cid:105) (cid:104) φ j | ,where ˆ u is hermitian and T r (ˆ u ˆ u ) = 1. The new energyfunctional G (cid:48) v [ γ ] = G v [ γ ] − κ ( γ u − p u ) will be station-ary with respect to an arbitrary variation in the ENRspace. Here, the Lagrange multiplier κ enforces the con-straint γ u = p u , where γ u = T r (ˆ γ ˆ u ) is the component ofˆ γ corresponding to ˆ u . The stationary condition δ G (cid:48) v = 0leads to the set of equations h ij − (cid:15) j δ ij − κu ij = 0 inthe basis of natural orbitals. The 1-matrix that satisfiesthese equations can be found by solving self-consistentlythe eigenvalue equation ˆ h | ξ i (cid:105) = ω i | ξ i (cid:105) together with γ ( x, x (cid:48) ) = (cid:80) ijk q i S ji S ∗ ki ξ j ( x ) ξ ∗ k ( x (cid:48) ), where S is the uni-tary matrix that diagonalizes the matrix u in the basis ofnatural orbitals, i.e., SuS † is diagonal. The energy of theself-consistent solution defines a function G (cid:48) v ( { q i } , p u ),whose minimum with respect to { q i } and p u is the groundstate energy. It may not be known in advance whetherthe response function of a given interacting system willhave null vectors. Therefore, it is helpful to understandhow null vectors occur.Null vectors of χ are connected with the so-callednonuniqueness problem in various extensions ofDFT. A system with the ground state Ψ is said to havea nonuniqueness problem if there is more than one ex-ternal potential for which Ψ is the ground state. TheSchr¨odinger equation defines a unique map from the ex-ternal potential to the ground state wavefunction (if it isnondegenerate), but when there is more than one exter-nal potential yielding the same ground state wavefunc-tion, the map cannot be inverted. In 1MFT the general-ity of the external potential (nonlocal in space and spincoordinates) allows greater scope for nonuniqueness thanin the other extensions of DFT. Of course, every degreeof nonuniqueness is a null vector of χ because if δ ˆ v does not change the ground state wavefunction, it does notchange the 1-matrix either, and hence it is a null vector.In fact, every null vector of χ is caused by nonuniqueness;the existence of a null vector δ ˆ v that induced a nonzero δ Ψ would contradict the one-to-one relationship γ ↔ Ψ proved by the extension of the HK theorem to 1MFT. It was mentioned above that there are null vectors asso-ciated with the pinned states. Suppose φ k is a naturalorbital with occupation number f k = 0 in the groundstate. The perturbation of the external potential δ ˆ v = λ | φ k (cid:105) (cid:104) φ k | does not change the ground state if the systemhas an energy gap between the ground state and excitedstates and λ is small enough because δ ˆ V Ψ = 0, where δ ˆ V = (cid:82) dxdx (cid:48) ˆ ψ † ( x ) δv ( x, x (cid:48) ) ˆ ψ ( x (cid:48) ) and ˆ ψ and ˆ ψ † are fieldoperators. If f k = 1, δ ˆ V Ψ = Ψ and the ground stateis again unchanged by the perturbation. The “vector” | φ k (cid:105) (cid:104) φ k | is therefore a null vector of χ if φ k is a pinnedstate. Another type of nonuniqueness, which has beencalled systematic nonuniqueness, is related to constantsof the motion. Suppose ˆ A = (cid:82) dxdx (cid:48) ˆ ψ † ( x ) a ( x, x (cid:48) ) ˆ ψ ( x (cid:48) ) isa constant of the motion. The ground state, if it is nonde-generate, is an eigenstate of ˆ A as constants of the motioncommute with the Hamiltonian. If the system has an en-ergy gap between the ground state and the first excitedstate, then a perturbation δ ˆ V = λ ˆ A will not change theground state wavefunction if λ is small enough. Thus, ˆ A is a null vector of χ .As in DFT, an exact and explicit expression for theuniversal energy functional in 1MFT is unknown in gen-eral. In actual calculations it is usually necessary to useapproximate functionals. Many of the approximate en-ergy functionals that have been introduced are expressedin terms of the natural orbitals and occupation numbers.Such functionals are valid 1-matrix energy functionals,but as the dependence on the 1-matrix is implicit ratherthan explicit, they have been called “implicit” function-als. Recently, the KS equations were derived for thiscase. It was found that the contribution to the KS po-tential from electron-electron interactions can be evalu-ated by applying the following chain rule to (17), w ( x, x (cid:48) ) = δWδγ ( x (cid:48) , x )= (cid:88) i (cid:90) dy δWδφ i ( y ) δφ i ( y ) δγ ( x (cid:48) , x )+ (cid:88) i (cid:90) dy δWδφ ∗ i ( y ) δφ ∗ i ( y ) δγ ( x (cid:48) , x )+ (cid:88) i ∂W∂f i δf i δγ ( x (cid:48) , x ) . (23) B. Iteration of the KS equations In this section we show that the “straightforward” pro-cedure for iterating the KS equations (5) and (16) is in-trinsically divergent.The KS equations are nonlinear because the KS Hamil-tonian itself depends on the 1-matrix. In favorable casessuch nonlinear equations can be solved by iteration.Given a good initial guess for the 1-matrix, iteration maylead to the self-consistent solution corresponding to theground state. In order to iterate (16), one needs an al-gorithm to define the 1-matrix of iteration step n + 1from the 1-matrix of step n , i.e., one needs to “close”the KS equations. In the previous section, we saw thatin the 1MFT-KS scheme the occupation numbers f i areheld fixed during the optimization of the natural orbitals.The following is a “straightforward” algorithm that op-timizes the natural orbitals: i) the KS Hamiltonian forstep n + 1 is defined byˆ h ( n +1) = ˆ h [ γ ( n ) ] , (24)where ˆ h is given by (15) and γ ( n ) is the 1-matrix of it-eration step n ; ii) the eigenstates of ˆ h ( n +1) are taken asthe natural orbitals of step n + 1; iii) the 1-matrix of step n +1 is constructed from the natural orbitals of step n +1by the expression γ ( n +1) ( x, x (cid:48) ) = (cid:88) i f i φ ( n +1) i ( x ) φ ∗ ( n +1) i ( x (cid:48) ) . (25)Let u i be the eigenstates of the KS Hamiltonian ˆ h ( n +1) .In operation (ii), the natural orbitals φ ( n +1) i are chosenfrom among the u i such that φ ( n +1) i has maximum over-lap with φ ( n ) i , i.e., T r ((ˆ γ ( n +1) − ˆ γ ( n ) ) ) is the minimumpossible. If this procedure converges to the stationary1-matrix γ min giving the lowest energy possible for thefixed set of occupation numbers, then it defines the func-tion G v ( { f i } ), introduced in the previous section, forwhich the ground state energy is the absolute minimum.Unfortunately, for any set of occupation numbers { f i } sufficiently close to the ground state occupation num-bers, this procedure does not converge to γ min . In otherwords, the “straightforward” algorithm defines an itera-tion map for which the ground state γ gs is an unstablefixed point. This is a consequence of the degeneracy ofthe KS spectrum at the ground state.The divergence of the iteration map is revealed by alinear analysis of the fixed point. Suppose the occupationnumbers are fixed to values perturbed from their groundstate values by δf i . Let us consider an iteration step n and ask whether the next iteration takes us closer to thestationary point γ min that gives the minimum energy forthe fixed occupation numbers. The linearization of theiteration map at the stationary point gives δ ˆ γ ( n +1) = ˆ γ ( n +1) − ˆ γ min ≈ ˆ χ s [ γ min ] (cid:0) ˆ v ( n +1) s − ˆ v mins (cid:1) = ˆ χ s [ γ min ] (cid:0) ˆ h ( n +1) − ˆ h min (cid:1) ≈ − ˆ χ s [ γ min ] ˆ χ − δ ˆ γ ( n ) , (26)where ˆ v s = ˆ v +ˆ v H +ˆ v xc is the KS potential. The responsefunction χ was defined in (22). The KS response function is χ s ( x, x (cid:48) ; y, y (cid:48) ) = δγ ( x, x (cid:48) ) δv s ( y, y (cid:48) )= (cid:88) i (cid:88) j (cid:54) = i f i − f j (cid:15) i − (cid:15) j × φ j ( x ) φ ∗ i ( x (cid:48) ) φ ∗ j ( y ) φ i ( y (cid:48) ) . (27)In the last line of (26), we have usedˆ h ( n +1) − ˆ h min ≈ − ˆ χ − δ ˆ γ ( n ) , (28)which can be established by the following arguments.First considerˆ h ( n +1) − ˆ h min ≈ δ ˆ h [ γ ] δγ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ min δ ˆ γ ( n ) , (29)which follows from (24). The KS Hamiltonian is an im-plicit functional of the 1-matrix, and (29) defines thefirst order change of the KS Hamiltonian with respectto a perturbation of that 1-matrix. Since the occupationnumbers are close to their ground state values, we maymake the replacement δ ˆ hδγ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ min → δ ˆ hδγ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) γ gs , (30)which is valid to O (max( | δf i | )). Thus, to establish (28)we need to show δ ˆ h/δγ = − ˆ χ − at the ground state 1-matrix γ gs . The KS Hamiltonian ˆ h is associated withthe original many-body Hamiltonian ˆ H , which has ex-ternal potential v ( x, x (cid:48) ). According to (20) , ˆ h [ γ gs ] = ˆ0.Consider now a Hamiltonian ˆ H (cid:48) with a slightly differ-ent external potential v (cid:48) ( x, x (cid:48) ) = v ( x, x (cid:48) ) + δv ( x, x (cid:48) ) suchthat its ground state 1-matrix is γ (cid:48) gs = γ gs + δγ . Theassociated KS Hamiltonian is ˆ h (cid:48) = ˆ h + ˆ v (cid:48) − ˆ v . At the newground state, ˆ h (cid:48) [ γ (cid:48) gs ] = ˆ0. This allows us to relate δ ˆ h to δ ˆ v as δ ˆ h = ˆ h [ γ (cid:48) gs ] − ˆ h [ γ gs ]= ˆ h [ γ (cid:48) gs ] + δ ˆ v − δ ˆ v = − δ ˆ v. (31)Finally, using δ ˆ v = ˆ χ − δ ˆ γ we obtain δ ˆ h/δγ = − ˆ χ − ,which verifies (28).Returning to the question of convergence, we see that(26) implies that the next iteration takes us fartherfrom the stationary point γ min . The reason is that (cid:12)(cid:12) det ˆ χ s ˆ χ − (cid:12)(cid:12) > γ is sufficiently close to the groundstate. According to (21) the KS response diverges as γ → γ gs because (cid:15) i − (cid:15) j ∼ O (max( | δf i | )). For a fixed setof occupation numbers sufficiently close to their groundstate values, the moduli of all eigenvalues of the opera-tor ˆ χ s ˆ χ − become greater than 1 (the null space of χ isassumed to be excluded). Therefore, a perturbation δγ from the ground state is amplified by iteration, and theground state is an unstable fixed point of the iterationmap. A fixed point is stable if and only if all eigenvaluesof the linearized iteration map have modulus less than 1. C. Level shifting method In the previous section, it was shown that the“straightforward” iteration of the KS equations is intrin-sically divergent. To obtain a practical KS scheme theiteration map must be modified. In this section, we con-sider the level shifting method and by linearizing themodified iteration map we obtain a criterion for conver-gence.Intrinsic divergent behavior can be encountered alsoin the Hartree-Fock approximation, and various modifi-cations of the iteration procedure have been introduced,for example, Hartree damping (also called configurationmixing) and level shifting. The level shifting method isparticularly attractive in 1MFT because it can preventthe collapse of eigenvalues that is the origin of divergentbehavior. Indeed, it has been shown that the level shift-ing method is capable of giving a convergent KS schemein 1MFT. In the “straightforward” iteration procedure, thechange of the orbitals, to first order, from iteration step n to iteration step n + 1 is δφ i ( x ) = φ ( n +1) i ( x ) − φ ( n ) i ( x )= (cid:88) j (cid:54) = i (cid:104) φ j | ˆ h ( n +1) − ˆ h ( n ) | φ i (cid:105) (cid:15) i − (cid:15) j φ j ( x ) , (32)where ˆ h ( n ) is the KS Hamiltonian for iteration step n and in the second line the orbitals and eigenvalues arefrom iteration step n . In the level shifting method, thefirst order change in the orbitals given by (32) is alteredby applying the shifts (cid:15) i → (cid:15) i + ζ i to the eigenvaluesin the denominator. To first order, this modification isequivalent to adding the term ˆ∆ = (cid:80) i ζ i (cid:12)(cid:12) φ ( n ) i (cid:11)(cid:10) φ ( n ) i (cid:12)(cid:12) tothe KS Hamiltonian for step n + 1. Let ˆ h ζ = ˆ h + ˆ∆define the level shifted Hamiltonian. Repeating the linearanalysis of the previous section for the iteration map forthis level shifted Hamiltonian, we findˆ δγ ( n +1) ≈ ˆ χ s [ γ min ] (cid:0) ˆ v ( n +1) s − ˆ v mins (cid:1) = ˆ χ s [ γ min ] (cid:0) ˆ h ( n +1) ζ − ˆ h minζ (cid:1) ≈ (cid:0) − ˆ χ s [ γ min ] ˆ χ − + ˆΩ (cid:1) ˆ δγ ( n ) , (33)where we have defined the operator ˆΩ with the kernelΩ( xx (cid:48) , yy (cid:48) ) = (cid:90) dzdz (cid:48) χ s ( xx (cid:48) , zz (cid:48) ) δ ∆( zz (cid:48) ) δγ ( yy (cid:48) )= (cid:88) i (cid:88) j (cid:54) = i φ j ( x ) φ ∗ i ( x (cid:48) ) φ ∗ j ( y ) φ i ( y (cid:48) ) . (34) From the last line of (33), we obtain a criterion for theconvergence of the iteration map. All of the eigenvaluesof the operator ˆ A = − ˆ χ s [ γ min ] ˆ χ − + ˆΩ (35)must have modulus less than 1. The dependence on thelevel shift parameters ζ i enters only through the shiftedeigenvalues in the denominator of χ s . The level shiftingmethod is effective because it prevents the divergence of χ s at the ground state and there is a cancellation betweenthe two terms in (35). Unfortunately, the convergencecriterion depends on χ , which is unknown at the outsetof a 1MFT calculation. In Sec. III, the level shiftingmethod is applied in an explicit example and the abovecriterion is verified. D. Properties of the KS system The distinguishing feature of the KS system in 1MFTis the degeneracy of the eigenvalue spectrum. This hassurprising consequences. It was shown in section II Athat the KS eigenvalue spectrum splits linearly as wemove away from the ground state 1-matrix. There-fore, the total KS energy changes linearly with respectto the displacement, i.e., E s [ γ ] − E s [ γ gs ] ∝ δγ , where E s [ γ ] = tr (ˆ h [ γ ] γ ) (for a specific example see Fig. 4in Sec. III B 2). This is surprising because such linearchanges do not occur for the energy functional E v (inthe VR space). The immediate implication is that E s [ γ ]is not stationary at the ground state. While this causesno difficultly in principle — we need only the functional E v to be stationary — it is intimately connected with thedivergence of the iteration map. Precisely at the groundstate E s [ γ gs ] = (cid:80) (cid:48) i (cid:15) i , where the prime indicates thatonly the pinned states with f i = 1 contribute to the sum.Away from the ground state the KS eigenvalue spectrumsplits, and E s [ γ ] is a multivalued functional due to thechoice implied in occupying the new KS levels. This is thesame choice encountered in the iteration of the KS equa-tions (see Sec. II B), where the natural orbitals φ ( n +1) i areselected from among the eigenstates of the KS Hamilto-nian. Near the self-consistent solution, there will be onesuch choice for which the resulting γ ( n +1) is very close to γ ( n ) .It was shown in the preceding section that the static re-sponse function of the KS system diverges at the groundstate. Thus, even an infinitesimal perturbation δ ˆ v s mayinduce a finite change of γ . At the ground state, all ofthe natural orbitals, except those which have an occu-pation number that is degenerate, are uniquely defined.The natural orbitals which belong to a degenerate occu-pation number are only defined modulo unitary rotationin the degenerate subspace. When a perturbation is in-troduced, the natural orbitals change discontinuously tothe eigenstates of the perturbed KS Hamiltonian ˆ h = δ ˆ v .These eigenstates may be any functions in the degenerateHilbert space because δ ˆ v is arbitrary. III. TWO-SITE HUBBARD MODEL The 1MFT-KS system has some unusual features, suchas the collapse of the KS eigenvalues at the ground state,so it is desirable to derive explicitly the KS equations for asimple model. The Hubbard model on two sites providesa convenient example because it is exactly solvable andespecially easy to interpret. Also, analytic expressions forthe 1-matrix energy functional and KS Hamiltonian canbe obtained. In the following sections, for the purposeof comparison, we find the ground state of the two-siteHubbard model by three methods — direct solution ofthe Schr¨odinger equation, 1MFT and DFT. A. Direct solution The Hamiltonian of the two-site Hubbard model isˆ H = ˆ T + ˆ U + ˆ V withˆ T = − (cid:88) σ (cid:16) t c † σ c σ + t c † σ c σ (cid:17) ˆ U = U (ˆ n ↑ ˆ n ↓ + ˆ n ↑ ˆ n ↓ )ˆ V = V 12 (ˆ n − ˆ n ) , (36)where t = t = t , c † iσ and c iσ are the creation andannihilation operators of an electron at site i with spin σ , and ˆ n i = (cid:80) σ c † iσ c iσ . We consider only the sector ofstates with N = 2 and S z = 0, i.e., a spin unpolarizedsystem. In this sector, the eigenstates of ˆ T + ˆ U areΦ = √ yxxy , Φ = √ − , Φ = √ − , Φ = √ x − y − yx , (37)in the site basis (cid:8) c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) (cid:9) . The following variables have been intro-duced x = cos( π/ − α / y = sin( π/ − α / α = U/ t with 0 ≤ α ≤ π/ 2. The eigenvalues ofˆ T + ˆ U for the states Φ i are λ = − By , λ = 0 λ = B ( x − y ) , λ = Bx , where B = √ U + T and T = 4 t . Φ , Φ , and Φ aresinglet states ( S = 0) and Φ is a triplet state with S = 1and S z = 0. We will omit Φ from consideration as it isnot coupled to the other states by the spin-independentexternal potential chosen in (36). The Hamiltonian maybe written as ˆ H = λ ˆ I + B ˆ K , where K = νy νy x νx νx (38) in the basis (Φ , Φ , Φ ). We have defined the dimension-less variable ν = V /B . The secular equation (cid:12)(cid:12) ˆ K − κ i ˆ I (cid:12)(cid:12) =0 is κ i − ( x + 1) κ i + ( x − ν ) κ i + ν y = 0 . (39)The normalized eigenvectors of ˆ H areΨ i = 1 η i ν xyνxκ i β i (40)and have energy E i = λ + Bκ i for i = 0 , , 3, where κ i isa root of the secular equation (39). We have also defined β i = κ i ( κ i − x ) − ν y (41)and η i = (cid:2) κ i ( ν (3 x − 1) + y ) + κ i ( x y (2 ν − ν y ( ν − y ) (cid:3) / . (42)The two dimensionless energy scales of the system arethe interaction strength U/T and the bias V /T . Thebehavior of the system with respect to these energy scalesis illustrated in Fig. 1. The quantity m = ( n − n ) / n i is the average ground state occupancy of site i , is plotted with respect to the external potential V forvarious values of the interaction strength U . For the FIG. 1: [color online] The density variable m = ( n − n ) / ν = V /B ( B = √ T + U ). Curves for U/T =(1 / , / , , 4) are shown as (solid [blue], dotted [green],dash-dotted [light blue], dashed [red]) curves, respectively. ground state, m = (cid:104) Ψ | ˆ m | Ψ (cid:105) = 2 νκ x η (cid:0) κ − x (cid:1) , (43)0where ˆ m = (ˆ n − ˆ n ) / 2. A weakly interacting system(e.g., the solid [blue] curve in Fig. 1) responds strongly tothe external potential. In contrast, a strongly interactingsystem (e.g., the dashed [red] curve) responds weakly upto a threshold V /B ∼ B ≈ U .) This behavior has a simple interpretation: inorder for the external bias to induce charge transfer, itmust overcome the on-site Hubbard interaction. In thelimit U → ∞ , the curve develops step-like behavior near V /B ∼ ± B. Solution by 1MFT In the first part of this section, we derive the en-ergy functional and KS Hamiltonian. In the second,we demonstrate the divergence of the iteration of theKS equations. In the third, we use the level shiftingmethod to obtain a convergent KS scheme. 1. Energy functional and KS Hamiltonian For lattice models such as the Hubbard model, the 1-matrix is defined as γ ( iσ, jτ ) = (cid:10) Ψ (cid:12)(cid:12) c † iσ c jτ (cid:12)(cid:12) Ψ (cid:11) . (44)One may ask whether the HK theorem (or Gilbert’s ex-tension in 1MFT) applies when the density (or 1-matrix)is defined over a discrete set of points, i.e., when thecontinuous density function n ( r ) is replaced by the siteoccupation numbers n i . This has been investigated, and it was found that the HK theorem remains valid.We consider here only spin unpolarized states ( S z = 0).Accordingly, we define the spatial γ ( ij ) = (cid:88) σ γ ( iσ, jσ ) . (45)The 1-matrix may be expressed as, cf. (5), γ ( ij ) = (cid:88) α f α φ α ( i ) φ ∗ α ( j ) , (46)where φ α are the spatial natural orbitals. As our sys-tem is spin unpolarized, the spin up and spin down spin-orbitals have the same spatial factors. Therefore, in (46)each spatial orbital φ α may be occupied twice (once by aspin up electron and once by a spin down electron), i.e.,0 ≤ f α ≤ 2. It is convenient to parametrize the naturalorbitals as φ a = (cid:18) cos( θ/ θ/ (cid:19) and φ b = (cid:18) sin( θ/ − cos( θ/ (cid:19) . (47)In terms of this parametrization, the 1-matrix in the sitebasis is γ = I + A (cos θσ z + sin θσ x )= I + (cid:126)γ · (cid:126)σ ; (cid:126)γ = ( γ x , γ z ) (48) where σ i are the Pauli matrices and A = ( f a − f b ) / α .For the two-site Hubbard model in the sector of sin-glet states with N = 2 and S z = 0, (44) may beinverted to express Ψ = Ψ [ γ ]. Explicitly, we findΨ = cos( α/ aa − sin( α/ bb , where Φ ii is the Slaterdeterminant composed of the natural spin orbitals φ i ↑ and φ i ↓ ( i = a, b ). The terms of the energy functional E [ γ ] = T [ γ ] + U [ γ ] + V [ γ ] are found to be T [ γ ] = (cid:68) Ψ (cid:12)(cid:12)(cid:12) ˆ T (cid:12)(cid:12)(cid:12) Ψ (cid:69) = − tA sin θU [ γ ] = (cid:68) Ψ (cid:12)(cid:12)(cid:12) ˆ U (cid:12)(cid:12)(cid:12) Ψ (cid:69) = U − U (cid:16) (cid:112) − A (cid:17) sin θV [ γ ] = (cid:68) Ψ (cid:12)(cid:12)(cid:12) ˆ V (cid:12)(cid:12)(cid:12) Ψ (cid:69) = V A cos θ. (49)The electron-electron interaction energy functional U [ γ ]agrees with the general exact result for 2-electron closedshell systems . We may partition U [ γ ] into theHartree energy E H [ γ ] = 12 (cid:88) ij n i n j U δ ij = U (cid:0) A cos θ (cid:1) (50)and the exchange-correlation energy E xc [ γ ] = U [ γ ] − E H [ γ ]= − U (cid:18) 12 + (cid:112) − A (cid:19) − U (cid:16) A + (cid:112) − A (cid:17) cos θ. (51)In Sec. II A the KS Hamiltonian was derived from thestationary principle for the energy. For the present modelthe KS Hamiltonian is a real 2 × h ( ij ) = (cid:10) (cid:12)(cid:12) c i ˆ hc † j (cid:12)(cid:12) (cid:11) . This matrixmay be expressed as h = (cid:126)h · (cid:126)σ with h x = − B (cid:18) cos α − sin α sin α cos α sin θ (cid:19) − B α (1 + sin α ) sin α cos α sin θ cos θh y = 0 h z = B α (1 + sin α ) sin α cos α sin θ cos θ + V . (52)In these expressions the variable α represents the depen-dence on the occupation numbers through the definition α = cos − A = cos − (( f a − f b ) / α = tan − ( U/ t ) isthe ground state value of α when V = 0, and θ representsthe dependence on the natural orbitals, c.f. (47). Let usverify (20) for the uniform case V = 0, for which theground state 1-matrix has θ = π/ α = α . At thesevalues h x = h y = h z = 0, which verifies the eigenvaluecollapse in this case.1 2. Iteration of the KS equations We demonstrate here the iteration of the KS equa-tions following the straightforward algorithm describedin Sec. II B. During the optimization of the orbitals theoccupation numbers (i.e. α ) are held fixed. Let us lookmore closely at each operation in the algorithm. In oper-ation (i), the KS Hamiltonian for step n + 1 is found byevaluating (52) at the 1-matrix γ ( n ) , i.e., at θ = θ n . Inoperation (ii), we find the eigenstates u i of ˆ h ( n +1) , whichwe parametrize in the form (47) with θ = θ n +1 . Theseeigenstates are taken as the natural orbitals φ ( n +1) i forstep n + 1. This implies setting each of the φ ( n +1) i equalto one of the u i . In the present case, the natural orbitalsare chosen such that θ n +1 is as close as possible to θ n . Inoperation (iii), γ ( n +1) is constructed from the φ ( n +1) i by(46). We may now condense these three operations intoa discrete iteration map on θ , i.e., a map θ n → θ n +1 . Itis defined bycos θ n +1 = sgn( A − A gs ) h z (cid:112) h x + h z (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) θ = θ n (53)for 0 < θ n < π and A > < α < π/ A gs isthe ground state value of A . An example of the iterationmap for t = 1, U = 1, V = 0, and A = A gs − . 02 is shownin Fig. 2. The solid [black] and dashed [red] curves arethe left and right hand sides of (53). The intersections ofthe two curves are fixed points of the iteration map. Theground state corresponds to the fixed point at θ = π/ FIG. 2: [color online] The iteration map (53) is shown for t = 1, U = 1, V = 0, and A = A gs − . 02. The solid [black]curve is the left hand side of (53). The dashed [red] curveis the right hand side. The dotted [blue] curve demonstratesthe first two iterations. The iteration map does not convergeto the ground state fixed point θ = π/ ternately drawing vertical lines from the solid curve tothe dashed curve and horizontal lines from the dashed curve to the solid curve. The dotted [blue] curve showsan example of the first two iterations beginning from aninitial guess θ . The next iterations θ and θ move far-ther away from the ground state, and the map does notconverge to the ground state fixed point θ = π/ V = 0, for which the ground state fixed point is θ = π/ 2. Linearization of (53) in terms of the variable m = ( n − n ) / A cos θ gives m n +1 ≈ sgn( A − A gs ) h z | h x |≈ − ξm n , (54)where ξ = (1 + sin α ) (cot α − cot α ) sin α cos α . (55)Suppose the occupation numbers are close to theirground state values, i.e., A = A gs + δA where δA isa small displacement. The leading approximation for ξ gives ξ ≈ − U ( U + B ) T B δA . (56)For any nonzero values of t and U , there is a thresh-old d > | δA | < d , | ξ | > 1. Therefore,the ground state is an unstable fixed point. In Sec. II B,the divergence of the iteration map was connected to thedivergence of the static KS response function. Let usverify (26) explicitly for the present case. As seen in(54), the linearized iteration map affects only the diag-onal elements of the 1-matrix, i.e., the density, which isdescribed by the variable m . Therefore, the relevant re-sponse functions are the density-density response for theKS system χ s = (cid:88) i (cid:88) j (cid:54) = i f i − f j (cid:15) i − (cid:15) j (cid:10) φ i (cid:12)(cid:12) ˆ m (cid:12)(cid:12) φ j (cid:11)(cid:10) φ j (cid:12)(cid:12) ˆ m (cid:12)(cid:12) φ i (cid:11) = 2 T U B δA (57)and the density-density response for the interacting sys-tem χ = (cid:88) k (cid:104) Ψ | ˆ m | Ψ k (cid:105) (cid:104) Ψ k | ˆ m | Ψ (cid:105) E − E k + c.c. = 2 U − BB ( B + U ) . (58)For the two-site Hubbard model, these response functionsare just constants. The KS response has a functional de-pendence on the 1-matrix. It diverges as the ground state2is approached, i.e., in the limit δA → 0. The linearizediteration map (26) is simply multiplication by a constant χ s χ − = − U ( U + B ) T B δA = ξ, (59)which agrees with the direct calculation (56).Of course, in actual calculations it is necessary to havea convergent iteration scheme. One possibility for ob-taining convergence is the level shifting method, whoseapplication in 1MFT was discussed in Sec. II C. In thelevel shifting method, one introduces artificial shifts ofthe KS eigenvalues in order to improve convergence. Ashift of the KS eigenvalue (cid:15) i by an amount ζ i is equiva-lent to adding a term ζ i (cid:12)(cid:12) φ i (cid:11)(cid:10) φ i (cid:12)(cid:12) to the KS Hamiltonian,where φ i is the orbital with eigenvalue (cid:15) i . The KS systemfor the two-site Hubbard model has two orbitals. As thedivergence of the iteration map is due to the degeneracyof the KS spectrum at the ground state, it seems sensi-ble to prevent degeneracy by introducing a separation 2 ζ between the levels. Thus, we add the following term tothe KS Hamiltonian at iteration step n − ζ (cid:12)(cid:12) φ a (cid:11)(cid:10) φ a (cid:12)(cid:12) + ζ (cid:12)(cid:12) φ b (cid:11)(cid:10) φ b (cid:12)(cid:12) = − ζ (sin θ n σ x + cos θ n σ z ) , (60)where φ a and φ b are evaluated at θ = θ n . An example ofthe effect of level shifting is shown in Fig. 3. Convergence FIG. 3: [color online] The iteration map for t = 1, U = 5, V = − . 5, and A = A gs − . 1. The solid [black] curve isthe left hand side of (53). The (dotted [blue], dash-dotted[green] , dashed [red]) curve is the right hand side with levelshift ζ = (0 , , ζ c ≈ . 07, which may be calculated with (35). is achieved when ζ exceeds a threshold, which may be cal-culated from the convergence criterion (35). The dashed[red] curve in Fig. 3 shows the iteration map with a levelshift value greater than the threshold. For the two-siteHubbard model, the criterion for convergence can be vi-sualized graphically as the condition that the magnitudeof the slope of the level shifted curve be less than theslope of the solid [black] curve at the fixed point.At each iteration step the KS system has an “instan-taneous” energy E s [ γ ] = tr (ˆ h [ γ ] γ ), which has, of course, no physical meaning when the KS system is not self-consistent. The KS energy is shown in Fig. 4 as a func-tion of the deviation δ(cid:126)γ = ( δγ x , δγ z ) of the 1-matrix(48) from the ground state 1-matrix. It is immediately FIG. 4: [color online] The KS energy for t = 1, U = 3 . V = 0 is shown as a function of the deviation ( δγ x , δγ z )from the ground state. The two surfaces represent the twobranches of the KS energy. The space curves show the energyas a function of θ for fixed occupation numbers. Optimizationof the orbitals corresponds to moving along one of these curvesto the stationary point. The ground state is the point of conicintersection at the origin. seen that the KS energy is not stationary at the groundstate 1-matrix, which is a cusp point where the energy E s changes linearly with respect to δγ . The KS energyis multivalued due to the choice implied in occupyingthe KS levels when the system is not self-consistent (seeSec. II D). The space curve in Fig. 4 shows the energy asa function of θ for fixed occupation numbers, i.e., for fixed A . The KS response is proportional to the inverse sepa-ration between the two branches of the space curve. Theseparation vanishes as the curve approaches the conicpoint, which is the origin of the divergent KS response. C. Solution by DFT The two-site Hubbard model with the local externalpotential chosen in (36) may be treated also with DFT.It is interesting to compare the DFT-KS scheme withthe 1MFT-KS scheme, especially with regard to theirconvergence behavior. The variational energy functionaland KS Hamiltonian may be constructed explicitly. Aninteresting result of the investigation is that the straight-forward iteration map is divergent when U > . t (for V = 0). We derive a general condition for the conver-gence of the DFT-KS equations.3 1. Energy functional The HK energy functional for a lattice is E [ n, v ] = (cid:88) i v ( i ) n i + F [ n ] , (61)where v ( i ) is the external potential at site i and F [ n ] is auniversal functional of the density (here, site occupancy)defined as F [ n ] = (cid:10) Ψ [ n ] (cid:12)(cid:12) ˆ T + ˆ U (cid:12)(cid:12) Ψ [ n ] (cid:11) , (62)where ˆ T is the kinetic energy operator and ˆ U is electron-electron interaction. In the following treatment of thetwo-site Hubbard model, we depart from standard prac-tice by enforcing the normalization condition n + n = N explicitly (i.e., through the parametrization), rather thanwith a Lagrange multiplier. Thus, we take as basic vari-able the single parameter m = ( n − n ) / V = v (1) − v (2). The functional F [ n ] in (62) is then justa function F ( m ), which may be constructed explicitly asfollows: i) a map m → Ψ is defined as the compositionof the maps m → V and V → Ψ and ii) the result-ing function Ψ ( m ) is used to evaluate (62). An explicitexpression for the map m → V can be found from the in-verse of (43). The second map v → Ψ was given in (40).The composition of these two maps gives the ground stateas a function of m , i.e., Ψ ( m ), with which the universalfunctional (62) may be evaluated. 2. KS Hamiltonian Following standard practice, the KS Hamiltonian takesover, unchanged, the kinetic energy operator from themany-body Hamiltonian. Thus, we consider the KSHamiltonian ˆ h = ˆ t + ˆ v s , (63)where ˆ t = − t ( c † c + c † c ) is the kinetic energy operatorand v s ( i ) is the KS potential at site i defined by v s ( i ) = ∂W∂n i , (64)where W [ n ] = E [ n, v ] − T s [ n ] contains the Hartree andexchange-correlation energy as well as the external poten-tial energy, and T s [ n ] is the kinetic energy of the KS sys-tem. We do not separate these contributions explicitly.The KS potential is spin independent because the groundstate density is spin unpolarized. Also, it is determinedonly to within an arbitrary additive constant, which wechoose such that v s (1) + v s (2) = 0. In the site basis, theKS Hamiltonian is a 2 × h = − tσ x + ( V s / σ z , where V s = v s (1) − v s (2). Thekinetic energy of the KS system is evaluated as T s = occ (cid:88) i f i (cid:10) φ i (cid:12)(cid:12) ˆ t (cid:12)(cid:12) φ i (cid:11) = 2 (cid:10) φ a (cid:12)(cid:12) ˆ t (cid:12)(cid:12) φ a (cid:11) = − t sin θ, (65)where φ a is the lowest energy eigenstate of (63) and istwice occupied (once by a spin up electron and once bya spin down electron.) It is parametrized as in (47) withtan θ = − t/V s . The density of the KS system is m s = occ (cid:88) i f i (cid:10) φ i (cid:12)(cid:12) ˆ m (cid:12)(cid:12) φ i (cid:11) = cos θ. (66)Thus, from (65) and (66) the kinetic energy T s is a knownfunction of m s . From (61), (64) and (65), the KS poten-tial is V s = ∂W∂m (cid:12)(cid:12)(cid:12)(cid:12) m = m s = ∂∂m ( E ( m, V ) − T s ( m )) (cid:12)(cid:12)(cid:12)(cid:12) m = m s = V + f ( m s ) − ∂T s ∂m (cid:12)(cid:12)(cid:12)(cid:12) m = m s (67)where V = v (1) − v (2) is the given external potential and f ( m s ) = ∂F∂m (cid:12)(cid:12)(cid:12)(cid:12) m = m s . (68)Eq. 67 is simply the familiar expression v s ( r ) = v ( r ) + v H ( r ) + v xc ( r ) with a different partitioning of the terms.It is seen that the terms f − ∂T s /∂m together correspondto the Hartree and exchange-correlation potentials. 3. Iteration of the KS equations Let us investigate the iteration of the KS equationsin the present context. The conventional iteration mapconsists of the following steps: i) the KS potential forstep n + 1 is determined from the density of step n using(67), i.e., V ( n +1) s = V s ( m ( n ) s ), ii) the eigenstates of ˆ h ( n +1) are found, and iii) the density of step n + 1 is calculatedwith (66).Consider step (i) in more detail. The KS potential isobtained from (67), V ( n +1) s = V + f ( m ( n ) s ) − ∂T s ∂m (cid:12)(cid:12)(cid:12)(cid:12) m = m ( n ) s . (69)The right hand side may be expressed differently by usingthe stationary conditions for the energy functional E [ n, v ]4and the KS energy E s = T s + (cid:80) i v s ( i ) n i . The stationarycondition ∂E/∂m = 0 applied to (61), gives f = − V ( m ),where V ( m ) is the external potential such that the inter-acting system has ground state density m . Similarly, thestationary condition applied to E s gives ∂T s /∂m = − V s .Substituting these relations in (69) yields V ( n +1) s = V − V ( m ( n ) s ) + V s ( m ( n ) s ) . (70)At self-consistency the V s terms cancel, and we obtain theexpected result V = V ( m gs ), where m gs is the groundstate density. For the present model, the ground statedensity could be found by solving V = V ( m ) as V ( m ) isknown exactly from (43). However, in general the groundstate must be found by iteration. Eq. 70 implies an it-eration map for the density, i.e., a map m ( n ) s → m ( n +1) s ,because V ( n +1) s determines m ( n +1) s . From (66) and thedefinition tan θ = − t/V s , we find the relationship V s = − t m s (cid:112) − m s . (71)The density may be iterated until self-consistency isreached. However, we encounter a technical difficultyfor the present model. In order to express explicitly theterm V ( m ( n ) s ) in (70), we must invert (43), which involvessolving a cubic equation. As the solutions are ratherunwieldy, we take here a different approach. We iterateinstead the external potential V ( m ). It may seem oddto iterate the external potential, which is given in thestatement of the problem. Nevertheless, the iterationmap for V provides an “image” of the iteration map for m s , by virtue of the HK theorem. Such an approachallows us to investigate certain features of the iterationmap, in particular its convergence behavior. In order toexpress (70) as an iteration map for V , we need to express V s as a function of V . In other words, we find the value of V s such that the KS system has density m s = m , where m is the density of the interacting system with V . Thecomposition of (71) and (43) yields the desired function˜ V s ( V ) = V s ( m ( V )) . (72)Using (72) in (70), we obtain the iteration map for theexternal potential˜ V s ( V ( n +1) ) = V − V ( n ) + ˜ V s ( V ( n ) ) , (73)which is expressed in implicit form.Examples of the iteration map for a uniform system( V = 0) are shown in Figs. 5 and 6, where the left andright hand sides of (73) are plotted. Suppose an initialvalue V (0) (cid:54) = 0 is chosen. For a system with V = 0, theground state has uniform density ( m = 0), but the initialdensity m (0) s associated with V (0) is not uniform. Uponiteration, we expect the KS system to relax to a uniformdensity, i.e., we expect the KS potential to be such as topush the system closer to uniform occupancy in the nextiteration. The solid [black] curves in Figs. 5 and 6 rep-resent the left hand side of (73), while the dashed [red] FIG. 5: The iteration map for the external potential V isshown for a weakly interacting system with U = t . The leftand right hand sides of (73) are shown as solid [black] anddashed [red] curves, respectively.FIG. 6: The iteration map for the external potential V isshown for a strongly interacting system with U = 4 t . Theleft and right hand sides of (73) are shown as solid [black]and dashed [red] curves, respectively. curves represent the right hand side. The iteration mapmay be demonstrated graphically by alternately drawingvertical lines from the solid curve to the dashed curve andhorizontal lines from the dashed curve to the solid curve.The map displays “charge oscillation.” The ground stateis a stable fixed point if the magnitude of the slope ofthe dashed curve at the origin is less than the slope of5the solid curve at the origin. For weakly interacting sys-tems the iteration map is convergent, while for stronglyinteracting systems it is nonconvergent. The thresholdfor convergence is U ≈ . t . 4. Linearization of the KS equations The nature of the fixed point and the origin of diver-ent behavior are revealed by linearization of the iterationmap. We linearize the map by expanding both sides of(70) with respect to δm s = m s − m gs , where m gs is theground state density. We find χ − s δm ( n +1) s ≈ χ − s δm ( n ) s − χ − δm ( n ) s δm ( n +1) ≈ χ s (cid:0) χ − s − χ − (cid:1) δm ( n ) s , (74)where χ s and χ are the density-density response func-tions defined in (57) and (58). The threshold for conver-gent behavior is (cid:12)(cid:12) − χ s χ − (cid:12)(cid:12) ≤ , (75)or equivalently, χ s χ − ≤ 2. Note the change from 1 for1MFT to 2 for DFT on the right hand side, cf. (26).Consider the case V = 0, which has uniform density inthe ground state ( m gs = 0). Using the (57) and (58) in(75), gives the threshold conditioncos( π/ − α / 2) = 4 (sin( π/ − α / . (76)The threshold is cos( π/ − α / ≈ . U ≈ . t . Let us consider the limit U → ∞ .The leading behavior of the KS response is independentof U , χ s ∼ T , (77)while the response of the interacting system vanishes as χ ∼ T U . (78)For sufficiently large U , the threshold (75) is crossed anddivergent behavior results. In DFT, as also in 1MFT,the source of divergent behavior is a KS response that istoo large in relation to the exact response. In 1MFT theimbalance results from a divergent KS response, whereasin DFT the KS response generally remains finite but theresponse of the interacting system becomes too small as U increases.In standard DFT (with continuous n ( r )), the analogof the linearized iteration map (74) may be written n ( n +1) ( r ) ≈ (cid:90) dr (cid:48) dr (cid:48)(cid:48) χ s ( r, r (cid:48) ) ( v c ( r (cid:48) , r (cid:48)(cid:48) ) + f xc ( r (cid:48) , r (cid:48)(cid:48) )) × n ( n ) ( r (cid:48)(cid:48) ) , (79)where n ( n ) ( r ) is the density of iteration step n , v c is thekernel of the Coulomb interaction, and f xc = δv xc /δn is the exchange-correlation kernel. The necessary andsufficient condition for convergence of the KS equationsis that all eigenvalues of the operatorˆ χ s (cid:0) ˆ v c + ˆ f xc (cid:1) (80)have modulus less than 1. IV. CONCLUSIONS The status of the KS system in 1MFT has been uncer-tain. Although Gilbert derived effective single-particleequations from the stationary conditions for the energyfunctional, the degeneracy of essentially all of the result-ing orbitals was thought to be paradoxical. We havepresented an alternative derivation of the KS equationsin which the degeneracy is lifted by constraining the oc-cupation numbers. Such a KS scheme is well-behaved inthe neighborhood of the ground state occupation num-bers. Therefore, the correct natural orbitals are obtainedin the limit that the ground state is approached. We haveconstructed explicitly the 1MFT-KS system for a simpletwo-site Hubbard model. While we find no paradoxi-cal results, the KS system has many striking features, inparticular the collapse of eigenvalues at the ground state.Although the KS eigenvalues do not have a physical in-terpretation as in DFT, the orbitals, which are callednatural orbitals, play an important role in the contextof configuration interaction, i.e., the expansion of thefull wavefunction as a sum of Slater determinants. Thismay be important in the search for approximate energyfunctionals.Beyond the question of the existence of the KS systemin 1MFT, there is the issue of its practicality. The KSsystem has been extremely useful in DFT calculations.Due to the implicit 1-matrix dependence of the single-particle potential, the KS equations are nonlinear. Suchequations are generally solved by iteration. As in DFT,there is a “straightforward” procedure for iteration. Incontrast to DFT, the “straightforward” procedure is al-ways divergent, in the sense that the ground state is anunstable fixed point. We have demonstrated the insta-bility of the ground state by linearization of the itera-tion map. The source of the instability is the divergenceof the KS static response function at the ground state,which in turn, is due to the degeneracy of the KS spec-trum. Degeneracy-driven instability is reminiscent of theJahn-Teller effect, and the connection is strengthened ifwe regard the implicit 1-matrix dependence of the KSHamiltonian as analogous to the parametric dependenceof the Born-Oppenheimer Hamiltonian on nuclear coordi-nates. In both cases, the energy spectrum splits linearlywith respect to displacement from the degeneracy point.Thus, the energy may always be lowered by displacement.For the 1MFT-KS system, this means that the KS energy tr (ˆ h ˆ γ ) may always be lowered by displacement from theground state, leading to an instability of the iterationprocedure. However, this is a fictitious energy and the6HK energy functional E v is of course always minimum atthe ground state. Acknowledgments We gratefully acknowledge helpful discussions withWei Ku. APPENDIX A: GROUND STATE NOTDETERMINED BY THE DENSITY We give here a simple example which shows that thedensity alone does not always uniquely determine theground state wavefunction if the external potential isnonlocal. Our example is the two-site Hubbard model,which was solved for the case of a local external potentialin Sec. III A. The Hamiltonian ˆ H = ˆ T + ˆ V + ˆ U is givenin (36). In such a lattice model, the hopping parameters t jk are real numbers that represent the kinetic energy. A“magnetic field” can be introduced by giving t jk a phase,i.e., by the transformation t jk → t jk exp (cid:0) i (cid:80) kn = j A ( n ) (cid:1) ,where A ( n ) is the “vector potential” at site n and the sumruns over a string of sites from site j to site k . For thetwo-site model this is just the transformation t → te iτ .We see that this magnetic field appears in the Hamil-tonian in exactly the same manner as a nonlocal exter-nal potential, such as v c † c + h.c. , because it modifiesthe nonlocal hopping terms. We can generate the abovephase transformation by the U (1) rotations c → c e iτ/ and c → c e − iτ/ . The eigenstates of the transformed Hamiltonian are readily generated from the eigenstatesof the original Hamiltonian by applying the same trans-formation. For example, without the magnetic field, theground state to first order in small V isΨ = Φ + (cid:10) Φ (cid:12)(cid:12) ˆ V (cid:12)(cid:12) Φ (cid:11) E − E Φ , (A1)where Φ i are given in (37). When the magnetic field isturned on, the Φ i change, e.g.,Φ → √ ye − iτ xxye iτ (A2)in the site basis (cid:8) c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) , c † ↑ c † ↓ | (cid:105) (cid:9) . Accordingly, the ground state acquiresa nontrivial dependence on the magnetic field ( τ -dependence). At the same time, the ground state 1-matrix is transformed to γ = γ + A (cid:18) cos θ sin θ e − iτ sin θ e iτ − cos θ (cid:19) , (A3)where A and θ are the ground state values (for τ = 0)of the variables defined in (47) and (48). The density isgiven by the diagonal elements, which are unaffected bythe transformation. Only the off-diagonal (nonlocal) ele-ments are sensitive to the magnetic field. Therefore, the1-matrix rather than the density is required to uniquelyspecify the ground state. ∗ Electronic address: [email protected] P. Hohenberg and W. Kohn, Phys. Rev. , B864 (1964). T. L. Gilbert, Phys. Rev. B , 2111 (1975). U. von Barth and L. Hedin, J. Phys C , 1629 (1972). O. Gunnarsson and B. I. Lundqvist, Phys. Rev. B , 4274(1976). G. Vignale and M. Rasolt, Phys. Rev. Lett. , 2360(1987). G. Vignale and M. Rasolt, Phys. Rev. B , 10685 (1988). O. Gritsenko, K. Pernal, and E. J. Baerends, J. Chem.Phys. , 204102 (2005). N. N. Lathiotakis, N. Helbig, and E. K. U. Gross, Phys.Rev. B , 195120 (2007). M. A. Buijse and E. J. Baerends, Mol. Phys. , 401(2002). N. Helbig, N. N. Latiotakis, M. Albrecht, and E. K. U.Gross, Euro. Phys. Lett. , 67003 (2007). A. M. K. M¨uller, Physics Letters , 446 (1984). W. Kohn and L. J. Sham, Phys. Rev. , A1133 (1965). S. M. Valone, J. Chem. Phys. , 1344 (1980). T. T. Nguyen-Dang, E. V. Ludena, and Y. Tal, J. Mol.Struct. , 247 (1985). A. Schindlmayr and R. W. Godby, Phys. Rev. B , 10427 (1995). N. Helbig, N. N. Latiotakis, M. Albrecht, and E. K. U.Gross, Phys. Rev. A , 030501(R) (2005). P. O. L¨owdin, Phys. Rev. , 1474 (1955). K. Pernal, Phys. Rev. Lett. , 233002 (2005). V. R. Saunders and I. H. Hillier, Int. J. Quant. Chem. ,699 (1973). F. Aryasetiawan, O. Gunnarsson, and A. Rubio, Europhys.Lett. (2002). R. M. Dreizler and E. K. U. Gross, Density functional the-ory (Springer-Verlag, Berlin, 1990). M. Levy, Proc. Natl. Acad. Sci. USA , 6062 (1979). A. J. Coleman, Rev. Mod. Phys. , 668 (1963). J. Cioslowski and K. Pernal, Chem. Phys. Lett. , 188(2006). H. Eschrig and W. E. Pickett, Solid State Commun. ,123 (2001). K. Capelle and G. Vignale, Phys. Rev. Lett. , 5546(2001). O. Gunnarsson and K. Sch¨onhammer, Phys. Rev. Lett. ,1968 (1986). G. Xianlong, et al. , Phys. Rev. B , 165120 (2006). H. Shull and P. O. L¨owdin, J. Chem. Phys. , 617 (1959). W. Kutzelnigg, Theor. Chim. Acta , 327 (1963). I. M. Gelfand and S. V. Fomin, Calculus of variations (Prentice-Hall, Englewood Cliffs, 1963). C. A. Ullrich and W. Kohn, Phys. Rev. Lett. , 093001(2001). W. Kohn, Phys. Rev. Lett. , 1596 (1983). The following theorem implies that the functional E v [ γ ]is stationary, i.e., δE v = 0 with respect to all variationsin the VR space. Theorem – A necessary condition for thedifferentiable functional J [ y ] to have an extremum (mini-mum) for y = y is that its first variation vanish for y = y ,i.e., that δJ [ h ] = 0 for y = y and all admissible variations h . Thus, the functional E v [ γ ] must be stationary because itis differentiable in the VR space and the quantum mechani-cal variation principle (Rayleigh-Ritz variational principle)implies that it is minimum at the ground state 1-matrix γ = γ gs . The differentiability of E v [ γ ] follows from the dif- ferentiability of the mapping γ → Ψ in the VR space, ifthe ground state Ψ is nondegenerate. If the DFT-KS system is degenerate, the occupation num-bers of the degenerate KS orbitals are not determined bythe aufbau principle. In this case, the KS system adoptsan ensemble state, and the occupation numbers of the de-generate orbitals are chosen such that the KS system isself-consistent and reproduces the density of the interact-ing system. We do not demonstrate this statement for χ as an operatoron an infinite dimensional space. However, it is valid in afinite dimensional basis (see Ref. 33 for a proof in DFT).37