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.

## I. INTRODUCTION

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 systems^{4} 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 DFT^{12} 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 $HN\u2261\u2227N[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 functional^{13} then follows immediately by determining the ground state energy of *H*

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 proposed^{17} to relax in (1) the set $PN$ to $EN$ by skipping the purity, leading to

with the ensemble functional $Fe(\gamma )\u2261minEN\u220b\Gamma \u21a6\gamma TrN[V\Gamma ]$ 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, $Fp\u2261Fe|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

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

## II. AN INSTRUCTIVE EXAMPLE: HUBBARD DIMER

The simplest way to refute the suggested equality $Fp\u2261Fe|PN1$ is to find one counterexample. A simple one is given by the *asymmetric* Hubbard dimer

a system of two electrons on two sites. Here, $ci\sigma \u2020$(*c*_{iσ}) denotes the creation(annihilation) operator of an electron at site *i* with spin *σ* and $ni\sigma =ci\sigma \u2020ci\sigma $ 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 exercise^{19,20} to determine the respective pure functional $Fp$

with *γ*_{ij} ≡ ⟨*i ↑* | *γ* | *j ↑*⟩ = ⟨*i ↓* | *γ* | *j ↓*⟩, *i*, *j* = 1, 2. $Fp(\gamma )$ is invariant under *γ*_{11} → (1 − *γ*_{11}) (particle-hole duality^{21}) 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 condition^{19,20} $(\gamma 11\u221212)2+(\gamma 122\u226414)$]. On the other hand, it is well-known^{4,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 $Fe\u2261Fp$ on $P21$.

## III. GEOMETRIC PICTURE OF LEVY’S CONSTRAINED SEARCH

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(\gamma )$ and $PN(\gamma )$ with respect to which the minimization (3) is carried out: While both sets $EN(\gamma )$ and $PN(\gamma )$ are also compact and $EN(\gamma )$ is convex for all *γ*, extremal states Γ of $EN(\gamma )$ 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(\gamma )$ despite the fact that it is not extremal within $EN$.

As already stated above, the ensemble functional $Fe(\gamma )$ follows for each $\gamma \u2208EN1$ by minimizing Tr_{N}[*V* Γ] over $EN(\gamma )$. Since Tr_{N}[*V* (·)] is linear and $EN(\gamma )$ is convex and compact, the minimum (i.e., $Fe(\gamma )$) is attained on the boundary of $EN(\gamma )$. This is a general (and rather obvious) fact from linear optimization: First, we observe that Tr_{N}[*V* Γ] is nothing else than the standard inner product on the Hilbert space of hermitian matrices, Tr_{N}[*V* Γ] ≡ ⟨*V*,Γ⟩_{N}. In that sense, there is given a notion of geometry on the space of density operators^{23–34} and thus *V* defines a direction in $EN(\gamma )$. The set of $\Gamma \u2208EN(\gamma )$ with a specific interaction energy $v$ = ⟨*V*,Γ⟩_{N} gives rise to a hyperplane, orthogonal to *V*. The minimum of Tr_{N}[*V* (·)] ≡ ⟨*V*,·⟩_{N} on $EN(\gamma )$ then follows by shifting the hyperplane along the direction −*V* (i.e., by reducing $v$) until an extremal point of $EN(\gamma )$ is reached. This final hyperplane is a so-called *supporting hyperplane*.^{35} By definition, this means that $EN(\gamma )$ is entirely contained in one of the two closed half-spaces bounded by that hyperplane and $EN(\gamma )$ 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 *exposed*^{35} 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 Tr_{N}[*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 $\gamma \u2208PN1$, the minimizer of Tr_{N}[*V* Γ] on $EN(\gamma )$ is denoted by $\Gamma \xaf$. Since $EN(\gamma )$ is convex and compact, $\Gamma \xaf$ can be expressed according to the Krein-Milman theorem^{36} as a convex combination of the extreme points of $EN(\gamma )$. This convex combination can be grouped into two parts, one ($wp\Gamma p$) arising from pure extremal states and one ($(1\u2212wp)\Gamma e$) arising from mixed states, i.e., $Fe(\gamma )=TrN[V\Gamma \xaf]=wpTrN[V\Gamma p]+(1\u2212wp)TrN[V\Gamma e]$ (with Γ_{p/e} normalized to unity). In a straightforward manner,^{18} this yields $Fe(\gamma )\u2265wpFp(\gamma )+(1\u2212wp)Fe(\gamma )$, implying that (if $wp$ > 0) $Fe(\gamma )\u2265Fp(\gamma )$. In combination with $Fe(\gamma )\u2264Fp(\gamma )$ [following from the definition of $Fp/e$ and $PN(\gamma )\u2282EN(\gamma )$], this eventually yields the suggested equality $Fe(\gamma )=Fp(\gamma )$ on $PN1$. It is exactly the hidden assumption $wp$ > 0 which is not justified: As explained above, the minimizer $\Gamma \xaf$ already lies on the boundary of $EN(\gamma )$. Even more importantly, according to a theorem from convex optimization,^{37} $\Gamma \xaf$ 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 $\Gamma \xaf$ is not *incidentally* a pure state and thus $Fe(\gamma )\u2260Fp(\gamma )$. 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=c1\u2191\u2020c1\u2193\u2020|0$, $|2=c2\u2191\u2020c2\u2193\u2020|0$, and $|3=[c1\u2191\u2020c2\u2193\u2020|0\u2212c1\u2193\u2020c1\u2191\u2020|0]/2$, where |0⟩ denotes the vacuum. Expressing any singlet state $\Gamma =\u2211i,j=13\Gamma 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)

The set $E2(\gamma )$ 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\Gamma ]=U(\Gamma 11+\Gamma 22)=12U1+\Gamma 11\u22122\gamma 11$ where Eq. (6) and the normalization of Γ have been used. For two exemplary $\gamma \u2208P21\u2261E21$, we illustrate in Fig. 3 the respective sets $E2(\gamma )$. 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(\gamma )$, shown as blue dots (in the right figure, one of them is shown in black since it coincides with the minimizer of Tr_{2}[*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(\gamma )$] at random, pure states appear with finite probability as minimizers of Tr_{2}[*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.

## IV. RELATING PURE AND ENSEMBLE FUNCTIONAL

Strongly inspired by Lieb’s seminal work^{38} 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:Rn\u2192R\u222a{\xb1\u221e}$ is defined as $f*(y)=supx\u2208Rn\u27e8y,x\u27e9\u2212f(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(\gamma )=\u221e$ and $Fe(\gamma )=\u221e$ for *γ* outside their original domains $PN1$ and $EN1$, respectively. By referring to those common definitions and identifying Tr_{1}[*hγ*] 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, $\mu \u2208R$) 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*,^{35} $conv(F)$, would not change the outcome of the minimization.

The second ingredient required for relating $Fp$ and $Fe$ is a theorem from convex analysis stating^{35} 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 implies^{35} lower semicontinuity (except for $\gamma \u2208\u2202EN1$). We assume in the following that $Fe$ is also lower semicontinuous on the boundary $\u2202EN1$ 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)

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(\gamma )$ follows as the minimisation of $\u2211iwiFp(\gamma i)$ with respect to all possible convex decompositions $\gamma =\u2211iwi\gamma i(0\u2264wi\u22641,\u2211iwi\u2009=\u20091)$ involving only 1RDMs *γ*_{i} from $PN1$, $Fe(\gamma )=min\u2009\u2211iwiFp(\gamma i)\u2211iwi\gamma i=\gamma ,\u2009\gamma i\u2208PN1$. 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.

## V. SUMMARY AND CONCLUSION

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 *F*_{p}’s domain $PN1$ is a proper subset of $EN1$. This latter point in conjunction with the refutation of the relation $Fp\u2261Fe|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.

## ACKNOWLEDGMENTS

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).