Reduced Density Matrix Functional Theory at Finite Temperature:Theoretical Foundations

# Reduced Density Matrix Functional Theory at Finite Temperature: Theoretical Foundations

## Abstract

We present an ab-initio approach for grand canonical ensembles in thermal equilibrium with local or nonlocal external potentials based on the one-reduced density matrix. We show that equilibrium properties of a grand canonical ensemble are determined uniquely by the eq-1RDM and establish a variational principle for the grand potential with respect to its one-reduced density matrix. We further prove the existence of a Kohn-Sham system capable of reproducing the one-reduced density matrix of an interacting system at finite temperature. Utilizing this Kohn-Sham system as an unperturbed system, we deduce a many-body approach to iteratively construct approximations to the correlation contribution of the grand potential.

###### pacs:
31.15.ec, 65.40.-b

## I Introduction

Based on the celebrated theorems of Hohenberg and KohnHohenberg and Kohn (1964), Kohn-Sham density functional theory (KS-DFT)Kohn and Sham (1965) is currently the method of choice for calculating groundstate (gs) properties of quantum systems.

There are, however, cases in which KS-DFT performs rather poorly. A prominent example is its failure in predicting the fundamental gap, in particular, of so-called Mott insulatorsGrüning et al. (2006); Giantomassi et al. (2011). KS-DFT with standard exchange-correlation approximations fails for this kind of strongly correlated system and typically yields a metallic ground state, while the true experimental low-temperature phase is that of an antiferromagnetic insulator. At finite temperature the description of strongly correlated systems is even more challenging. Genuine Mott insulators exhibit a characteristic feature: when heated up from their antiferromagnetic insulating gs, they stay insulating above the Néel temperature, i.e., in the absence of long-range magnetic order. Contrarily, weakly correlated insulators – so-called Slater insulators – become metallic at the Néel temperature.

A possible approach to tackle this challenge is to search for more accurate functionals in the framework of KS-DFT. Alternatively, one may look for other theoretical frameworks in which the treatment of strong correlation might be simpler. One candidate for such a framework is reduced density matrix functional theory (RDMFT). Through its more direct treatment of many-particle correlations it has promising potential for calculations of finiteBaerends (2001); Buijse and Baerends (2002); Helbig et al. (2007); Marques and Lathiotakis (2008); Piris et al. (2010); Helbig et al. (2009) as well as infinite systemsLathiotakis et al. (2010, 2009); Sharma et al. (2008). In particular, it was possible to predict insulating ground states for transition metal oxides without breaking the spin symmetrySharma et al. (2008).

Motivated by the success of RDMFT at zero temperature, the purpose of the present work is to lay the theoretical foundations for the finite-temperature version of RDMFT (FT-RDMFT). As a general ab-initio theory its applicability is not restricted to Mott insulators. There is a variety of physical phenomena, in particular in the warm dense matter regimeGraziani et al. (2014) which requires an accurate description of quantum effects at finite temperatureKnudson and Desjarlais (2009). These phenomena include temperature-driven magneticLevy (1979); Rosengaard and Johansson (1997) or superconductingProfeta et al. (2006); Cudazzo et al. (2008) phase transitions in solids, femtochemistry at surfaces of solidsGavnholt et al. (2009), properties of shock compressed noble gasesMilitzer (2006); Root et al. (2010), the properties of plasmasDharma-wardana and Perrot (1982); Perrot and Dharma-wardana (2000); Dharma-wardana and Murillo (2008), thermal conductivities of inertial confinement fusion capsulesAtzeni and Meyer-ter Vehn (2004), and planetary interiors and their formation processesNettelmann et al. (2012, 2013); Lorenzen et al. (2009, 2011); Redmer et al. (2011).

This paper is divided as follows: In Sec. II we derive and present the formalism of FT-RDMFT. First, in Sec. II.1 we introduce our notation. Note that we work in atomic units throughout where so that lenghts are expressed in Bohr radii, and energies in hartree. Then, in Sec. II.2 we lay the foundations of FT-RDMFT by showing that the grand potential of systems with generally nonlocal external potentials can be written as a functional of the one-reduced density matrix (1RDM). Next, in Sec. II.3 we show the existence of a KS system in FT-RDMFT and demonstrate how the KS Hamiltonian is explicitly constructed. Subsequently, in Sec. II.4 we derive the adiabatic connection formula which forms the basis for the construction of approximations to the correlation functional in FT-RDMFT. Finally, in Sec. II.5 the existence of a KS system and the adiabatic connection formula enable us to derive a methodology for iteratively constructing correlation functionals based on finite-temperature many-body perturbation theory (FT-MBPT). Furthermore, in Appendix A we give a detailed analysis of occupation numbers in interacting systems, in Appendix B we investigate the one-to-one mapping between the external potential and the wavefunction at zero temperature, in Appendix C we show that our iterative procedure for constructing functionals from FT-MBPT yields the finite-temperature Hartree-Fock functional as the first-order contribution, and in Appendix D we present the formulation of FT-RDMFT for a canonical ensemble.

## Ii Finite-temperature reduced density matrix functional theory

### ii.1 Background

The main thermodynamic variable in a grand canonical ensemble is the grand potential

 Ω[^D] =\textuptr{^D^Ω} (1)

given as a statistical average over the grand canonical operator

 ^Ω=^H+μ^N+^S/β, (2)

where , , and are the Hamiltonian, particle number, and entropy operators. In electronic structure theory the Hamiltonian is typically given by , where denotes the kinetic energy operator, the interelectronic repulsion in a Coulombic system, and represents a scalar external potential. The coupling to particle and heat baths is achieved via the Lagrangian multipliers denoting the chemical potential and denoting the temperature.

Statistical averages as in Eq. (1) are computed via the statistical density operator (SDO) which is defined as a weighted sum of projection operators on the underlying Hilbert space. The appropriate Hilbert space for grand canonical ensembles, where a change of particle number is allowed, is a direct sum of symmetrized tensor products of the one-particle Hilbert space – called the Fock space. Assuming that the system does not allow for mixing of states with different particle numbers, the set of all possible SDOs can be expressed just by projection operators on states with defined particle number :

 ^D=∑α,NwαN|ΨαN⟩⟨ΨαN|,wαN≥0,∑α,NwαN=1, (3)

where and are orthonormal -particle states and their corresponding weights.

The thermal equilibrium (eq) of a grand canonical ensemble is then defined as that SDO for which the grand potential is minimal. This definition leads to the finite-temperature Rayleigh-Ritz variational principleMermin (1963) which states that

 Ω[^D] >Ω[^Deq],^D≠^Deq, (4)

with

 ^Deq =e−β(^H−μ^N)/\textuptr{e−β(^H−μ^N)}. (5)

The 1RDM is defined by the SDO and the help of the common fermionic field operators as

 γ(x,x′) =\textuptr{^D^ψ+(x′)^ψ(x)}, (6)

where the variable denotes a combination of spin index and spatial coordinate r (). An integration over is therefore to be interpreted as an integration over r and a summation over . Since the 1RDM is hermitean by construction, it is commonly written in spectral representation as

 γ(x,x′) =∑iniϕ∗i(x′)ϕi(x) (7)

with real-valued eigenvalues and eigenfunctions , which are called occupation numbers and natural orbitalsLöwdin (1955). The neccessary and sufficient conditions for N-representabilityColeman (1963) of are that is a complete set and

 0≤ni≤1,∑ini=N . (8)

It is sometimes desirable to treat spin and spatial variables seperately. To this end we introduce a two-component (Pauli) spinor notation.

 Φi(r)=(ϕi1(r)ϕi2(r)), (9)

where are the orbitals of Eq. (7). The 1RDM can then be written as a matrix in spin space as

 γ(r,r′) =∑iniΦ†i(r′)⊗Φi(r) (10) =∑ini(ϕ∗i1(r′)ϕi1(r)ϕ∗i2(r′)ϕi1(r)ϕ∗i1(r′)ϕi2(r)ϕ∗i2(% r′)ϕi2(r)) . (11)

In the special case of collinear spin configuration different spin channels can be treated seperately. For these systems, the natural orbitals are so-called spin orbitals, i.e., spinors containing only one spin component, where

 Φi1(r)=(ϕi1(r% )0), Φi2(r)=(0ϕi2(r)) . (12)

This leads to a 1RDM which has only one nonvanishing entry in every 2x2 matrix of Eq. (11), either the 11 or the 22 element. Hence the complete 1RDM is diagonal with respect to the spin coordinate

 γσσ′(r,r′) =δσσ′∑iniσϕ∗iσ(r′)ϕiσ(r), (13)

where are the occupation numbers of the special spinors in Eq. (12). Spin spiral states are another special case where this seperation also appliesBaldsiefen et al. (2015).

### ii.2 Hohenberg-Kohn theorems for finite-temperature reduced density matrix functional theory

Mermin’s extension of the HK theorems to finite temperatureMermin (1965) immediately implies the one-to-one mappings

 γ(r,r′)⟶ρ(%r)\lx@stackrelHK⟷v(r)  and (14a) v(r)⟶Ψ(^D)⟶γ(r,r′), (14b)

i.e., between the 1RDM, the density, the eq-SDO, and the local external potential1.

However, in FT-RDMFT we need to go further than this and consider nonlocal external potentials, in which case the ground state is not uniquely determined by the density anymore, but by the 1RDMGilbert (1975). Going beyond local external potentials is necessary, because the KS potentials in FT-RDMFT are nonlocal in general, as we show in Sec. II.3.

Therefore, we lay the foundations of FT-RDMFT and establish the uniqueness of the KS scheme in FT-RDMFT by proving the one-to-one correspondences in Eqs. (14a) and (14b) for nonlocal external potentials at finite temperature. We also show that in this case the 1RDM is still sufficient to describe equilibrium properties. We divide this up into three steps, namely showing (i) that the map between and has to be invertible, implying the existence of a grand potential functional , (ii) the existence of a universal functional and (iii) that the minimization of leads to the eq-1RDM. Note that we consider only eq-1RDMs for the proof in step (i). However, we can relax this restriction in step (ii).

#### (i) Proof of ^Deq\lx@stackrel1−1⟷γeq(x,x′), i.e., one-to-one mapping between eq-SDO and eq-1RDM

We divide proving the existence of a one-to-one mapping between the eq-SDO and the eq-1RDM into two parts:

 ^Deq\lx@stackrel1−1⟷μ(x,x′)\lx@stackrel1−1⟷γeq(x,x′), (15)

where we prove (i.i) the one-to-one mapping between and the nonlocal chemical potential and (i.ii) the one-to-one mapping between and .

(i.i) Proof of , i.e., one-to-one mapping between eq-SDO and nonlocal chemical potential.    We show this with a proof by contradiction. Let and be two different Hamiltonians and assume they lead to the same SDO , where shall differ from only by a one-particle potential contribution . With Eq. (5) this reads

 e−β(^H−μ^N)/Z=e−β(^H+Δ^V−μ^N)/Z′, (16)

where and are the partition functions, e.g., . Solving Eq. (16) for yields

 Δ^V =∫dxdx′Δv(x′,x)^ψ+(x′)^ψ(x)=−1βlnZZ′. (17)

We now need to show that a one-particle potential fulfilling this equality cannot exist, thereby contradicting our initial assumption. To proceed, we assume three different Slater determinants , and in the basis . The potential in this basis is denoted by . Calculating the expectation value of both sides of Eq. (17) with respect to these three Slater determinants yields the following system of equations

 −1βlnZZ′ =Δv11=Δv22=Δv11+Δv22. (18)

This can only be fulfilled by and . A repetition of this argument for all possible bases then shows that only fulfills Eq. (17) which in turn proves the one-to-one correspondence between and , and hence the one-to-one correspondence between and .

This proof is valid for any finite temperature. It is based on the bijectivity of the exponential function which allows us to invert Eq. (16), leading to Eq. (18). At zero temperature, however, this bijectivity breaks down. Further elaborations on zero-temperature mappings between external potentials and wavefunctions are given in Appendix B.

(i.ii) Proof of , i.e., one-to-one mapping between nonlocal chemical potential and eq-1RDM.    In order to prove the one-to-one correspondence between and we use a proof by contradiction again. Consider two Hamiltonians and differing only in their external and chemical potentials. The corresponding grand potentials are given by

 Ω[^Deq] =\textuptr{^Deq(^H−μ^N+1/βln^Deq)}~{} ~{}and (19) Ω′[^D′eq] =\textuptr{^D′eq(^H′−μ′^N+1/βln^D′eq)}, (20)

where and are defined by Eq. (5). Using as we have proven in (i.i), the variational principle in Eq. (4) then leads to

 Ω[^Deq] <Ω[^D′eq] (21) =\textuptr{^D′eq(^H−μ^N+1/βln^D′eq)} (22) =Ω′[^D′eq]+\textuptr{^D′eq((^H−μ^N)−(^H′−μ′^N))} . (23)

By exchanging primed and unprimed quantities we obtain

 Ω[^Deq] <Ω′[^D′eq]+∫dxdx′(μ(x′,x)−μ′(x′,x))γ′(x,x′) (24) Ω′[^D′eq] <Ω[^Deq]+∫dxdx′(μ′(x′,x)−μ(x′,x))γ(x,x′) . (25)

 ∫dxdx′(μ′(x′,x)−μ(x′,x))(γ(x,x′)−γ′(x,x′)) >0. (26)

The existence of two different sets of external and chemical potentials yielding the same eq-1RDM lets the integral in Eq. (26) vanish which leads to a contradiction. Hence the initial assumption is falsified. This concludes our proof of Eq. (15).

Having established the existence of a one-to-one mapping between and we can now proceed and properly define the grand potential as a functional of the 1RDM as

 Ω[γeq] =\textuptr{^D[γeq](^H−μ^N+1/βln(^D[γeq]))} . (27)

#### (ii) Existence of a universal functional F[γeq]

In analogy to DFT, we define a universal functional by separating the external and the chemical potential contributions from Eq. (27):

 F[γeq] =\textuptr{^D[γeq](^T+^W+1/βln^D[γeq])} (28)

such that

 Ω[γeq] =∫dxdx′μ(x,x′)γeq(x,x′)+F[γeq] . (29)

Notice a subtlety involved with defining the universal functional in Eq. (28). In our proof we considered a restricted set of 1RDMs, namely those coming from eq-SDOs given by Eq. (5). However, the conditions to ensure that an arbitrary 1RDM comes from such a SDO are unknown. Nevertheless, following ideas of ValoneValone (1980) and LiebLieb (1983) we can resolve this subtle point and extend the domain of to the whole set of ensemble-N-representable 1RDMs. Accordingly we can now define the universal functional as

 F[γ] =inf^D→γ\textuptr{^D(^T+^W+1/βln^D)} (30)

such that

 Ω[γ] =∫dxdx′μ(x,x′)γ(x,x′)+F[γ]. (31)

#### (iii) Minimization of Ω[γ]

The variational principle in Eq. (4) now allows us to determine the equilibrium grand potential by

 Ωeq =minγ∈ΓNΩ[γ], (32)

a minimization over which is the set of all ensemble-N-representable 1RDMs. Additionally, we postulate

 δF[γ]δγ(x,x′)∣∣∣γeq+v(x,x′) =μ, (33)

the Euler-Lagrange equation for the eq-1RDM in FT-RDMFT2.

### ii.3 Kohn-Sham system for finite-temperature reduced density matrix functional theory

We have established the theoretical framework of FT-RDFMT by proving Hohenberg-Kohn-like theorems. The central problem for turning this theory into a practical scheme is finding approximations as a functional of the 1RDM. In analogy to DFT, one possible route for constructing such approximations requires us to introduce the KS scheme. Then we can exploit the existence of a KS system to derive a methodology for the iterative construction of functionals using methods from FT-MBPT.

Our starting point is an auxiliary system of noninteracting fermions described by the Hamiltonian

 ^H(1)S=^T+^VS,^H(1)S =∑iεi|ϕi⟩⟨ϕi| (34)

with eigenvalues and eigenfunctions , denoting the operator of the KS potential. Then we assume the existence of a nonlocal potential – the KS potential – which yields a ground-state or eq-1RDM that equals the true ground-state or eq-1RDM,

 γeq(x,x′)=γS,eq(x,x′) . (35)

Note that a KS system does not exist in RDMFT for Coulombic matter at zero temperature. The reason behind this is the presence of the electron-electron cusp emerging from the interelectronic repulsionKirzhnits (1957). Capturing this cusp requires a superposition of infinitely many Slater determinantsFriesecke (2003). Hence, the gs 1RDM for Coulombic systems has an infinite number of occupied orbitals, i.e., natural orbitals with occupation numbers .

In the following we show that in FT-RDMFT, however, such a KS system does indeed exist. For a grand canonical ensemble, the eq-1RDM is given by

 γS,eq(x,x′) =∑iniϕ∗i(x′)ϕi(x), (36)

where the natural orbitals of the 1RDM in the KS system are by definition identical to the natural orbitals of the true, interacting 1RDM defined in Eq. (7). The eigenvalues and the chemical potential completely determine the occupation numbers by the relation

 ni =11+eβ(εi−μ), (37)

which can be inverted to yield

 εi−μ =1βln(1−nini) . (38)

In contrast to the zero-temperature case, it is now possible to construct the KS Hamiltonian in the following way: The KS Hamiltonian is obtained via its spectral representation in Eq. (34), where its eigenvalues are determined from Eq. (38), while its eigenfunctions are given by the natural orbitals of the given 1RDM in Eq. (36). The occupation numbers of a KS system in thermal equilibrium at finite temperature cannot be 0 or 1, as can be seen from Eq. (37). This is also true for the interacting 1RDM of a grand canonical ensemble, as we show in Appendix A. Hence, it is ensured that the domain of the KS system includes the interacting system.

Furthermore, due to the variational principle the KS potential is generally nonlocalGilbert (1975). Its uniquess follows from the Hohenberg-Kohn-like theorems shown in Sec. II.2. It can be expressed explicitly as

 vS(x,x′) =∑i,j(δijεi−tij)ϕ∗i(x′)ϕj(x), (39)

where is the kinetic operator in the basis of natural orbitals. The requirement of locality can be imposed on the KS potential. This is computationally advantageous, but leaves the domain of an exact theory, because this requirement results in approximate natural orbitals that cannot be equal to the true natural orbitalsLathiotakis et al. (2014).

Having established the KS scheme, we can express the grand potential of the interacting system as

 Ω[γ] =F[γ]+V[γ]−μN[γ], (40)

where we express the universal functional in terms of common KS quantities as

 F[γ] =Ωk[γ]−S0[γ]β+ΩH[γ]+ΩX[γ]+ΩC[γ]. (41)

Here,

 Ωk[γ] =∫dx′limx→x′(−∇22)γ(x,x′), (42) V[γ] =∫dxdx′v(x,x′)γ(x′,x), (43) N[γ] =∫dxγ(x,x), (44) SS[γ] =−∑i(nilnni+(1−ni)ln(1−ni)), (45) ΩH[γ] =12∫dxdx′w(x,x′)γ(x,x)γ(x′,x′), and (46) ΩX[γ] =−12∫dxdx′w(x,x′)γ(x,x′)γ(x′,x) (47)

denote the functionals of kinetic energy, external potential, particle number, KS entropy, Hartree, and exchange energy, which are known explicitly3. The remaining term is the correlation contribution, but its exact form is not known explicitly.

### ii.4 Adiabatic connection formula in finite-temperature reduced density matrix functional theory

We derive the adiabatic connection formula in FT-RDMFT which allows us to connect the interacting system to the KS system with the same eq-1RDM and forms the basis for systematically constructing approximations to the correlation functional via FT-MBPT.

Closely following the standard zero-temperature DFT approachLevy and Perdew (1985); Görling and Levy (1993), we begin by introducing a coupling constant into the electronic Hamiltonian

 ^Hλ=^T+λ^W+^Vλ, (48)

where . The potential is chosen such that for any there is an associated eq-SDO that leaves the eq-1RDM invariant under a change of . Along with that we define an auxiliary Hamiltonian

 ^Hλa=^T+λ^W+^V, (49)

such that it agrees with Eq. (48) at full coupling-strength when , i.e., . Additionally, we also introduce an auxiliary potential such that . The grand potential for the auxiliary Hamiltonian becomes

 Ωλa[γ] =min^D→γ\textupTr{^D(^Hλa−μ^N+ln^D/β)} . (50)

With the aid of the auxiliary potential we obtain

 Ωλa[γ]=min^D→γ(\textupTr{^D(^Hλa+^Vλa−μ^N+ln^D/β)}−\textupTr{^D^Vλa}) . (51)

Since is a one-particle operator, we can take the last term out of the minimization and replace in this term by any . Then, the minimization in Eq. (51) yields the eq-SDO that is associated with yielding the eq-1RDM of the true interacting system that is invariant under a change of . Hence, the grand potential becomes

 Ωλa[γ]=\textupTr{^Dλ(^Hλa+^Vλa−μ^N+ln^Dλβ)}−\textupTr{^D^Vλa} . (52)

By definition , i.e., the auxiliary grand potential at full coupling strength is identical to the true interacting grand potential, therefore

 Ω[γ] =Ω0a+∫10dλdΩλa[γ]dλ . (53)

Taking the derivative with respect to the coupling constant is simplified by the fact that we consider a system in thermal equilibrium. Hence, only and contribute to the coupling-constant derivative in Eq. (53), yielding

 Ω[γ] =Ω0a+∫10dλ\textupTr{^Dλ^W} . (54)

Consider the grand potentials

 Ω[γ] =\textupTr{^D(^T+^V+^W−μ^N+ln^D/β)}, (55) Ω0a[γ] =\textupTr{^D0(^T+^V−μ^N+ln^D0/β)} (56)

and take into account that and yield the same eq-1RDM, hence, the same expectation values of one-particle operators, such as , , and .

Then, we can further reduce Eq. (54) and obtain the adiabatic connection formula for the entire interaction as

 ΩW[γ] =∫10dλWλ[γ], (57)

where we define , , and as the entropic correlation contribution.

Finally, by subtracting the Hartree and exchange contributions defined as we obtain the adiabatic connection formula for the correlation contribution

 ΩC[γ] =∫10dλλWλC[γ], (58)

where we define , , and .

In analogy to DFT, Eq. (58) allows us to express the correlation contribution to the KS system in FT-RDMFT as a contribution coming solely from the interaction potential. It is interesting to note another similarity between DFT and FT-RDMFT. In DFT, the adiabatic connection formula includes the kinetic correlation contribution, i.e., the difference between the kinetic energy of the interacting system and the KS system via the coupling constant integration. In FT-RDMFT, where there is no kinetic correlation contribution, the coupling-constant integral instead incorporates the entropic correlation contribution .

The adiabatic connection formulas derived in Eqs. (57) and (58) are central results, because they are our key for developing systematic approximations to and its correlation contribution based on FT-MBPT.

### ii.5 Constructing correlation functionals

With the aid of the adiabatic connection we can use methods from FT-MBPTFetter and Walecka (2003) to systematically construct approximations to the functionals and , where the KS system, defined by the Hamiltonian , serves as our reference system in a perturbative expansion. Our starting point from the perspective of FT-MBPT is to relate the temperature Green’s function to the adiabatic connection in Eq.(57). This relation is expressed as

 λWλ=12∫dxdx′limτ′→τ+(−δ(x−x′)∂∂τ−(−δ(x−x′)∇22+vS(x′,x)))Gλ(x,τ,x′,τ′)+12∫dxdx′(vS(x′,x)−vλ(x′,x))γ(x,x′), (59)

where denotes imaginary time and . The use of FT-MBPT in Eq. (59) is facilitated by the existence of the adiabatic connection which connects the true interacting system with the KS system, and hence, allows to express the resulting Feynman diagrams in terms of occupation numbers and natural orbitals of the 1RDM.

Well-known methods of FT-MBPT can now be applied. The unperturbed Hamiltonian is , whereas the peturbation consists of a two-particle interaction and a nonlocal one-particle potential . The proof of Wicks theorem is still applicable for this kind of perturbation and the same Feynman rules apply. We show our notation conventions in Table 1. {fmffile}gp.one \fmfsetdot_len1.5mm \fmfsetdash_len2.5mm,1mm \fmfsetarrow_len7 \fmfsetarrow_ang25

In particular, if the Hamiltonian is temperature-independent and the system is uniform, Eq. (54) can be written entirely in terms of Feynman diagrams as

 λWλunif (60)

where denotes the irreducible self-energy. Whereas in general, the irreducible self energy for the first-order contribution becomes

 Σ∗ (61)

Combining Eqs. (57), (54) and (61), we arrive at the first-order contribution to the interaction-induced grand potential functional in FT-RDMFT which is

 Ω(1)W[γ] =ΩH[γ]+ΩX[γ],~{}where~{} (62) ΩH[γ] =12∫dxdx′w(x,x′)γ(x,x)γ(x′,x′)~{}and~{} (63) ΩX[γ] =−12∫dxdx′w(x,x′)γ(x,x′)γ(x′,x). (64)

This justifies the definitions of the Hartree and exchange energies which we postulated in Eqs. (46) and (47). Note that the functional form of the first-order contributions are equivalent to the Hartree and exchange functionals in zero-temperature RDMFT4.

Approximations for the correlation functional can now be derived by expanding the Green’s function to higher orders and then solving Eqs. (57) and (59).

## Iii Summary and conclusions

In this work, we have derived and presented the foundations of FT-RDMFT. We have proven Hohenberg-Kohn-like theorems and shown that the equilibrium properties of a grand canonical ensemble with nonlocal external potential are determined uniquely by the eq-1RDM. This allows us to establish a functional theory for the grand potential in terms of the 1RDM and, in analogy to DFT, to define a universal functional. A minimization of that grand potential functional then yields the eq-1RDM.

Furthermore, we have shown that there exists a KS system in FT-RDMFT, in contrast to the zero-temperature case, and derived the adiabatic connection formula. Based on this, we have established an iterative procedure for constructing approximations to the correlation functional in FT-RDMFT by utilizing methods from FT-MBPT. We have further demonstrated that the minimization of the first-order functional in this perturbative scheme is equivalent to the solution of the finite-temperature Hartree-Fock equations.

The present work sparks the hope that FT-RDMFT might become the method of choice for quantum problems at finite temperature where the standard DFT approach fails and the thermal DFT approach has not been developed to satisfactionPribram-Jones et al. (2014).

The main task for the future is the development of correlation functionals for the grand potential and free energy in FT-RDMFT and the application to real systems. Some further developments, such as exchange-only functional for collinear and non-collinear spin, as well as correlation functionals, momentum distributions, and phase diagrams in the framework of FT-RDMFT will be presented in Ref. Baldsiefen et al. (2015).

## Appendix A Equilibrium occupation numbers in general systems

As we have pointed out in Section II.3, the eq-1RDM of a noninteracting system has occupation numbers strictly between 0 and 1. We now show that this is also true for the occupation numbers of eq-1RDMs of arbitrary systems, including interacting ones.

We start from the spectral representation of the eq-1RDM given by

 γ(x,x′) =∑iniϕ∗i(x′)ϕi(x) . (65)

The occupation number operator is now defined as

 ^ni =^c+i^ci, (66)

where creates and annihilates the natural orbital . An arbitrary occupation number of the eq-1RDM in grand canonical equilibrium can then be written as

 ni =\textuptr{^D^ni}=∑ewe⟨Ψe|^ni|Ψe⟩. (67)

The are eigenfunctions of the Hamiltonian and form a basis of the underlying Hilbert space. Another basis is formed by the Slater determinants which are constructed by the natural orbitals of the eq-1RDM. The transformation between these bases is governed by the expansion coefficients via

 Ψe =∑αceαΦα. (68)

Due to completeness and normalization of the and , the coefficients fulfill

 ∑e|ceα|2 =∑α|ceα|2=1. (69)

Expanding the in Eq. (67) in terms of the then leads to

 ni =∑ewe∑αβc∗eαceβ⟨Φα|^ni|Φβ⟩ . (70)

Since the Slater determinants are by definition eigenfunctions of the occupation number operator , this reduces to

 ni =∑α(∑ewe|ceα|2)fα⟨Φα|^ni|Φα⟩giα. (71)

Using Eq. (69) and the properties of the thermal weights, and , we see that

 fα >0 (72) ∑αfα =1. (73)

The factors are equal to 1, if the natural orbital appears in the Slater determinant . Otherwise vanishes. The summation over corresponds to a summation over a basis of the Hilbert space, which is the Fock space in the case of a grand canonical ensemble. Therefore, for a fixed , there will be at least one , such that and at least one for which . Combining this fact with Eqs. (72) and (73), we can rewrite Eq. (71) to yield the desired inequality

 0

## Appendix B Zero-temperature mapping between potentials and wavefunctions

Assume an arbitrary Hamiltonian

 ^H =^T+^V+^W (75)

with gs-1RDM

 γgs(x,x′) =∑iniϕ∗i(x′)ϕi(x) . (77)

Due to Gilbert’s theoremGilbert (1975), the wavefunction can be written as a functional of the 1RDM allowing us to define an energy functional

 E[γ] =⟨Ψ[γ]|^H|Ψ[γ]⟩ . (78)

From the variational principle we infer that this functional is minimal for the gs-1RDM. Therefore the following relations hold

 ∂E[γ]∂ni∣∣∣γgs =⎧⎪⎨⎪⎩ai>μ,ni=0μ,0

where we identify as the chemical potential of the system.

In the following we show that there is a one-to-one mapping between potential and gs-1RDM, if and only if there are no pinned occupation numbers, i.e., occupation numbers equal to 0 or 1. This is done in two steps: (i) For unpinned occupation numbers we show that the external potential is uniquely determined up to a constant; (ii) considering gs-1RDMs with pinned occupation numbers we show that one can explicitly construct infinitely many potentials which leave the gs-1RDM invariant.

(i) Unpinned states    The absence of pinned states allows us to use the following Euler-Lagrange equation

 δE[γ]δγ(x,x′)∣∣∣γgs =μδN[γ]δγ(x,x′)∣∣∣γgs=μδ(x,x′), (80)

where denotes the Dirac delta function. Adding an arbitrary potential contribution to the energy functional then yields

 δEu[γ]δγ(x,x′)∣∣∣γgs =μδ(x,x′)+u(x,x′). (81)

As every gs-1RDM has to fulfill Eq. (80) we deduce that the only choice of that leaves the gs-1RDM invariant is

 u(x,x′) =cδ(x,x′) (82)

with being an arbitrary constant.

(ii)Pinned states    For pinned occupation numbers the minimum of is at the boundary of the domain, and hence, we cannot use Eq. (80). It is possible to adjust the Euler-Lagrange equation by incorporating Kuhn-Tucker multipliersKuhn and Tucker (1951), but there is a simpler way, as described in the following.

We exploit the fact that the derivatives in Eq. (80) can be different from for pinned states and construct a one-particle potential which leaves the gs-1RDM invariant. This potential shall be governed by the generally nonlocal kernel . By choosing it diagonal in the natural orbital basis of the gs-1RDM we ensure that the orbitals do not change upon addition of the potential. For simplicity, we choose only one component to be non-vanishing, namely

 u(x,x′) =uϕ∗α(x′)ϕα(x) . (83)

We then define an energy functional

 Eα[γ] =E[γ]+∫dxdx′u(x,x′)γ(x′,x) (84)

for which the derivative with respect to the occupation numbers becomes