Repulsively diverging gradient of the density functional in the Reduced Density Matrix Functional Theory

The Reduced Density Matrix Functional Theory (RDMFT) is a remarkable tool for studying properties of ground states of strongly interacting quantum many body systems. As it gives access to the one-particle reduced density matrix of the ground state, it provides a perfectly tailored approach to studying the Bose-Einstein condensation or systems of strongly correlated electrons. In particular, for homogeneous Bose-Einstein condensates as well as for the Bose-Hubbard dimer it has been recently shown that the relevant density functional exhibits a repulsive gradient (called the Bose-Einstein condensation force) which diverges when the fraction of non-condensed bosons tends to zero. In this paper, we show that the existence of the Bose-Einstein condensation force is completely universal for any type of pair-interaction and also in the non-homogeneous gases. To this end, we construct a universal family of variational trial states which allows us to suitably approximate the relevant density functional in a finite region around the set of the completely condensed states. We also show the existence of an analogous repulsive gradient in the fermionic RDMFT for the $N$-fermion singlet sector in the vicinity of the set of the Hartree-Fock states. Finally, we show that our approximate functional may perform well in electron transfer calculations involving low numbers of electrons. This is demonstrated numerically in the Fermi-Hubbard model in the strongly correlated limit where some other approximate functionals are known to fail.


Introduction
The Bose-Einstein condensate (BEC) has been theoretically predicted by Bose and Einstein [1,2] in 1924. Seventy years passed before it was finally realised experimentally by Wieman and Cornell and independently by Ketterle in gases of ultracold atoms [3,4]. Nowadays, the BEC provides a widely-used experimental testing ground for quantum many-body theories, especially superfluidity, superconductivity [7,8,9] and the laser phenomenon [6,5]. In this work, we focus on the description of the BEC via the Density Functional Theories. Our developed methods for the BEC are subsequently extended to systems of strongly interacting electrons in the singlet sector.
The most successful and renowned functional theory is perhaps the one which uses the Hohenberg-Kohn density functional [13] (see also [14]). The natural variable of the Hohenberg-Kohn density functional is the single-particle density which unfortunately is not suitable for the purpose of this paper [15]. One of the reasons is that the singleparticle densities of the fully condensed quantum states can range through all possible densities. In other words, it is not possible to decide whether a given single-particle density describes a BEC condensate or not. Thus, one needs to use a more sophisticated generalization of the Hohenberg-Kohn density functional, namely the Reduced Density Matrix Functional Theory (RDMFT) due to Gilbert [16] and Levy [17]. RDMFT uses the single-particle reduced density matrix as its natural variable and hence is able to recover quantum correlations exactly. All density functional theories suffer the issue that their respective functionals are extremely challenging to find and hence they are usually not known. Remarkable exceptional cases in RDMFT where the functional is actually known analytically include the Bose-Hubbard dimer occupied by two bosons [18] and the singlet sector of the Fermi-Hubbard dimer for two spinful electrons [19,20]. The common practical strategy is thus to find suitable approximations to the relevant density functional. In RDMFT it is often the case that approximate functionals for some specific systems often lead to progress in finding more complicated approximate functionals for more general systems. For instance, in the fermionic RDMFT this has been the case with the well-known Hartree-Fock functional [28] which in turn led to the discovery of the Müller and Power functionals [29,30] which in turn inspired numerous more complicated hierarchies of functionals [31,32]. It is also a crucial task to construct approximate functionals that are simple enough to make explicit calculations tractable, but accurate enough to capture some desired information about the exact functional. Our presented work makes effort in this direction. Our results build on the recent tremendous progress in the foundations of both bosonic and fermionic RDMFT [22,21,18,23], especially in RDMFT for bosons and homogeneous Bose-Einstein condensates, where the remarkable concept of the repulsive BEC force has been introduced. In particular, it has been shown that the gradient of the relevant density functional in the space of the reduced density matrices diverges repulsively in the regime of the Bose-Einstein condensation as the inverse of the square root of the fraction of the non-condensed bosons. This provided an alternative and more fundamental explanation for the existence of the quantum depletion in the homogeneous condensate [23] (see also [25,26,27] for the underlying theory) and in the Bose-Hubbard dimer [18]. This is compared with the celebrated theories by Bogoliubov [10] and Gross-Pitaevskii [11] which work in the dilute or weakly interacting regime as well as in the the high density regime for charged bosons and that have been recently verified experimentally [12].
Our work shows that the repulsive BEC force is present also in non-homogeneous Bose-Einstein condensates, showing that the presence of the BEC force is an universal phenomenon making the above alternative explanation of the quantum depletion complete for any system of pairwise interacting bosons. Our arguments are remarkably versatile and they can be extended to other physical systems. In particular, similar methods are employed to prove the existence of an analogous phenomenon for the singlet fermionic systems. Our construction coincides with the exact functional in the dimer case (two fermions occupying two sites/orbitals) known from previous works [19]. We numerically test the accuracy of our approximate functional in the Fermi-Hubbard model. For two fermions (N = 2) we show that the approximate functional gives very accurate values of the ground state energy. For the Fermi-Hubbard chains, we employ the approximate functional for calculating the electron transfer. In particular, for the half-filled chain of length four, our functional recovers the exact electron transfer in the strongly correlated limit with a good accuracy. Remarkably, this is in contrast with the Müller and Power functionals [29,30] which are known to fail even in the dimer case [19]. This shows that our proposed functional or its generalizations (see Section 8) may be useful in electron transfer calculations, a task which has been pointed out as one of the current challenges of the RDMFT [33,32].
For the sake of clarity, we would like to make precise what is the form of the most general bosonic hamiltonian with pairwise interactions that we consider. We work in the discrete setting and consider quantum systems of N bosons occupying d sites interacting according to hamiltonians of the form whereĥ is the single-particle term which comprises of hopping and local potentials whilê W is an arbitrary pair interaction. Thus, the most general form ofĥ iŝ where t i,j are (possibly complex) hopping amplitudes and v i are on-site potentials. Recall that the creation and annihilation operators b i , b † j satisfy the bosonic commutation rules where Ω ijkl are complex parameters.

A brief recap of the RDMFT
The setting of the Reduced Density Matrix Functional Theory is to fix the interaction termŴ in Equation (1) and consider a family of hamiltoniansĤ(h) =ĥ +Ŵ . One of the key objects in RDMFT is the one-particle reduced density matrix (denoted here shortly by 1RDM) which in the case of N bosons occupying d sites is a d × d hermitian matrix assigned to any state |Ψ whose entries read For any one particle hamiltonianĥ where h is simply the matrix of coefficients [h i,j ]. The set of all 1RDMs coming from reductions of the pure states (also called the N -representable 1RDMs) will be denoted by Γ, i.e.
Note that not every γ ∈ Γ is the 1RDM of the ground state ofĤ(h) for some h. Thus, it is also necessary to consider the set of so-called v-representable 1RDMs associated witĥ W which is the set Γ v ⊂ Γ of 1RDMs stemming from ground states of hamiltonianŝ H(h). Following Levy [17], Gilbert [16] and Lieb [14] one defines the RDMFT functional F(γ) defined on Γ v as follows. Consider the ground state energy ofĤ(h) as a function of h, E gs (h). For any γ ∈ Γ v we define Then, by the variational principle we have E gs (h) = min γ∈Γv (Tr (hγ) + F(γ)) .
The following two notable problems arise: i) we do not know what the set Γ v is and ii) we do not know what F(γ) is. Levy has proposed to circumvent this problem by extending the domain of the RDMFT functional and the domain of the search (5) to all (possibly non-physical) 1RDMs from Γ [17] and define the extended functionalF as where the minimization is done over all pure states whose 1RDM is γ. FunctionalF has the following two properties: P1)F(γ) + Tr (hγ) ≥ E gs (h) for every h and γ ∈ Γ, P2) F(γ gs ) + Tr (hγ gs ) = E gs (h), if γ gs is the 1RDM of the ground state ofĤ(h). Properties P1 and P2 imply that for γ ∈ Γ v we haveF(γ) = F(γ). Thus, functionals F as well asF are universal, i.e. they are independent of h. The do however depend on the interactionŴ . Finally, let us point out that by Equation (6) the functionalF(γ) is upper bounded by Ψ|Ŵ |Ψ for any |Ψ whose 1RDM is γ. We will strongly rely on this elementary fact in the following parts of the paper.

The repulsive gradient in the Bose-Hubbard dimer
In this section we introduce the main concepts and build some key intuitions that we will develop further in the remaining parts of the paper. A simple and intuitive explanation of the relevant repulsive gradient in the the RDMFT functional that we are going to focus on in this section comes from the Bose-Hubbard dimer. This is a system of two sites labelled by the numbers i = 1, 2 with the corresponding creation operators b † 1 , b † 2 and the interaction given byŴ BH =n 2 1 +n 2 2 , wheren i is the particle number operator at the site i. To simplify things further, we will assume that the system at hand involves only N = 2 bosons and consider only real wave functions. Then, the corresponding two-boson states can be written in the position basis as the superposition where the basis states are 3) the 1RDM of |Ψ is a 2 × 2 real symmetric matrix whose coefficients are given by Remarkably, the above relations can be inverted in order to express the wave function's coefficients in terms of the variables γ i,j leading to an exact expression for the Levy's RDMFT functional according to Equation (6). This fact has been studied in the literature -see [18] for the complete derivation and [19] for a similar treatment of the two-electron Fermi-Hubbard dimer. The resulting expression for the exact RDMFT functional reads Obtaining exact expressions for the RDMFT functionals is generally not possible for larger systems as the multi-particle wave functions are parametrised by many more variables than the respective 1RDMs which makes realising the minimization in Equation (6) in an exact way very difficult. According to Equation (8) the domain of the RDMFT functional for the N = 2 dimer is the disk of unit radius γ 2 1,2 + (γ 1,1 − 1) 2 ≤ 1. This fact comes solely from the condition that the 1RDM γ has to be positive-semidefinite. Importantly, the 1RDMs which lie on the boundary of the disk have the property that γ 2 = 2γ, i.e. they have the eigenvalues ν 1 = 2, ν 2 = 0. In other words, such 1RDMs come from the states which are completely condensed in a single mode |ψ being some linear combination of the spatial modes |1 and |2 . In this paper, we will be mainly interested in the behaviour of the functionalF BH (γ) when γ is close to the set of 1RDMs which come from such completely condensed states. In contrast to the dimer, in higher dimensions the 1RDMs of the completely condensed states form only a subset of the boundary of the domain of the RDMFT functional. However, in the simple example at hand we can visualise the exact functional (Fig. 1a) and see what exactly is happening in the vicinity of its boundary (Fig. 1c). Namely, we see that the gradient ofF BH (γ) diverges at the boundary of its domain (this is often called a cusp singularity). More precisely, we introduce the distance and angle variables, D ∈ [0, 1] and φ ∈ [0, 2π[ (Fig. 1b), which parametrise the disk such that which diverges as − sin 2 φ √ 2 1 √ D in the leading order when D → 0. This repulsive divergence of the RDMFT functional at the boundary of its domain has been confirmed to exist in many more physical systems including the N -boson Bose-Hubbard dimer, homogeneous Bose-Einstein condensates (where it is called the Bose-Einstein condensation force [18,23,24]) and translationally invariant one-dimensional fermionic systems [21] (where it is called the exchange force). It has been conjectured that it is a property of generic interacting many-body systems. In this paper, we show that this is indeed the case. However, as we are not able to write down the exact functional for a general system and simply compute its derivative, we need to find an appropriate bound for the RDMFT functional whose diverging gradient will imply that the gradient of the exact RDMFT functional diverges as well when D → 0. To this end, for any 1RDM γ we construct an ansatz state |Ψ(γ) such that γ is the 1RDM of the state |Ψ(γ) . With such an ansatz at hand, we simply bound the exact functional bỹ and show that the bound diverges repulsively when the (appropriately defined) distance variable D tends to 0. Let us also remark that the distance variable D has the interpretation as the depletion parameter. To see this, we find the eigenvalues of γ to be ν 1 = 2 − D and ν 2 = D. Thus, D is the number of bosons that are outside of the highest occupied mode.
Let us next show how this strategy works out for the N -boson Bose-Hubbard dimer. We assume that the 1RDMs are normalised to tr γ = N . This implies that the set of 1RDMs is a disk of the radius N/2 given by the inequality (γ 1,1 − N/2) 2 + γ 2 1,2 ≤ N 2 /4. Thus, we use the parametrisation of γ(D, φ) which reads We will construct the ansatz state |Ψ(γ(D, φ)) ≡ |Ψ(D, φ) in the basis of eigenmodes The ansatz is the following superposition of the completely condensed state in the mode |φ with the doubly-excited state It is straightforward to verify that the state |Ψ(D, φ) gives the correct 1RDM (9) which is diagonal in the |φ , |φ ⊥ basis and whose diagonal entries in this basis are N − D and D. By evaluating the expression Ψ(D, φ)|Ŵ |Ψ(D, φ) we get the following upper bound for the RDMFT functional for any D ≤ 2 The upper bound is equal to the exact functional at D = 0 and its gradient diverges repulsively when D → 0 proportionally to −1/ √ D in the leading order which implies the desired result. It is however crucial to be assured that is generically strictly positive. Otherwise, the entire argument breaks down as there are no other terms in the bound that give the repulsive gradient. This is done by a straightforward but rather tedious calculation which relies on expressing the interactionŴ BH =n 2 1 +n 2 2 in terms of the creation and annihilation operators This way, we find that thus the above overlap vanishes only when φ ∈ {0, π}. A similar calculation allows us to obtain expressions for the remaining expectation values ofŴ BH . The resulting upper bound for the RDMFT functional is plotted for N = 4 and N = 10 in Fig. 2. It is given by the formulã Finally, the resulting bound for the gradient of the functional at D → 0 reads The above bound for the gradient can be compared with the results of [18] where the first term of the expansion of the gradient in the powers of D has been obtained. The results agree with Equation (11) up to some factors of N which come from the different choices of the normalisation of 1RDMs.
To end this section, let us remark that analytic bounds for the exact RDMFT functional are potentially very useful tools in approximately solving complex physical systems. This is because the functional is universal for all the single-particle termŝ h in Equation (1). For the Hubbard model the single-particle terms include the external on-site potentials and all the possible hopping terms. Thus, quite remarkably, an approximate functional can be applied to the Hubbard model of any geometry. Moreover, finding the ground state energy boils down to the easy task of finding the minimum of a surface (just like the surfaces on Fig. 2 and Fig. 1a, but generally of a higher dimension) which is tilted by the single-particle terms of the Hamiltonian. In particular, for the dimer this boils down to the almost trivial task of finding the minimum of a two-dimensional surface -independently of the number of bosons N .

The spectral simplex ∆ N,d
This section is devoted to explaining some key geometric properties of the set Γ defined in Equation (4). One of the main difficulties in RDMFT is the need of optimization over set Γ which is in general a complicated non-convex set. However, in contrast to its fermionic counterpart [37,40,41,42], set Γ for bosons is convex and has a tractable description in terms of a geometric object which we call the spectral simplex and denote by ∆ N,d . Before we define ∆ N,d , we need to take a closer look at a particular set of operations called the single-particle unitaries. The single-particle unitaries are unitary operators generated by single-particle hermitian operators, i.e. they are operations of Hence, with a given matrix A we have associated a single-particle unitaryÛ and a d × d unitary U . WhileÛ acts on N -boson wave functions, matrix U acts on 1RDMs via conjugation. Importantly, these two actions are equivalent, namely In particular, any 1RDM γ can be diagonalised by the above operations. This means that set Γ is parametrised by the spectra of 1RDMs coming from pure states of N bosons. Note that the spectrum of any 1RDM (n 1 , n 2 , . . . , n d ) has to satisfy 0 ≤ n i ≤ N and n 1 + . . . + n d = N.
These constraints define the spectral simplex ∆ N,d (see Fig 3). Simplex ∆ N,d is the is a vector of real numbers whose squares sum up to one. The 1RDM of |Ψ α is diagonal and it reads Thus, by changing the parameters α we can reach any convex combination of simplex's vertices. Stated another way, the above fact means that any d × d density matrix γ (i.e. a positive-definite hermitian matrix of trace N ) is the 1RDM of a pure bosonic state of the formÛ |Ψ α for some single-particleÛ and some α.
For a given typical d × d density matrix γ there exist many N -boson pure quantum states that have such a γ as their 1RDM. However, this changes diametrically when we consider a special subset of Γ which is defined as the set of 1RDMs of completely condensed states of N bosons. All such completely condensed states are of the formŝ U N (1) whereÛ is a single-particle unitary. By the property (12) the 1RDM of a completely condensed state is of the form U γ BEC U † , where

Thus, we define
Importantly, set Γ BEC is in a one-to-one correspondence with the set of completely condensed states and the correspondence is U γ BEC U † ↔Û N (1) . Hence, when computing Levy's functional (6) no minimization is required and we havẽ

The exposition of the main result
The proposed variational ansatz state is a combination of the completely condensed state N (1) and a weighted superposition of states with two bosons outside the condensate. Its precise normalised form reads where the variational parameters are ≥ 0 and σ = ±1. Numbers {α i } d−1 i=1 are arbitrary fixed numbers whose role is to probe all directions inside the simplex. They are all nonnegative and satisfy α 2 1 + . . . + α 2 d−1 = 1. Parameter determines the Euclidean distance of γ := γ (|Ψ ) from the 1RDM of the completely condensed state, γ BEC = diag (N, 0, . . . , 0). To see this, recall that the distance between two matrices is given by the Hilbert-Schmidt measure It is straightforward to see that γ is diagonal and its precise form reads are in a one-to-one correspondence with the possible directions from vertex γ BEC to the interior of the simplex ∆ N,d (see Fig. 4). A straightforward calculation shows that where β d := 2 1 + α 4 1 + . . . + α 4 d−1 . Note that the number of non-condensed bosons in γ is exactly N − N BEC = 2 2 /(1 + 2 ), thus D has the physical interpretation as the depletion parameter Note that distance D is invariant with respect to conjugating its arguments by a d × d unitary U . Hence, D is also equal to the distance from γ Û |Ψ = U γ U † to the 1RDM of the corresponding completely condensed state γ Û |N 1 = U γ BEC U † (see Fig. 5). The ansatz state (15) provides an upper bound for the RDMFT functional of the formF The above bound is defined for any 1RDM γ in the finite region around the set Γ BEC consisting of the 1RDMs that satisfy the condition N − N BEC ≤ 2. This is because any such γ can be written as U γ U † for some U and γ that are determined via the diagonalisation of γ. What is more, the matrix U γ BEC U † minimises the distance from γ to the set Γ BEC . The above upper bound has a particularly simple form, namely where constants c 1 , as well as the precise form ofŴ and the unitary U . The expressions for c 1 and c 2 can be viewed as the following combinations of the two-body integrals The two-body integrals can be obtained as closed-form expressions if d = 2, while for d > 2 they are computable as convergent series. Importantly, in the following sections we argue that the coefficient c 1 is generically nonzero. Using Equation (17) this yields the following upper bound for the functional in terms of the distance D < β d The final key observation is that the upper bound in Equation (19) is exact for = 0. Hence, while approaching the completely condensed state by taking << 1 (and hence D << 1), the gradient of the upper bound with respect to the distance D is also an upper bound for the gradient of the exact functional (see Fig. 6). Furthermore, the binary parameter σ can be chosen so that the gradient of the expression (19) is repulsive for small D and diverges to −∞ as − |c 1 | D in the leading order, a result anticipated in earlier works [18,23].

Calculating the coefficient c 1
In this section we provide arguments for the crucial fact that thanks to the specific form of our proposed ansatz state the coefficient c 1 in the bound (19) is generically non-zero.
Although our proof of the bound (19) is elementary and requires only minimal assumptions about the hamiltonian, it becomes heavy with calculations in its most general instance. However, most of the key arguments appear already when considering a dimer, d = 2. Hence, it will be most instructive to gradually increase the generality of our considerations, beginning with the Bose-Hubbard dimer, then moving to a dimer with an arbitrary pair interaction and finally formulating the proof in its full generality.

The Bose-Hubbard dimer
The interaction term for the Bose-Hubbard hamiltonian on d sites readŝ with V ∈ R being a real parameter. V > 0 corresponds to an on-site repulsion, while V < 0 corresponds to an on-site attraction. The most difficult part of the proof of bound (19) relies on computing U †Ŵ BH U . The general strategy is to express the interaction term in terms of generators of the unitary group and then use some wellknown facts concerning conjugation of the generators by unitary matrices (in particular, the correspondence with rotations in (d 2 − 1)-dimensions [34]). Let us next show how this strategy is realized in the case of the Hubbard dimer whereŴ BH = V (n 2 1 +n 2 2 ). Our first goal is to expressŴ BH in terms of the SU (2) generatorŝ In the following, we will often arrange the above generators into a vectorR = R 1 ,R 2 ,R 3 , whereR 1 =X,R 2 =Ŷ ,R 3 =Ẑ. In the above uniform notation, the generators satisfy the standard SU (2) commutation relations, i.e. [R i ,R j ] = 2i ijkRk , where ijk is the Levi-Civita symbol (not to be confused with the real parameter ).
With the above notation established, the Bose-Hubbard interaction term can be written asŴ whereN is the total particle number operator. For d > 2, we analogously definê Z ij =n i −n j to obtain Let us next take a slight detour and revisit some algebraic properties of conjugation in d = 2. By definition, any unitary operator for d = 2 can be written in terms of generators (20) asÛ = e iαr·R , where r is a vector in three dimensions of unit length. Recall the following key formula for the conjugation in d = 2 Matrix [C i,j ](α, r) is real and is well-known to describe a rotation by angle 2α about the axis r [34]. The exact expressions for its entries can be found e.g. in [35].
In particular, for the Bose-Hubbard dimer, we have The last line of the above equation follows from the fact that Ψ ,σ |R jRi |Ψ ,σ = 0 whenever j = 3 and i = 3 or i = 3 and i = 3. This is because |Ψ ,σ is a superposition of |N 1 and b 2 0 b † 1 2 |N 1 , so there have to be two different annihilations inR jRi to get a nonzero expectation value. Moreover, |Ψ ,σ is a real vector, thus Ψ ,σ | (XŶ +ŶX) |Ψ ,σ = 0. Next, we check the following facts by a straightforward calculation where we used the shorthand notation |K, L for the normalised state |vac . The above calculation is key for our proof. Namely, computing expectation values ofR 2 i generates terms proportional to due to the fact that squares of the form √ N ± σ 2(N − 1) 2 appear on the way. Collecting all contributions proportional to and 2 , we find the following constants from Equation (18).

The general case
Our previous considerations concerning an interacting dimer generalise fairly straightforwardly to d > 2. Firstly, we have more generators of the single-particle unitary operations. Their forms are completely analogous to the dimer case, namely for every pair 1 ≤ i < j ≤ d we define the following operators forming the SU (d) algebrâ Thus, we have 3d(d − 1)/2 generators. However, the number of linearly independent generators is d 2 − 1 due to the relationsẐ i,j +Ẑ j,k =Ẑ i,k . We will also arrange all generators in a linear order and use a single greek symbol to enumerate the generators, so that the set of linearly independent generators can be written in a uniform way as {R µ } d 2 −1 µ=1 . Thus, in full analogy to the dimer case, we can write a general pair interaction asŴ = µ≤νŴ µ,ν , wherê with ω µ,ν and τ µ,ν being some given real constants. ConjugatingR µ by a single-particle unitary of the formÛ = e iÂ withÂ = µ a µRµ results with the transformation where coefficients C µ,ν (A) are real and can be computed by taking the exponent of the so-called adjoint action ad ] describes a rotation in d 2 − 1 dimensions [34]. Thus, we get whereω γ,δ = µ≤ν ω µ,ν (C γ,µ C δ,ν + C γ,ν C δ,µ ) , Let us remark that the coefficients in Equation (26) and Equation (27) depend only on the Lie algebra element A, not on the particular representation. In particular, the same equations hold for the fermionic systems described in in Section 7. The rest of the proof relies on keeping track of contributions proportional to and 2 in expressions of the type Ψ ,σ |R γRδ |Ψ ,σ . The calculation is tedious, but straightforward, thus we move it to Appendix A.1. For the sake of completeness, below we write down the explicit expression for c 1 from Equation (18). As it turns out, the only contributions to c 1 come from expressions Ψ ,σ |R 2 γ |Ψ ,σ , whereR γ =X 1,k orR γ =Ŷ 1,k for some k. Thus, the result is completely analogous to the formula describing c 1 in the dimer case. Denote by ω where {α k } d−1 k=1 are parameters of the ansatz state (15).

Repulsive gradient in the fermionic RDMFT for the singlet sector
In this section we show that an analogous repulsive gradient force is present in systems of interacting spin-1/2 fermions for hamiltonians preserving the total spin which we assume to be zero by selecting the singlet sector. The fermions occupy d sites/orbitals with the corresponding annihilation operators {a i↑ , a i↓ } d i=1 . We start by defining the auxiliary operators f i,j which will be helpful in carrying the calculations in an almost complete analogy to the bosonic case.
The one-particle reduced density matrix is given by The single-particle operators form an SU (d) algebra and are represented on the singlet sector aŝ where i < j. As before, for any single-particle unitaryÛ the one-particle reduced density matrices transform equivariantly according to the formula γ Û |Ψ = U γ (|Ψ ) U † . We also consider the spectral polytope P M,d which is the polytope formed by the spectra of the one-particle reduced density matrices of the singlet states of N = 2M fermions occupying d orbitals (or sites). Importantly, the polytope P M,d has a tractable description [37] it is defined only by the Pauli constraints where n i = n i↑ +n i↓ is the number of fermions occupying the ith orbital. Clearly, P M,d is not a simplex when d > 3 or M > 1. It is a convex hull of d M vertices whose occupation numbers n i ∈ {0, 2} and sum up to 2M (see Fig. 7).
For 1 ≤ i 1 < i 2 < . . . < i M ≤ d we define the corresponding singlet Hartree-Fock state as The one-particle reduced density matrix corresponding to the above state is easily seen to be diagonal and to correspond to the convex combination of the polytope's vertices with the coefficients |α i 1 ,...,i M | 2 . Any hamiltonian conserving the total spin of 2M fermions occupying d orbitals when truncated to the singlet sector can be written aŝ where the interaction termŴ is an arbitrary linear combination of products of operators (29) with real coefficients (i.e. a real polynomial inX f i,j ,Ŷ f i,j andẐ f i,j ). In particular, the Fermi-Hubbard interaction is obtained asŴ For the attractive interaction (V i,j < 0) the Hamiltonian (30) with the Bose-Hubbard interaction and only real hopping terms (y i,j = 0) extended to the entire Fock space is known to have a unique singlet ground state [36]. For the repulsive interaction, the ground state is a singlet state only at the half-filling (N = d) and only when the graph defined by the hopping geometry is bipartite with the two disjoint sets of vertices having the same cardinality. This is satisfied for instance for the chain-geometry. In what follows, we will assume thatŴ is a pair interaction, i.e. it only consists of terms of degree 2 in the operatorsX f i,j ,Ŷ f i,j andẐ f i,j . We consider the following variational ansatz state.
where σ ∈ {−1, 1} and α i,j are real parameters such that M i=1 d j=M +1 α 2 i,j = 1. The one-particle reduced density matrix of the above ansatz state is diagonal and its entries γ (|Ψ ,σ ) k,k = Ψ ,σ |n k↓ +n k↑ |Ψ ,σ read Clearly, for = 0 the ansatz state is a Hartree-Fock state and γ Ψ f 0,σ = diag (2, . . . , 2, 0, . . . , 0) =: γ HF . If > 0, the parameters α i,j swipe through all possible directions coming out of the vertex γ HF and going inside the polytope P M,d (see Fig. 7). Consequently, for a fixed choice of the coefficients α i,j and for any single-particle unitarŷ U the parameter parametrizes the shortest path from the point U γ Ψ f ,σ U † to the set of the Hartree-Fock states Alternatively, one may think of the ansatz stateÛ Ψ f ,σ as having the form (31) where the annihilation operators {a i↑ , a i↓ } d i=1 are replaced by the natural orbital basis annihilators {ã i↑ ,ã i↓ } d i=1 ordered decreasingly in terms of their corresponding natural occupation numbers. The relevant distance measure is the Hilbert-Schmidt distance from Equation (16). In particular, for any single-particle unitaryÛ we have . By evaluating the sum of the d − M lowest natural occupation numbers of the stateÛ Ψ f ,σ , we obtain N − N HF = d l=M +1 γ (|Ψ ,σ ) l,l . By the number N HF we understand the total number of fermions occupying the M most occupied natural orbitals. Using the fact that the squares of the coefficients α i,j sum up to one, we get that N − N HF = 2 2 /(1 + 2 ) and thus the distance D has the interpretation as Our goal here is to show that the corresponding fermionic RDMFT functional is upper bounded in the vicinity of the set Γ HF bỹ By the reasoning presented in Section 5 this implies that the gradient of the RDMFT functional with respect to the distance D diverges repulsively as − D in the vicinity of the Hartree-Fock point U γ HF U † . To this end, we show that the following holds for any pair-interactionŴ . As before, the parameters c f 1 and c f 2 are linear combinations of the appropriate two-body integrals. For the most general interaction given by Equation (25)  . Then, where α k,l are the parameters of the ansatz state (31). Thus, the coefficient c f 1 is generically nonzero.

The form of the approximate functional
For any 1RDM γ the approximate functional is computed from the ansatz state (31) by using the one-particle basis of the natural orbitals corresponding to the decreasingly ordered spectrum (i.e. the natural occupation numbers) of γ. Denote by 2 ≥ ν 1 ≥ ν 2 ≥ . . . ≥ ν d ≥ 0 the decreasingly ordered natural occupation numbers of γ. Then, γ is a d × d matrix γ = d i=1 ν i |φ i φ i | and we define the operators {ã i↓ ,ã i↑ } d i=1 so that they refer to the annihilation operators corresponding to the natural orbitals {|φ i } d i=1 . According to the Equation (31), after slight modifications, the corresponding ansatz state reads where |HF (γ) : where the two-fermion integrals read Thus, the computation of G(γ) requires performing a convex minimization over the set ∆(ν). Importantly, the region of applicability of the approximate functional is given by the condition N HF ≥ N − 2. This is always satisfied when N ≤ 4, however when N > 4 the approximate functional is defined only on a subset of the set of all the Nrepresentable 1RDMs. Thus, we expect the approximate functional to perform best for low particle numbers.
In fact, this is the exact RDMFT functional whenever γ = diag (1,1). This is because any singlet state of two spin-1/2 fermions with d = 2 whose one-body reduced density matrix is not the identity matrix necessarily has the form (35) up to a relative phase, when written in the natural orbital basis. The expression (38) is equivalent to the expressions for the exact density functional studied in other works [19,18,43]. Another simplification arises when N = 2, regardless of d. Then, the coefficients {α 1,j } d j=2 follow immediately from the Equation (36) as

Electron transfer in the Fermi-Hubbard model
The repulsive Fermi-Hubbard chain is a good testing ground for our proposed approximate functional. The ground state of N electrons on the Fermi-Hubbard chain of the length d = 2N (half-filling) is known to be a singlet state [36]. The Hamiltonian readŝ In the calculations to follow, we will assume an open chain with i) uniform hopping terms, i.e. t i,i+1 ≡ −t with t ≥ 0 for all i, ii) uniform interaction strength, i.e. V i ≡ V for all i and iii) the on-site potentials v i that change by even . Then, up to a polynomial term in the total number of electrons, the hamiltonianĤ F H is easily seen to have the formĤ The exact RDMFT functional in the dimer case (d = 2) has been studied in [19] where its performance has been compared with other approximate density functionals. In particular, the authors of [19] studied the electron transfer problem which asks about how the diagonal entries of the 1RDM of the ground state for the Hamiltonian (39) depend on the on-site potential difference v between the endpoints of the chain. It has been shown that the Müller and Power approximate density functionals [29,30] give qualitatively incorrect results for the electron transfer, where the failure is most visible in the strong interaction limit (t >> V ). Other works show that approximate molecular density functionals generally fail to give the qualitatively correct description of the electron transfer and point out that this is an important challenge to improve the performance of the approximate functionals in this area [33,32]. In the following part of this section, we show that out proposed approximate density functional does indeed perform well in the electron transfer calculations of small Hubbard chains. In particular, we show that our functional gives accurate values of the ground state energy for the two-electron singlet systems and that show that it describes the electron transfer with a remarkable accuracy for the four-electron chain with four sites.
By applying the functional (37) to the two-electron (N = 2) Fermi-Hubbard chain, we observed numerically that for the uniform interaction and for all the sampled values of the parameters t ∈ [0, 1] and v ∈ [−2, 2] of the Hamiltonian (39), the approximate functional G agreed with the exact functional within the numerical accuracy. In other words, the ground state energy error |G(γ gs ) + tr hγ gs − E gs | was of the mean order 10 −14.5 (see Fig. 8a and Fig. 8b) which is the same as the error of the LU diagonalization algorithm which was used in the calculation of E gs . Thus, the approximate functional coincides with the exact functional for the Hubbard chain with the Hamiltonian (39) within the numerical accuracy range. We also tested other geometries of the Hubbard model. In particular, for the Fermi-Hubbard model on d = 6 sites with a uniform repulsive interaction we sampled 10 3 values of the hopping and the on-site potentials. There, the hopping and the on-site potentials were allowed to be arbitrary and they did not correspond to the chain geometry. The above parameters were sampled from the flat distributions t i,j ∈ [−1, 1] and v i ∈ [−2, 2]. For each set of the sampled parameters a minimization of the approximate density functional G(γ) has been performed. This showed the mean energy error of the order 10 −(13±3) (see Fig. 8c). This number is to be compared with the error of the LU diagonalization algorithm which was used in the calculation of E gs which was of the order 10 −15 or with the typical difference between the lowest and the highest eigenvalue of the random Hamiltonian which was of the order 10 1 .  (39). b) The absolute errors between E min and E gs from the plot a). c) The distribution of the ground state energy errors for the Fermi-Hubbard model with d = 6 sites with all the hopping amplitudes nonezro. The hopping amplitudes and the on-site potentials were randomly sampled 10 3 times. The resulting mean error was of the order 10 −(13±3) . The maximum error for the sample was of the order 10 −3 , but for 77% of the sample the error was of the order 10 −14 or less.
The above good results for the ground state energy of the two-electron systems are due to the fact that the excitations of the tested systems exhibit localization in the natural orbital basis making the superposition of doubly-excited Slater determinants in the ansatz state (31) an accurate guess for the ground state of the system. The ground states of systems of N > 2 electrons generally involve higher excitations and thus the correlation energy is expected to be recovered by our proposed functional only to a small extent. This can be seen already in the half-filled chain of length four (N = d = 4) studied in the remaining part of this section. However, as the numerical results show, the doubly-excited ansatz is enough to recover the electron transfer in the ground state of this system.
Thus, the minimization in Equation (37) is done simply over the interval containing the coefficient α 1,3 . The results for the ground state energy and the electron transfer are plotted in Fig. 9 and Fig. 10, where the interaction strength has been fixed to V = 1 due to the fact that the physically meaningful quantity is the ratio t/V . As we mentioned before, our approximate functional gives a qualitatively correct electron transfer with a remarkably good accuracy, especially in the limit of the strong interaction. On the other hand, the correlation energy is recovered only to a small extent as can be seen in Fig. 10.
Technical note: the numerical minimization of G(γ) over the set of 1RDMs has been done using the adaptive Nelder-Mead method [38,39] where the convergence treshold was the difference in function's arguments between iterations which was set to 10 −8 . Typically, the convergence was attained after about 10 4 iterations.   (39) computed using the exact diagonalization (red dots), the Hartree-Fock method (black pluses) and the approximate functional (blue crosses). The approximate functional generally does not recover the correlation energy. This is most visible in the strongly correlated limit (t = 0.05) for the values of v close to 0. There, the fraction of the recovered correlation energy is only about 70%.

Summary and discussion
Using the ansatz states (15) and (31) we have constructed universal upper bounds (19) and (32) for the bosonic and fermionic RDMFT functionals in a finite region around of the set of fully condensed (respectively the Hartree-Fock) states for any pair interaction. The upper bounds depend on the (generalised) depletion parameter which is the number of non-condensed bosons, N −N BEC , or the number of fermions which do not occupy the N/2 highest occupied natural orbitals, N − N HF . Because the upper bound coincides with the exact RDMFT functional when N BEC = N (all the bosons condensed) or N HF = N (system in a Hartree-Fock state), the gradient of the upper bound is also an upper bound for the gradient of the exact functional with respect to the depletion parameter. The gradient is shown to diverge repulsively when N BEC → N (N HF → N respectively). This provides an alternative explanation for the existence of the quantum depletion effect which uses only the very fundamental properties of the quantum theory: the variational principle and the geometry of the set of one-particle reduced density matrices.
We also tested our proposed approximate functional by computing the electron transfer in small Hubbard chains involving N = 2 and N = 4 electrons. We showed that in these systems the functional performs remarkably better that some other widely used approximate functionals, and thus our results shed some light on the possible ways of improving the performance of the density matrix functional methods in the electron transfer computations. In particular, the ansatz states (15) and (31) take into account only double excitations in the natural orbitals. Such ansatz states are expected to work well only in small systems (like the ones tested in this paper) or in extremely inhomogeneus ones where the excitations are effectively localised in a small part of the natural orbital basis. Our functional seems to perform well only in the specific task of computing the electron transfer, as it generally does not recover the correlation energy in the strong interaction limit, as seen in Fig. 10. A possible way of constructing more accurate functionals would be to superpose Slater determinants with higher numbers of excitations, as one expects that the number of excitations grows linearly with the size of the system. Finally, another practical potential use for our results could be a testing criterion for finding accurate density functionals -an accurate approximate density functional should reproduce the repulsive gradient.
Let us remark that one of the main obstructions in finding the exact RDMFT functional is that it requires a detailed knowledge of the geometry of the set of singleparticle reduced density matrices. The description of this set has been long-recognised as one of the most important fundamental problems in quantum chemistry [44]. The solution of this problem in the late 2000s [37] provided new insights to the RDMFT and it may lead to describing the counterparts of the quantum depletion effect for more general systems of fermions or distinguishable particles using the generalised doublyexcited states [45]. Finally, note that throughout the paper we are implicitly assuming that the Levy's RDMFT functionalF is actually differentiable at points from the set Γ BEC (Γ HF ). Earlier works studying the exact RDMFT functional for Hubbard dimers [18,19] suggest that this is a reasonable assumption. However, we are not aware of any mathematically rigorous proof of the fact that the extended functional is differentiable at Γ BEC (Γ HF ) for any pair interactionŴ . Perhaps so-called symplectic slice techniques used in [45] may shed some light on a solution of this problem.

Acknowledgments
I would like to thank Christian Schilling for introducing me to the subject of the Reduced Density Matrix Functional Theory and for the feedback at the early stages of the manuscript. I am also very grateful to Jonathan Robbins for providing careful feedback about the manuscript and for many useful discussions.
Appendix A. Calculation of c 1 for d ≥ 2 Appendix A. 1

. The bosonic density functional
Recall the definition of our ansatz state where ≥ 0, α 2 1 + . . . + α 2 d−1 = 1, σ = ±1 and the normalization constant reads N = 1/ √ 1 + 2 . The basis states are defined in terms of bosonic creation and annihilation operators We Recall also the definitions of generators of the single-particle unitaries.
In particular, we want to extract all terms which yield the constant c 1 from Equation (18), i.e. terms which are proportional to N 2 in Equation (A.3). To this end, we focus on the expression b † l b k ± b † k b l |Ψ ,σ , k ≤ l, as its knowledge determines all the vectorŝ R µ |Ψ ,δ and thus all the expectation values in Equation (A.3).