Functional Theory for Bose-Einstein Condensates

One-particle reduced density matrix functional theory would potentially be the ideal approach for describing Bose-Einstein condensates. It namely replaces the macroscopically complex wavefunction by the simple one-particle reduced density matrix, therefore provides direct access to the degree of condensation and still recovers quantum correlations in an exact manner. We eventually initiate and establish this novel theory by deriving the respective universal functional $\mathcal{F}$ for general homogeneous Bose-Einstein condensates with arbitrary pair interaction. Most importantly, the successful derivation necessitates a particle-number conserving modification of Bogoliubov theory and a solution of the common phase dilemma of functional theories. We then illustrate this novel approach in several bosonic systems such as homogeneous Bose gases and the Bose-Hubbard model. Remarkably, the general form of $\mathcal{F}$ reveals the existence of a universal Bose-Einstein condensation force which provides an alternative and more fundamental explanation for quantum depletion.


I. INTRODUCTION
Bose-Einstein condensation (BEC) is one of the most fascinating quantum phenomena. While its theoretical prediction by Einstein [1] based on Bose's work [2] dates back to 1925, the realization of BEC for ultracold atoms in 1995 [3][4][5] has led to a renewed interest. The development of the respective field of ultracold gases has opened new research avenues and revealed new phenomena such as the crossover from BEC-superfluidity to BCSsuperconductivity [6][7][8][9]. The general need to describe bosonic quantum systems within and also beyond the ordinary BEC regime has urged us very recently to put forward a novel physical theory for describing interacting bosonic quantum systems [10]. This bosonic reduced density matrix functional theory (RDMFT) is based on a generalization of the Hohenberg-Kohn theorem, abandons the complex N -boson wavefunction but still recovers quantum correlations in an exact way. Since it involves the one-particle reduced density matrix (1RDM) γ as the natural variable it would be particularly wellsuited for the accurate description of Bose-Einstein condensates. Indeed, according to the Penrose-Onsager criterion [11], BEC is present whenever the largest eigenvalue of the 1RDM, n max = max ϕ ϕ|γ|ϕ , is proportional to the total particle number N . As a matter of fact, n max quantifies the number of condensed bosons, without requiring any preceding information about the form of the maximally populated one-particle state |ϕ max .
While bosonic RDMFT would potentially be the ideal theory for describing BECs (including the regime of fractional BEC as well as quasicondensation [12]), RDMFT of course does not trivialize the ground state problem. Instead, it is the fundamental challenge in RDMFT to * c.schilling@physik.uni-muenchen.de construct reliable approximations of the universal interaction functional F(γ), determine its leading order behavior in certain physically regimes or its exact form for simplified model systems. Results along any of those lines are typically quite rare, however, and their significance for the general development of RDMFT could hardly be overestimated. The latter is due to the fact that improved functional approximations do often build upon previous ones (see, e.g., [13][14][15] and references therein). In fermionic RDMFT, the elementary Hartree-Fock functional [16] can be seen as the first level of the hierarchy of functional approximations. It has directly led to the celebrated Müller functional [17,18] which in turn inspired more elaborated functional approximations [13,15]. In bosonic RDMFT even the analog of the Hartree-Fock functional has not been established yet. It is therefore the main goal of our work to initiate and establish this novel bosonic RDMFT by deriving such a first-level functional in a comprehensive way. Due to the significance of BEC, we identify systems of interacting bosons in the BEC regime as the starting point for the hierarchy of functional approximations. It is worth noticing that this regime as described by the Bogoliubov theory covers a large range of systems, including in particular the experimentally realized dilute ultracold Bose gases as well as charged bosons in the high density regime. The respective first-level functional may then not only serve as a starting point for the development of further functional approximations but its concrete form will also reveal a remarkable new physical concept. The gradient will namely be found to diverge repulsively in the regime of almost complete BEC, preventing quantum systems of interacting bosons from ever reaching complete condensation. This BEC force will thus provide an alternative explanation for quantum depletion which is most fundamental : It emerges from the geometry of density matrices and the properties of the partial trace, independently of the pairinteraction between the bosons and other system-specific features.
The paper is structured as follows. In Sec. II we discuss the foundation of bosonic RDMFT, recall conventional Bogoliubov theory and explain why the latter is incompatible with RDMFT from a conceptual point of view. Then, in Sec. III we present a particle-number conserving modification of Bogoliubov's theory which eventually allows us to derive the universal functional within the BEC regime. We then illustrate in Sec. IV how bosonic RDMFT is applied and present functionals for a number of different systems. In Sec. V we establish and illustrate the novel concept of a BEC force. In the Summary and Conclusion, we provide also a general idea for constructing higher order functional approximations based on a perturbational theoretical generalization of Bogoliubov theory.

II. NOTATION AND CONCEPTS
We outline in this section foundational aspects of RDMFT and its application to homogeneous bosonic systems. We also introduce the most relevant concepts which will be used in subsequent sections and briefly recap conventional Bogoliubov theory. Since our work involves rather different concepts some of which are not broadly known yet this section shall be rather comprehensive to make our paper self-contained.
A. RDMFT The one-particle reduced density matrix (1RDM)γ of an N -fermion/boson quantum stateΓ is obtained by tracing out all except one particlê Equivalently, it can be characterized as the mathematically most primitive object which still determines the expectation values of all one-particle observablesĥ, By exploiting the latter, Gilbert in 1975 [19] has proven for quantum systems of identical fermions/bosons with HamiltoniansĤ the existence of a universal interaction functional FŴ (γ) of the 1RDM: The energy and 1RDM of the ground state ofĤ(ĥ) for any one-particle Hamiltonianĥ can be determined by minimizing the total energy functional The significance of this reduced density matrix functional theory (RDMFT) is based on the fact that the interaction functional FŴ (γ) does not depend on the choice of the one-particle Hamiltonianĥ but only on the in-teractionŴ . Since the latter is typically fixed in each scientific field (we therefore drop the indexŴ in the following), RDMFT is a particularly economic approach for addressing the ground state problem. Indeed, any effort to approximate F(γ) contributes to the solution of the ground state problem ofĤ(ĥ) for allĥ simultaneously. This shall be contrasted with wavefunction-based methods whose application toĤ(ĥ) does in general not provide any simplifying information towards solving other systemsĤ(ĥ ). In Gilbert's RDMFT the domain of the universal functional comprises exactly those 1RDMsγ which follow from ground states, i.e., there exists anĥ withĤ(ĥ) → |Ψ →γ. To circumvent the problem of describing this complicated set the constrained search formalism has been established which is based on the following consideration [20][21][22] .
The variational principle in the first line can refer to either pure or ensemble N -particle quantum statesΓ. Depending on that choice, the constrained search formalism leads to the pure/ensemble 1RDM-functional F with a domain given by all pure/ensemble N -representable 1RDMs. In the context of fermionic quantum systems, the complexity of the pure one-body N -representability conditions (generalized Pauli constraints) will thus hamper the calculation of either the functional's domain or the functional itself [23].

B. Bosonic RDMFT for homogeneous systems
While RDMFT has been developed and applied so far only for fermionic quantum systems, various concepts are applicable to bosonic quantum systems as well, yet with one crucial simplification. Since in case of bosons any 1RDM,γ ≡ α λ α |α α|, is pure N -representable (e.g. to |Φ = 1/ √ N α √ λ α |α, . . . , α ) the respective bosonic RDMFT is not hampered by one-body Nrepresentability constraints. Another crucial reason for us for putting forward bosonic RDMFT [10] was that the 1RDM is the crucial entity for the description and quantification of BEC [11] which is one of the most fascinating quantum phenomena. In that respect, it is worth noticing that the widely used density functional theory [24][25][26] fails to provide a direct description of BEC since its natural variable is the too rudimentary particle density.
Since our work is concerned with homogeneous BECs, we consider only one-particle Hamiltonians which are diagonal in the momentum representation, i.e., there is only a kinetic energy operatorst contributing toĥ but no external potential,ĥ ≡t. Implementing this within the constrained search formalism (5) identifies the momentum occupation numbers n ≡ (n p ) as the natural variables and the pure functional follows as While we are focussing in the following on the pure functional, it is worth recalling that the corresponding ensemble functional would follow as the lower convex envelop of the pure functional F [23]. Also their two domains coincide. To describe , let us first use the normalization constraint to get rid of the entry n 0 = N − p =0 n p . Then the functional's domain follows as In case of finite lattice models there are finitely many momenta p (forming a discrete Brillouin zone), while in case of continuous systems or infinite lattices, n will have infinitely many entries. It will be instructive to also understand the functional's domain from a geometric point of view. Apparently, is a convex set which after all takes the form of a simplex with vertices 0 and v p = N e p , where e p has only one non-vanishing entry 1 at position p. In the following, we are mainly interested in the regime of BEC which is characterized by an occupation number n 0 close to N . This corresponds in the simplex to the neighborhood of the vertex 0, which can equivalently be characterized by the simultaneous saturation of the constraints n p ≥ 0 for all p = 0.
Last but not least, we would like to stress that the parity-symmetry of common physical spaces implies the additional symmetry n p = n −p for all momenta p. This does not really change the geometric form of the functional's domain but just allows us to skip in the definition (7) for every pair (p, −p) of momenta one of the two occupation numbers n ±p . In the context of RDMFT, respecting this common symmetry would mean to restrict the kinetic energy operatorst ≡ p ε pnp to those with ε p = ε −p .

C. Recap of conventional Bogoliubov theory
In this section be recap the most important aspects of Bogoliubov's [27] well-known and experimentally confirmed [28] theory to describe BEC in homogenous bosonic quantum systems and the effect of depletion of the condensate as a result of the interaction.
The Hamiltonian describing a homogenous system of N interacting spinless bosons in first quantization ( ≡ 1) is given bŷ Its second quantized form in momentum representation for particles in a large box of volume V = L 3 and size L with periodic boundary conditions then readŝ where W p is the Fourier transform of W (·). In case of an isotropic pair interactions, W in Eq. (8) would depend only on the modulus of the distance between the particles i and j which in turn would imply W p ≡ W |p| . The most crucial feature of the HamiltonianĤ and the pair interactionŴ is that they are conserving the particle number as well as the total momentum. Assuming a BEC at T = 0, the standard approach to determine the ground state energy (and the low lying excited states) of the Hamiltonian (9) is the Bogoliubov approximation [27]. It is based on the fact that in the regime of BEC the zero-momentum mode is macroscopically occupied and interactions between non-condensed bosons can be neglected due to the conservation of momentum: Since application of a creation/annihilation operator a ( †) 0 to the BEC ground state leads to macroscopically large prefactors of the order √ N , terms in the expansion (9) ofŴ involving less then two 0-indices are dropped. The resulting quartic interaction is further simplified by replacing the condensate operatorŝ a 0 ,â † 0 → √ n 0 ≈ √ N by a c-number. This eventually leads to the quadratic Bogoliubov Hamiltonian which involves (besides the kinetic energyt and some trivial contributions) for each pair (p, −p) an anomalous term of the formâ † pâ † −p +â pâ−p . The Bogoliubov Hamiltonian can then easily be diagonalized by a Bogoliubov The respective ground state follows as |Ψ =Û B |N , where |N ≡ (N !) −1/2 (â † 0 ) N |0 is the ground state of the non-interacting system and |0 the vacuum state. The phases θ p are chosen such that the anomalous terms in the Hamiltonian, containing either two quasiparticle annihilation (b p ≡Û † Bâ pÛB ) or creation operators (b † p ) vanish to eventually obtain a diagonal quadratic form inb p (see also textbook [29] for more details). Bogoliubov's approach can also be interpreted as the variational minimization of the Bogoliubov Hamiltonian over all trial states of the formÛ B |N .
D. Incompatibility of RDMFT and conventional Bogoliubov theory As explained in the previous section, Bogoliubov's approximation results in a Hamiltonian which is not particle-number conserving anymore. At the same time, RDMFT defines a universal functional F(n) (or more generally F(γ)) by minimizing the interaction Hamiltonian according to (6) with respect to quantum states with a fixed total particle number N and fixed momentum occupation numbers n. Replacing in Eq. (6)Ŵ by Bogoliubov's approximated Hamiltonian would therefore erroneously ignore the important anomalous termŝ a † pâ † −p +â pâ−p . At first sight, this incompatibility of Bogoliubov's conventional approximation and RDMFT seems to be paradoxical. Yet, it is worth recalling that the merits of the unitary Bogoliubov transformation lie in the simple calculation of the (low-lying) energy spectrum while its violation of particle-number conservation can lead to conceptual difficulties beyond RDMFT as well. At the same time, since RDMFT has the distinctive goal to (partly) solve the ground state problem forĤ(ĥ) for allĥ simultaneously, it requires apparently a mathematically more rigid and well-defined framework than the one provided by conventional Bogoliubov theory.
Before we discuss in the following section such a welldefined mathematical framework for realizing Bogoliubov's ideas within RDMFT, we briefly comment on an alternative natural idea for circumventing the outlined difficulties. Instead of applying the constrained search formalism to a fixed particle number sector, one could also extend (6) to the entire Fock space. This would result in a Fock space RDMFT and the anomalous terms would contribute to the functional. Yet, there would be a crucial drawback. The respective functional would namely allow one for any Hamiltonian (3) to only cal-culate the overall ground state on the Fock space. For instance, for specific kinetic energy operators or pair interactions, this overall minimum may lie in the sector of zero or infinitely many bosons. Also adding a chemical potential term µN for steering the particle number to a preferred one would only work in case the Fock space functional was convex in the total particle number.

III. DERIVATION OF THE UNIVERSAL FUNCTIONAL
A. Particle-number conserving Bogoliubov theory As discussed in the context of Eq. (5) and motivated in Sec. II D, the derivation of the universal functional for BECs requires a particle-number conserving variant of conventional Bogoliubov theory. Exactly such a modification has been provided by Girardeau [30] (see also [31][32][33][34]) in the context of pair theory. In the following, we will outline and then apply this theory which in particular improves upon Bogoliubov theory by including more terms of the Hamiltonian. The idea behind pair theory is that in the regime of BEC, excitations of pairs (p, −p) from the condensate dominate and thus the interacting ground state is well approximated by a state with a corresponding pairing structure [30,32]. Restricting the original HamiltonianĤ to the space of such pair excitation states and assuming that the zero-momentum state is macroscopically occupied means to effectively deal with a modified interactionŴ P of pair excitation type [30], The terms in the first line of Eq. (10) give rise to the Bogoliubov Hamiltonian (after the replacementâ 0 ,â † 0 → √ N ) while those in the second line improve upon Bogoliubov theory. For the sake of simplicity, we omit in the following derivations the constant term N (N −1)W0
To determine a variational ground state energy ofĤ(ĥ), Girardeau's idea was then to employ a particle-number conserving analogue of Bogoliubov trial stateÛ B |N . For this, one first introduces the operatorŝ which annihilate/create a boson in the condensate, yet without changing the normalization of the respective quantum state. Girardeau's N -boson trial states of pair excitation form are defined by the following operatorŝ with θ p ∈ R and θ p = θ −p . The operatorsÛ G are particle-number conserving as desired, [Û G ,N ] = 0, which is due to the additional operatorsβ 0 andβ † 0 . Since its exponent is antihermitian,Û G is still unitary (asÛ B ). The price one has to pay for the more complicated exponent, however, is that no compact exact expression can be found for the quasiparticle operatorsÛ † Gâ pÛG anymore. Instead the result known from Bogoliubov theory holds only approximately, where φ p ≡ tanh(θ p ). A careful mathematical estimate of the difference between left and right side of (14) has been provided in [33] (yet involving a slightly different but conceptually similar definition ofÛ G ). It effectively allows us to treat (14) and the implied Eq. (15) as exact relations for our further derivation. The particle number expectation value of the momentum mode p = 0 in the interacting ground state |Ψ (12) then follows as

B. Calculation of the Functional
Relation (15) is the crucial ingredient for our derivation of the universal functional F in the regime of BEC. This connection between the family of variational trial states of fixed particle number and the momentum occupation numbers n will drastically simplify the constrained search (6) and will allow us eventually to determine the explicit form of F. For this, we observe that relation (15) can be inverted up to binary degrees of free- This sign ambiguity is conceptually very similar to the so-called phase dilemma in fermionic RDMFT [35]. The latter resembles the fact that general phase changes of the natural orbitals (eigenstates of the 1RDM) affect Ψ|Ŵ |Ψ in (5) via the N -fermion wavefunction |Ψ while keeping the 1RDM invariant. In contrast to fermionic RDMFT, however, the minimizing signs {σ p } can be found in our case of bosons in the BEC regime.
We combine now various concepts and ideas to determine the universal functional F(n) for BECs. According to the constrained search formalism (6) we need to minimize for any vector n the expectation value of the interac-tionŴ over all N -boson quantum states with momentum occupation numbers n. Our focus on the regime of BECs then allows us to restrict this to Girardeau's N -boson trial states (12) with the additional effect thatŴ simplifies toŴ P in Eq. (10), i.e., Ψ|Ŵ |Ψ = Ψ|Ŵ P |Ψ = N |Û † G W PÛG |N . The operatorÛ † G W PÛG should then be expressed in terms of the quasiparticle operatorsξ p given by Eq. (14), allowing us to eventually calculate its action on the state |N . Since the trial states |Ψ are almost uniquely determined by n according to (15) we are only left with a minimization over all possible combinations of signs σ p . Keeping only terms which do not vanish in the thermodynamic limit N → ∞, V → ∞ and n = N/V = cst. yields then the final result for the Girardeau approximated functional where For general n ∈ one cannot overcome the common phase dilemma and in particular the minimizing sign factors σ p = ±1 in (17) depend on n. This in turn leads to a partitioning of the functional's domain into cells characterized by different signs {σ p }, similarly to the Ising cells corresponding to different spin configurations (see also Fig. 2 for an illustration). As already explained in Sec. III A, Girardeau's approach based on pair theory goes beyond Bogoliubov theory by including additional terms of the original Hamiltonian (see also second line of (10)). Yet, since those involve fewer creation/annihilation operators a ( †) 0 and since the Girardeau approach uses at the end (almost) the same trial states (12) as Bogoliubov, we expect that the additional terms I 1 , I 2 in (17) have only a minor quantitative rather than a significant qualitative effect on the description of BECs. Whether this changes beyond the regime of BEC is not clear since one still restricts to the common BEC trial states (12).
In the regime of BEC the two terms in Eq. (17) involving I 1 and I 2 , respectively, are significantly smaller than the term proportional to n 0 . Accordingly, in the regime of interest the minimization of various σ p can be executed analytically, leading to Also, the possible approximation n 0 ≈ N would be of the same order as neglecting the less significant Girardeau terms I 1 and I 2 . Eventually, implementing those two last approximations leads to one of our key results, the Bogoliubov approximated functional (n ≡ N/V ) The distinctive form of the Bogoliubov functional F B resembles clearly the decoupling of various momentum pairs (p, −p) from each other within Bogoliubov theory.
Remarkably, the Bogoliubov approximated functional F B is convex, in contrast to common pure functionals in fermionic RDMFT. The pure functional F B therefore coincides with the corresponding ensemble functional since the latter is given by the lower convex envelop of the former [23]. We also would like to reiterate that due to the general significance of BECs, the functional (20) can be seen as the first-level approximation of the universal functional in bosonic RDMFT. In analogy to the Hartree-Fock [16] and the Müller functional [17,18] in fermionic RDMFT and the local density approximation in density functional theory [36], F B and F G will represent a promising starting point for the construction of more elaborated functional approximations. In that sense, we expect that our key results (17) and (20) will initiate and establish eventually bosonic RDMFT. In the following, we simplify our notation by skipping the index B, G of the functional, also since both functionals (almost) coincide in the relevant regime of BEC.

A. Dilute Bose gas in 3D
We now apply the concepts of RDMFT to the homogenous dilute Bose gas, the system for which Bogoliubov's theory [27] was originally developed. This will also allow us to demonstrate how the well-known expression for the ground state energy of a dilute Bose gas [37] can be obtained using RDMFT.
Let us introduce for the following considerations the degree D of quantum depletion (N BEC ≡ n 0 ) From a geometric point of view, D is nothing else than the l 1 -distance of n in the simplex (7) to the vertex 0 corresponding to complete BEC. We also recall that the ground state energy ofĤ(t) =t +Ŵ follows in RDMFT by minimizing the respective energy functional over the space of occupation number vectors n, wheret ≡ p ε pnp , assuming w.l.o.g. ε 0 = 0, and ε·n ≡ p =0 ε p n p . To calculate for the realistic dilute Bose gas the ground state energy and the degree of condensation, we would need to plug in for the kinetic energy in (22) the specific dispersion relation of free particles, i.e., ε p = p 2 /2m (ignoring boundary effects). It is worth reiterating that in principle systems with any kinetic energyt could be considered in RDMFT. From an experimental point of view, one could indeed imagine a modified dispersion relation due to a specific background medium and in case of lattice models both the rate and range of the hopping can be actually varied (see, e.g., Refs. [38,39]). Because of this, we are for the moment still considering a general t and ε, respectively.
Finding the minimizer n of the energy functional then means to solve Using the explicit form of the Bogoliubov functional (20) then leads to This is nothing else than the well-known result for the momentum occupation numbers [29].
Considering now the specific case of a realistic dilute Bose gas then allow us to determine the ground state energy explicitly. For this, we first evaluate the universal functional at the minimum n, leading to (see Appendix A) It depends only on the first two terms of the Born series of the s-wave scattering length a [37], As it is shown in Appendix A, adding the kinetic energy and reintroducing the omitted constant term W 0 N (N − 1)/2V ≈ W 0 N n/2 leads to the well-known ground state energy [37]: is infinite dimensional in case of a continuous system in a box it is difficult to graphically illustrate the Bogoliubov functional. Yet, to visualize at least some of its most crucial features we define two paths within , both starting at the physical point n (corresponding to ε p = −p 2 /2m) and terminating at the vertex 0 describing complete BEC. The first one is just the straight path s between those two points, parameterized by t ∈ [0, 1], The l 1 -distance D(t) of n(t) to 0 follows directly as and the functional's concrete values F[n(t)] along that path can easily be calculated by exact numerical means.
As second path κ, we consider the one experimentally realized in Ref. [28] by continuously reducing the coupling strength κ of the pair interaction from one to zero. Since the interaction HamiltonianŴ in RDMFT is fixed, this path has to be realized equivalently by increasing the strength of the kinetic energy according to p 2 /2mκ. The respective distance D of n(κ) to 0 follows as and the functional F[n(κ)] along that path is given by In Eq. (31) one may replace W 0 by a 0 according to (26). The result for F as a function of the fraction D of noncondensed bosons along the two paths s and κ is shown in Fig. 1 for the parameters n = 10 −3 , W 0 = m = 1 and a 1 = −0.01. This choice of parameters (recall that we set several physical constants to one) corresponds to realistic dilute Bose gases as our following results of the small degrees of depletion will confirm. We observe that the Bogoliubov functional F goes to zero for D = 0 which corresponds to complete BEC. Also, it can be seen that the gradient of F increases for smaller distances D. In Sec. V we will show that the gradient of F actually diverges in the limit D → 0 and provide a detailed discussion of this remarkable and far-reaching observation.
A generalization of F given by Eq. (20) to dimensions d = 3 within the s-wave scattering approximation is possible if the Bogoliubov approximation for the given set of parameters, i.e. the interaction strength and the density, is valid. Two-dimensional dilute systems are weakly interacting if the condition n|a| 2 2D 1 [40,41] is satisfied where a 2D is now the respective two-dimensional s-wave scattering length. In contrast to higher dimensional systems, a one-dimensional Bose gas is weakly interacting in the limit of high densities and the validity of the Bogoliubov approximation in that limit was shown in Ref. [42]. Due to their distinctive role, we will study a one-dimensional model in Sec. IV C.

B. Charged Bose gas in 3D
In contrast to the dilute Bose gas discussed in the previous section, the scattering of charged bosons cannot be described within the s-wave scattering approximation anymore. This is due to the infinite range of the Coulomb interaction W (r) ∝ 1/r. The respective Fourier coefficients W p can still be determined analytically though. In case of an additional uniform background they follow as For charged bosons the weak interaction regime corresponds to the high density limit [43][44][45]. This regime to which Bogoliubov's approximation refers to is characterized by a small "gas parameter", r s ≡ (3/4π) 1/3 me 2 n −1/3 1.
To illustrate again how RDMFT works, we calculate the energy and momentum occupation numbers n p of the ground state for the most realistic case of a kinetic energy given byt = p p 2 2mn p . For this, we add the exact kinetic energy functional Tr 1 [t(·)] to the universal interaction functional F with Fourier coefficients W p given by Eq. (32). Then, we minimize the total energy functional with respect to all n ∈ , leading to the minimizer n which is given by Eq. (24). Evaluating then the functional at n is straightforward (in contrast to the dilute neutral Bose gas) and one finds (recall n ≡ N/V ) Eq. (34) verifies that the depletion of the condensate decreases with increasing density n. Adding the kinetic energy to Eq. (33) leads to the known result for the ground state energy ( = 4π 0 = 1) [46]: As a consistency check, this confirms the correctness of Eq. (33). Next, in analogy to Sec. IV A we consider again the straight path s and the curved path κ. The latter is again defined as the curve n(κ) obtained by reducing by factor κ ∈ [0, 1] the coupling strength of the Hamiltonian above (which led to the results Eqs. (33), (34) and (35)). Evaluating the distance D along the path κ yields D(κ) = Dκ 3/4 (36) and the functional F takes the values For the path s we have D(t) = (1 − t)D and the concrete values of the functional F along that path can be evaluated by exact numerical means. The right panel of Fig. 1 shows F along the two paths s and κ. The curves have qualitatively similar shapes to those for the dilute neutral Bose gas shown on the left of Fig. 1. This is not surprising because both setups correspond to a weakly interacting system in which the Bogoliubov approximated functional Eq. (20) is valid. However, we will see in Sec. V that the momentum dependence of W p can alter the behaviour of the gradient of F.

C. Bose-Hubbard model for five lattice sites
As a third example, we discuss in this section the onedimensional Bose-Hubbard model. For illustrative purposes, we consider the specific case of just L = 5 lattice sites and N = 100 bosons since this allows us to visualize the functional and its gradient on the entire domain . Indeed, for L = 5 there are only two independent momentum occupation numbers due to the general parity symmetry n p = n −p and normalization n 0 = N − p =0 n p .
We start by discussing a few conceptual aspects which are valid for any number L of sites (assuming for simplicity L odd). The one-dimensional Brillouin zone comprises momenta p = 2πν/L where ν takes integer values in the interval described by |ν| ≤ (L − 1)/2. The bosons interact via Bose-Hubbard on-site interaction as described by the operator U 2 L j=1n j (n j − 1). It is worth recalling that the universal functional depends on the entire interaction Hamiltonian, i.e., it includes a priori the coupling constant U as well. Yet, due to the linear structure of the constrained search (5), (6) any non-negative prefactor could be separated from the interaction Hamil-tonianŴ and added instead in front of the respective universal functional, It is crucial to observe that the same does not apply to possible sign factors since otherwise this would mean to change the minimization in (5), (6) to a maximization. Because of this, we consider in the following the interaction HamiltonianŴ ≡ sgn(U ) 2 L j=1n j (n j − 1) and add eventually the coupling constant |U | in front of the respective universal functional FŴ . To proceed, it is then an elementary exercise to determine the corresponding Fourier coefficients W p = sgn(U ) which are in particular independent of the (one-dimensional) momentum p. The universal functional in the BEC regime is obtained by plugging the concrete result for the Fourier coefficients W p into the general formula for the Bogoliubov functional (20) or its extension (17) based on Girardeau's approach. Just to reiterate, the respective functionals are valid in the regime of BEC, i.e., for weak interactions. In contrast to their higher dimensional counterparts, onedimensional systems require high densities n = N/L 1 to be weakly interacting [42].
From a general point of view, the context of lattice models emphasizes very well the conceptual advantages of RDMFT relative to wavefunction based methods. After having determined the universal interaction functional FŴ (or decent approximations thereof) the ground state energy of every HamiltonianĤ(t) =t + |U |Ŵ can be calculated with relatively little computational effort by minimizing the total energy functional with respect to all n ∈ . In that sense, RDMFT represents a highly economic approach for solving simultaneously the ground state problem for the entire class {Ĥ(t)} of Hamiltonians. For continuous systems the benefits of this are less obvious since there is essentially one particularly relevant choice for the kinetic energy operatort. This is quite different for lattice models since both the rate and the range of the hopping can be varied in experiments (see, e.g., Refs. [38,39]). Nonetheless, we focus in the following on hopping just between neighbouring sites at a rate t ≥ 0, i.e., we chooset = −2t p cos(p) − 1 n p and w.l.o.g. fix |U | ≡ 1. In analogy to the most realistic dilute and charged Bose gas as discussed in Sec. IV A and Sec. IV B, respectively, we pick t = U = 1 as a reference point for further investigations and illustrations .
To illustrate and compare the Bogoliubov-(20) and the Girardeau-approximated functionals (17) for the specific case of L = 5 sites we first need to execute the minimization of the sign factors in (17). As it is shown in Appendix C, this can be done analytically due to the specific Fourier coefficients and leads to a partitioning of the functional's domain into three regions. Just for illustrative purposes, we present in Fig. 2 the entire domain of the functional (17) (recall its validity refers to the regime of BEC only) and the three cells which are characterized by different minimizing sign configurations (σ p1 , σ p2 ) in (17). As the two independent occupation numbers we choose here the momenta p 1 = 2π/5 and p 2 = 4π/5 which can take values n pj ∈ [0, 50]. The vector (0, 0) corresponds to complete BEC, i.e., N BEC ≡ n 0 = N = 100 and its vicinity represents the BEC-regime to which our functionals  In Fig. 3 we present now the Bogoliubov functional (20) and its extension (17) based on Girardeau's approach in the form of a contour plot in the BEC-regime of not too large quantum depletion. The results for the two functionals are in quite good agreement for small degrees D of depletion. The occupation number vector n = (0.91, 0.44) obtained from minimizing the total energy functional for the reference point (t, U ) = (1, 1) is shown in Fig. 3 as well. The corresponding degree of depletion, D = 2.7%, justifies in retrospective the treatment of the interactionŴ within the Bogoliubov theory and the usage of the functionals (20) and (17), respectively. For stronger quantum depletion the two functionals in Fig. 3 begin to differ also qualitatively. Their (small) deviation already in the regime of BEC with a degree of depletion around 2% emphasizes the quantitative significance of the additional terms I 1 and I 2 and the usage of the exact value n 0 rather than its approximation to N in Eq. (17).
For the discussion in Sec. V of the new concept of a BEC force, we define in Fig. 3 five qualitatively different paths towards the polytope boundary ∂ , all starting from the point n. The path denoted by s corresponds to a straight path towards complete BEC and κ denotes the path where the interaction strength of the model is reduced by increasing the kinetic energy by a factor 1/κ with κ ∈ [0, 1]. The path denoted by a runs perpendicular towards the hyperplane defined by n p1 + n p2 = 0. Consequently, it corresponds to the path with the fastest increase of the condensate fraction (yet it will not reach complete BEC). In the cases b and c one occupation number is fixed while the other one is continuously decreased to zero. In Fig. 4 we present the functional F as a function of D = 1 − N BEC /N along those five paths. The black dots emphasize that the value of F at the boundary ∂ remains finite (quite in contrast to its derivative as shown and discussed in the subsequent section). The  Fig. 3 convexity of the curves in Fig. 4 corresponding to the four straight paths a, b, c, s just reflects the local convexity of the exact universal functional in the regime of not too large quantum depletion. In this context, we would like to reiterate that this convex behavior and the repulsive character of the functional's gradient close to the boundary emerges from the minimization of the sign factors in (17) leading to (19).

V. REPULSIVE BOSE-EINSTEIN CONDENSATION FORCE
In this section we explore in more detail the behavior of the functional and its gradient close to the boundary of their domain in the regime of BEC. This will eventually allow us to reveal and establish the novel concept of a BEC force. Due to its potentially far-reaching consequences for our understanding of bosonic quantum systems, we will calculate and illustrate the BEC force in Sec. V B for the three different systems studied in Sec. IV.

A. General results
The functional (20) which is based on the Bogoliubov approximation is convex on its entire domain (7). Since this approximate functional is exact in leading order in the regime of BEC with not too large quantum depletion all conclusions drawn from it are valid for the exact universal functional of (8) as well. The illustrations in the previous section for dilute and charged Bose gases in 3D and the Bose-Hubbard model have also confirmed the distinctive convex behavior of the universal functional in the BEC-regime. While the functional itself remains finite even at the point 0 ∈ of complete condensation, the same will not be true anymore for the functional's gradient. This can easily be deduced from the form of the Bogoliubov functional (20). To be more specific, approaching the vertex 0 of the simplex means to simultaneously send all momentum occupation numbers n p with p = 0 to zero. Taking then the derivative of (20) (or of its extension (17)) with respect to n p sufficiently close to 0 yields in leading order It is worth noticing that the divergence of this derivative for n p → 0 is always repulsive for any interactionŴ . This remarkable feature follows directly from the minimization of the sign factors in (17), leading to (19). The repulsive nature of the diverging gradient also proves universally that occupation numbers in interacting bosonic quantum systems can never attain the exact mathematical value 0. Although our work refers to homogeneous systems in their BEC regime only, we have little doubt that this conclusion is also valid for any generic nonhomogeneous interacting bosonic quantum system, also beyond the BEC regime. The general result (39) implies that the point 0 of complete condensation can never be reached, independent of the path towards 0 that is envisaged. Since the functional F is finite this seems to be paradoxical as far as the energy is concerned. Yet, the reader shall note that it is the kinetic energy which will need to diverges according (23) to enforce such a path towards 0.
We also would like to emphasize that the divergence of the gradient of F along a straight path is always proportional to 1/ √ D and its prefactor depends on the direction of the path, i.e., the angle at which 0 is approached. To confirm this in a quantitative way, let us consider a general straight path from a starting point n in the regime of BEC to 0, linearly parameterized by t ∈ [0, 1], The degree D (21) of quantum depletion along that path reduces according to The gradient of F projected onto that path is then nothing else than the weighted sum of individual contribu-tions (39) from every p, This second key result of our work establishes the new concept of a BEC force which prevents interacting bosonic quantum systems from ever exhibiting complete BEC. This novel concept is conceptually very similar to the fermionic exchange force that we have recently revealed and established in fermionic lattice models [47].

B. Examples
In this section we illustrate the novel concept of a BEC force (42) for various systems introduced in Sec. IV.

Bose gases in 3D
We revisit the 3D Bose gas for neutral atoms in the low density and for charged atoms in the high density regime. The aim is to calculate for those concrete systems the explicit values of the BEC force (42). For both systems, the derivative of F with respect to the degree D of quantum depletion along the path s defined by (40) is given by Eq. (42). The summation over p = 0 can be converted into an integral in the thermodynamic limit where N → ∞, V → ∞ and n = N/V = cst.. This eventually allows us (see Appendix B) to obtain a compact analytic expression for the BEC force, where η(a 0 , n, m) is a positive constant and a 0 and a 1 are the first two terms in the Born series for the scattering length a. It is worth reiterating that according to the general result (42) the BEC force is always repulsive. Since only the second term in Eq. (43) is negative (recall a 1 < 0) this imposes in turn a bound on the maximal valid distance D of the starting point n ∈ to the regime of complete BEC.
For the charged Bose gas the last expression in the second line in Eq. (42) can only be calculated by exact numerical means. Nonetheless, this also allows us to confirm the square root dependence of the divergence. In general, the functional's gradient diverges as 1/ dist(n, ∂ ) along straight paths reaching any arbitrary point on the boundary ∂ in the regime of BEC.
Moreover, we determine for both systems the BEC force along the curved path n(κ) which is defined by reducing an additional coupling constant κ in front of W from one to zero. Since exactly this path has been implemented in a very recent experiment [28] this may suggest a first experimental setup for realizing and visualizing our novel concept of a BEC force. The explicit calculation of the BEC force along the κ-path follows directly from differentiation of the expressions in (31) and (37), respectively, leading to and dF charged dD

Bose-Hubbard model
We illustrate the BEC force and the diverging behaviour of the functional's gradient close to the boundary ∂ of its domain in general for the Bose-Hubbard model. For this we consider again as in Sec. IV C the case of N = 100 bosons on L = 5 sites. We then determine the directional derivative of the functional along the five paths which were defined in Fig. 3. Since for all five paths the distance D of the occupation number vector n to 0 is monotonously decreasing we can parameterize the functional's derivative along each path by D. The respective results are depicted in Fig. 6. There, the vertical solid lines correspond to the values of D at which the respective paths reach the boundary of (see also Figs. 3,4). We first observe that for all five paths −∂F/∂D is diverging at the end point of each path on the boundary ∂ . As a rather elementary analysis reveals (not shown here) this divergence is always proportional to 1/ dist(n, ∂ ). As far as the four straight paths a, b, c, s are concerned, this was expected given the general results of Sec. V A. In contrast to the two continuous Bose gases, however, the same applies in the Bose-Hubbard model also for the curved κ-path which is obtained by just reducing the coupling strength.

VI. SUMMARY AND CONCLUSION
By using a particle-number conserving modification of the conventional Bogoliubov theory we succeeded in deriving the respective universal interaction functional F(n) for homogeneous Bose-Einstein condensates (BECs). The crucial ingredient of our derivation has been the strong connection (15) between the vectors n of momentum occupation numbers and the Bogoliubov trial states. This has drastically simplified the Levy-Lieb constrained search and the minimization (17) with respect to the remaining independent degrees of freedom (sign factors) could be executed as well, eventually leading to our key result (20). The universal functional F(n) comprises separate contributions of various momentum pairs (p, −p) (recall n p = n −p ). This of course resembles the fact that Bogoliubov theory describes independent pair excitations of condensed bosons, (0, 0) → (p, −p).
While typical experimental realizations of BEC involve ultracold gases in their dilute regime, it is worth recalling that the scope of Bogoliubov theory is wider. As a matter of fact, it applies to any homogeneous system in arbitrary dimensionality as long as it exhibits BEC with not too large depletion. Thus, to identify the range of applicability of our derived first-level bosonic RDMFT functional one just needs to understand the conditions under which different physical systems exhibit BEC. For instance, for neutral bosons in dimensions d ≥ 2 this would be the dilute regime (which implies effectively a weak coupling) while for d = 1 and also for charged bosons in d = 3 it would be the high density regime. It will be one of the future challenges to determine also the functional (or at least approximations thereof) in other physical regimes beyond BEC. Similarly to the Hartree-Fock functional in fermionic RDMFT [16] and the local density approximation in density functional theory [26], the Bogoliubov functional could then serve as a starting point for such approximations. That's particularly promising, since in contrast to the Hartree-Fock functional our functional already involves some quantum correlations while the former always leads to occupation numbers identical to just one and zero [16].
The second key result of our work is the emergence of a universal BEC force which prevents general quantum systems of interacting bosons from ever reaching complete BEC. To be more specific, the gradient of the universal functional has been found to repulsively diverge as 1/ 1 − N BEC /N in the regime of almost complete BEC. This BEC force can be seen as a collective force since it comprises individual diverging terms −∂F/∂n p ∼ n|W p |/2 √ n p from each momentum p. It is worth emphasizing that it has been exactly the minimization (17) of the phase factors within the constrained search formalism which eventually made this BEC force and more generally ∇ n F close to the boundary of its domain repulsive regardless of the signs of various Fourier coefficients W p . The BEC force provides an alternative explanation for quantum depletion which is most fundamental in the sense that it is merely based on the geometry of density matrices and the properties of the partial trace. Indeed, the type of interaction between the bosons only affects the prefactor of the BEC force, while its diverging behaviour proportional to 1/ 1 − N BEC /N is universal. Due to this universal behaviour one may expect that the respective prefactor could provide important insights into system-specific properties, in some analogy to the so-called Tan's contact [48][49][50]. Moreover, the BEC force is nothing else than the bosonic analogue of the recently discovered fermionic exchange force [47], showing exactly the same diverging behaviour close to the boundary of the fermionic functional's domain. Last but not least, we comment on important followup challenges besides the possible experimental realization of the BEC force in ultracold gases and the construction of functionals with a larger range of validity. One such direction would be to quantify in a mathematically rigorous way the deviation of the Bogoliubov functional from the exact functional FŴ of the full interaction HamiltonianŴ in (8). In particular, we are wondering whether the mathematical techniques used in Ref. [33] to estimate the accuracy of Bogoliubov's theory rigorously can be adapted to the context of RDMFT. In that respect, it is worth noticing that our result (17) with phases σ p = sgn(W p ) (or any other choice of phases) represents a universal upper bound to the functional of the full pair interactionŴ in the entire domain. This is due to the following two reasons. First, (17) is found through a variational ansatz and second this ansatz of paring states has only an overlap with the pairing interaction Hamiltonian W P (recall Eq. (10)). Proving in a mathematically rigorous way the BEC force might be particularly challeng-ing. Even if a lower and upper bound can be proven for the difference between the Bogoliubov and the full functional FŴ , the gradient cannot easily be controlled since FŴ may (at least in principle) strongly oscillated between both bounds. A more promising approach would be to exploit the one-to-one correspondence established through the constrained search formalism between n and the Bogoliubov trial state, n ↔Û n |N . The idea for obtaining higher order terms of (20) would be to conjugate in a first stepŴ ,Û nŴÛ As a consistency test we start now from Eq. (A2) and add the kinetic energy. The second term in Eq. (A2) can be split into two parts such that it cancels the divergence in the integral for the kinetic energy as follows: Inserting a 0 and a 1 leads to the ground state energy which is in agreement with Ref. [37].
Appendix B: BEC force for the dilute Bose gas We calculate the derivative of F with respect to the distance along a straight path denoted by s towards complete BEC starting at the occupation number vector n. Then, n p (t) = n p (1 − t) and for t ≈ 1 or equivalently D(t) 1 we can approximate dF(n) dD s = 1 D(t) p =0 nW p n p (t) 1 − 2n p (t) + 1 2 n p (t)(n p (t) + 1) The summation in Eq. (B1) can be replaced by an integral ( p → V such that the integral over p is converging after replacing W p by the constant value W 0 . The first two terms in the Born series for the scattering length a for identical particles are given by Eq. (26). Thus, the summation in the second line of Eq. (B2) can be identified with a 1 and the result of the integration in the first line will depend on W 0 which can be replaced by a 0 through Eq. (26).