A crucial theorem in Reduced Density Matrix Functional Theory (RDMFT) suggests that the universal pure and ensemble functionals coincide on their common domain of pure N-representable one-matrices. We refute this by a comprehensive analysis of the geometric picture underlying Levy’s constrained search. Moreover, we then show that the ensemble functional follows instead as the lower convex envelop of the pure functional. It is particularly remarkable that the pure functional determines the ensemble functional even outside its own domain of pure N-representable one-matrices. From a general perspective, this demonstrates that relaxing pure RDMFT to ensemble RDMFT does not necessarily circumvent the complexity of the one-body pure N-representability conditions (generalized Pauli constraints). Instead, the complexity may simply be transferred from the underlying space of pure N-representable one-matrices to the structure of the universal one-matrix functional.

Reduced density matrix functional theory (RDMFT)1–5 extends the widely used density functional theory (DFT)6–9 by involving the full one-particle reduced density matrix (1RDM) γ rather than just the spatial density. This therefore facilitates the exact description of the energy of any one-particle Hamiltonian h (including, e.g., the kinetic energy or a non-local external potential). Furthermore, RDMFT allows explicitly for fractional occupation numbers as it is required in the description of strongly correlated systems4 and thus offers promising prospects of overcoming the fundamental limitations of DFT.10,11 At the same time, involving the full 1RDM leads also to drawbacks relative to DFT: The complexity of the N-representability problem, e.g., is not only hidden in the structure of the universal functional as in DFT12 but also the space of underlying 1RDMs is already non-trivial. To explain this aspect crucial to our work, we consider Hamiltonians of the form H = h + V on the N-fermion Hilbert space HNN[H1], where h is a one-particle Hamiltonian and V is some interaction (e.g., Coulomb pair interaction) which is fixed for the following. Moreover, we assume a finite-dimensional one-particle Hilbert space H1 and denote the convex set of N-fermion density operators Γ by EN and the subset of pure states by PN. A general expression for the universal functional13 then follows immediately by determining the ground state energy of H

E(h)=minΓPNTrN[(h+V)Γ]=minγPN1Tr1[hγ]+minPNΓγTrN[VΓ]minγPN1Tr1[hγ]+Fp(γ).
(1)

In the second line, we have introduced the set PN1 of pure N-representable 1RDMs γ and the last line gives rise to the universal pure functional Fp defined on PN1. The crucial point is now that PN1 is not only constrained by the simple Pauli exclusion principle, 0 ≤ γ ≤ 1 but also there are rather involved additional one-body pure N-representability conditions (generalized Pauli constraints), linear conditions on the eigenvalues of the 1RDM.14–16 To circumvent at first sight the complexity of those generalized Pauli constraints, Valone proposed17 to relax in (1) the set PN to EN by skipping the purity, leading to

E(h)=minγEN1Tr1[hγ]+Fe(γ)
(2)

with the ensemble functional Fe(γ)minENΓγTrN[VΓ] defined on the convex set EN1 of ensemble N-representable 1RDMs. One may now expect that the complexity of the one-body pure N-representability conditions is simply transferred within ensemble RDMFT from the underlying set of 1RDMs to the structure of the exact functional Fe. This, however, seems not to happen according to Ref. 18, suggesting and proving that Fe and Fp coincide on their common domain PN1 of pure N-representable 1RDMs, FpFe|PN1. In our work, we refute this fundamental theorem in RDMFT and show that the ensemble functional follows instead as the lower convex envelop of the pure functional. For this, we first need to develop a better understanding for the space of N-fermion density matrices exploited in Levy’s constrained search,13 i.e., the sets

PN(γ){ΓPN|Γγ},EN(γ){ΓEN|Γγ}
(3)

of N-fermion density operators Γ mapping to a given 1RDM γ.

The simplest way to refute the suggested equality FpFe|PN1 is to find one counterexample. A simple one is given by the asymmetric Hubbard dimer

H=tσc1σc2σ+c2σc1σ+σϵ1n1σ+ϵ2n2σ+Un1n1+n2n2,
(4)

a system of two electrons on two sites. Here, ciσ(c) denotes the creation(annihilation) operator of an electron at site i with spin σ and niσ=ciσciσ is the corresponding occupation number operator. The first two terms in Eq. (4) represent the kinetic and external potential energy, while the last one describes the interaction (V) between the electrons. We restrict H to the three-dimensional singlet space which contains the ground state. It is an elementary exercise19,20 to determine the respective pure functional Fp

Fp(γ)=U12γ122114γ1224(γ1112)2+(γ1112)2γ122+(γ1112)2,
(5)

with γij ≡ ⟨i ↑ | γ | j ↑⟩ = ⟨i ↓ | γ | j ↓⟩, i, j = 1, 2. Fp(γ) is invariant under γ11 → (1 − γ11) (particle-hole duality21) and γ12 → −γ12.

On the one hand, result (5) and its graphical illustration in Fig. 1 reveal that the pure functional Fp for the Hubbard dimer (4) is not convex on the set P21=E21 [described by the condition19,20(γ1112)2+(γ12214)]. On the other hand, it is well-known4,22 and rather elementary to verify that the ensemble functional Fe is always convex. As a consequence, the Hubbard dimer already refutes the suggested equality FeFp on P21.

It is instructive to understand the geometric picture of density matrices Γ underlying Levy’s constrained search (1) and (2). This will, in particular, reveal the loophole in the derivation in Ref. 18. Let us first recall that the set EN of N-fermion density matrices is convex and also compact as a subset of the space of hermitian matrices with fixed trace (i.e., it is bounded and closed). Its extremal points are given by the pure states, forming the compact but non-convex set PN. These are the idempotent matrices, Γ = Γ2 (i.e., their eigenvalues all vanish except one). It is worth noticing that a “point” Γ in EN lies on the boundary if and only if Γ is not strictly positive, i.e., at least one of its eigenvalues vanishes. As a consequence, most boundary points are not extremal points. It is one of the crucial insights of our work that this changes considerably if we restrict this consideration to the subsets EN(γ) and PN(γ) with respect to which the minimization (3) is carried out: While both sets EN(γ) and PN(γ) are also compact and EN(γ) is convex for all γ, extremal states Γ of EN(γ) are not necessarily pure anymore. The general reason for this is that a convex decompositions of Γ (e.g., the spectral decomposition into pure states) involves states whose 1RDMs typically differ from γ. Thus, a mixed (i.e., non-pure) Γ might be extremal within EN(γ) despite the fact that it is not extremal within EN.

As already stated above, the ensemble functional Fe(γ) follows for each γEN1 by minimizing TrN[V Γ] over EN(γ). Since TrN[V (·)] is linear and EN(γ) is convex and compact, the minimum (i.e., Fe(γ)) is attained on the boundary of EN(γ). This is a general (and rather obvious) fact from linear optimization: First, we observe that TrN[V Γ] is nothing else than the standard inner product on the Hilbert space of hermitian matrices, TrN[V Γ] ≡ ⟨V,Γ⟩N. In that sense, there is given a notion of geometry on the space of density operators23–34 and thus V defines a direction in EN(γ). The set of ΓEN(γ) with a specific interaction energy v = ⟨V,Γ⟩N gives rise to a hyperplane, orthogonal to V. The minimum of TrN[V (·)] ≡ ⟨V,·⟩N on EN(γ) then follows by shifting the hyperplane along the direction −V (i.e., by reducing v) until an extremal point of EN(γ) is reached. This final hyperplane is a so-called supporting hyperplane.35 By definition, this means that EN(γ) is entirely contained in one of the two closed half-spaces bounded by that hyperplane and EN(γ) has at least one boundary point on the hyperplane. This geometric picture underlying Levy’s constrained search is illustrated in Fig. 2 for different interactions (“directions”) V. There are three conceptually different boundary points which can be characterized by referring to two distinctive features: On the one hand, points ΓA and ΓB have a unique supporting hyperplane (unique “normal” vector V), in contrast to ΓC supported by infinitely many hyperplanes. On the other hand, points ΓA and ΓC are exposed35 in contrast to ΓB, i.e., they are supported by hyperplanes which do not contain any further boundary points. In other words, ΓA and ΓC can be obtained as unique minimizers of TrN[V Γ] for some V.

After having explained the geometric picture underlying Levy’s constrained search, we can now identify the loophole of the proof in Ref. 18 which we briefly recap: For γPN1, the minimizer of TrN[V Γ] on EN(γ) is denoted by Γ¯. Since EN(γ) is convex and compact, Γ¯ can be expressed according to the Krein-Milman theorem36 as a convex combination of the extreme points of EN(γ). This convex combination can be grouped into two parts, one (wpΓp) arising from pure extremal states and one ((1wp)Γe) arising from mixed states, i.e., Fe(γ)=TrN[VΓ¯]=wpTrN[VΓp]+(1wp)TrN[VΓe] (with Γp/e normalized to unity). In a straightforward manner,18 this yields Fe(γ)wpFp(γ)+(1wp)Fe(γ), implying that (if wp > 0) Fe(γ)Fp(γ). In combination with Fe(γ)Fp(γ) [following from the definition of Fp/e and PN(γ)EN(γ)], this eventually yields the suggested equality Fe(γ)=Fp(γ) on PN1. It is exactly the hidden assumption wp > 0 which is not justified: As explained above, the minimizer Γ¯ already lies on the boundary of EN(γ). Even more importantly, according to a theorem from convex optimization,37Γ¯ is with probability one (i.e., for generic V) already extremal and even exposed. The application of Krein-Milman’s theorem is therefore rather meaningless, and the assumption wp > 0 is violated as long as the minimizer Γ¯ is not incidentally a pure state and thus Fe(γ)Fp(γ). To illustrate all those general aspects, we revisit in the following the Hubbard dimer.

As an orthonormal reference basis for the singlet spin sector underlying the Hubbard dimer, we choose |1=c1c1|0, |2=c2c2|0, and |3=[c1c2|0c1c1|0]/2, where |0⟩ denotes the vacuum. Expressing any singlet state Γ=i,j=13Γij|ij| with respect to that basis and restricting it (as usually in quantum chemistry) to real values, the 1RDM in spatial representation follows as (recall γij ≡ ⟨i ↑ | γ | j ↑⟩ = ⟨i ↓ | γ | j ↓⟩, i, j = 1, 2)

γ11=1γ22=Γ11+12Γ33,γ12=γ21=12(Γ13+Γ23).
(6)

The set E2(γ) can thus be parameterized by three independent real variables. We choose (Γ11, Γ12, Γ13) and find for the expectation value of the Hubbard interaction Tr2[VΓ]=U(Γ11+Γ22)=12U1+Γ112γ11 where Eq. (6) and the normalization of Γ have been used. For two exemplary γP21E21, we illustrate in Fig. 3 the respective sets E2(γ). Levy’s minimization of the Hubbard interaction V is illustrated as a set of black hyperplanes with the black normal vector corresponding to −V. For generic γ, there are only two pure states on the boundary of E2(γ), shown as blue dots (in the right figure, one of them is shown in black since it coincides with the minimizer of Tr2[V (·)]). All other points on the boundary turn out to be mixed states (see also the color scheme representing the purity 1 − Tr[Γ2]). Since almost all boundary points are exposed and describe mixed states (thus violating the assumption wp > 0 in Ref. 18) one may now even wonder why the functionals Fp and Fe do not differ almost everywhere on P21. The answer to this is the following: By choosing an interaction V [i.e., a “direction” in E2(γ)] at random, pure states appear with finite probability as minimizers of Tr2[V (·)]. This is due to the fact (see Fig. 3) that each pure state has a whole range of supporting hyperplanes (see also point ΓC in Fig. 2), whose normal vectors cover a finite angular range.

Strongly inspired by Lieb’s seminal work38 on DFT for Coulomb systems, we resort to convex analysis, particularly to the concept of convex conjugation, to relate pure (Fp) and ensemble functionals (Fe) for arbitrary interaction V. The conjugate f* (also called Legendre-Fenchel transform) of a function f:RnR{±} is defined as f*(y)=supxRny,xf(x). Allowing f to take infinite values “has the advantage that technical nuisances about effective domains can be suppressed almost entirely”35 and we therefore extend Fp and Fe to the respective Euclidean space of hermitian matrices by defining Fp(γ)= and Fe(γ)= for γ outside their original domains PN1 and EN1, respectively. By referring to those common definitions and identifying Tr1[] as the inner product ⟨h, γ⟩ on the Euclidean space of hermitian matrices, we make the crucial observation that the energy E(h) [recall Eqs. (1) and (2)] is nothing else than the conjugate of the universal functionals Fp and Fe, respectively (up to an overall minus sign and a reflection h ↦ −h). The conjugation and thus the minimizations in (1) and (2) have a clear geometric meaning as it is illustrated in Fig. 4: For any fixed “normal vector” h, one considers the respective hyperplanes in the Euclidean space of vectors (γ, μ) (γ is a hermitian matrix, μR) defined by μ = ⟨h, γ⟩ + u and determines the largest u such that the upper closed halfspace of the respective hyperplane still contains the entire graph of Fp/e. E(h) then follows as the intercept of that hyperplane with the F-axis, i.e., the maximal u. This interpretation of the conjugation, in particular, explains in a geometric way why some 1RDMs γ (such as those on the line segment between γA and γB in Fig. 4) are not pure v-representable.7,8,20,39 Moreover, it shows that replacing Fp in (1) by the lower convex envelop,35conv(F), would not change the outcome of the minimization.

The second ingredient required for relating Fp and Fe is a theorem from convex analysis stating35 that the biconjugate f** coincides with f whenever f is convex and lower semicontinuous (a weaker form of continuity). Moreover, for arbitrary f, f** is (the closure of) the lower convex envelop of f. It is straightforward to apply those mathematical results to the functionals Fp, Fe, and E: First, since Fe is convex, it is continuous in the interior of EN1. This immediately implies35 lower semicontinuity (except for γEN1). We assume in the following that Fe is also lower semicontinuous on the boundary EN1 of EN1. The latter seems to be particularly difficult to verify (also since the interaction V is arbitrary). In case this assumption turns out to be wrong, our final result (7) will be valid on the interior of EN1 only (which does not reduce at all the significance and scope from any practical point of view). According to the theorem mentioned above, Fe therefore coincides with its biconjugation. Furthermore, the biconjugate of Fp coincides with its lower convex envelop conv(Fp). Yet, since both functionals Fe=Fe** and conv(Fp)=Fp** follow as the conjugate of the same functional, namely, the energy E (up to the minus signs), we eventually obtain (see also Fig. 5)

Feconv(Fp).
(7)

It is particularly remarkable that the pure functional Fp determines the ensemble functional Fe on its whole domain EN1, despite the fact that Fp’s effective domain PN1 is a proper subset of EN1. To be more specific, (7), namely, states that Fe(γ) follows as the minimisation of iwiFp(γi) with respect to all possible convex decompositions γ=iwiγi(0wi1,iwi=1) involving only 1RDMs γi from PN1, Fe(γ)=miniwiFp(γi)iwiγi=γ,γiPN1. This is also illustrated on the right panel of Fig. 5 for a γ outside PN1, also emphasizing the important fact that the extremal points of PN1 and EN1 coincide.

A fundamental theorem in RDMFT suggested that the pure (Fp) and ensemble functionals (Fe) would coincide on their common domain PN1 of pure N-representable 1RDMs. Based on a comprehensive study of the geometric picture of density matrices underlying Levy’s constrained search, we have refuted this crucial theorem. By exploiting concepts from convex analysis, we have then shown that Fe follows instead as the lower convex envelop of Fp. This relation [see Eq. (7)] which holds for any interaction V is particularly remarkable: The pure functional Fp together with PN1 determines the ensemble functional Fe on its whole domain EN1, despite the fact that Fp’s domain PN1 is a proper subset of EN1. This latter point in conjunction with the refutation of the relation FpFe|PN1 demonstrates that relaxing pure RDMFT to ensemble RDMFT does not necessarily circumvent the complexity of the one-body pure N-representability conditions. Instead, it may simply be transferred from the underlying space of pure N-representable one-matrices into the structure of the universal one-matrix functional Fe. In that case, an additional conceptual insight would follow: Approximating the universal functional would have at least the same computational complexity as the problem of determining all generalized Pauli constraints. Moreover, taking the generalized Pauli constraints into account may facilitate the development of more accurate functionals.

We are very grateful to E. J. Baerends, O. Gritsenko, D. Kooi, N. N. Lathiotakis, M. Piris, and particularly also K. J. H. Giesbertz for inspiring and helpful discussions. C.S. acknowledges financial support from the UK Engineering and Physical Sciences Research Council (Grant No. EP/P007155/1).

1.
T. L.
Gilbert
,
Phys. Rev. B
12
,
2111
(
1975
).
2.
J.
Cioslowski
,
Many-Electron Densities and Reduced Density Matrices
(
Springer Science & Business Media
,
2000
).
3.
M.
Piris
, “
Natural orbital functional theory
,” in
Reduced-Density-Matrix Mechanics: With Application to Many-Electron Atoms and Molecules
, edited by
D. A.
Mazziotti
(
Wiley-Blackwell
,
2007
), Chap. 14, p.
387
.
4.
K.
Pernal
and
K. J. H.
Giesbertz
, “
Reduced density matrix functional theory (RDMFT) and linear response time-dependent rdmft (TD-RDMFT)
,” in
Density-Functional Methods for Excited States
, edited by
N.
Ferré
,
M.
Filatov
, and
M.
Huix-Rotllant
(
Springer International Publishing
,
Cham
,
2016
), p.
125
.
5.
R.
Schade
,
E.
Kamil
, and
P.
Blöchl
,
Eur. Phys. J.: Spec. Top.
226
,
2677
(
2017
).
6.
P.
Hohenberg
and
W.
Kohn
,
Phys. Rev.
136
,
B864
(
1964
).
7.
R. G.
Parr
and
W.
Yang
,
Annu. Rev. Phys. Chem.
46
,
701
(
1995
).
8.
E.
Gross
and
R.
Dreizler
,
Density Functional Theory
(
Springer Science & Business Media
,
2013
), Vol. 337.
9.
R. O.
Jones
,
Rev. Mod. Phys.
87
,
897
(
2015
).
10.
N. N.
Lathiotakis
,
N.
Helbig
, and
E. K. U.
Gross
,
Phys. Rev. B
75
,
195120
(
2007
).
11.
N. N.
Lathiotakis
and
M. A. L.
Marques
,
J. Chem. Phys.
128
,
184103
(
2008
).
12.
N.
Schuch
and
F.
Verstraete
,
Nat. Phys.
5
,
732
(
2009
).
13.
M.
Levy
,
Proc. Natl. Acad. Sci. U. S. A.
76
,
6062
(
1979
).
14.
A.
Klyachko
,
J. Phys.: Conf. Ser.
36
,
72
(
2006
).
15.
M.
Altunbulak
and
A.
Klyachko
,
Commun. Math. Phys.
282
,
287
(
2008
).
16.
A.
Klyachko
, e-print arXiv:0904.2009 (
2009
).
17.
S. M.
Valone
,
J. Chem. Phys.
73
,
1344
(
1980
).
18.
T. T.
Nguyen-Dang
,
E. V.
Ludena
, and
Y.
Tal
,
Comput. Theor. Chem.
120
,
247
(
1985
).
19.
M.
Saubanère
and
G. M.
Pastor
,
Phys. Rev. B
84
,
035111
(
2011
).
20.
A. J.
Cohen
and
P.
Mori-Sánchez
,
Phys. Rev. A
93
,
042511
(
2016
).
21.
K.
Yasuda
,
Phys. Rev. A
63
,
032517
(
2001
).
22.
G.
Zumbach
and
K.
Maschke
,
J. Chem. Phys.
82
,
5604
(
1985
).
23.
A. J.
Coleman
,
Rev. Mod. Phys.
35
,
668
(
1963
).
24.
H.
Kummer
,
J. Math. Phys.
8
,
2063
(
1967
).
25.
R. M.
Erdahl
,
J. Math. Phys.
13
,
1608
(
1972
).
26.
E.
Davidson
,
Reduced Density Matrices in Quantum Chemistry
(
Elsevier
,
2012
), Vol. 6.
27.
J. E.
Harriman
,
Phys. Rev. A
17
,
1249
(
1978
).
28.
J. E.
Harriman
,
Phys. Rev. A
17
,
1257
(
1978
).
29.
M.
Hübner
,
Phys. Lett. A
163
,
239
(
1992
).
30.
D.
Petz
and
C.
Sudár
,
J. Math. Phys.
37
,
2662
(
1996
).
31.
D. C.
Brody
,
J. Phys. A: Math. Theor.
44
,
252002
(
2011
).
32.
S. A.
Ocko
,
X.
Chen
,
B.
Zeng
,
B.
Yoshida
,
Z.
Ji
,
M. B.
Ruskai
, and
I. L.
Chuang
,
Phys. Rev. Lett.
106
,
110501
(
2011
).
33.
D. A.
Mazziotti
,
Phys. Rev. Lett.
108
,
263002
(
2012
).
34.
J.
Chen
,
Z.
Ji
,
M. B.
Ruskai
,
B.
Zeng
, and
D.-L.
Zhou
,
J. Math. Phys.
53
,
072203
(
2012
).
35.
R. T.
Rockafellar
,
Convex Analysis
(
Princeton University Press
,
1997
).
36.
M.
Krein
and
D.
Milman
,
Stud. Math.
9
,
133
(
1940
).
37.
J.
Bolte
,
A.
Daniilidis
, and
A. S.
Lewis
,
Math. Oper. Res.
36
,
55
(
2011
).
38.
E. H.
Lieb
,
Int. J. Quantum Chem.
24
,
243
(
1983
).
39.
D.
Van Neck
,
M.
Waroquier
,
K.
Peirs
,
V.
Van Speybroeck
, and
Y.
Dewulf
,
Phys. Rev. A
64
,
042512
(
2001
).