The Kohn-Sham system in one-matrix functional theory

# The Kohn-Sham system in one-matrix functional theory

## Abstract

A system of electrons in a local or nonlocal external potential can be studied with 1-matrix functional 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 12, 2111 (1975)] effective single-particle equations analogous to the Kohn-Sham (KS) equations in DFT. The self-consistent solution of these 1MFT-KS equations reproduces not only the density of the original electron system but also its 1-matrix. While in DFT it is usually possible to reproduce the density using KS orbitals with integer (0 or 1) occupancy, in 1MFT reproducing the 1-matrix requires in general fractional occupancies. The variational principle implies that the KS eigenvalues of all fractionally occupied orbitals must collapse at self-consistency to a single level, equal to the chemical potential. We show that as a consequence of the degeneracy the 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 introduce an alternative derivation of the 1MFT-KS equations that allows control of the eigenvalue collapse by constraining the occupancies. As an explicit example, we apply the 1MFT-KS scheme to calculate the ground state 1-matrix of an exactly solvable two-site Hubbard model.

###### pacs:
71.15.Mb,31.15.xr

## I Introduction

Density functional theory (DFT) benefits from operating with the electron density, which as a function of just three coordinates is much easier to work with than the full many-body wavefunction. According to the Hohenberg-Kohn (HK) theorem, Hohenberg and Kohn (1964) the density of an electron system in a local external potential may be found by minimizing a universal energy functional , whose basic variable is the density. Remarkably, the density uniquely determines the ground state wavefunction (if it is nondegenerate), i.e., there can be only one ground state wavefunction yielding a given density, no matter what is. However, if the external potential is nonlocal, then the density alone is generally not sufficient to uniquely determine the ground state (see Appendix A for a simple example). Gilbert Gilbert (1975) extended the HK theorem to systems with nonlocal and spin dependent external potential , where . It was proved that i) the ground state wavefunction is uniquely determined by the ground state 1-matrix (one-particle reduced density matrix) and ii) there is a universal energy functional of the 1-matrix, which attains its minimum at the ground state 1-matrix. The 1-matrix is defined as

 γ(x,x′)=N∫dx2…dxNρ(x,x2,…xN;x′,x2,…xN), (1)

where and is the full -electron density matrix with ensemble weights such that . An external potential may be nonlocal with respect to the space coordinates and/or the spin coordinates. For example, pseudopotentials are nonlocal in space, and Zeeman coupling , where is the vector of Pauli matrices, is nonlocal in spin space. The coupling of electron motion and an external vector potential, , may also be treated as a nonlocal potential because is a differential operator. It is rather intuitive that for such external potentials, which couple to the system in more complex ways than the local potential , it is necessary — in order to permit statements analogous to the HK theorem — to refine the basic variable accordingly. Hence spin-DFT, von Barth and Hedin (1972); Gunnarsson and Lundqvist (1976) whose basic variables are the density and the magnetization density, applies to systems with Zeeman coupling. Current-DFT, Vignale and Rasolt (1987, 1988) whose basic variables are the density and the paramagnetic current density, has the scope to treat systems in which the current is coupled to an external magnetic field. Generally, if one considers an external potential that is nonlocal in space and spin, the necessary basic variable is the one-matrix,Gilbert (1975) which contains all of the single-particle information of the system, including the density, magnetization density and paramagnetic current density.

The DFT-type approach that takes the 1-matrix as basic variable will be referred to here as 1-matrix functional theory (1MFT). As in DFT, an exact and explicit energy functional is generally unknown. An important difference between 1MFT and DFT is that the kinetic energy is a simple linear functional of the 1-matrix, while it is not a known functional of the density. Thus, in 1MFT the only part of the energy not known explicitly is the electron-electron interaction energy . Several approximate 1-matrix energy functionals have been proposed and tested recently (see Refs. Gritsenko et al., 2005 and Lathiotakis et al., 2007 and references therein.) Notably, the so-called BBC approximations,Gritsenko et al. (2005) which are modifications of the Buijse-Baerends functional,Buijse and Baerends (2002) have given fairly accurate results for the potential energy curves of diatomic moleculesGritsenko et al. (2005) and the momentum distribution and correlation energy of the homogeneous electron gas.Lathiotakis et al. (2007) In Ref. Lathiotakis et al., 2007, a density dependent fitting parameter was introduced into the BBC1 functional such that the resulting functional yields the correct correlation energy of the homogeneous electron gas at all values of density. There is also the prospect of using 1MFT to obtain accurate estimates for the band gaps of non-highly correlated insulators.Helbig et al. (2007) Many of the approximate functionals that have been proposed are similar to an early approximation by Müller.Müller (1984)

Actual calculations in 1MFT are more difficult than in DFT. The energy functional must be minimized in a space of higher dimension because the 1-matrix is a more complex quantity than the density. In the calculations cited above, the energy has been minimized directly by standard methods, e.g., the conjugate gradient method. In DFT the energy is generally not minimized by such direct methods. Instead, the Kohn-Sham (KS) schemeKohn and Sham (1965) provides an efficient way to find the ground state density. In this scheme, one introduces an auxiliary system of noninteracting electrons, called the KS system, which experiences an effective local potential . This effective potential is a functional of the density such that the self-consistent ground state of the KS system reproduces the ground state density of the interacting system. It is interesting to ask whether there is also a KS scheme in 1MFT. The question may be stated as follows: does there exist a 1-matrix dependent effective potential such that, at self-consistency, a system of noninteracting electrons experiencing this potential reproduces the exact ground state 1-matrix of the interacting system? Although Gilbert derived such an effective potential,Gilbert (1975) the implications were thought to be “paradoxical” because the KS system was found to have a high (probably infinite) degree of degeneracy. Evidently, the KS eigenvalues in 1MFT do not have the meaning of approximate single-particle energy levels, in contrast to DFT and other self-consistent-field theories, where the eigenvalues may often be interpreted as the negative of ionization energies, owing to Koopmans’ theorem. The status of the KS scheme in 1MFT appears to have remained unresolved, Valone (1980); Nguyen-Dang et al. (1985) and recently it has been argued that the KS scheme does not exist in 1MFT. Schindlmayr and Godby (1995); Helbig et al. (2005, 2007); Lathiotakis et al. (2007) Gilbert derived the KS equations from the stationary principle for the energy. The KS potential was found to be

 vs(x,x′)=v(x,x′)+δW/δγ(x′,x). (2)

In this article, we propose an alternative derivation of the KS equations, which, in our view, gives insight into the nature of the “paradoxical” degeneracy of the KS system.

One-matrix energy functionals are often expressed in terms of the so-called natural orbitals and occupation numbers. Löwdin (1955) This makes them similar to “orbital dependent” functionals in DFT. The natural orbitals are the eigenfunctions of the 1-matrix, and the occupation numbers are the corresponding eigenvalues.Löwdin (1955) These quantities play a central role in 1MFT. Recently, it was shown that when a given energy functional is expressed in terms of the natural orbitals and occupation numbers, the KS potential can be found by using a chain rule to evaluate the functional derivative in Eq. 2.Pernal (2005)

Although the concept of the KS system can indeed be extended to 1MFT, it has in this setting some very unusual properties. In particular, the KS orbitals must be fractionally occupied, for otherwise the KS system could not reproduce the 1-matrix of the interacting system, which always has noninteger eigenvalues (occupation numbers). This is different from the situation in DFT, where it is usually possible to reproduce the density using only integer (0 or 1) occupation numbers, or in any case, only a finite number of fractionally occupied states. Due to the necessity of fractional occupation numbers, the 1MFT-KS system cannot be described by a single Slater determinant. However, we find that it can be described by an ensemble of Slater determinants, i.e., a mixed state. In order that the variational principle is not violated, all the states that comprise the ensemble must be degenerate. This implies that the eigenvalues of all fractionally occupied orbitals collapse to a single level, equal to the chemical potential. The degeneracy has important consequences for the solution of the KS equations by iteration. We prove that the iteration of the KS equations is intrinsically divergent because the KS system has a divergent response function at the ground state. Fortunately, convergence can always be obtained with the level shifting method.Saunders and Hillier (1973) To illustrate explicitly the unique properties of the 1MFT-KS system, we apply it to a simple Hubbard model with two sites. The model describes approximately systems which have two localized orbitals with a strong on-site interaction, e.g., the hydrogen molecule with large internuclear separation. Aryasetiawan et al. (2002) The Schrödinger equation for this model is exactly solvable, and we find that the KS equations in 1MFT and in DFT can be derived analytically. It is interesting to compare 1MFT and DFT in this context. We demonstrate that divergent behavior will appear also in DFT when the operator , where and are the density response functions of the interacting and KS systems, respectively, has any eigenvalue with modulus greater than 1. In this expression the null space of is assumed to be excluded.

This article is organized as follows. In Sec. II, we derive the KS equations in 1MFT and discuss how to solve them self-consistently by iteration. In Sec. III, we compare three approaches to ground state quantum mechanics — direct solution of the Schrödinger equation, 1MFT and DFT — by using them to solve the two-site Hubbard model.

## Ii Kohn-Sham system in 1MFT

It is not obvious that a KS-type scheme exists in 1MFT for the following reason. Recall that in DFT the KS system consists of noninteracting particles and reproduces the density of the interacting system. The density of the KS system, if it is nondegenerate, is the sum of contributions of the lowest energy occupied orbitals

 n(→r)=occ∑i|ϕi(→r)|2. (3)

On the other hand, in 1MFT the KS system should reproduce the 1-matrix of the interacting system. The eigenfunctions of the 1-matrix are the so-called natural orbitals, and the eigenvalues are the corresponding occupation numbers. Löwdin (1955) Occupying the lowest energy orbitals in analogy to (3), one obtains

 γ(x,x′)=occ∑iϕi(x)ϕ∗i(x′). (4)

Such an expression, in which the orbitals have only integer (0 or 1) occupation, cannot reproduce the 1-matrix of an interacting system because the orbitals of an interacting system have generally fractional occupation (see the discussion in the following section.) The difference between the 1-matrix in (4) and the 1-matrix of an interacting system is clearly demonstrated by the so-called idempotency property. The 1-matrix in (4) is idempotent, i.e., , while the 1-matrix of an interacting system is never idempotent. However, if the KS system is degenerate and its ground state is an ensemble state, the 1-matrix becomes

 γ(x,x′)=∑ifiϕi(x)ϕ∗i(x′) (5)

with fractional occupation numbers . The -particle ground state density matrix of the KS system is , where the are Slater determinants each formed from degenerate KS orbitals. The occupation numbers are related to the ensemble weights by

 fi=∑jwjΘji, (6)

where equals if is one of the orbitals in the determinant and 0 otherwise. Dreizler and Gross (1990)

### ii.1 Derivation of the 1MFT Kohn-Sham equations

In this section, we discuss Gilbert’s derivationGilbert (1975) of the KS equations in 1MFT and propose an alternative derivation. We begin by reviewing the definition of the universal 1-matrix energy functional .

One-matrix functional theory describes the ground state of a system of electrons with the Hamiltonian , where is the kinetic energy operator, is the local or nonlocal external potential operator, and is the electron-electron interaction (in atomic units ). The ground state 1-matrix and ground state energy can be found by minimizing the functional

 Ev[γ]=Tr((^t+^v)γ)+W[γ], (7)

where

 W[γ]=⟨Ψ0∣∣^W∣∣Ψ0⟩. (8)

By extending the HK theorem, Gilbert provedGilbert (1975) that a nondegenerate ground state wavefunction, , is uniquely determined by the ground state 1-matrix, i.e., is a functional of . For this reason the interaction energy, as defined in (8), is a functional of . It is apparent that (8) defines only for that are the ground state 1-matrices of some system (with Hamiltonian ). In this article, a 1-matrix is said to be -representable (VR) if it is the ground state 1-matrix of some system with local or nonlocal external potential. Gilbert remarked (see the discussion between equations (2.24) and (2.25) in Ref. Gilbert, 1975) that, in principle, the domain of can be extended to the space of ensemble -representable (ENR) 1-matrices. A 1-matrix is said to be ENR if it can be constructed via (1) from some -particle density matrix , which is not required to be a ground state ensemble. One possible extension to the ENR space is provided by the so-called constrained search functionalLevy (1979); Valone (1980)

 W[γ]=min^ρ→γTr(^W^ρ), (9)

where the interaction energy is minimized in the space of -particle density matrices that yield via (1). The definition (9) is a natural extension to the ENR space because when it is adopted (7) may be expressed as

 Ev[γ]=min^ρ→γTr(^H^ρ). (10)

This is a variational functional which attains its minimum at the ground state 1-matrix, as seen from

 minγEv[γ]=min^ρTr(^H^ρ)=E0, (11)

where is the ground state energy. The extension to the ENR domain is significant, especially for applications of the variational principle, because the conditions a 1-matrix must satisfy to be ENR are known and simple to impose on a trial 1-matrix, while the conditions for -representability are unknown in general. The necessary and sufficient conditions Coleman (1963) for a 1-matrix to be ENR are i) must be Hermitean, ii) , and iii) all eigenvalues of (occupation numbers) must lie in the interval . The third condition is a consequence of the Pauli exclusion principle.

The 1MFT-KS equations were derivedGilbert (1975) from the stationary conditions for the energy with respect to arbitrary independent variations of the natural orbitals and angle variables chosen to parametrize the occupation numbers according to (). For the purpose of describing a variation in the ENR space, this set of variables, namely , is redundant. An arbitrary set of such variations may or may not correspond to an ENR variation, and when it does the variations will not be linearly independent. This causes no difficulty, of course, because it is always possible to formulate stationary conditions in a space whose dimension is higher than necessary, provided the appropriate constraints are enforced with Lagrange multipliers. Accordingly, the Lagrange multiplier terms which maintain the orthogonality of the orbitals and the term which maintains the total particle number were introduced. The KS equations were found to be

 (^t+^vs)|ϕi⟩=ϵi|ϕi⟩, (12)

where the kernel of the effective potential is , if the functional derivative exists. The stationary conditions imply that all fractionally occupied KS orbitals have the same eigenvalue . Gilbert described this result as “paradoxical” because in interacting systems essentially all orbitals are fractionally occupied.

The above stationary conditions assume to be stationary with respect to variations in the ENR space (except variations of occupation numbers equal to exactly 0 or 1 in the ground state, which are excluded by the parametrization). However, it is not known, in general, whether is stationary in the ENR space; the minimum property (11) ensures only that it is variational. Recall that the ENR space consists of all that can be constructed from an ensemble, and the energy of an ensemble is not stationary with respect to variations of the many-body density matrix . Therefore, the stationary conditions applied in the ENR space may be too strong. On the other hand, the quantum mechanical variational principle guarantees that is stationary with respect to variations in the VR space,1 but it is not known how to determine whether a given is VR. Hence, it is not known how to constrain the variations of to the VR space. Nevertheless, it may be that in some systems the entire neighborhood of in the ENR space is also VR. In such cases, is stationary in the ENR space, and the stationary conditions applied in Ref. Gilbert, 1975 are satisfied at the ground state.

We find it helpful to construct an alternative derivation of the KS equations. Consider the energy functional

 Gv[γ]=Ev[γ]−∑jϵj(fj−qj), (13)

where are Lagrange multipliers that constrain the occupation numbers of the natural orbitals to chosen values , which satisfy and . These Lagrange multipliers allow us to investigate the degeneracy of the KS eigenvalues, which leads to the “paradox” described by Gilbert. We have omitted the Lagrange multipliers used in Gilbert’s derivation. They are not necessary in our derivation because we formulate the stationary conditions with respect to variations of the 1-matrix instead of the orbitals. We adopt the definition (9) for , so the domain of is the ENR space. A variation will be said to be admissible if is ENR. For convenience we assume that the static response function for the interacting system under consideration has no null vectors apart from the null vectors associated with a) a constant shift of the potential (which is a null vector also in DFT) and b) integer occupied orbitals, i.e., orbitals with occupation numbers exactly 0 or 1. If has additional null vectors, the following derivation must be modified; the necessary modifications are discussed below. Granting the above assumption, is guaranteed to be stationary, and the KS equations can be derived from the stationary condition with respect to an arbitrary admissible variation of . The first variation of is

 δGv = Tr((^t+^v)δ^γ)+Tr(^wδ^γ)−∑jϵjδfj (14) = ∑ij⟨ϕi∣∣(^t+^v+^w)∣∣ϕj⟩⟨ϕj∣∣δ^γ∣∣ϕi⟩ −∑ijϵj⟨ϕi|ϕj⟩⟨ϕj∣∣δ^γ∣∣ϕi⟩ = ∑ij(hij−ϵjδij)δγji,

where the variation of the 1-matrix is expressed as in the basis of the ground state natural orbitals, and the relation defines a single-particle operator . In (14), we have also introduced the definition of the Hermitean operator

 ^h=^t+^v+^w, (15)

which will be seen to be the KS Hamiltonian. If the last line of (14) is to be zero for an arbitrary Hermitian matrix , then we must have for all and . ENR condition (i) has been maintained explicitly by requiring the variation to be Hermitian. ENR conditions (ii) and (iii) do not impose any constraint on the space of admissible variations as they are maintained by the Lagrange multipliers . The matrix elements are functionals of the 1-matrix, and the 1-matrix that satisfies the stationary conditions can be found by solving self-consistently the single-particle equations

 ^h∣∣ϕi⟩=ϵi∣∣ϕi⟩ (16)

together with (5). These are the KS equations in 1MFT. If they are solved self-consistently with the occupation numbers fixed to the values , they give the orbitals which minimize subject to the constraints . The KS potential is . The term is the effective contribution of the electron-electron interaction to the KS potential. In coordinate space, its kernel is

 w(x,x′) = ⟨x∣∣^w∣∣x′⟩ (17) = δWδγ(x′,x),

which recovers Gilbert’s result. The kernel of the KS Hamiltonian may be written in the familiar form

 h(x,x′) = δ(x−x′)(−12∇2r)+v(x,x′) (18) +δ(x−x′)vH(x)+vxc(x,x′),

where is the external potential and has been divided into the Hartree and exchange-correlation potentials. In 1MFT, the exchange-correlation potential is nonlocal.

The 1MFT-KS scheme optimizes the orbitals for a chosen set of occupation numbers, but it does not itself provide a rule for choosing the occupation numbers. On this point, it is different from the DFT-KS scheme, where the occupation numbers are usually uniquely determined by the aufbau principle ( Fermi statistics).2 In 1MFT, the KS equations have a self-consistent solution for any chosen set of occupation numbers that satisfy and . The Lagrange multipliers , which are seen to be the KS eigenvalues, adopt values such that the minimum of occurs for a 1-matrix whose occupation numbers are precisely the set . Therefore, the unconstrained minimum of coincides with the minimum of subject to the constraints . To find the ground state occupation numbers, , one can search for the minimum of the function , where the minimization of can be performed by the KS scheme.

What can be said about the KS eigenvalues ? As the minimum of is a stationary point, we have

 0=∂Gv∂fj∣∣∣γmin=∂Ev∂fj∣∣∣γmin−ϵj({qi}) (19)

for all . is a variational functional in the ENR space. It attains its minimum at the ground state 1-matrix , where must vanish for all fractionally occupied () orbitals, for otherwise the energy could be lowered. When for all , , and (19) implies for all fractionally occupied orbitals. Thus, we find that the KS eigenvalues of all orbitals that are fractionally occupied in the ground state must collapse to a single level when the chosen set of occupation numbers approach their ground state values, i.e., as . In Gilbert’s derivation of the 1MFT-KS equations, all fractionally occupied KS orbitals were found to have the eigenvalue , where is the chemical potential. As we have not introduced the chemical potential in our derivation (we consider a system with a fixed number of electrons), the eigenvalues collapse to instead of . The above arguments do not apply to orbitals with occupation numbers exactly or because these values lie on the boundary of the allowed interval specified by ENR condition (iii). All that can be concluded from the fact that has a minimum in the ENR space is for orbitals with and for orbitals with . States with occupation numbers exactly or have been called “pinned states.”Helbig et al. (2007); Lathiotakis et al. (2007) Instances of such states in real systems have been reported,Cioslowski and Pernal (2006) though their occurrence is generally considered to be exceptional.Helbig et al. (2007); Lathiotakis et al. (2007)

Due to the collapse of the eigenvalues at the ground state, the KS Hamiltonian becomes the null operator

 ^h[γgs]=^0 (20)

in the subspace of fractionally occupied orbitals. This is of course analogous to the familiar condition for the extremum of a function . Gilbert describedGilbert (1975) a similar result (with replacing ) as “paradoxical,” a statement that has been repeated.Valone (1980); Nguyen-Dang et al. (1985) The problem with (20) is that while we expect the KS Hamiltonian to define the natural orbitals, any state is an eigenstate of the null operator. However, the KS Hamiltonian is a functional of the 1-matrix, and when the occupation numbers are perturbed from their ground state values, the degeneracy is lifted and the KS Hamiltonian does define unique orbitals. In the KS scheme outlined above, this corresponds to the optimization of the orbitals with occupation numbers fixed to values , perturbed from the ground states values. In the limit that the occupation numbers approach their ground state values, the optimal orbitals approach the ground state natural orbitals. The degenerate eigenvalues generally split linearly with respect to perturbations away from the ground state. In particular,

 ∂ϵi∂qj = ∫dydy′⟨ϕi|δ^hδγ(y,y′)∂γ(y,y′)∂qj|ϕi⟩ (21) = −∫dxdx′dydy′ϕ∗i(x)ϕ∗j(y′)χ−1(xx′,yy′) ×ϕi(x′)ϕj(y).

Here, is the static response function defined as

 χ(x,x′;y,y′)=δγ(x,x′)δv(y,y′). (22)

The relation used in (21) is derived in the following section. If has a null space, its inverse is defined only on a restricted space. For example, (21) does not apply to pinned states as there is a null vector associated with each pinned state (see below).

Our derivation of the 1MFT-KS equations in fact assumes that the static response function of the interacting system has no null vectors except for those associated with pinned states and a constant shift of the potential. We now show that this guarantees is stationary. If the interacting system has any other null vectors we can no longer be certain that is stationary and the derivation should be modified as described below. We have remarked already (see Ref. 34) that is stationary in the VR space, i.e., it satisfies the stationary condition with respect to an arbitrary variation of the 1-matrix in the VR space. However, our derivation of the KS equations requires to be stationary in the ENR space. As the VR space is a subspace of the ENR space, this is a stronger condition. The assumption that has no null vectors (apart from those associated with pinned states) is equivalent to assuming that any ENR variation (apart from variations of the pinned occupation numbers) is also a VR variation. For if has no null vectors, then it is invertible and any ENR variation can be induced by the perturbation .3 Hence, with the above assumption, is guaranteed to be stationary with respect to any ENR variation. The above arguments do not apply to variations of pinned occupation numbers because there are null vectors associated with such variations; nevertheless, is stationary with respect to such variations as this is maintained by the Lagrange multipliers . It is now clear how to modify the derivation of the KS equations when has additional null vectors. By introducing new Lagrange multipliers, the additional null vectors can be treated in analogy with the pinned states. For example, suppose has one additional null vector , where is hermitian and . The new energy functional will be stationary with respect to an arbitrary variation in the ENR space. Here, the Lagrange multiplier enforces the constraint , where is the component of corresponding to . The stationary condition leads to the set of equations in the basis of natural orbitals. The 1-matrix that satisfies these equations can be found by solving self-consistently the eigenvalue equation together with , where is the unitary matrix that diagonalizes the matrix in the basis of natural orbitals, i.e., is diagonal. The energy of the self-consistent solution defines a function , whose minimum with respect to and is the ground state energy. It may not be known in advance whether the response function of a given interacting system will have null vectors. Therefore, it is helpful to understand how null vectors occur.

Null vectors of are connected with the so-called nonuniqueness problemvon Barth and Hedin (1972); Eschrig and Pickett (2001); Capelle and Vignale (2001) in various extensions of DFT. A system with the ground state is said to have a nonuniqueness problem if there is more than one external potential for which is the ground state. The Schrödinger equation defines a unique map from the external potential to the ground state wavefunction (if it is nondegenerate), but when there is more than one external potential yielding the same ground state wavefunction, the map cannot be inverted. In 1MFT the generality of the external potential (nonlocal in space and spin coordinates) allows greater scope for nonuniqueness than in the other extensions of DFT. Of course, every degree of nonuniqueness is a null vector of because if does not change the ground state wavefunction, it does not change 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 that induced a nonzero would contradict the one-to-one relationship proved by the extension of the HK theorem to 1MFT.Gilbert (1975) It was mentioned above that there are null vectors associated with the pinned states. Suppose is a natural orbital with occupation number in the ground state. The perturbation of the external potential does not change the ground state if the system has an energy gap between the ground state and excited states and is small enough because , where and and are field operators. If , and the ground state is again unchanged by the perturbation. The “vector” is therefore a null vector of if is a pinned state. Another type of nonuniqueness, which has been called systematic nonuniqueness,Capelle and Vignale (2001) is related to constants of the motion. Suppose is a constant of the motion. The ground state, if it is nondegenerate, is an eigenstate of as constants of the motion commute with the Hamiltonian. If the system has an energy gap between the ground state and the first excited state, then a perturbation will not change the ground state wavefunction if is small enough. Thus, is a null vector of .

As in DFT, an exact and explicit expression for the universal energy functional in 1MFT is unknown in general. In actual calculations it is usually necessary to use approximate functionals. Many of the approximate energy functionals that have been introduced are expressed in 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 rather than explicit, they have been called “implicit” functionals. Recently, the KS equations were derived for this case.Pernal (2005) It was found that the contribution to the KS potential from electron-electron interactions can be evaluated by applying the following chain rule to (17),

 w(x,x′) = δWδγ(x′,x) (23) = ∑i∫dyδWδϕi(y)δϕi(y)δγ(x′,x) + ∑i∫dyδWδϕ∗i(y)δϕ∗i(y)δγ(x′,x) + ∑i∂W∂fiδfiδγ(x′,x).

### ii.2 Iteration of the KS equations

In this section we show that the “straightforward” procedure for iterating the KS equations (5) and (16) is intrinsically divergent.

The KS equations are nonlinear because the KS Hamiltonian itself depends on the 1-matrix. In favorable cases such nonlinear equations can be solved by iteration. Given a good initial guess for the 1-matrix, iteration may lead to the self-consistent solution corresponding to the ground state. In order to iterate (16), one needs an algorithm to define the 1-matrix of iteration step from the 1-matrix of step , i.e., one needs to “close” the KS equations. In the previous section, we saw that in the 1MFT-KS scheme the occupation numbers are held fixed during the optimization of the natural orbitals. The following is a “straightforward” algorithm that optimizes the natural orbitals: i) the KS Hamiltonian for step is defined by

 ^h(n+1)=^h[γ(n)], (24)

where is given by (15) and is the 1-matrix of iteration step ; ii) the eigenstates of are taken as the natural orbitals of step ; iii) the 1-matrix of step is constructed from the natural orbitals of step by the expression

 γ(n+1)(x,x′)=∑ifiϕ(n+1)i(x)ϕ∗(n+1)i(x′). (25)

Let be the eigenstates of the KS Hamiltonian . In operation (ii), the natural orbitals are chosen from among the such that has maximum overlap with , i.e., is the minimum possible. If this procedure converges to the stationary 1-matrix giving the lowest energy possible for the fixed set of occupation numbers, then it defines the function , introduced in the previous section, for which the ground state energy is the absolute minimum. Unfortunately, for any set of occupation numbers sufficiently close to the ground state occupation numbers, this procedure does not converge to . In other words, the “straightforward” algorithm defines an iteration map for which the ground state is an unstable fixed point. This is a consequence of the degeneracy of the KS spectrum at the ground state.

The divergence of the iteration map is revealed by a linear analysis of the fixed point. Suppose the occupation numbers are fixed to values perturbed from their ground state values by . Let us consider an iteration step and ask whether the next iteration takes us closer to the stationary point that gives the minimum energy for the fixed occupation numbers. The linearization of the iteration map at the stationary point gives

 δ^γ(n+1) = ^γ(n+1)−^γmin (26) ≈ ^χs[γmin](^v(n+1)s−^vmins) = ^χs[γmin](^h(n+1)−^hmin) ≈ −^χs[γmin]^χ−1δ^γ(n),

where is the KS potential. The response function was defined in (22). The KS response function is

 χs(x,x′;y,y′) = δγ(x,x′)δvs(y,y′) (27) = ∑i∑j≠ifi−fjϵi−ϵj × ϕj(x)ϕ∗i(x′)ϕ∗j(y)ϕi(y′).

In the last line of (26), we have used

 ^h(n+1)−^hmin ≈ −^χ−1δ^γ(n), (28)

which can be established by the following arguments. First consider

 ^h(n+1)−^hmin ≈ δ^h[γ]δγ∣∣ ∣∣γminδ^γ(n), (29)

which follows from (24). The KS Hamiltonian is an implicit functional of the 1-matrix, and (29) defines the first order change of the KS Hamiltonian with respect to a perturbation of that 1-matrix. Since the occupation numbers are close to their ground state values, we may make the replacement

 (30)

which is valid to . Thus, to establish (28) we need to show at the ground state 1-matrix . The KS Hamiltonian is associated with the original many-body Hamiltonian , which has external potential . According to (20)4, . Consider now a Hamiltonian with a slightly different external potential such that its ground state 1-matrix is . The associated KS Hamiltonian is . At the new ground state, . This allows us to relate to as

 δ^h = ^h[γ′gs]−^h[γgs] (31) = ^h[γ′gs]+δ^v−δ^v = −δ^v.

Finally, using we obtain , which verifies (28).

Returning to the question of convergence, we see that (26) implies that the next iteration takes us farther from the stationary point . The reason is that if is sufficiently close to the ground state. According to (21) the KS response diverges as because . For a fixed set of occupation numbers sufficiently close to their ground state values, the moduli of all eigenvalues of the operator become greater than 1 (the null space of is assumed to be excluded). Therefore, a perturbation from the ground state is amplified by iteration, and the ground state is an unstable fixed point of the iteration map. A fixed point is stable if and only if all eigenvalues of the linearized iteration map have modulus less than 1.

### ii.3 Level shifting method

In the previous section, it was shown that the “straightforward” iteration of the KS equations is intrinsically divergent. To obtain a practical KS scheme the iteration map must be modified. In this section, we consider the level shifting methodSaunders and Hillier (1973) and by linearizing the modified iteration map we obtain a criterion for convergence.

Intrinsic divergent behavior can be encountered also in the Hartree-Fock approximation, and various modifications of the iteration procedure have been introduced, for example, Hartree damping (also called configuration mixing) and level shifting. The level shifting method is particularly attractive in 1MFT because it can prevent the collapse of eigenvalues that is the origin of divergent behavior. Indeed, it has been shown that the level shifting method is capable of giving a convergent KS scheme in 1MFT.Pernal (2005)

In the “straightforward” iteration procedure, the change of the orbitals, to first order, from iteration step to iteration step is

 δϕi(x) = ϕ(n+1)i(x)−ϕ(n)i(x) (32) = ∑j≠i⟨ϕj∣∣^h(n+1)−^h(n)|ϕi⟩ϵi−ϵjϕj(x),

where is the KS Hamiltonian for iteration step and in the second line the orbitals and eigenvalues are from iteration step . In the level shifting method, the first order change in the orbitals given by (32) is altered by applying the shifts to the eigenvalues in the denominator. To first order, this modification is equivalent to adding the term to the KS Hamiltonian for step . Let define the level shifted Hamiltonian. Repeating the linear analysis of the previous section for the iteration map for this level shifted Hamiltonian, we find

 ^δγ(n+1) ≈ ^χs[γmin](^v(n+1)s−^vmins) (33) = ^χs[γmin](^h(n+1)ζ−^hminζ) ≈ (−^χs[γmin]^χ−1+^Ω)^δγ(n),

where we have defined the operator with the kernel

 Ω(xx′,yy′) = ∫dzdz′χs(xx′,zz′)δΔ(zz′)δγ(yy′) (34) = ∑i∑j≠iϕj(x)ϕ∗i(x′)ϕ∗j(y)ϕi(y′).

From the last line of (33), we obtain a criterion for the convergence of the iteration map. All of the eigenvalues of the operator

 ^A=−^χs[γmin]^χ−1+^Ω (35)

must have modulus less than 1. The dependence on the level shift parameters enters only through the shifted eigenvalues in the denominator of . The level shifting method is effective because it prevents the divergence of at the ground state and there is a cancellation between the two terms in (35). Unfortunately, the convergence criterion depends on , which is unknown at the outset of a 1MFT calculation. In Sec. III, the level shifting method is applied in an explicit example and the above criterion is verified.

### ii.4 Properties of the KS system

The distinguishing feature of the KS system in 1MFT is the degeneracy of the eigenvalue spectrum. This has surprising consequences. It was shown in section II.1 that the KS eigenvalue spectrum splits linearly as we move away from the ground state 1-matrix. Therefore, the total KS energy changes linearly with respect to the displacement, i.e., , where (for a specific example see Fig. 4 in Sec. III.2.2). This is surprising because such linear changes do not occur for the energy functional (in the VR space). The immediate implication is that is not stationary at the ground state. While this causes no difficultly in principle — we need only the functional to be stationary — it is intimately connected with the divergence of the iteration map. Precisely at the ground state , where the prime indicates that only the pinned states with contribute to the sum. Away from the ground state the KS eigenvalue spectrum splits, and is a multivalued functional due to the choice implied in occupying the new KS levels. This is the same choice encountered in the iteration of the KS equations (see Sec. II.2), where the natural orbitals are selected from among the eigenstates of the KS Hamiltonian. Near the self-consistent solution, there will be one such choice for which the resulting is very close to .

It was shown in the preceding section that the static response function of the KS system diverges at the ground state. Thus, even an infinitesimal perturbation may induce a finite change of . At the ground state, all of the natural orbitals, except those which have an occupation number that is degenerate, are uniquely defined. The natural orbitals which belong to a degenerate occupation number are only defined modulo unitary rotation in the degenerate subspace. When a perturbation is introduced, the natural orbitals change discontinuously to the eigenstates of the perturbed KS Hamiltonian . These eigenstates may be any functions in the degenerate Hilbert space because is arbitrary.

## Iii Two-site Hubbard model

The 1MFT-KS system has some unusual features, such as the collapse of the KS eigenvalues at the ground state, so it is desirable to derive explicitly the KS equations for a simple model. The Hubbard model on two sites provides a convenient example because it is exactly solvable and especially easy to interpret. Also, analytic expressions for the 1-matrix energy functional and KS Hamiltonian can be obtained. In the following sections, for the purpose of comparison, we find the ground state of the two-site Hubbard model by three methods — direct solution of the Schrödinger equation, 1MFT and DFT.

### iii.1 Direct solution

The Hamiltonian of the two-site Hubbard model is with

 ^T = −∑σ(t12c†1σc2σ+t21c†2σc1σ) ^U = U(^n1↑^n1↓+^n2↑^n2↓) ^V = V12(^n1−^n2), (36)

where , and are the creation and annihilation operators of an electron at site with spin , and . We consider only the sector of states with and , i.e., a spin unpolarized system. In this sector, the eigenstates of are

 Φ0=1√2⎛⎜ ⎜ ⎜⎝yxxy⎞⎟ ⎟ ⎟⎠,Φ1=1√2⎛⎜ ⎜ ⎜⎝01−10⎞⎟ ⎟ ⎟⎠,Φ2=1√2⎛⎜ ⎜ ⎜⎝100−1⎞⎟ ⎟ ⎟⎠,Φ3=1√2⎛⎜ ⎜ ⎜⎝x−y−yx⎞⎟ ⎟ ⎟⎠, (37)

in the site basis , , , . The following variables have been introduced , , and with . The eigenvalues of for the states are

 λ0=−By2,λ1=0 λ2=B(x2−y2),λ3=Bx2,

where and . , , and are singlet states () and is a triplet state with and . We will omit from consideration as it is not coupled to the other states by the spin-independent external potential chosen in (36). The Hamiltonian may be written as , where

 K=⎛⎜⎝0νy0νyx2νx0νx1⎞⎟⎠ (38)

in the basis . We have defined the dimensionless variable . The secular equation is

 κ3i−(x2+1)κ2i+(x2−ν2)κi+ν2y2=0. (39)

The normalized eigenvectors of are

 Ψi=1ηi⎛⎜⎝ν2xyνxκiβi⎞⎟⎠ (40)

and have energy for , where is a root of the secular equation (39). We have also defined

 βi=κi(κi−x2)−ν2y2 (41)

and

 ηi = [κ2i(ν2(3x2−1)+y2)+κi(x2y2(2ν2−1)) (42) +ν2y2(ν2−y2)]1/2.

The two dimensionless energy scales of the system are the interaction strength and the bias . The behavior of the system with respect to these energy scales is illustrated in Fig. 1. The quantity , where is the average ground state occupancy of site , is plotted with respect to the external potential for various values of the interaction strength .