2.4 Mean-field approximation
3.1.3 Mapping to an effective noninteracting system
matrix which yields the minimum in Eq. (3.22) for this hopping matrix t must be the pure-state density ˆρ0 = |Ψ0[γ]ihΨ0[γ]|constructed from the ground state |Ψ0[γ]i associated withγ. Therefore,We[γ] = TrN{ρˆ0Wˆ} = hΨ0[γ]|Wˆ|Ψ0[γ]i must hold, which coincides with Eq. (3.14). This concludes the proof of statement(a).
The constrained-search formulation (3.19) extends the interaction-energy func-tional (3.14) to the domain Γe(N) of ensemble N-representable SPDMs, which are easy to characterize by the condition (3.18). This opens the way to practical imple-mentations of the variational principle for the energy functional (3.21), from which the ground-state energy andgs-SPDM, and in principle all ground-state observables can be obtained by virtue of Corollary3.2. The next section addresses the problem of implementing the variational principle for the ground-state energy in practice. It will be shown that the interacting many-particle system can be mapped to a sys-tem ofN noninteracting particles with an effective hopping matrixts[γ]which is a functional of theSPDMγ itself. In this way, the ground-state problem forN interact-ing particles can be formally reduced to the selfconsistent solution of sinteract-ingle-particle equations.
For theSPDMγwhich minimizesEs[γ], i. e., for thegs-SPDMassociated with the non-interacting system characterized by the hopping matrixts, the corresponding Euler-Lagrange functional
Ls = Es[γ] −X
kσ
λkσ
X
i
uikσ∗ uikσ −1
−µ
X
kσ
ηkσ −N
(3.26)
must be stationary. Lagrange multipliersλkσ andµ have been introduced in order to enforce the normalization of the eigenvectorsukσ and to satisfy the second condition in Eq. (3.18) which accounts for the number of particles. At a stationary point of Ls
its derivative with respect tou∗jkσ must vanish, and by using Eqs. (3.23) and (3.24) we thus obtain the following eigenvalue equation for the vectorsukσ
ηkσ
X
i
tijσs uikσ =λkσujkσ ∀jkσ. (3.27) Forηkσ = 0 theSPDM(3.24) is independent of the eigenvectorukσ and therefore it can be chosen as an arbitrary normalized vector. Forηkσ > 0 we can recast Eq. (3.27) in the form
tσs>
ukσ =εkσukσ ∀kσ, (3.28)
with the spin-resolved hopping matrixtsσ andεkσ =λkσ/ηkσ. Equation (3.28) demon-strates that the eigenvectorsukσ of the gs-SPDM associated with the given nonin-teracting system are common eigenvectors of the corresponding transposed3 spin-resolved hopping matrix tσs>. The fact thattσs is hermitian ensures that the eigen-vectorsukσ of thegs-SPDMcan be chosen to be mutually orthogonal.
Having identified the eigenvectorsukσ of thegs-SPDMas the eigenvectors of the transposed hopping matrix tσs>, it remains to determine the corresponding occupa-tion numbersηkσ. To this aim, we combine Eqs. (3.23), (3.24), and (3.28), and obtain the single-particle energy functional in diagonal form
Es[γ]= X
kσ
εkσηkσ. (3.29)
Thus, the occupation numbers 0 ≤ ηkσ ≤ 1 which correspond to the ground state of the given noninteracting system are immediately identified as
ηkσ =1 for εkσ < µ ηkσ =0 for εkσ > µ 0 ≤ηkσ ≤ 1 for εkσ =µ,
(3.30)
3Notice, that the hopping matrixtσs is hermitian, such that transposition is equivalent to complex conjugation and therefore this operation can be discarded in the usual case of a real hopping matrix.
where the chemical potentialµmust be chosen such that the conditionP
kσηkσ = N is satisfied. If the Fermi levelεkσ = µ is nondegenerate and the numberNσ = P
kηkσ of particles with spin polarizationσ is integer, we haveηkσ ∈ {0,1} for allk and the correspondingSPDMis thus idempotent, i. e.,γσ2 =γσ. The ground state |Ψ0iof the given noninteracting system is then a single Slater determinant of the form
|Ψ0i=Y
kσ
bˆkσ† ηkσ
|vaci, (3.31)
where
bˆkσ† =X
i
uikσ∗ cˆiσ† (3.32)
are the creation operators associated to single-particle states which we will refer to as natural orbitals. Conversely it is clear, that the only kind of many-body state which can give rise to an idempotentSPDMis a single Slater determinant of the form (3.31).
Notice, however, that the ground state of an interacting system can not be repre-sented by a single Slater determinant, and that the correspondinggs-SPDMγ is con-sequently not idempotent. This meansγ must always involve fractional occupation numbers 0<ηkσ < 1 if it is associated with the ground state of an interacting many-particle system. Therefore, if we try to map a given interacting system to an auxiliary noninteracting system which is supposed to yield the samegs-SPDM, we will face the difficulty that the Fermi levelεkσ = µσ of the auxiliary noninteracting system must be degenerate in order to produce fractional ground-state occupation numbers. This means that thegs-SPDMof the auxiliary noninteracting system can not be unique, since variations of the occupation numbersηkσ at the degenerate Fermi level have no impact on the energy as long as the number of particles is kept fixed. In order to resolve these difficulties we will further below propose to map the given interacting system to an auxiliary noninteracting system in equilibrium at a fictitious finite tem-peratureTa > 0, since the occupation numbers in equilibrium at a finite temperature are always fractional.
Before doing so, let us return to the given interacting system and minimize the energy functional (3.21) with respect to all ensemble N-representable SPDMsγ ∈ Γe(N). To this aim, we seek for the stationary points of the corresponding Euler-Lagrange functional
L =Ee[γ] −X
kσ
λkσ
X
i
uikσ∗ uikσ −1
−µ
X
kσ
ηkσ −N
. (3.33)
Again, we have introduced Lagrange multipliersλkσ and µ in order to enforce the normalization of the eigenvectorsukσ of theSPDMand to obtain the desired particle
number. At a stationary point, the derivative ofLwith respect tou∗jkσ vanishes, such that we obtain the following eigenvalue equation for the eigenvectorsukσ
ηkσ
X
i
tijσ + ∂We[γ]
∂γijσ
uikσ =λkσujkσ ∀jkσ. (3.34) This eigenvalue equation is formally identical to the corresponding equation (3.27) for the noninteracting system if we choose the hopping integrals of the auxiliary system as
tijσs [γ]=tijσ + ∂We[γ]
∂γijσ . (3.35)
Therefore, thegs-SPDMγ of the interacting system is also ags-SPDM of the auxil-iary noninteracting system with the hopping integralstijσs [γ]. Equations (3.24), (3.28), (3.30), and (3.35) make up an iterative scheme, by which the gs-SPDMof a given in-teracting system can be determined from the solution of an effective single-particle problem. These equations must be solved in a selfconsistent manner, since the ef-fective hopping integrals (3.35) depend on the SPDMγ itself. As already discussed above, a difficulty arises from the fact that fractional occupation numbers 0<ηkσ < 1, which are characteristic for interacting ground states, are obtained only if the auxil-iary noninteracting system has a degenerate ground state. In this case, the fractional occupation numbersηkσ at the Fermi level are not unique and only constrained by the conditionN =P
kσηkσ. This renders the practical implementation of the iterative scheme based on the ground state of the auxiliary noninteracting system impractical.
One possibility to resolve the issues arising from the requirement of fractional occu-pation numbers was proposed by Saubanére, Lepetit, and Pastor [87]. They suggested to consider an auxiliary noninteracting system at a fictitious finite temperatureTa > 0, such that the correspondingeq-SPDMequals the one in the ground state of the given interacting system. The advantage of an auxiliary noninteracting system at a finite temperature is that the occupation numbers in equilibrium are fractional. In order to introduce the method ofSaubanéreet al., let us consider the eigenvectorsukσ of theSPDM as fixed for the moment, such that the energy (3.21) can be regarded as a functional of the occupation numbersηkσ alone, i. e.,E = Ee[η]. In order to seek for an extremum of the energy functionalEe[η]under the constraint of a given particle numberN, we consider the corresponding Euler-Lagrange functional
L =Ee[η] −µ
X
kσ
ηkσ −N
. (3.36)
At a stationary point ofL its variation with respect to allηkσ vanishes, such that
∂Ee[η]
∂ηkσ = µ ∀kσ. (3.37)
Let us now consider an auxiliary noninteracting system with spin-resolved hopping matrixtσs which is chosen such that Eq. (3.28) is fulfilled for the given set of eigenvec-torsukσ and some set of corresponding eigenvaluesεkσ. The energy of this noninter-acting system is then given by
Es = X
kσ
εkσηkσ, (3.38)
whereηkσ are the occupation numbers of the natural orbitals associated with the cre-ation operators (3.32). Working in a grand-canonical ensemble at a temperatureTa > 0 and chemical potentialµ, the equilibrium occupation numbers of the noninteracting system follow a Fermi-Dirac distribution
ηkσ = 1
1+e(εkσ−µ)/kBTa . (3.39) We require that the equilibrium occupation numbers of the auxiliary noninteracting system coincide with the eigenvalues of thegs-SPDMof the given interacting system, such that they lead to a stationary point ofL. This means, the equilibrium occupa-tion numbers of the auxiliary noninteracting system must satisfy the condioccupa-tions (3.37) and (3.39) together. The energy levels of the auxiliary noninteracting system are there-fore obtained as
εkσ[η]= ∂Ee[η]
∂ηkσ +kBTa log1−ηkσ
ηkσ
. (3.40)
Since these energy levels depend on the occupation numbers ηkσ themselves, Eqs. (3.39) and (3.40) must be solved in a selfconsistent manner. The full iterative procedure for the minimization of the energy functionalEe[γ]can be implemented as follows:
1. Start from an arbitrary ensembleN-representableSPDMγ, compute its eigen-valuesηkσ and the effective hopping integrals (3.35). Solve the eigenvalue equa-tion (3.28) for the eigenvectorsukσ.
2. Optimize the occupation numbersηkσ by keeping the currentukσ fixed. To this aim, select a fixed auxiliary temperatureTa > 0, compute the effective energy levels (3.40) and use them to updateηkσ via Eq. (3.39). Iterate until convergence is achieved.
3. Compute the new SPDM γijσ = P
kuikσηkσujkσ∗ from the current eigenvec-torsukσ and the optimized occupation numbersηkσ. Use Eq. (3.35) to determine new effective hopping integralstijσs .
4. Solve the eigenvalue equation (3.28) with the updated hopping integralstijσs in order to obtain new eigenvectorsukσ0 .
5. Exit if convergence is achieved, i. e.,uikσ0 = uikσ ∀ikσ. Otherwise, return to step2after having updatedukσ withukσ0 .
The rate of convergence in the presented iterative procedure is affected by the choice of the initial SPDMγ and the auxiliary temperatureTa > 0. Practical im-plementations for model Hamiltonians, such as the Hubbard model, have shown that a good choice for the initialSPDMγ is obtained from the ground state of the corre-sponding noninteracting system which is obtained by setting ˆW = 0. If the noninter-acting ground state is degenerate, one can freely choose the occupation numbers at the Fermi level such that the conditionP
kσηkσ =N is satisfied. Concerning the auxil-iary temperatureTa, it is clear that extreme values should be avoided. Good values for the auxiliary temperature are found in the order of the hopping integralskBTa ' |tijσ|. Moreover,Ta may be tuned to some extend, depending on the specific lattice model and interaction strength under consideration, in order to further improve the rate of convergence.