$\hat{t}\equiv\sum_{\boldsymbol{p}}\varepsilon_{\boldsymbol{p}}\hat{n}_{\boldsymbol{p}}$,
assuming w.l.o.g. $\varepsilon_{\mathbf{0}}=0$, and
$\boldsymbol{\varepsilon}\cdot\boldsymbol{n}\equiv\sum_{\boldsymbol{p}\neq\mathbf{0}}\varepsilon_{\boldsymbol{p}}n_{\boldsymbol{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
(3.34) the specific dispersion relation of free particles, i.e.,
$\varepsilon_{\boldsymbol{p}}=p^{2}/2m$ (ignoring boundary effects). It is
worth reiterating that in principle systems with any kinetic energy $\hat{t}$
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. [77, 78]). Because of this, we are
for the moment still considering a general $\hat{t}$ and
$\boldsymbol{\varepsilon}$, respectively.

Finding the minimizer $\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}$ of the
energy functional then means to solve

$\varepsilon_{\boldsymbol{p}}=-\frac{\partial\mathcal{F}}{\partial
n_{\boldsymbol{p}}}(\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}})\,,\quad\forall\boldsymbol{p}\,.$
(3.35)

Using the explicit form of the Bogoliubov functional (3.32) then leads to

$\accentset{\rule{3.5pt}{0.5pt}}{n}_{\boldsymbol{p}}=\frac{1}{2}\left(\frac{\varepsilon_{\boldsymbol{p}}+nW_{\boldsymbol{p}}}{\sqrt{\varepsilon_{\boldsymbol{p}}(\varepsilon_{\boldsymbol{p}}+2nW_{\boldsymbol{p}}})}-1\right)\,.$
(3.36)

This is nothing else than the well-known result for the momentum occupation
numbers [63].

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
$\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}$, leading to (see Appendix B)

$\mathcal{F}(\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}})=\frac{128\sqrt{\pi}}{3m}Na_{0}^{5/2}n^{3/2}+\frac{4\pi
Na_{1}n}{m}\,.$ (3.37)

It depends only on the first two terms of the Born series of the s-wave
scattering length $a$ [76],

$a_{0}=\frac{mW_{\mathbf{0}}}{4\pi}\,,\quad a_{1}=-\frac{1}{4\pi
V}\sum_{\boldsymbol{p}\neq\boldsymbol{0}}\frac{W_{\boldsymbol{p}}^{2}m^{2}}{p^{2}}\,.$
(3.38)

As it is shown in Appendix B, adding the kinetic energy and reintroducing the
omitted constant term $W_{\mathbf{0}}N(N-1)/2V\approx W_{\mathbf{0}}Nn/2$
leads to the well-known ground state energy [76]:

$\displaystyle E_{0}$ $\displaystyle=$ $\displaystyle\frac{2\pi
Nn}{m}\left(a_{0}+a_{1}+\frac{128}{15\sqrt{\pi}}a_{0}(na_{0}^{3})^{1/2}\right)\,.$
(3.39)

This can be recast by using the scattering length $a$ which eventually leads
(up to higher order terms) to the compact expression [63] $E_{0}=\frac{2\pi
Nna}{m}\left(1+\frac{128}{15\sqrt{\pi}}(na^{3})^{1/2}\right)$.

Since the underlying domain $\triangle$ 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 $\triangle$, both starting at the physical
point $\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}$ (corresponding to
$\varepsilon_{\boldsymbol{p}}=\boldsymbol{p}^{2}/2m$) and terminating at the
vertex $\mathbf{0}$ describing complete BEC. The first one is just the
straight path $s$ between those two points, parameterized by $t\in[0,1]$,

$\boldsymbol{n}(t)=\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}-t\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}\,.$
(3.40)

The $l_{1}$-distance $D(t)$ of $\boldsymbol{n}(t)$ to $\mathbf{0}$ follows
directly as

$D(t)=\frac{\sqrt{n}}{3\pi^{2}}(mW_{\mathbf{0}})^{3/2}(1-t)$ (3.41)

and the functional’s concrete values $\mathcal{F}[\boldsymbol{n}(t)]$ along
that path can easily be calculated by exact numerical means. As second path
$\kappa$, we consider the one experimentally realized in Ref. [68] by
continuously reducing the coupling strength $\kappa$ of the pair interaction
from one to zero. Since the interaction Hamiltonian $\hat{W}$ in RDMFT is
fixed, this path has to be realized equivalently by increasing the strength of
the kinetic energy according to $p^{2}/2m\kappa$. The respective distance $D$
of $\boldsymbol{n}(\kappa)$ to $\mathbf{0}$ follows as

$D(\kappa)=\frac{\sqrt{n}}{3\pi^{2}}(mW_{\mathbf{0}}\kappa)^{3/2}$ (3.42)

and the functional $\mathcal{F}[\boldsymbol{n}(\kappa)]$ along that path is
given by

$\begin{split}\mathcal{F}[\boldsymbol{n}(\kappa)]&=4nNW_{\mathbf{0}}D(\kappa)\\\
\
&\quad+\frac{3^{2/3}4\pi^{7/3}Nn^{2/3}a_{1}}{m^{2}W_{\mathbf{0}}}D^{2/3}(\kappa)\,.\end{split}$
(3.43)

In Eq. (3.43) one may replace $W_{\mathbf{0}}$ by $a_{0}$ according to (3.38).

The result for $\mathcal{F}$ as a function of the fraction $D$ of non-
condensed bosons along the two paths $s$ and $\kappa$ is shown in Fig. 3.2 for
the parameters $n=10^{-3}$, $W_{\mathbf{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 $\mathcal{F}$ goes to zero for $D=0$ which corresponds to complete
BEC. Also, it can be seen that the gradient of $\mathcal{F}$ increases for
smaller distances $D$. In Sec. 3.7 we will show that the gradient of
$\mathcal{F}$ actually diverges in the limit $D\to 0$ and provide a detailed
discussion of this remarkable and far-reaching observation.

Figure 3.2: Bogoliubov functional $\mathcal{F}$ as a function of the relative
depletion $D$ along the straight path $s$ (3.40) and the curved path $\kappa$.
Dilute Bose gas in 3D for $n=10^{-3}$, $W_{\mathbf{0}}=m=1$ and $a_{1}=-0.01$
(left), charged Bose gas in 3D for $2m=e^{2}/2=1$ and $n=100$ (right).

A generalization of $\mathcal{F}$ given by Eq. (3.32) to dimensions $d\neq 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 (c.f. Sec. 3.1.4). Two-dimensional dilute systems
are weakly interacting if the condition $n|a|_{2\mathrm{D}}^{2}\ll 1$ [79, 80]
is satisfied where $a_{2\mathrm{D}}$ 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.
[81]. Due to their distinctive role, we will study a one-dimensional model in
Sec. 3.6.3.

####  3.6.2 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)\propto 1/r$. The respective Fourier coefficients
$W_{\boldsymbol{p}}$ can still be determined analytically though. In case of
an additional uniform background they follow as (see Appendix D)

$W_{\textbf{0}}=0\,,\quad W_{\boldsymbol{p}}=\frac{4\pi
e^{2}}{p^{2}}\,,\,\forall\boldsymbol{p}\neq\mathbf{0}\,.$ (3.44)

For charged bosons the weak interaction regime corresponds to the high density
limit [82, 83, 67]. This regime to which Bogoliubov’s approximation refers to
is characterized by a small “gas parameter”,
$r_{s}\equiv(3/4\pi)^{1/3}me^{2}n^{-1/3}\ll 1$. To illustrate again how RDMFT
works, we calculate the energy and momentum occupation numbers
$n_{\boldsymbol{p}}$ of the ground state for the most realistic case of a
kinetic energy given by
$\hat{t}=\sum_{\boldsymbol{p}}\frac{p^{2}}{2m}\hat{n}_{\boldsymbol{p}}$. For
this, we add the exact kinetic energy functional
$\mbox{Tr}_{1}[\hat{t}(\cdot)]$ to the universal interaction functional
$\mathcal{F}$ with Fourier coefficients $W_{\boldsymbol{p}}$ given by Eq.
(3.44). Then, we minimize the total energy functional with respect to all
$\boldsymbol{n}\in\triangle$, leading to the minimizer
$\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}$ which is given by Eq.
(3.36). Evaluating then the functional at
$\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}$ is straightforward (in
contrast to the dilute neutral Bose gas) and one finds (recall $n\equiv N/V$)

$\mathcal{F}(\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}})=\frac{2\Gamma\left(-\frac{1}{4}\right)\Gamma\left(\frac{7}{4}\right)Nn^{1/4}e^{5/2}m^{1/4}}{3\pi^{5/4}}$
(3.45)

and the respective fraction of non-condensed bosons $D=1-N_{\mathrm{BEC}}/N$
follows as

$\accentset{\rule{3.5pt}{0.5pt}}{D}\equiv
D(\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}})=-\frac{\Gamma\left(-\frac{3}{4}\right)\Gamma\left(\frac{5}{4}\right)m^{3/4}e^{3/2}}{4\pi^{7/4}n^{1/4}}\,.$
(3.46)

Eq. (3.46) verifies that the depletion of the condensate decreases with
increasing density $n$. Adding the kinetic energy to Eq. (3.45) leads to the
known result for the ground state energy ($\hbar=4\pi\epsilon_{0}=1$) [84]:

$E_{0}=-\frac{4\Gamma\left[-\frac{5}{4}\right]\Gamma\left[\frac{7}{4}\right]Nn^{1/4}m^{1/4}e^{5/2}}{3\pi^{5/4}}\,.$
(3.47)

As a consistency check, this confirms the correctness of Eq. (3.45).

Next, in analogy to Sec. 3.6.1 we consider again the straight path $s$ and the
curved path $\kappa$. The latter is again defined as the curve
$\boldsymbol{n}(\kappa)$ obtained by reducing by factor $\kappa\in[0,1]$ the
coupling strength of the Hamiltonian above (which led to the results Eqs.
(3.45), (3.46) and (3.47)). Evaluating the distance $D$ along the path
$\kappa$ yields

$D(\kappa)=\accentset{\rule{3.5pt}{0.5pt}}{D}\kappa^{3/4}$ (3.48)

and the functional $\mathcal{F}$ takes the values

$\mathcal{F}[\boldsymbol{n}(\kappa)]=qD^{1/3}(\kappa)\,,$ (3.49)

$q\equiv
2^{5/3}\Gamma\left(-\frac{1}{4}\right)\Gamma\left(\frac{7}{4}\right)Nn^{1/3}e^{2}/3\left[-\Gamma\left(-\frac{3}{4}\right)\Gamma\left(\frac{5}{4}\right)\right]^{1/3}\\!\pi^{2/3}$.
For the path $s$ we have $D(t)=(1-t)\accentset{\rule{3.5pt}{0.5pt}}{D}$ and
the concrete values of the functional $\mathcal{F}$ along that path can be
evaluated by exact numerical means. The right panel of Fig. 3.2 shows
$\mathcal{F}$ along the two paths $s$ and $\kappa$. The curves have
qualitatively similar shapes to those for the dilute neutral Bose gas shown on
the left of Fig. 3.2. This is not surprising because both setups correspond to
a weakly interacting system in which the Bogoliubov approximated functional
Eq. (3.32) is valid. However, we will see in Sec. 3.7 that the momentum
dependence of $W_{\boldsymbol{p}}$ can alter the behaviour of the gradient of
$\mathcal{F}$.

####  3.6.3 Bose-Hubbard model for five lattice sites

As a third example, we discuss in this section the one-dimensional 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 $\triangle$. Indeed, for
$L=5$ there are only two independent momentum occupation numbers due to the
general parity symmetry $n_{\boldsymbol{p}}=n_{-\boldsymbol{p}}$ and
normalization
$n_{\boldsymbol{0}}=N-\sum_{\boldsymbol{p}\neq\boldsymbol{0}}n_{\boldsymbol{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\pi\nu/L$ where $\nu$ takes integer values in the
interval described by $|\nu|\leq(L-1)/2$. The bosons interact via Bose-Hubbard
on-site interaction as described by the operator
$\frac{U}{2}\sum_{j=1}^{L}\hat{n}_{j}\left(\hat{n}_{j}-1\right)$. 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 (2.5), (2.65) any
non-negative prefactor could be separated from the interaction Hamiltonian
$\hat{W}$ and added instead in front of the respective universal functional,

$\mathcal{F}_{U\hat{W}}=|U|\,\mathcal{F}_{\mathrm{sgn}(U)\hat{W}}\,.$ (3.50)

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 (2.5), (2.65) to
a maximization. Because of this, we consider in the following the interaction
Hamiltonian
$\hat{W}\equiv\frac{\mathrm{sgn}(U)}{2}\sum_{j=1}^{L}\hat{n}_{j}\left(\hat{n}_{j}-1\right)$
and add eventually the coupling constant $|U|$ in front of the respective
universal functional $\mathcal{F}_{\hat{W}}$. To proceed, it is then an
elementary exercise to determine the corresponding Fourier coefficients
$W_{p}=\mathrm{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 (3.32) or its extension
(3.28) 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, one-dimensional systems
require high densities $n=N/L\gg 1$ to be weakly interacting [81].

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
$\mathcal{F}_{\hat{W}}$ (or decent approximations thereof) the ground state
energy of every Hamiltonian $\hat{H}(\hat{t})=\hat{t}+|U|\,\hat{W}$ can be
calculated with relatively little computational effort by minimizing the total
energy functional with respect to all $\boldsymbol{n}\in\triangle$. In that
sense, RDMFT represents a highly economic approach for solving simultaneously
the ground state problem for the entire class $\\{\hat{H}(\hat{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 operator $\hat{t}$. 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. [77, 78]). Nonetheless, we focus in the following on hopping just
between neighbouring sites at a rate $t\geq 0$, i.e., we choose
$\hat{t}=-2t\sum_{p}\big{(}\\!\cos(p)-1\big{)}\hat{n}_{p}$ and w.l.o.g. fix
$|U|\equiv 1$. In analogy to the most realistic dilute and charged Bose gas as
discussed in Sec. 3.6.1 and Sec. 3.6.2, respectively, we pick $t=U=1$ as a
reference point for further investigations and illustrations .

Figure 3.3: Domain $\triangle$ of the universal functional is shown for $L=5$
sites. Minimization of signs $(\sigma_{p_{1}},\sigma_{p_{2}})$ in (3.28)
partitions $\triangle$ into three cells (see text for details).

To illustrate and compare the Bogoliubov- (3.32) and the Girardeau-
approximated functionals (3.28) for the specific case of $L=5$ sites we first
need to execute the minimization of the sign factors in (3.28). As it is shown
in Appendix E, this can be done analytically due to the specific Fourier
coefficients and leads to a partitioning of the functional’s domain
$\triangle$ into three regions. Just for illustrative purposes, we present in
Fig. 3.3 the _entire_ domain $\triangle$ of the functional (3.28) (recall its
validity refers to the regime of BEC only) and the three cells which are
characterized by different minimizing sign configurations
$(\sigma_{p_{1}},\sigma_{p_{2}})$ in (3.28). As the two independent occupation
numbers we choose here the momenta $p_{1}=2\pi/5$ and $p_{2}=4\pi/5$ which can
take values $n_{p_{j}}\in[0,50]$. The vector $(0,0)$ corresponds to complete
BEC, i.e., $N_{\mathrm{BEC}}\equiv n_{0}=N=100$ and its vicinity represents
the BEC-regime to which our functionals refer to.

Figure 3.4: Left: Contour plot of the Bogoliubov-approximated functional
(3.32) for the Bose-Hubbard model with $N=100$ bosons on $L=5$ sites in the
BEC-regime of not too large depletion. Right: Girardeau’s extension (3.28) for
the same system.

In Fig. 3.4 we present now the Bogoliubov functional (3.32) and its extension
(3.28) 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 $\accentset{\rule{3.5pt}{0.5pt}}{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.4 as well. The corresponding degree of depletion,
$D=2.7\%$, justifies in retrospective the treatment of the interaction
$\hat{W}$ within the Bogoliubov theory and the usage of the functionals (3.32)
and (3.28), respectively. For stronger quantum depletion the two functionals
in Fig. 3.4 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_{\mathbf{0}}$ rather than its
approximation to $N$ in Eq. (3.28).

For the discussion in Sec. 3.7 of the new concept of a BEC force, we define in
Fig. 3.4 five qualitatively different paths towards the polytope boundary
$\partial\triangle$, all starting from the point
$\accentset{\rule{3.5pt}{0.5pt}}{n}$. The path denoted by $s$ corresponds to a
straight path towards complete BEC and $\kappa$ denotes the path where the
interaction strength of the model is reduced by increasing the kinetic energy
by a factor $1/\kappa$ with $\kappa\in[0,1]$. The path denoted by $a$ runs
perpendicular towards the hyperplane defined by $n_{p_{1}}+n_{p_{2}}=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. 3.5 we present the functional $\mathcal{F}$ as a
function of $D=1-N_{\mathrm{BEC}}/N$ along those five paths. The black dots
emphasize that the value of $\mathcal{F}$ at the boundary $\partial\triangle$
remains finite (quite in contrast to its derivative as shown and discussed in
the subsequent section).

Figure 3.5: Universal functional $\mathcal{F}$ for the Bose-Hubbard model
along the five paths defined in Fig. 3.4

The convexity of the curves in Fig. 3.5 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 behaviour and the repulsive
character of the functional’s gradient close to the boundary emerges from the
minimization of the sign factors in (3.28) leading to (3.31).

###  3.7 Bose-Einstein force

In this section we explore in more detail the behavior of the functional and
its gradient close to the boundary of their domain $\triangle$ 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. 3.7.2 for the three different systems studied in Sec.
3.6.

####  3.7.1 General results

The functional (3.32) which is based on the Bogoliubov approximation is convex
on its entire domain $\triangle$ (2.66). 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 (3.14) 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 behaviour of the universal functional in the
BEC-regime. While the functional itself remains finite even at the point
$\mathbf{0}\in\triangle$ 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 (3.32). To be more specific, approaching the
vertex $\mathbf{0}$ of the simplex $\triangle$ means to simultaneously send
all momentum occupation numbers $n_{\boldsymbol{p}}$ with
$\boldsymbol{p}\neq\mathbf{0}$ to zero. Taking then the derivative of (3.32)
(or of its extension (3.28)) with respect to $n_{\boldsymbol{p}}$ sufficiently
close to $\mathbf{0}$ yields in leading order

$\frac{\partial\mathcal{F}}{\partial
n_{\boldsymbol{p}}}(\boldsymbol{n})\sim-\frac{n|W_{\boldsymbol{p}}|}{2}\frac{1}{\sqrt{n_{\boldsymbol{p}}}}\,.$
(3.51)

It is worth noticing that the divergence of this derivative for
$n_{\boldsymbol{p}}\rightarrow 0$ is always repulsive for _any_ interaction
$\hat{W}$. This remarkable feature follows directly from the minimization of
the sign factors in (3.28), leading to (3.31). 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 this chapter 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 (3.51) implies that the point $\mathbf{0}$ of complete
condensation can never be reached, independent of the path towards
$\mathbf{0}$ that is envisaged. Since the functional $\mathcal{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 diverge
according (3.35) to enforce such a path towards $\mathbf{0}$.

We also would like to emphasize that the divergence of the gradient of
$\mathcal{F}$ along a straight path is always proportional to $1/\sqrt{D}$ and
its prefactor depends on the direction of the path, i.e., the angle at which
$\mathbf{0}$ is approached. To confirm this in a quantitative way, let us
consider a general straight path from a starting point
$\accentset{\rule{3.5pt}{0.5pt}}{n}$ in the regime of BEC to $\mathbf{0}$,
linearly parameterized by $t\in[0,1]$,

$\boldsymbol{n}(t)=(1-t)\,\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}\,.$
(3.52)

The degree $D$ (3.33) of quantum depletion along that path reduces according
to

$D(t)=(1-t)\,\accentset{\rule{3.5pt}{0.5pt}}{D}\equiv(1-t)\,\frac{1}{N}\sum_{\boldsymbol{p}\neq\mathbf{0}}\accentset{\rule{3.5pt}{0.5pt}}{n}_{\boldsymbol{p}}\,.$
(3.53)

The gradient of $\mathcal{F}$ projected onto that path is then nothing else
than the weighted sum of individual contributions (3.51) from every
$\boldsymbol{p}$,

$\displaystyle\frac{\partial\mathcal{F}}{\partial D}\Big{|}_{\mathrm{path}}$
$\displaystyle=$
$\displaystyle\nabla_{\\!\boldsymbol{n}}\mathcal{F}\cdot\frac{\partial\boldsymbol{n}}{\partial
D}\Big{|}_{\mathrm{path}}$ $\displaystyle=$
$\displaystyle\nabla_{\\!\boldsymbol{n}}\mathcal{F}\cdot\frac{\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}}}{\accentset{\rule{3.5pt}{0.5pt}}{D}}\sim-\frac{n}{2}\sum_{\boldsymbol{p}\neq\mathbf{0}}|W_{\boldsymbol{p}}|\sqrt{\frac{\accentset{\rule{3.5pt}{0.5pt}}{n}_{p}}{\accentset{\rule{3.5pt}{0.5pt}}{D}}}\,\frac{1}{\sqrt{D}}\,.$

This second key result in this chapter 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 [85].

####  3.7.2 Examples

In this section we illustrate the novel concept of a BEC force (3.7.1) for
various systems introduced in Sec. 3.6.

##### 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 (3.7.1). For both
systems, the derivative of $\mathcal{F}$ with respect to the degree $D$ of
quantum depletion along the path $s$ defined by (3.52) is given by Eq.
(3.7.1). The summation over $\boldsymbol{p}\neq\boldsymbol{0}$ can be
converted into an integral in the thermodynamic limit where $N\to\infty$,
$V\to\infty$ and $n=N/V=\mathrm{cst.}$. This eventually allows us (see
Appendix C) to obtain a compact analytic expression for the BEC force,

$\left.\frac{\mathrm{d}\mathcal{F}(\accentset{\rule{3.5pt}{0.5pt}}{\boldsymbol{n}})}{\mathrm{d}D}\right|_{s}\sim
N\left[\eta(a_{0},n,m)+\frac{2\pi
na_{1}}{m\sqrt{\accentset{\rule{3.5pt}{0.5pt}}{D}}}\right]\frac{1}{\sqrt{D}}$
(3.55)

where $\eta(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 (3.7.1) the BEC force is
always repulsive. Since only the second term in Eq. (3.55) is negative (recall
$a_{1}<0$) this imposes in turn a bound on the maximal valid distance
$\accentset{\rule{3.5pt}{0.5pt}}{D}$ of the starting point
$\accentset{\rule{3.5pt}{0.5pt}}{n}\in\triangle$ to the regime of complete
BEC.

For the charged Bose gas the last expression in the second line in Eq. (3.7.1)
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/\sqrt{\mbox{dist}(\boldsymbol{n},\partial\triangle)}$ along straight paths
reaching any _arbitrary_ point on the boundary $\partial\triangle$ in the
regime of BEC.

Moreover, we determine for both systems the BEC force along the curved path
$\boldsymbol{n}(\kappa)$ which is defined by reducing an additional coupling
constant $\kappa$ in front of $\hat{W}$ from one to zero. Since exactly this
path has been implemented in a very recent experiment [68] 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 $\kappa$-path
follows directly from differentiation of the expressions in (3.43) and (3.49),
respectively, leading to

$\left.\frac{\mathrm{d}\mathcal{F}_{\mathrm{dilute}}}{\mathrm{d}D}\right|_{\kappa}\propto-\frac{1}{D^{1/3}}$
(3.56)

and

$\left.\frac{\mathrm{d}\mathcal{F}_{\mathrm{charged}}}{\mathrm{d}D}\right|_{\kappa}\propto-\frac{1}{D^{2/3}}\,.$
(3.57)

Fig. 3.6 displays the BEC force along the straight $s$-path and the curved
$\kappa$-path for both ultracold gas systems. The linear behaviour shown in
this log-log plot confirms the algebraic dependence of the BEC force on the
degree $D$ of quantum depletion along both paths. For the dilute Bose gas, the
gradient of $\mathcal{F}$ according to Eq. (3.55) and Eq. (3.56) diverges
faster along the path $s$ than along the path $\kappa$. For the charged Bose
gas we observe the opposite behaviour.

Figure 3.6: BEC force $|\mathrm{d}\mathcal{F}/\mathrm{d}D|$ along the straight
path $s$ (blue) and the curved path $\kappa$ (red) is shown for the dilute
Bose gas in 3D with $n=10^{-3}$, $W_{\mathbf{0}}=m=1$ and $a_{1}=-0.01$ (left)
and for the charged Bose gas in 3D with $2m=e^{2}/2=1$ and $n=100$ (right).

##### Bose-Hubbard model

We illustrate the BEC force and the diverging behaviour of the functional’s
gradient close to the boundary $\partial\triangle$ of its domain in general
for the Bose-Hubbard model. For this we consider again as in Sec. 3.6.3 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.4. Since for all five paths the distance $D$ of the occupation number vector
$\boldsymbol{n}$ to $\mathbf{0}$ is monotonously decreasing we can parametrize
the functional’s derivative along each path by $D$. The respective results are
depicted in Fig. 3.7. There, the vertical solid lines correspond to the values
of $D$ at which the respective paths reach the boundary of $\triangle$ (see
also Figs. 3.4, 3.5).

Figure 3.7: Gradient of the universal functional $\mathcal{F}$ for the Bose-
Hubbard model along the five paths defined in Fig. 3.4. The results for
$\kappa$ and $s$ almost coincide.

We first observe that for all five paths $-\partial\mathcal{F}/\partial D$ is
diverging at the end point of each path on the boundary $\partial\triangle$.
As a rather elementary analysis reveals (not shown here) this divergence is
always proportional to
$1/\sqrt{\mbox{dist}(\boldsymbol{n},\partial\triangle)}$. As far as the four
straight paths $a,b,c,s$ are concerned, this was expected given the general
results of Sec. 3.7.1. In contrast to the two continuous Bose gases, however,
the same applies in the Bose-Hubbard model also for the curved $\kappa$-path
which is obtained by just reducing the coupling strength.

##  Chapter 4 Excited State RDMFT

In the previous two chapters, we focused on the foundations of ground state
RDMFT for bosons and successfully applied it to BECs. In this chapter, we
introduce a novel method, namely a bosonic ensemble RDMFT for excited states.
Just like ground state RDMFT, the development of excited state RDMFT requires
a solid mathematical foundation. Therefore, we first extend in Sec. 4.1 the
collection of mathematical tools, mainly from convex analysis, started in Sec.
2.1. Based on these concepts, we introduce the new
$\boldsymbol{\omega}$-ensemble RDMFT in Sec. 4.2. Furthermore, we apply a
relaxation scheme to turn $\boldsymbol{\omega}$-ensemble RDMFT into a viable
method in Sec. 4.3. The main task is to characterize the domain of the relaxed
$\boldsymbol{\omega}$-ensemble universal functional in Sec. 4.4. Since Sec.
4.4 is quite technical, we proceed in Sec. 4.5 by discussing several examples
and illustrations. Applying the excited state RDMFT to homogeneous Bose gases
allows us to derive an excited state universal functional for BECs in Sec.
4.6. Moreover, we illustrate the $\boldsymbol{\omega}$-ensemble functional and
its domain for the Bose-Hubbard Dimer in Sec. 4.7.

###  4.1 Mathematical Preliminaries

####  4.1.1 Vector majorization

For any $\boldsymbol{x}=(x_{1},x_{2},...,x_{d})\in\mathbb{R}^{d}$, we denote
the vector of decreasingly ordered entries by $\boldsymbol{x}^{\downarrow}$
such that $x_{1}^{\downarrow}\geq x_{2}^{\downarrow}\geq...\geq
x_{d}^{\downarrow}$. The superscript $\downarrow$ might be omitted in the
following sections if it is clear from the context which of the two vectors
$\boldsymbol{x}$, $\boldsymbol{x}^{\downarrow}$ is meant.

A vector $\boldsymbol{x}\in\mathbb{R}^{d}$ is said to be majorized by
$\boldsymbol{y}\in\mathbb{R}^{d}$, denoted by
$\boldsymbol{x}\prec\boldsymbol{y}$, if and only if the two conditions

$\displaystyle\sum_{i=1}^{k}x_{i}^{\downarrow}\leq\sum_{i=1}^{k}y_{i}^{\downarrow}\,,\quad
1\leq k\leq d-1$ (4.1)
$\displaystyle\sum_{i=1}^{d}x_{i}^{\downarrow}=\sum_{i=1}^{d}y_{i}^{\downarrow}$
(4.2)

hold. Replacing the equality in (4.2) by an inequality with $\leq$ leads to
the concept of weak majorization,
$\boldsymbol{x}\prec_{\mathrm{w}}\boldsymbol{y}$. Clearly, the original order
of the entries in $\boldsymbol{x}$, $\boldsymbol{y}$ does not enter the
majorization conditions and therefore it does not follow from
$\boldsymbol{x}\prec\boldsymbol{y}$ and $\boldsymbol{y}\prec\boldsymbol{x}$
that $\boldsymbol{x}=\boldsymbol{y}$ because the two vectors can have the same
entries but in a different order.

Hardy, Littlewood, and Pólya [86] provided an alternative definition of
majorization, than through the partial sums in Eq. (4.1) and Eq. (4.2), which
might be more practical for some purposes (see also Sec. 4.4.3). We first
recall that a matrix $P\in\mathbb{R}^{d\times d}$ is called doubly stochastic
if all its entries are nonnegative and the entries in each row and column sum
up to one. Then, there is a strong connection between doubly stochastic
matrices and the majorization of vectors provided by the following theorem
[86]: For any $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$ we have
$\boldsymbol{x}\prec\boldsymbol{y}$ if and only if
$\boldsymbol{x}=P\boldsymbol{y}$ for some doubly stochastic matrix
$P\in\mathbb{R}^{d\times d}$.

Since we intent to use the concept of majorization in the context of RDMFT and
thus apply it to quantum systems, we require a connection to quantum states
and in particular density matrices. Alberti and Uhlmann [87] proved that for
two density matrices
$\hat{\Gamma},\hat{\Gamma}^{\prime}:\mathcal{H}\to\mathcal{H}$ acting on a
$D$-dimensional Hilbert space $\mathcal{H}$, the density operator
$\hat{\Gamma}^{\prime}$ can be decomposed into a convex combination

$\hat{\Gamma}^{\prime}=\sum_{i=1}^{D}q_{i}\hat{U}_{i}^{\phantom{\dagger}}\hat{\Gamma}\hat{U}_{i}^{\dagger}\,,$
(4.3)

where $\hat{U}_{i}$ are unitary operators, $0\leq q_{i}\leq 1$ and
$\sum_{i=1}^{D}q_{i}=1$, if and only if the spectrum of
$\hat{\Gamma}^{\prime}$ is majorized by the spectrum of $\hat{\Gamma}$,
$\mathrm{spec}(\hat{\Gamma}^{\prime})\prec\mathrm{spec}(\hat{\Gamma})$.

A special subclass of doubly stochastic matrices are the permutation matrices
which are square matrices with exactly one entry equal to one in each row and
column, and zeros in all other entries. Birkhoff’s theorem [88, 89] then
states that the doubly stochastic matrices represent the convex hull of the
permutation matrices. Further, the permutation matrices constitute the
extremal points of the convex set of doubly stochastic matrices. For
$P\in\mathbb{R}^{d\times d}$, there exist $d!$ many permutation matrices
$\pi_{i},i=1,...,d!$.

####  4.1.2 Permutohedra and Rado’s theorem

A polyhedron is a subset $P\subset\mathbb{R}^{d}$ which is given by the
intersection of finitely many halfspaces and thus the set of solutions
$\boldsymbol{x}$ to finitely many inequalities
$A\boldsymbol{x}\leq\boldsymbol{y}$ where $A$ is a $m\times d$ matrix and
$\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$. Note that following this
definition any polyhedron is convex. A bounded polyhedron is called a
polytope. This leads to the so-called hyperplane representation
($H$-representation) of a polytope. In general, one finds arbitrary many
different $H$-representations for one polytope but there always exists a
unique minimal $H$-representation up to scaling. The inequalities belonging to
the minimal $H$-representation are called facet-defining. Moreover, every
$H$-representation can be turned into a vertex representation
($V$-representation) in which the polytope is given by the convex hull of a
finite set of extremal points called vertices. A special case of polytopes are
simplices and another important subclass of polytopes is the permutohedron
which is defined as

$P_{\boldsymbol{x}}\equiv\mathrm{conv}(\\{\pi(\boldsymbol{x})\,|\,\pi\in\mathcal{S}^{d}\\})$
(4.4)

for a vector $\boldsymbol{x}\in\mathbb{R}^{d}$ and $\mathcal{S}^{d}$ denotes
the group of permutations $\pi$ of $d$ elements. Thus, $P_{\boldsymbol{x}}$ is
given by the convex hull of all possible permutations of the entries of
$\boldsymbol{x}$. Using Eq. (4.4), we obtain the following theorem by Rado
[90] relating permutohedra to the majorization of vectors: Let
$\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$ be two vectors in
$\mathbb{R}^{d}$. Then, $\boldsymbol{x}\prec\boldsymbol{y}$ if and only if
$\boldsymbol{x}\in P_{\boldsymbol{y}}$. To prove Rado’s theorem we first show
that for $\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d}$ the majorization
$\boldsymbol{x}\prec\boldsymbol{y}$ implies $\boldsymbol{x}\in
P_{\boldsymbol{y}}$. It follows directly from combining the theorem by Hardy,
Littlewood, Pólya with Birkhoff’s theorem, that

$\begin{split}\boldsymbol{x}&=P\boldsymbol{y}=\sum_{i=1}^{d!}p_{i}(\pi_{i}\boldsymbol{y})\,,\end{split}$
(4.5)

where $\pi_{1},...,\pi_{d!}$ are the permutation matrices, $p_{i}\geq 0$ and
$\sum_{i=1}^{d!}p_{i}=1$. Therefore, $\boldsymbol{y}$ lies in the convex hull
of all permutations of the entries of $\boldsymbol{y}$ generating the set
$P_{\boldsymbol{y}}$. To prove the opposite direction $\boldsymbol{x}\in
P_{\boldsymbol{y}}\Rightarrow\boldsymbol{x}\prec\boldsymbol{y}$ one simply has
to reverse the above argument.

####  4.1.3 Generalization of Rayleigh-Ritz variational principle

One of the most well-known methods to derive the ground state energy of a
quantum system is the Rayleigh-Ritz variational principle, which we already
explained in context of the energy minimization in Eq. (2.25). Moreover, this
variational principle is the underlying concept of Levy’s constrained search
discussed in Sec. 2.5 providing a more viable formulation of RDMFT for ground
states than Gilbert’s theorem. Also, the RDMFT for excited states, which we
will develop in this chapter, is based on a variational principle and a
constrained search formalism. In 1988, Gross, Oliviera, and Kohn [27, 28, 29]
provided such a generalization of the Rayleigh-Ritz variational method for
DFT. In the following, we refer to this generalization as the GOK variational
principle. The GOK variational principle holds for observables $\hat{H}$
acting on a $D$-dimensional Hilbert space $\mathcal{H}$. We denote the
increasingly ordered eigenvalues by $E_{1}\leq E_{2}\leq...\leq E_{D}$. Let
$\boldsymbol{\omega}\in\mathbb{R}^{D}$ be a vector with decreasingly arranged
entries $\omega_{1}\geq\omega_{2}\geq...\geq\omega_{D}\geq 0$ and
$\sum_{j}\omega_{j}=1$. Then, the weighted sum of the eigenvalues follows from
[27, 28, 29]:

$E_{\boldsymbol{\omega}}\equiv\sum_{j=1}^{D}\omega_{j}E_{j}=\min_{\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega)}}\mbox{Tr}\left[\hat{\Gamma}\hat{H}\right]\,.$
(4.6)

The minimizer state in Eq. (4.6) is given by

$\hat{\Gamma}_{\boldsymbol{\omega}}=\sum_{j=1}^{D}\omega_{j}\mbox{$|\Psi_{j}\rangle$}\\!\mbox{$\langle\Psi_{j}|$}$
(4.7)

with $\hat{H}\mbox{$|\Psi_{j}\rangle$}=E_{j}\mbox{$|\Psi_{j}\rangle$}$, and
$\mathcal{E}^{N}(\boldsymbol{\omega)}$ denotes the set of all $N$-particle
density operators with spectrum $\boldsymbol{\omega}$. In the case of
degenerate eigenvalues, one has to assign an arbitrary ordering to the
orthonormal states in the corresponding degenerate subspace of states and keep
the labels fixed afterwards [27].

In the case that $\hat{H}$ is a Hamiltonian, the $E_{j}$’s are the
corresponding eigenenergies, as already anticipated by the notation. In
physics, one is often interested in the first few low-lying excited state
energies and thus chooses only finitely many weights to be non-zero. We denote
the number of non-vanishing weights by $r$. Then, knowing the expression for
$E_{\boldsymbol{\omega}}$ in Eq. (4.6) allows extracting the eigenvalues
$E_{j}$ for $j\leq r$ by taking appropriate derivatives. Alternatively, one
can evaluate $E_{\boldsymbol{\omega}}$ for different $\boldsymbol{\omega}$ and
apply appropriate gradient triangles to determine the eigenenergies. For
example, $r=2$ is already sufficient to calculate the energy gap $\Delta E$
between the ground state and the first excited state. In that case, we have
$\boldsymbol{\omega}=(\omega,1-\omega,0,...)$ and evaluating
$E_{\boldsymbol{\omega}}$ for $\omega=1$ and a second
$1/2\leq\omega^{\prime}<1$ yields

$\Delta E\equiv
E_{2}-E_{1}=\frac{E_{(\omega^{\prime},1-\omega^{\prime},0,...)}-E_{(1,0,...)}}{1-\omega^{\prime}}\,.$
(4.8)

###  4.2 Introduction of $\boldsymbol{\omega}$-ensemble RDMFT

In the following derivation of an RDMFT for excited states, we adopt the GOK
variational principle established for DFT and combine it with a constrained
search formalism in the spirit of Levy’s [53] and Valone’s [55] ideas.
However, this procedure is far from being trivial and several obstacles occur.
For fermions, such a theory was just recently proposed and explained in Ref.
[40]. In the following, we focus mainly on the bosonic version which uses the
same underlying concepts and ideas up to modifications due to the bosonic
statistics.

In Eq. (4.6), we already introduced the set
$\mathcal{E}^{N}(\boldsymbol{\omega)}$ of all $\boldsymbol{\omega}$-ensemble
$N$-particle density operators $\hat{\Gamma}$ which is defined by

$\mathcal{E}^{N}(\boldsymbol{\omega)}\equiv\\{\hat{\Gamma}\in\mathcal{E}^{N}\,|\,\mathrm{spec}^{\downarrow}(\hat{\Gamma})=\boldsymbol{\omega}\\}\,.$
(4.9)

In analogy to Eq. (2.31) and Eq. (2.32), we obtain by tracing out $N-1$
particles

$\mathcal{E}_{N}^{1}(\boldsymbol{\omega})=N\mbox{Tr}_{N-1}[\mathcal{E}^{N}(\boldsymbol{\omega})]\,.$
(4.10)

The 1RMDs $\hat{\gamma}\in\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$ are then
called $\boldsymbol{\omega}$-ensemble $N$-representable in contrast to the
ensemble $N$-representable 1RDMs $\hat{\gamma}\in\mathcal{E}_{N}^{1}$
introduced in Sec. 2.2. Thus, a practical description of the set
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$ requires the knowledge of these
additional constraints. For one non-vanishing weight, i.e.
$\boldsymbol{\omega}_{0}=(1,0,...)$, there exists only one minimizer state in
Eq. (4.6) which is pure. Then, Eq. (4.6) reduces to the minimization performed
to obtain the ground state energy and ground state 1RDM. From the definition
of $\mathcal{E}^{N}(\boldsymbol{\omega)}$ in Eq. (4.9) follows that for
$\boldsymbol{\omega}_{0}$, the set of all $\boldsymbol{\omega}$-ensemble
$N$-particle density operators reduces to

$\mathcal{E}^{N}(\boldsymbol{\omega}_{0})=\mathcal{P}^{N}\,.$ (4.11)

Moreover, by tracing out all except one particle we obtain

$\mathcal{E}_{N}^{1}(\boldsymbol{\omega}_{0})=\mathcal{P}_{N}^{1}\,,$ (4.12)

and thus the set of all pure state $N$-representable 1RDMs
$\hat{\gamma}\in\mathcal{P}_{N}^{1}$. Recall that for fermions, the set
$\mathcal{P}_{N}^{1}$ is usually unknown due to the generalized Pauli
constraints. This hints that for general $r$, the additional
$\boldsymbol{\omega}$-ensemble $N$-representability constraints increase the
complexity of the problem drastically. From a conceptual point of view, this
has an even more remarkable effect in case of bosons: Solving the pure state
$N$-representability problem for bosons is trivial because the ground state
occupation numbers are only restricted through positivity and normalization of
the 1RDM. However, for a bosonic
$\hat{\gamma}\in\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$, the description of
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$ leads to additional constraints on
the natural occupation numbers interpreted as generalized exclusion
constraints for bosons. Altogether, we conclude that finding the solution to
the $\boldsymbol{\omega}$ ensemble $N$-representability problem is an
extremely difficult task and tackle this problem in Sec. 4.3.

Next, we focus on some important properties of the two sets in Eq. (4.9) and
Eq. (4.10). First, note that the restriction of the set $\mathcal{E}^{N}$ of
all ensemble $N$-particle density operators to the subset
$\mathcal{E}^{N}(\boldsymbol{\omega})$ of all $N$-particle density operators
with fixed spectrum $\boldsymbol{\omega}$ is non-linear. As a result, the two
sets $\mathcal{E}^{N}(\boldsymbol{\omega})$ and
$\mathcal{E}_{N}^{1}(\boldsymbol{\omega})$ are in general not convex,
contrarily to $\mathcal{E}^{N}$ and $\mathcal{E}_{N}^{1}$. Furthermore, we can
exploit the invariance of the spectrum of a density operator under unitary
transformations to simplify the description of the set
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$. We observe that if
$\hat{\gamma}\in\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$ it follows that also
the unitary transformation $\hat{u}:\mathcal{H}_{1}\to\mathcal{H}_{1}$ of
$\hat{\gamma}$ satisfies
$\hat{\gamma}^{\prime}\equiv\hat{u}\hat{\gamma}\hat{u}^{\dagger}\in\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$.
To proof this statement we first lift the unitary transformation of
$\hat{\gamma}$ to the $N$-particle level through $\hat{\gamma}^{\prime}\equiv
N\mbox{Tr}_{N-1}[\hat{u}^{\otimes^{N}}\hat{\Gamma}(\hat{u}^{\dagger})^{\otimes^{N}}]$.
If $\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})$ is fulfilled, it
follows that
$\hat{\Gamma}^{\prime}\equiv\hat{u}^{\otimes^{N}}\hat{\Gamma}(\hat{u}^{\dagger})^{\otimes^{N}}\in\mathcal{E}^{N}(\boldsymbol{\omega})$
because the set $\mathcal{E}^{N}(\boldsymbol{\omega})$ is defined only through
the spectral constraint
$\mathrm{spec}^{\downarrow}(\hat{\Gamma})=\boldsymbol{\omega}$. Therefore, it
is invariant under unitary transformations of its elements. The conclusion
that $\hat{\gamma}^{\prime}\in\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$
follows directly from tracing out $N-1$ particles of $\hat{\Gamma}^{\prime}$.
Thus, to determine whether a 1RDM $\hat{\gamma}$ belongs to the set
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$ or not, it is sufficient to know
its spectrum. We will return to the consequences of this unitary invariance in
Sec. 4.4 and explain how it facilitates the description of
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$.

Inspired by Levy’s constrained search and the GOK variational principle, we
now introduce a constrained search on the set
$\mathcal{E}^{N}(\boldsymbol{\omega})$ such that

$\begin{split}E_{\boldsymbol{\omega}}(\hat{h})&=\min_{\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})}\mbox{Tr}_{N}[(\hat{h}+\hat{W})\hat{\Gamma}]\\\
\
&=\min_{\hat{\gamma}\in\mathcal{E}_{N}^{1}(\boldsymbol{\omega})}\Big{[}\min_{\mathcal{E}^{N}(\boldsymbol{\omega})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[(\hat{h}+\hat{W})\hat{\Gamma}]\Big{]}\\\
\
&=\min_{\hat{\gamma}\in\mathcal{E}_{N}^{1}(\boldsymbol{\omega})}\Big{[}\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]+\underbracket{\min_{\mathcal{E}^{N}(\boldsymbol{\omega})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}]}_{\equiv\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})}\Big{]}\,.\end{split}$
(4.13)

The above equation provides a definition of the new universal functional
$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$ in the same spirit as in Eq.
(2.5). Note that Eq. (4.13) does not require a generalization of the Gilbert
theorem to excited states to prove the existence of a
$\boldsymbol{\omega}$-ensemble functional. In that case, the domain would be
defined by all 1RDMs which follow as $\boldsymbol{\omega}$-minimizers for a
particular Hamiltonian $\hat{H}(\hat{h})$, leading again to a sort of modified
$v$-representability problem. However, the constrained search formalism in Eq.
(4.13) allows us to circumvent this $v$-representability problem from the very
beginning and avoid difficulties due to possible degenerate states by
extending the domain of $\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$ to
the set $\mathcal{E}_{N}^{1}(\boldsymbol{\omega})$ of all
$\boldsymbol{\omega}$-ensemble $N$-representable 1RDMs.

Due to the nonlinear restriction of $\mathcal{E}^{N}$ to
$\mathcal{E}^{N}(\boldsymbol{\omega})$, the $\boldsymbol{\omega}$-ensemble
functional $\mathcal{F}_{\boldsymbol{\omega}}$ is in general not convex on the
entire domain $\mathcal{E}_{N}^{1}(\boldsymbol{\omega})$. This shall be
contrasted with the ensemble ground state functional
$\mathcal{F}_{e}(\hat{\gamma})$ in Eq. (2.59). As a consistency check, we
observe that for $r=1$, the $\boldsymbol{\omega}$-ensemble functional
$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$ reduces to the pure state
universal functional

$\mathcal{F}_{p}(\hat{\gamma})=\mathcal{F}_{\boldsymbol{\omega}_{0}}(\hat{\gamma})\,.$
(4.14)

In the more interesting case $\boldsymbol{\omega}\neq\boldsymbol{\omega}_{0}$,
the minimization in the second line of Eq. (4.13) is performed over all
ensemble $N$-particle density operators
$\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})$ and thus the ensemble
functional $\mathcal{F}_{e}(\hat{\gamma})$ is related to
$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$ through

$\begin{split}\mathcal{F}_{e}(\hat{\gamma})&=\min_{\mathcal{E}^{N}\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}]\\\
\
&=\min_{\boldsymbol{\omega}}\min_{\mathcal{E}^{N}(\boldsymbol{\omega})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}]\\\
\
&\equiv\min_{\boldsymbol{\omega}}\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})\,.\end{split}$
(4.15)

In Eq. (4.15), we extend the domain of $\mathcal{F}_{\boldsymbol{\omega}}$ to
the full set $\mathcal{E}_{N}^{1}$ by defining

$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})\equiv\begin{cases}\min_{\mathcal{E}^{N}(\boldsymbol{\omega})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}],&\forall\hat{\gamma}\in\mathcal{E}_{N}^{1}(\boldsymbol{\omega})\,,\\\
\infty,&\forall\hat{\gamma}\notin\mathcal{E}_{N}^{1}(\boldsymbol{\omega})\,.\end{cases}$
(4.16)

Thus, it follows that the extension of the domain of
$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$ to the entire set of all
ensemble $N$-representable 1RDMs does not affect any minimization process to
calculate $E_{\boldsymbol{\omega}}(\hat{h})$ because the minimum will always
be attained in $\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$. It is worth
noticing that whenever a 1RDM $\hat{\gamma}$ is ground state
$v$-representable, the minimizer $\boldsymbol{\omega}_{\mathrm{min}}$ in Eq.
(4.15) is given by
$\boldsymbol{\omega}_{\mathrm{min}}=\boldsymbol{\omega}_{0}$.

###  4.3 Relaxation of $\boldsymbol{\omega}$-ensemble RDMFT

The $\boldsymbol{\omega}$-ensemble RDMFT introduced in the above section is
not practically feasible due to the involved $\boldsymbol{\omega}$-ensemble
$N$-representability constraints. For this purpose, we resort in this section
to the well-known convex relaxation method to circumvent this problem. It is
based on the fact, that it is always possible to replace a non-convex
minimization problem by the corresponding convex minimization problem (c.f.
Sec. 2.1.1 and Sec. 2.5). Applied to the constrained search formalism in Eq.
(4.13), this means replacing the $\boldsymbol{\omega}$-ensemble functional
$\mathcal{F}_{\\!\boldsymbol{w}}(\hat{\gamma})$ by its lower convex envelope

$\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})\equiv\mathrm{conv}\left(\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})\right)\,.$
(4.17)

In addition, the domain of the relaxed $\boldsymbol{\omega}$-ensemble
functional is given by the convex hull of
$\mathcal{E}^{1}_{N}(\boldsymbol{\omega})$, and thus

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\equiv\mathrm{conv}\left(\mathcal{E}_{N}^{1}(\boldsymbol{\omega})\right)\,.$
(4.18)

We will refer to a 1RDM $\hat{\gamma}\in\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$ in the following as being relaxed
$\boldsymbol{\omega}$-ensemble $N$-representable. Then, the weighted sum of
lowest eigenenergies, $E_{\boldsymbol{\omega}}$, follows from minimizing the
new energy functional $\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]+\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ as

$E_{\boldsymbol{\omega}}=\min_{\hat{\gamma}\in\accentset{\rule{3.5pt}{0.5pt}}{\mathcal{E}}_{N}^{1}(\boldsymbol{\omega})}\left[\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]+\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})\right]\,.$ (4.19)

The above expression yields not only $E_{\boldsymbol{\omega}}$ but also the
minimizer 1RDM
$\hat{\gamma}=N\mbox{Tr}_{N-1}[\hat{\Gamma}_{\boldsymbol{\omega}}]$
corresponding to the $\boldsymbol{\omega}$-minimizer
$\hat{\Gamma}_{\boldsymbol{\omega}}$ in Eq. (4.7). The convex relaxation has
two main advantages: First, every local minimum of $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ is also a global one facilitating
the minimization of the energy functional. However, even more important is
that we can find a concrete description of the compact convex set
$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ which is not
hampered by the $\boldsymbol{\omega}$-ensemble $N$-representability problem.
We present a concrete strategy to obtain $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$ in Sec. 4.4.

Alternatively, we can introduce the relaxed $\boldsymbol{\omega}$-ensemble
RDMFT following Valone’s ideas [55] of ensemble RDMFT by extending on the
$N$-particle level the set $\mathcal{E}^{N}(\boldsymbol{\omega})$ to its
convex hull,

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{N}\\!(\boldsymbol{w})\equiv\mathrm{conv}\left(\mathcal{E}^{N}(\boldsymbol{\omega})\right)\,.$
(4.20)

Before applying the constrained search formalism to $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{N}\\!(\boldsymbol{w})$, we further investigate the properties of
$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{N}\\!(\boldsymbol{w})$ as well as their
important consequences. Since according to Eq. (4.9), the non-convex set
$\mathcal{E}^{N}\\!(\boldsymbol{w})$ is fully characterized through the
spectral constraint
$\mathrm{spec}(\hat{\Gamma}^{\downarrow})=\boldsymbol{\omega}$, and the
spectrum is invariant under unitary transformations, we can apply Uhlmann’s
theorem in Sec. 4.1.1 to show that

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{N}\\!(\boldsymbol{w})=\bigcup\limits_{\boldsymbol{\omega}^{\prime}\prec\,\boldsymbol{\omega}}\mathcal{E}^{N}(\boldsymbol{\omega}^{\prime})\equiv\\{\hat{\Gamma}\in\mathcal{E}^{N}\,|\,\mathrm{spec}(\hat{\Gamma})\prec\boldsymbol{\omega}\\}\,.$
(4.21)

Thus, the set $\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{N}\\!(\boldsymbol{w})$ constitutes of all
$N$-particle density operators $\hat{\Gamma}$ whose spectrum is majorized by
$\boldsymbol{\omega}$. Since the partial trace map $\mathrm{Tr}_{N-1}[\cdot]$
is linear, we are allowed to change its order with the convex hull operation
$\mathrm{conv}(\cdot)$. From this consideration it follows immediately that

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})=N\mathrm{Tr}_{N-1}[\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{N}\\!(\boldsymbol{w})]=\bigcup\limits_{\boldsymbol{\omega}^{\prime}\prec\,\boldsymbol{\omega}}\mathcal{E}^{1}_{N}(\boldsymbol{\omega}^{\prime})\,.$
(4.22)

Moreover, the last equality in Eq. (4.22) has another striking consequence

$\boldsymbol{\omega}^{\prime}\prec\boldsymbol{\omega}\,\,\Leftrightarrow\,\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{\omega}^{\prime})\subset\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\,.$ (4.23)

Thus, the set of all relaxed $\boldsymbol{\omega}$-ensemble $N$-representable
1RDMs becomes smaller if a vector $\boldsymbol{\omega}^{\prime}$ is majorized
by another vector $\boldsymbol{\omega}$. We illustrate this inclusion relation
of the set $\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ for the Bose-Hubbard
dimer in Sec. 4.7. Next, recall from our discussion in context of Eq. (4.11)
and (4.12) that we recover the ground state RDMFT for
$\boldsymbol{\omega}=\boldsymbol{\omega}_{0}\equiv(1,0,...)$. Then, it follows
that

$\mathcal{E}_{N}^{1}=\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}_{N}^{1}(\boldsymbol{\omega}_{0})\equiv\mathrm{conv}(\mathcal{E}_{N}^{1}(\boldsymbol{\omega}_{0}))=\mathrm{conv}(\mathcal{P}_{N}^{1})\,.$
(4.24)

We prove the statement in Eq. (4.24) in a general way without distinguishing
between bosons and fermions. First, we show that
$\mathcal{E}_{N}^{1}\subset\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}_{N}^{1}(\boldsymbol{\omega}_{0})$. For every
$\hat{\gamma}\in\mathcal{E}_{N}^{1}$ we consider a $N$-particle operator
$\hat{\Gamma}$ with $\mathcal{E}^{N}\ni\hat{\Gamma}\mapsto\hat{\gamma}$. Due
to the Krein-Milman theorem every element in the convex and compact set
$\mathcal{E}^{N}$ can be expressed as a convex combination of its extremal
points which are given by pure states $\hat{\Gamma}_{i}\in\mathcal{P}^{N}$
such that $\hat{\Gamma}=\sum_{j}p_{j}\hat{\Gamma}_{j}$. Since $\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}^{N}(\boldsymbol{\omega}_{0})=\mathrm{conv}(\mathcal{E}^{N}(\boldsymbol{\omega}_{0}))=\mathrm{conv}(\mathcal{P}^{N})$,
it follows that $\hat{\Gamma}\in\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}^{N}(\boldsymbol{\omega}_{0})$ and thus $\hat{\gamma}\in\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{\omega}_{0})$ proving that
$\mathcal{E}_{N}^{1}\subset\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}_{N}^{1}(\boldsymbol{\omega}_{0})$. Next, we have to show that
$\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}_{N}^{1}(\boldsymbol{\omega}_{0})\subset\mathcal{E}_{N}^{1}$. Every
$\hat{\Gamma}\in\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}^{N}(\boldsymbol{\omega}_{0})$ with
$\hat{\Gamma}\mapsto\hat{\gamma}\in\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{E}$\kern
0.0pt}}}_{N}^{1}(\boldsymbol{\omega}_{0})$ can be written as a convex
combination $\hat{\Gamma}=\sum_{j}p_{j}\hat{\Gamma}_{j}$ with
$\hat{\Gamma}_{j}\in\mathcal{E}^{N}(\boldsymbol{\omega}_{0})=\mathcal{P}^{N}$.
This implies that $\hat{\Gamma}\in\mathcal{E}^{N}$ and therefore
$\hat{\gamma}\in\mathcal{E}_{N}^{1}$ which finishes the proof of Eq. (4.24).
Note that this elaborated proof is redundant for bosons since Eq. (4.24)
follows directly from $\mathcal{P}_{N}^{1}=\mathcal{E}_{N}^{1}$ derived in
Sec. 2.3.2.

In analogy to the definition of the ensemble functional
$\mathcal{F}_{e}(\hat{\gamma})$ in Eq. (2.59), we obtain the relaxed
functional $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ through the constrained search
formalism

$\begin{split}\min_{\,\scriptsize{\hbox{\vbox{ \hrule height=0.6pt\kern
0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{N}(\boldsymbol{w})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}]&=\min_{\begin{subarray}{c}\sum_{j}p_{j}\hat{\Gamma}_{j}\mapsto\hat{\gamma},\\\
\hat{\Gamma}_{j}\in\mathcal{E}^{N}(\boldsymbol{\omega})\end{subarray}}\sum_{j}p_{j}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma_{j}}]\\\
\ &=\min_{\begin{subarray}{c}\sum_{j}p_{j}\hat{\gamma}_{j}=\hat{\gamma},\\\
\hat{\gamma}_{j}\in\mathcal{E}_{N}^{1}(\boldsymbol{\omega})\end{subarray}}\,\,\min_{\mathcal{E}^{N}(\boldsymbol{\omega})\ni\hat{\Gamma}_{j}\mapsto\hat{\gamma}_{j}}\sum_{j}p_{j}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma_{j}}]\\\
\
&\equiv\min_{\begin{subarray}{c}\sum_{j}p_{j}\hat{\gamma}_{j}=\hat{\gamma},\\\
\hat{\gamma}_{j}\in\mathcal{E}_{N}^{1}(\boldsymbol{\omega})\end{subarray}}\sum_{j}p_{j}\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma}_{j})\\\
\
&=\mathrm{conv}\left(\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})\right)\\\
&\equiv\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})\,.\end{split}$ (4.25)

In the first line, we use that every $N$-particle density operator
$\hat{\Gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{N}\\!(\boldsymbol{w})$ can be written as a
convex combination $\hat{\Gamma}=\sum_{j}p_{j}\hat{\Gamma}_{j}$ of
$\hat{\Gamma}_{j}\in\mathcal{E}_{N}(\boldsymbol{\omega})$. In the third line,
we insert the definition of $\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$
leading directly to the expression for the lower convex envelope of
$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})$, defined as $\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ in agreement with Eq. (4.17). To
summarize, the constrained search formalism provides a concrete approach to
derive the relaxed functional resulting into a convex minimization problem
with

$\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})\equiv\min_{\,\scriptsize{\hbox{\vbox{
\hrule height=0.6pt\kern 0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{N}(\boldsymbol{w})\ni\hat{\Gamma}\mapsto\hat{\gamma}}\mbox{Tr}_{N}[\hat{W}\hat{\Gamma}]\,.$
(4.26)

Moreover, the domain of $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ is the compact, convex set
$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ of all relaxed
$\boldsymbol{\omega}$-ensemble $N$-representable 1RDMs. It is thus our next
task to determine a conclusive description of $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$, which provides a convenient procedure to
check whether a given 1RDM belongs to the set or not.

###  4.4 Characterization of $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$

In this section, we derive a systematic approach to characterize the set
$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$. To be more
specific, we will show that it is sufficient to describe a spectral polytope
which is only determined through the spectra of all 1RDMs
$\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$. We will then obtain
a concrete description of this spectral polytope in its vertex representation
in Sec. 4.4.1. Unfortunately, it turns out that the vertex representation is
in many cases still not sufficient for practical purposes. Thus, we need to
translate the vertex representation of the spectral polytope into its
halfspace representation in Sec. 4.4.2. Although both representations are
equivalent, as discussed in Sec. 4.1.2, it can be quite complicated to
transform the two representations into each other.

####  4.4.1 General procedure

Let us start with the general procedure to determine the set $\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$ which is based on various mathematical
concepts explained in Sec. 2.1 and Sec. 4.1.

##### Duality principle

The duality correspondence for convex sets, discussed in Sec. 2.1.3, states
that $\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ can be characterized
through its support function $\sigma_{\,\scriptsize{\hbox{\vbox{ \hrule
height=0.6pt\kern 0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{1}_{N}(\boldsymbol{w})}(h)$. We illustrate the duality
correspondence for the compact convex set $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$ in Fig. 4.1. Therefore, Eq. (2.10) shows
that the convex set $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$
follows from minimizing $\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]$ for all possible
hermitian one-particle Hamiltonians $\hat{h}$ and determining the
corresponding minimizers $\hat{\gamma}_{\hat{h}}$. Since
$\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]\equiv\langle\hat{h},\hat{\gamma}\rangle_{1}$,
the inner product $\langle\hat{h},\hat{\gamma}\rangle_{1}$ on the Euclidean
space of hermitian matrices defines a hyperplane in the set of 1RDMs
$\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$. Moreover, the one-
particle Hamiltonian $\hat{h}$ determines the normal vector of the hyperplane.
Therefore, minimizing $\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]$ over all
$\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ corresponds to
shifting the hyperplane in direction $-\hat{h}$ until it touches the boundary
of the compact convex set $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$. Taking the convex hull of all minimizers
$\hat{\gamma}_{\hat{h}}$ eventually leads to

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})=\mathrm{conv}\left(\left\\{\underset{\hat{\gamma}\in\,\scriptsize{\hbox{\vbox{
\hrule height=0.6pt\kern 0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{1}_{N}(\boldsymbol{w})}{\mathrm{arg\,min}}\,\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]\Big{|}\,\hat{h}\text{
herm.}\right\\}\right)\,.$ (4.27)

The extremal points of the set $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$
are given by those 1RDMs $\hat{\gamma}\in\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$, which follow as unique minimizers for a
particular choice of $\hat{h}$. If there exists more than one
$\hat{\gamma}_{\hat{h}}$ for a given $\hat{h}$, the boundary points are not
extremal and can be expressed as convex combinations of extremal elements in
$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$.

Figure 4.1: Left: Schematic illustration of the minimization over a convex set
$S$ as described in the text. Right: Performing the minimization for all
possible directions determines the boundary of $S$.

##### Lifting minimization to $N$-particle level

Since the set $\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{N}\\!(\boldsymbol{w})$ contains all
information about the spectral vector $\boldsymbol{\omega}$ on the
$N$-particle level, we lift the minimization to the $N$-boson level. This is
implemented through

$\min_{\hat{\gamma}\in\,\scriptsize{\hbox{\vbox{ \hrule height=0.6pt\kern
0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{1}_{N}(\boldsymbol{w})}\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]=\min_{\hat{\Gamma}\in\,\scriptsize{\hbox{\vbox{
\hrule height=0.6pt\kern 0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{N}(\boldsymbol{w})}\mbox{Tr}_{N}[\hat{h}\hat{\Gamma}]\,.$
(4.28)

For simplicity, we drop the identities attached to $\hat{h}$ due to the
lifting process to the $N$-particle level on the right hand side. Further, the
extremal points of $\mathcal{E}^{N}(\boldsymbol{\omega})$ and $\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{N}\\!(\boldsymbol{w})$ coincide, leading to

$\min_{\hat{\Gamma}\in\,\scriptsize{\hbox{\vbox{ \hrule height=0.6pt\kern
0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{N}(\boldsymbol{w})}\mbox{Tr}_{N}[\hat{h}\hat{\Gamma}]=\min_{\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})}\mbox{Tr}_{N}[\hat{h}\hat{\Gamma}]\,.$
(4.29)

The above equality is based on the observation that relaxing the minimization
from $\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})$ to
$\hat{\Gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{N}\\!(\boldsymbol{w})$ does not change the
outcome of the minimization. We then apply the GOK variational principle in
Eq. (4.6) to obtain

$\begin{split}\min_{\hat{\gamma}\in\,\scriptsize{\hbox{\vbox{ \hrule
height=0.6pt\kern 0.90417pt\hbox{\kern-1.39998pt$\mathcal{E}$\kern
0.0pt}}}}\normalsize^{1}_{N}(\boldsymbol{w})}\mbox{Tr}_{1}[\hat{h}\hat{\gamma}]&=\min_{\hat{\Gamma}\in\mathcal{E}^{N}(\boldsymbol{\omega})}\mbox{Tr}_{N}[\hat{h}\hat{\Gamma}]\\\
\ &=\sum_{j}\omega_{j}E_{j}\,.\end{split}$ (4.30)

Thus, the minimization on the one-particle level on the left hand-side of the
above equation, is determined through the weighted sum of the eigenenergies of
the one-particle Hamiltonian $\hat{h}$ on the $N$-particle level. It is thus
the next step to determine the eigenstates of $\hat{h}$ on the $N$-boson
Hilbert space $\mathcal{H}_{N}$ defined in Eq. (2.16). This also requires that
we restrict solely to bosonic quantum systems in the following.

##### Configuration states

Using an orthonormal basis $\\{\mbox{$|i\rangle$}\\}_{i=1}^{d}$ with
$d=\mathrm{dim}(\mathcal{H}_{1})$, the eigenstates of $\hat{h}$ on the
$N$-boson Hilbert space are given by the configuration states

$\mbox{$|\boldsymbol{i}\rangle$}\equiv\mbox{$|i_{1},...,i_{N}\rangle$}\equiv\hat{a}_{i_{1}}^{\dagger}...\hat{a}_{i_{N}}^{\dagger}\mbox{$|0\rangle$}\,,$
(4.31)

where $|0\rangle$ denotes the vacuum state and $\hat{a}_{i}^{\dagger}$ is a
bosonic creation operator creating a boson in state $|i\rangle$. The one-
particle Hamiltonian $\hat{h}$ on $\mathcal{H}_{1}$ in its spectral
decomposition reads

$\hat{h}\equiv\sum_{i=1}^{d}h_{i}\mbox{$|i\rangle$}\\!\mbox{$\langle i|$}\,.$
(4.32)

Since bosons are allowed to occupy the same quantum state, the set of all
possible configurations in the $N$-boson Hilbert space follows as

$\mathcal{I}_{N,d}\equiv\\{\boldsymbol{i}\equiv(i_{1},...,i_{N})|\,1\leq
i_{1}\leq i_{2}\leq...\leq i_{N}\leq d\\}\,.$ (4.33)

Note that referring to a $N$-boson configuration state
$|\boldsymbol{i}\rangle$ is always associated with a respective reference
basis on the one-particle Hilbert space $\mathcal{H}_{1}$. Further, in the
case of unique minimizers in Eq. (4.27), the eigenstates of $\hat{\gamma}_{h}$
and $\hat{h}$ have to be equal.

##### Spectral polytope

Next, we choose a fixed set of natural orbitals (NO) by resorting to the
unitary invariance

$\hat{u}\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\hat{u}^{\dagger}=\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\,,$ (4.34)

which holds for any unitary operator
$\hat{u}:\mathcal{H}_{1}\to\mathcal{H}_{1}$. Let $\\{\mathcal{B}_{1}\\}$
denote the set of all possible ordered orthonormal bases for
$\mathcal{H}_{1}$. Then, we split $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$
into subsets for each basis $\mathcal{B}_{1}\in\\{\mathcal{B}_{1}\\}$. As a
result, $\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ is given by

$\begin{split}\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})&\equiv\bigcup\limits_{\mathcal{B}_{1}}\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\Big{|}_{\mathrm{NOs}=\mathcal{B}_{1}}\\\ \
&=\mathrm{spec}^{\downarrow}\left(\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\right)\times\\{\mathcal{B}_{1}\\}\,,\end{split}$
(4.35)

where $\mathrm{spec}^{\downarrow}(\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w}))$ denotes the set of all possible spectra of
1RDMs $\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$
with decreasingly ordered entries. The larger set of all possible (unordered)
spectra is denoted by $\mathrm{spec}(\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w}))$ and we define the spectral polytopes

$\begin{split}\Sigma(\boldsymbol{\omega})&\equiv\mathrm{spec}\left(\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\right)\,,\\\ \
\Sigma^{\downarrow}(\boldsymbol{\omega})&\equiv\mathrm{spec}^{\downarrow}\left(\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})\right)\,.\end{split}$ (4.36)

The sets $\Sigma(\boldsymbol{\omega})$ and
$\Sigma^{\downarrow}(\boldsymbol{\omega})$ are called spectral polytopes
because they only depend on spectral constraints on the 1RDMs
$\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$. We will show below
that they are indeed polytopes.

To characterize the set $\,\hbox{\vbox{ \hrule height=0.6pt\kern
1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$
it is sufficient to determine the spectral polytope
$\Sigma(\boldsymbol{\omega})$. It is important to note that this finding
simplifies our task to determine the domain $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$ of the relaxed
$\boldsymbol{\omega}$-ensemble functional $\,\hbox{\vbox{ \hrule
height=0.6pt\kern 1.29167pt\hbox{\kern-1.99997pt$\mathcal{F}$\kern
0.0pt}}}_{\\!\boldsymbol{w}}(\hat{\gamma})$ tremendously because we can
restrict to a fixed choice of NOs. Next, we introduce the set

$\Delta=\\{\boldsymbol{\lambda}\in\mathbb{R}^{d}\,|\,N\geq\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{d}\geq
0\\}$ (4.37)

of all possible 1RDMs with decreasingly ordered natural occupation numbers to
relate the two sets $\Sigma(\boldsymbol{\omega})$ and
$\Sigma^{\downarrow}(\boldsymbol{\omega})$. Then, it follows that

$\Sigma^{\downarrow}(\boldsymbol{\omega})=\Sigma(\boldsymbol{\omega})\cap\Delta\,.$
(4.38)

Therefore, it is sufficient to determine the natural occupation number (NON)
vector $\boldsymbol{\lambda}$ for all $\boldsymbol{\omega}$-minimizers
$\hat{\Gamma}_{\boldsymbol{\omega}}$ for all possible choices of $\hat{h}$ in
order to obtain the set $\Sigma^{\downarrow}(\boldsymbol{\omega})$. In the
next step we show that it is indeed sufficient to determine
$\Sigma^{\downarrow}(\boldsymbol{\omega})$ instead of the full spectral
polytope $\Sigma(\boldsymbol{\omega})$.

##### Generating vertices

The class of all one-particle Hamiltonians $\hat{h}$ considered in the
minimization (4.27) can be restricted to those with arbitrary, but in the
following fixed, eigenbasis
$\mathcal{B}_{1}\equiv\\{\mbox{$|i\rangle$}\\}_{i=1}^{d}$ and
$h_{1}<h_{2}<...<h_{d}$. It follows that for a finite number $r$ of non-
vanishing weights $\omega_{j}$ there will be only a finite number $R$ of NON
vectors $\boldsymbol{v}^{(j)}$ with $j=1,2,...,R<\infty$. The entries of the
NON vectors $\boldsymbol{v}^{(j)}$ are by definition ordered decreasingly due
to $h_{i}<h_{i+1}$ and $\omega_{i}>\omega_{i+1}$. Thus, the spectral sets
$\Sigma(\boldsymbol{\omega})$ and $\Sigma^{\downarrow}(\boldsymbol{\omega})$
take indeed the form of convex polytopes and $\Sigma(\boldsymbol{\omega})$
follows as the convex hull of all possible permutations of the vector entries
of all $\boldsymbol{v}^{(j)}$, $j=1,...,R$:

$\Sigma(\boldsymbol{\omega})=\mathrm{conv}\Big{(}\big{\\{}\pi(\boldsymbol{v}^{(j)})\,\big{|}\,j=1,...,R,\pi\in\mathcal{S}^{d}\big{\\}}\Big{)}\,.$
(4.39)

Since the spectral polytope $\Sigma(\boldsymbol{\omega})\subset\mathbb{R}^{d}$
is invariant under permutations, we can from now on restrict to the set
$\Sigma^{\downarrow}(\boldsymbol{\omega})$, which contains all generating
vertices $\boldsymbol{v}^{(j)}$.

##### Partial ordering of configurations

The next step is to determine in a systematic way the excitation spectrum for
different choices of $N,d$ and $r$. For a system with $N$ bosons, all of them
can occupy the lowest orbital of a fixed eigenbasis
$\mathcal{B}_{1}\equiv\\{|i\rangle\\}_{i=1}^{d}$ and $h_{1}<h_{2}<...<h_{d}$.
Clearly, the lowest configuration for an arbitrary number of bosons is given
by $(1,...,1)$ followed by $(1,...,1,2)$. For more than two non-vanishing
weights we have to consider all possible excitations while taking the
indistinguishability of the bosons into account. To achieve this we first
introducing a total ordering on $\mathcal{I}_{N,d}$ through

$\boldsymbol{i}\leq_{\boldsymbol{h}}\boldsymbol{j}:\,\,\sum_{k=1}^{N}h_{i_{k}}\leq\sum_{k=1}^{N}h_{j_{k}}\,,$
(4.40)

where $\boldsymbol{h}\equiv\boldsymbol{h}^{\uparrow}\equiv(h_{1},...,h_{d})$,
$h_{1}<h_{2}<...<h_{d}$. Performing this ordering for all configurations in a
sequence of length $r$ leads for a fixed $\boldsymbol{h}$ to

$\boldsymbol{i}_{1}\leq_{\boldsymbol{h}}\boldsymbol{i}_{2}\leq_{\boldsymbol{h}}...\leq_{\boldsymbol{h}}\boldsymbol{i}_{r}\,.$
(4.41)

This ordering has to be performed for all possible, and thus in general
infinitely many $\boldsymbol{h}$. However, depending on the number $r$ of non-
vanishing weights, many one-particle Hamiltonians lead to the same sequence
and thus are considered as equivalent. As a result, we only obtain a finite
number $R$ of distinct sequences and to each sequence we assign a vector
$\boldsymbol{v}$ according to

$\boldsymbol{v}=\sum_{j=1}^{r}\omega_{j}\boldsymbol{n}_{\boldsymbol{i}_{j}}\,.$
(4.42)

This NON vector $v$ contains all the information about the occupancies in the
consecutive configurations of the sequence. Further, the $k$-th entry of the
occupation number vector $\boldsymbol{n}_{\boldsymbol{i}}$ of the
configuration $\boldsymbol{i}$ can take values $0\leq
n_{\boldsymbol{i}}^{(k)}\leq N$. Thus, the lowest configuration is given by
$\boldsymbol{n}_{(1,...,1)}=(N,0,...)$.

To determine all possible sequences for a given number $r$ of non-vanishing
weights, we introduce the following partial ordering on the set
$\mathcal{I}_{N,d}$

$\boldsymbol{i}\leq\boldsymbol{j}:\,\,\Leftrightarrow\boldsymbol{i}\leq_{\boldsymbol{h}}\boldsymbol{j}\,,\quad\forall\boldsymbol{h}\,.$
(4.43)

From the ordering $h_{1}<h_{2}<...<h_{d}$ and Eq. (4.40) it follows
immediately that $\boldsymbol{i}\leq\boldsymbol{j}$ is fulfilled if and only
if $i_{k}\leq j_{k}$ for all $k=1,2,...,N$. Moreover, the condition $i_{k}\leq
j_{k}$ for all $k=1,2,...,N$ holds if and only if
$\sum_{k=1}^{N}h_{i_{k}}\leq\sum_{k=1}^{N}h_{j_{k}}$ for all $\boldsymbol{h}$.
Thus, the order described by Eq. (4.43) is determined by the eigenenergies of
the one-particle Hamiltonian. The so-called Gale poset [91] in Eq. (4.43)
results into lineups of length $r$. The number of lineups is equal to the
number $R$ of generating vertices defined in Eq. (4.42) which, in turn,
determine the vertices of the spectral polytope $\Sigma(\boldsymbol{\omega})$
according to Eq. (4.39).

Let us now summarize in detail the procedure to construct the sequence

$\boldsymbol{h}\mapsto\hat{\Gamma}\mapsto\hat{\gamma}\mapsto\boldsymbol{v}$
(4.44)

needed to obtain the spectral polytope $\Sigma(\boldsymbol{\omega})$ from the
partial ordering in Eq. (4.43) in a systematic way. Every vector
$\boldsymbol{h}$ determines a lineup

$\boldsymbol{i}_{1}\to\boldsymbol{i}_{2}\to\ldots\to\boldsymbol{i}_{r}$ (4.45)

of configurations. It is important to notice that, even if there exist
uncountably many $\boldsymbol{h}$, we only obtain a finite number $R$ of
lineups since many $\boldsymbol{h}$ are equivalent concerning the linear order
introduced in Eq. (4.40). For every lineup, the $N$-boson density operator
$\hat{\Gamma}$ follows as

$\hat{\Gamma}=\sum_{j=1}^{r}w_{j}\mbox{$|\boldsymbol{i}_{j}\rangle$}\\!\mbox{$\langle\boldsymbol{i}_{j}|$}\,,$
(4.46)

and tracing out $N-1$ particles yields the corresponding $\hat{\gamma}$. Since
the occupation number vectors $\boldsymbol{n}_{\boldsymbol{i}}$ introduced in
Eq. (4.42) are nothing else than the spectrum of $\hat{\gamma}$, we eventually
obtain

$\boldsymbol{v}=\sum_{j=1}^{r}w_{j}\mathrm{spec}\left(N\mbox{Tr}_{N-1}\left[\mbox{$|\boldsymbol{i}_{j}\rangle$}\\!\mbox{$\langle\boldsymbol{i}_{j}|$}\right]\right)\,.$
(4.47)

To illustrate the partial ordering in Eq. (4.43), we show in Fig. 4.2 the
excitation spectrum for $N=3$ bosons. The first configuration corresponds to
the ground state where all three bosons occupy the same state. The
corresponding lineup consists of one configuration, namely $(1,1,1)$. Also the
excitation pattern for only one excitation in the system is unique and is
represented by the lineup $(1,1,1)\to(1,1,2)$. Thus, for $r=1$ and $r=2$ we
obtain in each case one NON vector according to Eq. (4.42). For a larger
number of excitations, it depends on the values of the eigenenergies $h_{i}$
which configuration follows. For $r=3$, this consideration results into two
possible lineups

$\displaystyle(1,1,1)\to(1,1,2)\to(1,1,3)\,,$ (4.48)
$\displaystyle(1,1,1)\to(1,1,2)\to(1,2,2)\,.$ (4.49)

This also implies that in every excitation pattern with $r\geq 4$, both
configurations $(1,1,3)$ and $(1,2,2)$ must appear before $(1,2,3)$, but the
exact sequence depends on the particular choice of $\boldsymbol{h}$. We
discuss all those examples in more detail in Sec. 4.5.

Figure 4.2: Illustration of the excitation spectrum and linear ordering for
$N=3$ bosons.

It is worth noticing, that the corresponding excitation spectrum for fermions
shown in Ref. [92] can be obtained from the bosonic excitation spectrum in
Fig. 4.2 by adding the tuple $(0,1,2)$ to each bosonic configuration. For
example, the lowest configuration for fermions is given by $(1,2,3)$ followed
by $(1,2,4)$.

####  4.4.2 Hyperplane representation of spectral polytopes

In the above section, we succeeded in determining all generating vertices
$\boldsymbol{v}^{(j)}$, $j=1,...,R$, and thus the spectral set
$\Sigma(\boldsymbol{\omega})$ through Eq. (4.39) leading to $\,\hbox{\vbox{
\hrule height=0.6pt\kern 1.29167pt\hbox{\kern 0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})$. However, in practice, this vertex
representation of the spectral polytope is not very convenient since it is
extremely complicated to check whether a NON vector $\boldsymbol{\lambda}$
lies inside the polytope or not. Therefore, the vertex description is of
little use in practice. In contrast, the hyperplane representation of the
permutation-invariant polytope $\Sigma(\boldsymbol{\omega})$ provides a more
convenient way to test a given NON vector because one has to check only
finitely many linear hyperplane conditions $D_{k}(\boldsymbol{\lambda})\geq
0$. Thus, in the next step, we need to turn the vertex representation of the
spectral polytope $\Sigma(\boldsymbol{\omega})$ into a hyperplane
representation.

The procedure for $r=1,2$ is trivial because in each case there is only one
lineup, as discussed in the section above. Consequently, we also have only one
generating vertex $\boldsymbol{v}$. The spectral polytope
$\Sigma(\boldsymbol{\omega})$ is then obtained by considering all possible
permutations of the entries of $\boldsymbol{v}$. Note that even for $R=1$, the
total number of vertices in $\Sigma(\boldsymbol{\omega})$ can become
arbitrarily large because it depends on the dimension $d$ of the one-boson
Hilbert space. As a result, already $r=1,2$ demonstrates clearly that the
vertex representation is not efficient to test whether a given occupation
number vector $\boldsymbol{\lambda}$ lies inside the permutohedron
$\Sigma(\boldsymbol{\omega})$ or not. However, for maximally two non-vanishing
weights, we can directly apply Rado’s theorem [90] explained in Sec. 4.1.2 to
obtain

$\Sigma(\boldsymbol{\omega})=\\{\boldsymbol{\lambda}\in\mathbb{R}^{d}\,|\,\boldsymbol{\lambda}\prec\boldsymbol{v}\\}\,.$
(4.50)

The spectral polytope $\Sigma(\boldsymbol{\omega})$ is then determined by the
$d$ inequalities following from the majorization condition
$\boldsymbol{\lambda}\prec\boldsymbol{v}$. For more than one lineup, the
spectral polytope $\Sigma(\boldsymbol{\omega})$ defined in Eq. (4.39) is not a
permutohedron anymore, and thus Rado’s theorem does not apply. However,
$\Sigma(\boldsymbol{\omega})$ is still permutation-invariant. To translate its
vertex into a hyperplane representation, we thus require a modification of
Rado’s theorem adjusted to the new situation. Such a crucial generalization of
Rado’s theorem, introduced and proven in Ref. [92], states that for finitely
many vectors $\boldsymbol{v}^{(1)},...,\boldsymbol{v}^{(R)}\in\mathbb{R}^{d}$,
the polytope

$\mathcal{P}=\mathrm{conv}\left(\left\\{\pi(\boldsymbol{v}^{(j)})\,\Big{|}\,j=1,...,R,\pi\in\mathcal{S}^{d}\right\\}\right)$
(4.51)

is equivalent to

$\mathcal{P}=\Big{\\{}\boldsymbol{\lambda}\,\Big{|}\,\exists\,\,\text{conv.
comb.
}\sum_{j=1}^{R}p_{j}\boldsymbol{v}^{(j)}\equiv\boldsymbol{v}:\boldsymbol{\lambda}\prec\boldsymbol{v}\Big{\\}}\,.$
(4.52)

In Sec. 4.5.3, we show for the concrete example of $N=3$ and $r=3$ how the
generalized Rado theorem in Eq. (4.52) can be used to determine a minimal
hyperplane description up to scaling. The minimal hyperplane representations
for $r>3$ require a mathematically more complex derivation presented in Ref.
[93], introducing a mathematically rigorous derivation of the unique minimal
hyperplane representation for arbitrary values of $N$ and $d$, based on the
concept of the normal fan of a polytope and its extremal rays. However, $r\leq
2$ already incorporates the ground state energy and its gap, and large values
of $r$ are often less interesting from the physical point of view.

Furthermore, the system stabilizes in the bosonic case for $N\geq r-1$ and
$d\geq r$. This means, that for sufficiently large $N$ and $d$, increasing the
particle number $N$ or dimension $d$ does not result into new inequalities
[93]. We comment on the relation between the spectral potytopes for different
settings $(N,d)$ in more detail in Sec. 4.4.3, and several examples are also
provided in Sec. 4.5.

Remarkably, increasing $r$ leads to a hierarchy of inequalities [93] which
means that for increasing $r$, all inequalities obtained for a smaller $r$ are
still facet-defining and only finitely many additional inequalities occur. It
is worth noticing that these additional constraints are entirely independent
of the interaction itself and arise directly from the geometry of the
underlying set of density matrices. This shall be contrasted with the quantum
depletion due to interactions in Ch. 3. Thus, these additional constraints are
interpreted as a generalized Pauli exclusion principle for bosons. The number
of generating vertices $\boldsymbol{v}^{(j)}$ and the corresponding number of
inequalities for $r=1,...,12$ [93] are presented in Tab. 4.1.

$r$ | $1$ | $2$ | $3$ | $4$ | $5$ | $6$ | $7$ | $8$ | $9$ | $10$ | $11$ | $12$
---|---|---|---|---|---|---|---|---|---|---|---|---
$\\#\boldsymbol{v}^{(j)}$ | $1$ | $1$ | $2$ | $4$ | $8$ | $17$ | $37$ | $82$ | $184$ | $418$ | $967$ | $2278$
$\\#\text{ineq}$ | $1$ | $2$ | $3$ | $5$ | $8$ | $13$ | $22$ | $36$ | $59$ | $99$ | $171$ | $299$
Table 4.1: Number of generating vertices $\boldsymbol{v}^{(j)}$ and number of
facet-defining inequalities of $\Sigma^{\downarrow}(\boldsymbol{\omega})$ for
$r\leq 12$.

####  4.4.3 Generalization to larger particle numbers

In this section, we derive a concrete, general relation between the spectral
polytopes for different particle numbers. This will further explain why no new
inequalities appear in the minimal hyperplane representation for increasing
particle number $N$ (and $d$). We keep $d$ fixed from now on and only vary
$N$. Moreover, we restrict our discussion to large enough $N$ and $d$ such
that the number of configurations, and thus the number of lineups, are
independent of $N$ and $d$. Since in the lowest configuration all bosons
occupy the same orbital, the highest occupied orbital can be $i_{N}=r$ for a
fixed number of non-vanishing weights $r$ (c.f. Fig. 4.2). Thus, for the
lineups to be independent of $N$ and $d$, they must fulfil $N\geq r-1$ and
$d\geq r$.

For $N\geq r-1$, increasing the particle number to $N^{\prime}>N$ amounts to
placing more bosons in the lowest orbital. Thus, it does not change the number
of lineups, and for each of them it only changes the first entry of the
corresponding NON vector $\boldsymbol{v}^{(i)}$ which depends explicitly on
$N$. This also implies that for fixed $d\geq r$ the generating vertices
$\boldsymbol{v}^{(i)}$ for $N$ and $N^{\prime}>N$ are related through

$\boldsymbol{v}^{(i)}_{N}=\boldsymbol{v}^{(i)}_{N^{\prime}}-\delta\boldsymbol{e}_{1}\,,$
(4.53)

where $\delta=N^{\prime}-N$. Thus, generalizing the vertex representation of
the spectral polytope $\Sigma(\boldsymbol{\omega})$ from $N$ to $N^{\prime}$
is trivial. A higher dimensional setting with $d^{\prime}>d$ can be obtained
from $\Sigma_{N}(\boldsymbol{\omega})$ for N particles by first extending the
$d$-dimensional vector $\boldsymbol{v}^{(i)}_{N}$ to a
$d^{\prime}$-dimensional vector by adding zero entries and only afterwards
using the relation in Eq. (4.53).

Recall that according to Eq. (4.39), we obtain the spectral polytopes
$\Sigma(\boldsymbol{\omega})$ by considering all possible permutations of the
entries of all generating vertices. The key result, which we will prove in
this section, is that the spectral polytope
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ with $\delta=N^{\prime}-N>0$ is the
Minkowski sum of $\Sigma_{N}(\boldsymbol{\omega})$ and a permutation-invariant
rescaled simplex

$\Sigma_{N^{\prime}}(\boldsymbol{\omega})\equiv\Sigma_{N}(\boldsymbol{\omega})+\mathcal{C}\,,$
(4.54)

where $\mathcal{C}$ is a rescaled simplex given by

$\mathcal{C}\equiv\mathrm{conv}(\\{\pi(\delta\boldsymbol{e}_{1})\,|\,\pi\in\mathcal{S}^{d}\\})\,.$
(4.55)

To prove Eq. (4.54), we start by deriving the relation between the polytope
$\Sigma_{N}(\boldsymbol{\omega})$ and
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ for only one generating vertex
$\boldsymbol{v}$ which corresponds to $N$ bosons. Due to Eq. (4.53), we also
have only one generating vertex $\boldsymbol{v}+\delta\boldsymbol{e}_{1}$ for
$N^{\prime}$ bosons. According to a theorem by Hardy, Littlewood and Pólya
[86] presented in Sec. 4.1.1, we have
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})\ni\boldsymbol{\mu}\prec\boldsymbol{v}+\delta\boldsymbol{e}_{1}$
if and only if $\boldsymbol{\mu}=P(\boldsymbol{v}+\delta\boldsymbol{e}_{1})$
for some double stochastic matrix $P\in\mathbb{R}^{d\times d}$. Combining this
statement with Birkhoff’s theorem (see Sec. 4.1.1) leads to

$\Sigma_{N^{\prime}}(\boldsymbol{\omega})\ni\boldsymbol{\mu}=\sum_{i=1}^{d!}p_{i}\left(\pi_{i}(\boldsymbol{v}+\delta\boldsymbol{e}_{1})\right)=\sum_{i=1}^{d!}p_{i}(\pi_{i}\boldsymbol{v})+\delta\sum_{i=1}^{d!}p_{i}(\pi_{i}\boldsymbol{e_{1}})\,,$
(4.56)

where $\pi_{1},...,\pi_{d!}$ are permutation matrices, $p_{i}\geq 0$ and
$\sum_{i=1}^{d!}p_{i}=1$. It follows that
$\boldsymbol{\mu}\in\Sigma_{N}(\boldsymbol{\omega})+\mathrm{conv}(\\{\pi(\delta\boldsymbol{e}_{1})|\pi\in\mathcal{S}^{d}\\}$
where the set

$\mathcal{C}\equiv\mathrm{conv}(\\{\pi(\delta\boldsymbol{e}_{1})\,|\,\pi\in\mathcal{S}^{d}\\})$
(4.57)

is a rescaled simplex with edge length $\delta$. Then,
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})=\Sigma_{N}(\boldsymbol{\omega})+\mathcal{C}=\\{\boldsymbol{\lambda}+\boldsymbol{c}\,|\,\boldsymbol{\lambda}\in\Sigma_{N}(\boldsymbol{\omega}),\boldsymbol{c}\in\mathcal{C}\\}$
is the Minkowski sum of two permutation invariant polytopes. Clearly, the sum
of two convex sets is also convex. Also, every
$\boldsymbol{\mu}\in\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ is correctly
normalized since all $\boldsymbol{\lambda}\in\Sigma_{N}(\boldsymbol{\omega})$
are normalized to $N$ and all $\boldsymbol{c}\in\mathcal{C}$ are normalized to
$\delta=N^{\prime}-N$.

Next, we consider the more general case of more than one generating vertex,
i.e. $R>1$. Using the generalized Rado theorem in Eq. (4.52), we show that
vector $\boldsymbol{\mu}\in\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ can be
related to a vector $\boldsymbol{\lambda}\in\Sigma_{N}(\boldsymbol{\omega})$,
where

$\begin{split}\Sigma_{N}(\boldsymbol{\omega})&=\Big{\\{}\boldsymbol{\lambda}\,\Big{|}\,\exists\,\,\text{conv.
comb.
}\sum_{j=1}^{R}p_{j}\boldsymbol{v}^{(j)}\equiv\boldsymbol{v}:\boldsymbol{\lambda}\prec\boldsymbol{v}\Big{\\}}\,,\\\
\
\Sigma_{N^{\prime}}(\boldsymbol{\omega})&=\Big{\\{}\boldsymbol{\mu}\,\Big{|}\,\exists\,\,\text{conv.
comb.
}\sum_{j=1}^{R}p_{j}(\boldsymbol{v}^{(j)}+\delta\boldsymbol{e}_{1})\equiv\boldsymbol{u}:\boldsymbol{\mu}\prec\boldsymbol{u}\Big{\\}}\,.\end{split}$
(4.58)

From Eq. (4.52) follows that every
$\boldsymbol{\mu}\in\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ can be written
as a convex combination

$\begin{split}\Sigma_{N^{\prime}}(\boldsymbol{\omega})\ni\boldsymbol{\mu}&=\sum_{i=1}^{R}\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}\pi(\boldsymbol{v}^{(i)}+\delta\boldsymbol{e}_{1})=\sum_{i=1}^{R}\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}\left(\pi(\boldsymbol{v}^{(i)})+\pi(\delta\boldsymbol{e}_{1})\right)\\\
\
&\equiv\sum_{i=1}^{R}p_{i}\boldsymbol{\lambda}_{i}+\delta\sum_{i=1}^{R}\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}\pi(\boldsymbol{e}_{1})\\\
\
&=\boldsymbol{\lambda}+\delta\sum_{i=1}^{R}\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}\pi(\boldsymbol{e}_{1})\,,\end{split}$
(4.59)

where $\boldsymbol{\lambda}\in\Sigma_{N}(\boldsymbol{\omega})$, and we defined
$p_{i}\equiv\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}$ and
$\boldsymbol{\lambda}_{i}\equiv\sum_{\pi\in\mathcal{S}^{d}}q_{i,\pi}\pi(\boldsymbol{v}^{(i)})/p_{i}$.
Then, we eventually obtain
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})\equiv\Sigma_{N}(\boldsymbol{\omega})+\mathcal{C}$,
where $\mathcal{C}$ is the rescaled simplex given by Eq. (4.57) with edge
length $\delta$. Hence, the relation between the spectral polytopes belonging
to different total particle numbers stated in (4.54) holds for any $R>1$. Note
that this is in agreement with the observation that every extremal point in
$\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ can be obtained from adding a
vertex in $\Sigma_{N}(\boldsymbol{\omega})$ and a vertex in $\mathcal{C}$.

Eq. (4.54) can now be used to explain how the hyperplane representation for
$N^{\prime}>N$ follows from the inequalities for the smaller setting $N$.
Combining $\boldsymbol{\mu}\prec\boldsymbol{\mu}^{\downarrow}$ with

$\begin{split}&\sum_{j=1}^{k}\mu_{j}^{\downarrow}=\sum_{j=1}^{k}(\lambda_{j}+c_{j})^{\downarrow}\leq\sum_{j=1}^{k}\lambda_{j}^{\downarrow}+c_{j}^{\downarrow}\,,\quad
1\leq k\leq d-1\,,\\\ \
&\sum_{j=1}^{d}\mu_{j}^{\downarrow}=\sum_{j=1}^{d}\lambda_{j}^{\downarrow}+c_{j}^{\downarrow}\end{split}$
(4.60)

leads to the majorization condition
$\boldsymbol{\mu}\prec\boldsymbol{\lambda}^{\downarrow}+\boldsymbol{c}^{\downarrow}$.
Since in addition $0\leq c_{j}\leq\delta$ for all $j=1,...,d$, we have

$\boldsymbol{\mu}\prec\boldsymbol{\lambda}^{\downarrow}+\delta\boldsymbol{e}_{1}$
(4.61)

and it follows immediately that

$\lambda_{1}^{\downarrow}\geq\mu_{1}^{\downarrow}-\delta\,.$ (4.62)

The above equation can then be used together with Eq. (4.61) for the other
vector entries to replace $\boldsymbol{\lambda}^{\downarrow}$ in the linear
constraints $D_{k}(\boldsymbol{\lambda}^{\downarrow})\geq 0$ by
$\boldsymbol{\mu}^{\downarrow}-\delta\boldsymbol{e}_{1}$ to obtain the set of
linear constraints for the larger setting with $N^{\prime}$ bosons. Since
increasing the particle number does not lead to new generating vertices and
those are related for different particles numbers through Eq. (4.53), also no
additional inequalities appear as already explained in Sec. 4.4.2. For a
mathematically rigorous approch including the extension of the one-particle
Hilbert space to larger dimensions, we refer the reader to Ref. [93].

Next, we present a few examples illustrating how the minimal hyperplane
representation for a larger particle number $N^{\prime}$ can be obtained from
the inequalities for a smaller particle number $N$ ($d$ is still fixed).
First, we use the generalized Rado theorem (4.52) together with Eq. (4.53) to
show that a vector
$\boldsymbol{\mu}\in\Sigma_{N^{\prime}}(\boldsymbol{\omega})$ (see Eq. (4.58))
satisfies $\boldsymbol{\mu}\prec\boldsymbol{u}+\delta\boldsymbol{e}_{1}$ for a
convex combination
$\boldsymbol{u}=\sum_{i=1}^{R}p_{i}\boldsymbol{v}_{N}^{(i)}$. Hence, it
follows for the ordered NON vector
$\boldsymbol{\mu}^{\downarrow}\in\Sigma^{\downarrow}_{N}(\boldsymbol{\omega})$
that

$\tilde{\lambda}_{1}\equiv\mu_{1}^{\downarrow}-\delta\leq
u_{1}^{\downarrow}\,.$ (4.63)

Thus, the new vector $\tilde{\boldsymbol{\lambda}}$ with the entries
$\tilde{\lambda}_{1}\equiv\mu_{1}-\delta$ and $\tilde{\lambda}_{i}=\mu_{i}$,
$1<i\leq d$ is normalized to $N$, but its entries are not necessarily ordered
decreasingly anymore. To illustrate that this statement is in agreement with
Eq. (4.62), we start by discussing the two simplest examples $R=1,2$. Recall
that it is sufficient to restrict to the ordered case
$\lambda_{1}\geq\lambda_{2}\geq...\geq\lambda_{d}$ since the polytopes are
permutation invariant. For $r=1$ and the example of $N=3$ and $N^{\prime}=5$
we get the following two conditions for the ordered NON vectors
$\boldsymbol{\lambda}^{\downarrow}$ ($N=3$) and
$\boldsymbol{\mu}^{\downarrow}$ ($N=5$)

$\begin{split}\lambda_{1}^{\downarrow}&\leq 3\,,\quad\mu_{1}^{\downarrow}\leq
5\end{split}$ (4.64)

which are equivalent to the linear hyperplane conditions

$D(\boldsymbol{\lambda}^{\downarrow})\equiv 3-\lambda_{1}^{\downarrow}\geq
0\,,\quad D(\boldsymbol{\mu}^{\downarrow})\equiv 5-\mu_{1}^{\downarrow}\geq
0\,.$ (4.65)

Thus, $D(\boldsymbol{\mu}^{\downarrow})$ can be obtained from
$D(\boldsymbol{\lambda}^{\downarrow})$ by replacing
$\boldsymbol{\lambda}^{\downarrow}\to\boldsymbol{\mu}^{\downarrow}-\delta\boldsymbol{e}_{1}$.
This statement holds obviously also for $r=2$ where there is also only one
possible sequence of configurations. Next, we consider the case of multiple
vectors $\boldsymbol{v}^{(i)},\,\,i=1,...,R$, and denote again by
$\boldsymbol{\lambda}$ the NON vector in the polytope of $N$ particles. The
NON vector $\boldsymbol{\mu}$ corresponds to a larger particle number
$N^{\prime}>N$. Tab. 4.2 compares the linear inequalities obtained from the
majorization condition in the generalized Rado theorem for $r=3$, where we
have two generating vertices, $R=2$, by using the relation of the vectors
$\boldsymbol{v}^{(i)}$ for different total particle numbers in Eq. (4.53).

$N$ |  $N^{\prime}>N$
---|---
$\lambda_{1}^{\downarrow}\leq qv_{1}^{(1)}+(1-q)v_{1}^{(2)}$ |  $\mu_{1}^{\downarrow}\leq qv_{1}^{(1)}+(1-q)v_{1}^{(2)}+\delta$
$\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}\leq q(v_{1}^{(1)}+v_{2}^{(1)})+(1-q)(v_{1}^{(2)}+v_{2}^{(2)})$ |  $\mu_{1}^{\downarrow}+\mu_{2}^{\downarrow}\leq q(v_{1}^{(1)}+v_{2}^{(1)})+(1-q)(v_{1}^{(2)}+v_{2}^{(2)})+\delta$
$\sum_{i=1}^{d-1}\lambda_{i}^{\downarrow}\leq N$ |  $\sum_{i=1}^{d-1}\mu_{i}^{\downarrow}\leq N^{\prime}$

Table 4.2: Comparison of the linear inequalities obtained from the
majorization condition in the generalized Rado theorem for $r=3$ and two
different boson numbers $N$ and $N^{\prime}$.

In accordance with Tab. 4.1, we obtain three inequalities for the minimal
hyperplane representation of $\Sigma^{\downarrow}(\boldsymbol{\omega})$. We
then compare the linear conditions,

$\begin{split}D_{1}(\boldsymbol{\lambda}^{\downarrow})&\equiv
qv_{1}^{(1)}+(1-q)v_{1}^{(2)}-\lambda_{1}^{\downarrow}\geq 0\,,\\\ \
D_{2}(\boldsymbol{\lambda}^{\downarrow})&\equiv
q\left(v_{1}^{(1)}+v_{2}^{(1)}\right)+(1-q)\left(v_{1}^{(2)}+v_{2}^{(2)}\right)-\lambda_{1}^{\downarrow}-\lambda_{2}^{\downarrow}\geq
0\,,\\\ \ D_{r}(\boldsymbol{\lambda}^{\downarrow})&\equiv
N-\sum_{i=1}^{d-1}\lambda_{i}^{\downarrow}\geq 0\,,\end{split}$ (4.66)

corresponding to the smaller particle number $N$ in the left column of Tab.
4.2, to the linear hyperplane conditions

$\begin{split}D_{1}(\boldsymbol{\mu}^{\downarrow})&\equiv
qv_{1}^{(1)}+(1-q)v_{1}^{(2)}+\delta-\mu_{1}^{\downarrow}\geq 0\,,\\\ \
D_{2}(\boldsymbol{\mu}^{\downarrow})&\equiv
q\left(v_{1}^{(1)}+v_{2}^{(1)}\right)+(1-q)\left(v_{1}^{(2)}+v_{2}^{(2)}\right)+\delta-\mu_{1}^{\downarrow}-\mu_{2}^{\downarrow}\geq
0\,,\\\ \ D_{r}(\boldsymbol{\mu}^{\downarrow})&\equiv
N^{\prime}-\sum_{i=1}^{d-1}\mu_{i}^{\downarrow}\geq 0\ \end{split}$ (4.67)

for $N^{\prime}$ bosons in the right column. Then, (4.67) can be obtained from
Eq. (4.66) by replacing $\lambda_{1}^{\downarrow}$ by
$\mu_{1}^{\downarrow}-\delta$. It can be easily shown that the same
observation holds for arbitrary $R$ using $\sum_{i=1}^{R}q_{i}=1$ and the
generalized Rado theorem in Eq. (4.52). To obtain the minimal hyperplane
representation one needs to eliminate the $q_{j}$’s arising from the convex
combinations of generating vertices which is only feasible for $r\leq 3$ as
shown in Sec. 4.5.3. However, the purpose of these examples for small $R$ is
to illustrate Eq. (4.62) and to show how the inequalities for a larger
particle number can be obtained from the inequalities for a smaller particle
number. A mathematically more elaborate treatment is provided in Ref. [93],
where the facet-defining inequalities are derived for arbitrary $N$ and $d$
right from the beginning.

###  4.5 Examples and illustrations

In this section, we discuss the description of the spectral set
$\Sigma(\boldsymbol{\omega})$ proposed in Sec. 4.4 for several small settings.
Moreover, for $r\leq 3$, we illustrate how the redundant inequalities in the
hyperplane representation are identified and how the facet-defining
inequalities can be generalized to arbitrary particle number $N$ and dimension
$d$ of the one-particle Hilbert space.

####  4.5.1 Illustration of the spectral polytopes for two particles and three
non-vanishing weights

The main purpose of our first example is to illustrate the concepts of lineups
and partial ordering introduced in Sec. 4.4 for a small system size which does
not satisfy the stabilization conditions $N\geq r-1$ and $d\geq r$. Therefore,
the number of inequalities does not necessarily match Tab. 4.1.

The excitation spectrum with partial ordering for $N=3$ bosons and
$d=\mathrm{dim}(\mathcal{H}_{1})=2$ is shown in Fig. 4.3. Since both
configurations, $(1,3)$ and $(2,2)$, are lower than $(2,3)$ according to Eq.
(4.43), we obtain the following two possible sequences of length $6$:

$\begin{split}&(1,1)\to(1,2)\to(1,3)\to(2,2)\to(2,3)\to(3,3)\,,\\\
&(1,1)\to(1,2)\to(2,2)\to(1,3)\to(2,3)\to(3,3)\,.\end{split}$ (4.68)

For $r\geq 3$ non-vanishing weights $\omega_{j}$, the two sequences are
distinct, leading to two NON vectors $\boldsymbol{v}^{(1)}$ and
$\boldsymbol{v}^{(2)}$, whereas for $r\leq 2$ there is only one such NON
vector. For the maximal number of non-vanishing weights, $r=6$, they are given
by

$\begin{split}\boldsymbol{v}^{(1)}&=(2\omega_{1}+\omega_{2}+\omega_{3},\omega_{2}+2\omega_{4}+\omega_{5},\omega_{3}+\omega_{5}+2\omega_{6})\,,\\\
\
\boldsymbol{v}^{(2)}&=(2\omega_{1}+\omega_{2}+\omega_{4},\omega_{2}+2\omega_{3}+\omega_{5},\omega_{4}+\omega_{5}+2\omega_{6})\,.\end{split}$
(4.69)

Fig. 4.4 shows the spectral polytopes $\Sigma(\boldsymbol{\omega})$ and
$\Sigma^{\downarrow}(\boldsymbol{\omega})=\Sigma(\boldsymbol{\omega})\cap\Delta$
for $r=6$ non-vanishing weights on the left and $r=2$ on the right hand side.
Using the normalization of $\lambda_{1}+\lambda_{2}+\lambda_{3}$ we can omit
the third dimension and substitute $\lambda_{3}=2-\lambda_{1}-\lambda_{2}$ to
obtain the simplex $\Delta$ defined in Eq. (4.37). The two NON vectors
$\boldsymbol{v}^{(1)}$ and $\boldsymbol{v}^{(2)}$ for $r=6$ are marked by two
red dots in the left panel of Fig. 4.4. According to Eq. (4.39), the extremal
points of the spectral polytope $\Sigma(\omega)$ follow from all possible
permutations of the vector entries of $\boldsymbol{v}^{(1)}$ and
$\boldsymbol{v}^{(2)}$. In contrast, for $r=2$ there is only one NON vector
$\boldsymbol{v}$ whose permutations lead to the full spectral polytope
$\Sigma(\omega)$ since we have just a single corresponding sequence in Fig.
4.3.

Figure 4.3: Excitation spectrum with partial ordering for $d=3$ and $N=2$
bosons.

Figure 4.4: Illustration of the spectral polytopes
$\Sigma(\boldsymbol{\omega})$ and
$\Sigma^{\downarrow}(\boldsymbol{\omega})=\Sigma(\boldsymbol{\omega})\cap\Delta$
for $d=3$ and $N=2$ bosons for two different choices of $\boldsymbol{\omega}$,
and thus $r$.

The facet-defining inequalities for $r=6$ are depicted by the three dashed
lines in the left subfigure of Fig. 4.4. To derive them, we first observe
that, according to the generalization of Rado’s theorem in Eq. (4.52), we need
to determine the set of all NON vectors $\boldsymbol{\lambda}$ for which there
exists at least one convex combination
$\boldsymbol{u}=q\boldsymbol{v}^{(1)}+(1-q)\boldsymbol{v}^{(2)}$ such that
$\boldsymbol{\lambda}\prec\boldsymbol{u}$. This leads to the following set of
constraints on the NONs $\lambda_{j}^{\downarrow}$,

$\begin{split}\lambda_{1}^{\downarrow}&\leq
2\omega_{1}+\omega_{2}+\omega_{4}+q(\omega_{3}-\omega_{4})\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}&\leq
2(\omega_{1}+\omega_{2}+\omega_{3})+\omega_{4}+\omega_{5}-q(\omega_{3}-\omega_{4})\,,\\\
\
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
2\,,\end{split}$ (4.70)

where the last inequality is automatically fulfilled by the normalization.
Since $q\in[0,1]$, the upper bound on the largest NON
$\lambda_{1}^{\downarrow}$ can vary between
$2\omega_{1}+\omega_{2}+\omega_{4}$ and $2\omega_{1}+\omega_{2}+\omega_{3}$
where the latter bound must always be satisfied and the value of
$\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}$ can never exceed
$2(\omega_{1}+\omega_{2}+\omega_{3})+\omega_{4}+\omega_{5}$. Combining these
bounds with the inequality arising from the elimination of the parameter $q$
leads to

$\begin{split}\lambda_{1}^{\downarrow}&\leq
2\omega_{1}+\omega_{2}+\omega_{3}\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}&\leq
2\omega_{1}+2\omega_{2}+2\omega_{3}+\omega_{4}+\omega_{5}\,,\\\ \
2\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}&\leq
4\omega_{1}+3\omega_{2}+2\omega_{3}+2\omega_{4}+\omega_{5}\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
2\,,\end{split}$ (4.71)

where the last inequality is automatically fulfilled by the normalization of
$\boldsymbol{\lambda}$.

For $r=2$ shown in the right panel of Fig. 4.4, we only have one NON vector
$\boldsymbol{v}$ and therefore the first two inequalities (with
$\omega_{3}=\omega_{4}=\omega_{5}=0$) in Eq. (4.71) are sufficient to describe
$\Sigma^{\downarrow}(\boldsymbol{\omega})$. The two dashed lines then
correspond to the saturation of the bounds on $\lambda_{1}^{\downarrow}$ and
$\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}$.

####  4.5.2 Two non-vanishing weights

For the case of $r=1$ and $N$ bosons, we recover ground state RDMFT with only
one configuration given by $(1,1,...,1)$, i.e. all bosons in the ground state,
and $\omega$-minimizer
$\hat{\Gamma}_{\boldsymbol{\omega}_{0}}=\mbox{$|1,1,...,1\rangle$}\\!\mbox{$\langle
1,1,...,1|$}$. The corresponding NON vector reads $\boldsymbol{v}=(N,0,...)$.

The only sequence for $N=3$ particles and two non-vanishing weights is given
by $(1,1,1)\to(1,1,2)$, yielding the corresponding minimizer
$\hat{\Gamma}=\omega_{1}\mbox{$|1,1,1\rangle$}\\!\mbox{$\langle
1,1,1|$}+\omega_{2}\mbox{$|1,1,2\rangle$}\\!\mbox{$\langle 1,1,2|$}$ and

$\boldsymbol{v}=(3\omega_{1}+2\omega_{2},\omega_{2},0,...)\,,$ (4.72)

with $\omega_{1}+\omega_{2}=1$. Then, applying Rado’s theorem, introduced in
Sec. 4.1.2, leads to

$\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern
0.0pt}}}^{1}_{N}(\boldsymbol{w})=\\{\hat{\gamma}\in\mathcal{E}^{1}_{N}\,|\,\mathrm{spec}(\hat{\gamma})\prec\boldsymbol{v}\\}\,.$
(4.73)

The constraints on the natural occupation numbers of
$\hat{\gamma}\in\,\hbox{\vbox{ \hrule height=0.6pt\kern 1.29167pt\hbox{\kern
0.0pt$\mathcal{E}$\kern 0.0pt}}}^{1}_{N}(\boldsymbol{w})$ for arbitrary
$(N,d)$ with $d=\mathrm{dim}(\mathcal{H}_{1})$ follow as

$\begin{split}&\sum_{i=1}^{d-1}\lambda_{i}^{\downarrow}\leq N\,,\\\ \
&\lambda_{1}^{\downarrow}\leq N-1+\omega_{1}\,.\end{split}$ (4.74)

All further possible constraints are automatically fulfilled by the
normalization of $\hat{\gamma}$ to the total particle number. Note that we
obtain only one additional constraint in the last line of Eq. (4.74) compared
to $r=1$ in agreement with Tab. 4.1. Thus, the facet-defining inequalities for
$r=1$ and $r=2$ represent the first two levels of the hierarchy of facet-
defining inequalities introduced in Sec. 4.4.2.

####  4.5.3 Three non-vanishing weights

For $N=3$, we have two sequences

$\begin{split}&(1):\,\,(1,1,1)\to(1,1,2)\to(1,1,3)\,,\\\ \
&(2):\,\,(1,1,1)\to(1,1,2)\to(1,2,2)\,,\end{split}$ (4.75)

and therefore the two NON vectors

$\begin{split}&\boldsymbol{v}^{(1)}=(2+\omega_{1},\omega_{2},\omega_{3},0,...)\,,\\\
\
&\boldsymbol{v}^{(2)}=(1+2\omega_{1}+\omega_{2},2-2\omega_{1}-\omega_{2},0,...)\,.\end{split}$
(4.76)

The majorization condition $\mathrm{spec}(\hat{\gamma})\prec\boldsymbol{u}$
from the generalized Rado theorem (see Eq. (4.52)) for a convex combination
$\boldsymbol{u}=q\boldsymbol{v}^{(1)}+(1-q)\boldsymbol{v}^{(2)}$ requires that

$\begin{split}\lambda_{1}^{\downarrow}&\leq
1+2\omega_{1}+\omega_{2}+q(1-\omega_{1}-\omega_{2})\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}&\leq
3-q(1-\omega_{1}-\omega_{2})\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
3\,.\end{split}$ (4.77)

To derive the facet-defining inequalities of
$\Sigma^{\downarrow}(\boldsymbol{\omega})$, we need to eliminate the parameter
$q$ in the next step. Therefore, we first notice that the upper bound on
$\lambda_{1}^{\downarrow}$ can vary between $1+2\omega_{1}+\omega_{2}$ and
$2+\omega_{1}$ since $q$ is restricted to $q\in[0,1]$. Moreover, the upper
bound on $\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}$ can vary between
$3$ and $2+\omega_{1}+\omega_{2}$. Thus, it must always hold that
$\lambda_{1}^{\downarrow}\leq 2+\omega_{1}$ which requires to adjust the value
of $q$ according to

$q\geq\frac{\lambda_{1}^{\downarrow}-1-2\omega_{1}-\omega_{2}}{1-\omega_{1}-\omega_{2}}\,.$
(4.78)

The above inequality is then used to obtain the facet-defining upper bound on
$\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}$. Thus, we obtain for the
minimal hyperplane representation of
$\Sigma^{\downarrow}(\boldsymbol{\omega})$,

$\begin{split}\lambda_{1}^{\downarrow}&\leq 2+\omega_{1}\,,\\\ \
2\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}&\leq
4+2\omega_{1}+\omega_{2}\,,\\\ \
\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
3\,,\end{split}$ (4.79)

where the last inequality is automatically fulfilled by the normalization of
the 1RDM. Generalizing the inequalities to arbitrary $N$ and $d$ leads to

$\begin{split}\lambda_{1}^{\downarrow}&\leq N-1+\omega_{1}\\\ \
\lambda_{1}^{\downarrow}-\sum_{j=3}^{d}\lambda_{j}^{\downarrow}&\leq
N-2+2\omega_{1}+\omega_{2}\\\ \ \sum_{j=1}^{d-1}\lambda_{j}^{\downarrow}&\leq
N\,.\end{split}$ (4.80)

Thus, the second inequality is the only additional one compared to $r=2$
illustrating again the hierarchy of a generalized exclusion principle for
bosons explained in Sec. 4.4.2. In the next example, we thus only show the
additional inequalities for $r=4$ and for the remaining inequalities refer to
the minimal hyperplane representation for smaller values of $r$.

####  4.5.4 Four non-vanishing weights

For four non-vanishing weights and $N=3$ particles we obtain four sequences,

$\begin{split}&(1,1,1)\to(1,1,2)\to(1,1,3)\to(1,1,4)\,,\\\ \
&(1,1,1)\to(1,1,2)\to(1,1,3)\to(1,2,2)\,,\\\ \
&(1,1,1)\to(1,1,2)\to(1,2,2)\to(1,1,3)\,,\\\ \
&(1,1,1)\to(1,1,2)\to(1,2,2)\to(2,2,2)\,.\end{split}$ (4.81)

The four corresponding generating vertices for arbitrary particle number
$N\geq 3$ are summarized in Tab. 4.3. Thus, different particle numbers $N$
only affect the first entry $v_{1}$ by adding different constants depending on
$N$ to the weights $\omega_{j}$. This is generally true for $r$ non-vanishing
weights, $N\geq r-1$ particles and $d\geq r$.

| $\boldsymbol{v}^{(1)}$ | $\boldsymbol{v}^{(2)}$ | $\boldsymbol{v}^{(3)}$ | $\boldsymbol{v}^{(4)}$
---|---|---|---|---
$v_{1}$ | $\omega_{1}+N-1$ | $2\omega_{1}+\omega_{2}+\omega_{3}+N-2$ | $2\omega_{1}+\omega_{2}+\omega_{4}+N-2$ | $3\omega_{1}+2\omega_{2}+\omega_{3}+N-3$
$v_{2}$ | $\omega_{2}$ | $\omega_{2}+2\omega_{4}$ | $\omega_{2}+2\omega_{3}$ | $\omega_{2}+2\omega_{3}+3\omega_{4}$
$v_{3}$ | $\omega_{3}$ | $\omega_{3}$ | $\omega_{4}$ | $0$
$v_{4}$ | $\omega_{4}$ | $0$ | $0$ | $0$

Table 4.3: Entries of the NON vectors $\boldsymbol{v}^{(i)}$ for $r=4$ and
arbitrary $N\geq r-1$ and $d\geq r$.

Due to the larger number of generating vertices $\boldsymbol{v}^{(i)}$ than
for $r\leq 3$, turning the vertex representation into a halfspace
representation requires a mathematically more complex framework than in the
previous examples. We present the formalism in Ref. [93] and provide only the
inequalities in the following. For the setting $(N,d)=(3,4)$, the minimal
hyperplane representation consists of the three facet-defining inequalities
for $r=3$, as well as

$\begin{split}2\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
4+2\omega_{1}+\omega_{2}+\omega_{3}\,,\\\ \
3\lambda_{1}^{\downarrow}+2\lambda_{2}^{\downarrow}&\leq
6+3\omega_{1}+2\omega_{2}+\omega_{3}\,.\end{split}$ (4.82)

Thus, we obtain in total five constraints on the NON vector
$\boldsymbol{\lambda}^{\downarrow}$ and two new inequalities in Eq. (4.82) in
agreement with Tab. 4.1. For arbitrary $N$ and $d$, one obtains the following
inequalities [93],

$\begin{split}2\lambda_{1}^{\downarrow}+\lambda_{2}^{\downarrow}+\lambda_{3}^{\downarrow}&\leq
2N-2+2\omega_{1}+\omega_{2}+\omega_{3}\,,\\\ \
3\lambda_{1}^{\downarrow}+2\lambda_{2}^{\downarrow}&\leq
3N-3+3\omega_{1}+2\omega_{2}+\omega_{3}\,.\end{split}$ (4.83)

Recall that the minimal hyperplane representation for
$\Sigma(\boldsymbol{\omega})$ follows from the facet-defining inequalities for
$\boldsymbol{\lambda}^{\downarrow}$ by permuting the coefficients in the
inequalities in all possible ways and replacing $\lambda_{i}^{\downarrow}$ by
$\lambda_{i}$.

###  4.6 Application to Bose-Einstein condensation

In this section, we apply the excited state RDMFT to BEC’s using the
Bogoliubov approximated interaction (see Eq. (3.16))

$\hat{W}_{\mathrm{B}}=\frac{N(N-1)W_{\boldsymbol{0}}}{2V}+\frac{1}{2V}\sum_{\boldsymbol{p}\neq
0}W_{\boldsymbol{p}}\left[2\hat{n}_{\boldsymbol{0}}\hat{n}_{\boldsymbol{p}}+\hat{a}_{\boldsymbol{p}}^{\dagger}\hat{a}_{-\boldsymbol{p}}^{\dagger}\hat{a}_{\boldsymbol{0}}^{2}+\big{(}\hat{a}_{\boldsymbol{0}}^{\dagger}\big{)}^{2}\hat{a}_{\boldsymbol{p}}\hat{a}_{-\boldsymbol{p}}\right]\,.$
(4.84)

The goal is then, similar to Ch. 3, to derive a first-level functional for
excited states valid in the regime close to complete condensation. The
interacting ground state of a BEC at $T=0$ has the form
$\mbox{$|\Psi_{0}\rangle$}=\hat{U}\mbox{$|N\rangle$}$ where
$\mbox{$|N\rangle$}=(N!)^{-1/2}(\hat{a}_{\boldsymbol{0}}^{\dagger})^{N}\mbox{$|0\rangle$}$
and $\hat{U}$ is a unitary operator. The GOK variational principle requires
that $E_{0}\leq E_{1}\leq...$, where
$\hat{H}\mbox{$|\Psi_{0}\rangle$}=E_{0}\mbox{$|\Psi_{0}\rangle$}$ is the
eigenvalue equation for the ground state and ground state energy. Thus, the
$\boldsymbol{\omega}$-minimizer for $r=1$ would be
$\hat{\Gamma}_{\boldsymbol{\omega}}=\mbox{$|\Psi_{0}\rangle$}\\!\mbox{$\langle\Psi_{0}|$}$.
Following this consideration, the first excited state is the excitation of one
quasiparticle. Just as Levy’s constrained search for ground state RDMFT, the
constrained search formalism for $\boldsymbol{\omega}$-ensemble RDMFT requires
fixed total particle numbers $N$. Thus, we have to use particle number
conserving quasiparticle operators introduced by Girardeau [69] (see also Sec.
3.4).

####  4.6.1 Derivation of the $\boldsymbol{\omega}$-ensemble functional

The excitation spectrum of a homogeneous BEC in 3D in the thermodynamic limit
is gapless. However, for finite systems there will be a gap between the ground
state and the first excited state because the momenta $\boldsymbol{p}$ are
discrete in that case. We will therefore restrict to finite systems to derive
the excited state functional for Bogoliubov approximated systems in the
following.

The particle number conserving quasiparticles operators, appearing in the
diagonal Bogoliubov Hamiltonian, read

$\hat{c}_{\boldsymbol{q}}^{\dagger}=\hat{U}\hat{a}_{\boldsymbol{q}}^{\dagger}\hat{U}^{\dagger}\hat{\beta}_{0}\,.$
(4.85)

As discussed in Sec. 3.3, RDMFT requires that the unitary operator $\hat{U}$
lives on the $N$-boson Hilbert space and conserves the total particle number.
The Bogoliubov approximated Hamiltonian $\hat{H}_{\mathrm{B}}$ consists of a
constant term yielding the ground state energy plus a second term which is
diagonal in the Bogoliubov quasiparticle operators with dispersion
$\omega_{\boldsymbol{q}}$ [72]. The quasiparticle dispersion is given by

$\omega_{\boldsymbol{q}}=\sqrt{\varepsilon_{\boldsymbol{q}}(\varepsilon_{\boldsymbol{q}}+2nW_{\boldsymbol{q}})}\,,\quad\boldsymbol{q}\neq\boldsymbol{0}\,.$
(4.86)

Thus, the momentum $\boldsymbol{q}$ leading to the lowest
$\omega_{\boldsymbol{p}}$ defines the first excited state $|\Psi_{1}\rangle$
such that
$\hat{H}_{\mathrm{B}}\mbox{$|\Psi_{1}\rangle$}=E_{1}\mbox{$|\Psi_{1}\rangle$}$
with $E_{0}<E_{1}$. Moreover, this reveals a fundamental obstacle in RDMFT
because a possible crossing of energy levels splits the domain of the
$\boldsymbol{\omega}$-ensemble functional into different cells, each
corresponding to a different first excited state and thus a different
functional $\mathcal{F}_{\boldsymbol{\omega}}^{\boldsymbol{q}}(\hat{\gamma})$.
This immediately implies that we cannot obtain a closed expression for the
universal functional but rather have

$\mathcal{F}_{\boldsymbol{\omega}}(\hat{\gamma})=\min_{\boldsymbol{q}}\mathcal{F}_{\boldsymbol{\omega}}^{\boldsymbol{q}}(\hat{\gamma})\,.$
(4.87)

Thus, while deriving the expression for
${F}_{\boldsymbol{\omega}}^{\boldsymbol{q}}(\hat{\gamma})$ in the following,
we implicitly assume that we already made the correct choice of the momentum
$\boldsymbol{q}$ corresponding to the first exited state eigenenergy $E_{1}$.
Since the Hamiltonian is invariant under inversion
$\boldsymbol{p}\to-\boldsymbol{p}$, also the excited states
$\hat{c}_{\boldsymbol{q}}^{\dagger}\mbox{$|\Psi_{0}\rangle$}$ and
$\hat{c}_{-\boldsymbol{q}}^{\dagger}\mbox{$|\Psi_{0}\rangle$}$ correspond to
the same energy. Thus, this holds in particular for any superposition

$\mbox{$|\Psi_{1}\rangle$}_{x}=\left(x\hat{c}_{\boldsymbol{q}}^{\dagger}+\sqrt{1-|x|^{2}}\hat{c}_{-\boldsymbol{q}}^{\dagger}\right)\mbox{$|\Psi_{0}\rangle$}=\hat{U}\left(x\hat{a}_{\boldsymbol{q}}^{\dagger}+\sqrt{1-|x|^{2}}\hat{a}_{-\boldsymbol{q}}^{\dagger}\right)\mbox{$|N-1\rangle$}$
(4.88)

with $x\in\mathbb{C}$ and the quasiparticle operators given by Eq. (4.85). The
eigenvalues of the inversion operator can take the values $\pm 1$ and in the
following we thus use the excited state

$\mbox{$|\Psi_{1}\rangle$}=\frac{1}{\sqrt{2}}\left(\hat{c}_{\boldsymbol{q}}^{\dagger}+\hat{c}_{-\boldsymbol{q}}^{\dagger}\right)\mbox{$|\Psi_{0}\rangle$}=\frac{1}{\sqrt{2}}\hat{U}\left(\hat{a}_{\boldsymbol{q}}^{\dagger}+\hat{a}_{-\boldsymbol{q}}^{\dagger}\right)\mbox{$|N-1\rangle$}$
(4.89)

which corresponds to the symmetric case. The expectation value of
$\hat{H}_{\mathrm{B}}$ in the state $|\Psi_{1}\rangle$ reads

$\begin{split}E_{1}&=\mbox{$\langle\Psi_{1}|$}\hat{H}_{B}\mbox{$|\Psi_{1}\rangle$}=E_{0}+\frac{1}{2}\left(\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}+\sqrt{\varepsilon(-\boldsymbol{q})(\varepsilon(-\boldsymbol{q})+2nW_{-\boldsymbol{q}})}\right)\\\
\
&=E_{0}+\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}\,,\end{split}$
(4.90)

where we used the degeneracy in the last line. To derive the excited state
functional, $\boldsymbol{q}$ and $-\boldsymbol{q}$ can be either treated
together or independently. Both approaches work, but it is important to
remember which convention has been used to be consistent. This will be
particularly important when we derive $\mathcal{F}_{\boldsymbol{\omega}}$
using the Legendre-Fenchel transform because, if the last line in Eq. (4.90)
is used, also the kinetic energy must be treated as
$\langle\hat{t}\rangle=2\sum_{\boldsymbol{p}>\boldsymbol{0}}\varepsilon(\boldsymbol{p})n_{\boldsymbol{p}}$.
For the weighted sum of the two lowest energies $E_{0}$ and $E_{1}$ we
eventually obtain from Eq. (4.90)

$\begin{split}E_{\boldsymbol{\omega}}&=\omega
E_{0}+(1-\omega)\left(E_{0}+\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}\right)\\\
\
&=E_{0}+(1-\omega)\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}\,,\end{split}$
(4.91)

where

$E_{0}=-\sum_{\boldsymbol{p}>\boldsymbol{0}}\left(\varepsilon(\boldsymbol{p})+nW_{\boldsymbol{p}}-\sqrt{\varepsilon(\boldsymbol{p})(\varepsilon(\boldsymbol{p})+2nW_{\boldsymbol{p}})}\right)$
(4.92)

is the well-known result for the ground state energy [63, 72].

In the following, we present three different, but equivalent, approaches to
derive the excited state functional
$\mathcal{F}_{\boldsymbol{\omega}}^{\boldsymbol{q}}(\hat{\gamma})$ which lead
to the same functional, as required. We start by performing the Levy’s
constrained search. Then, we derive
$\mathcal{F}_{\boldsymbol{\omega}}^{\boldsymbol{q}}(\hat{\gamma})$ by
calculating the Legendre-Fenchel transform of (4.91). As a third option to
calculate the functional we use the solution for the variational parameters
$\phi_{\boldsymbol{p}}$ from the ground state problem. We believe that this
extended discussion of several approaches to obtain the excited state
universal functional illustrate different perspectives on how RDMFT works and
thus promotes a more comprehensive understanding of the method.

##### Constrained search

For two non-vanishing weights, $r=2$, the minimization in Eq. (4.13) is
performed over all $N$-boson density operators

$\hat{\Gamma}_{\boldsymbol{\omega}}=\omega\mbox{$|\Psi_{0}\rangle$}\\!\mbox{$\langle\Psi_{0}|$}+(1-\omega)\mbox{$|\Psi_{1}\rangle$}\\!\mbox{$\langle\Psi_{1}|$}\,.$
(4.93)

Then, the average numbers of particles with different momenta $\boldsymbol{p}$
representing the diagonal elements of the 1RDM $\hat{\gamma}$, follow as

$n_{\boldsymbol{p}}=\mbox{Tr}[\hat{\Gamma}_{\boldsymbol{\omega}}\hat{n}_{\boldsymbol{p}}]=\begin{cases}\frac{\phi_{\boldsymbol{p}}^{2}}{1-\phi_{\boldsymbol{p}}^{2}}\,,&\quad\text{if
}\boldsymbol{p}\neq\pm\boldsymbol{q}\,,\\\
\omega\frac{\phi_{\boldsymbol{p}}^{2}}{1-\phi_{\boldsymbol{p}}^{2}}+\frac{1}{2}(1-\omega)\frac{1+3\phi_{\boldsymbol{p}}^{2}}{1-\phi_{\boldsymbol{p}}^{2}}\,,&\quad\text{if
}\boldsymbol{p}=\pm\boldsymbol{q}\,,\end{cases}$ (4.94)

where we used the expression for the unitary operator in Eq. (3.23) as well as
Eq. (3.25), determining the connection between the phases
$\phi_{\boldsymbol{p}}$ and the variational parameters in the unitary operator
$\hat{U}$. Inverting the above expression yields

$\phi_{\boldsymbol{p}}=\begin{cases}\sigma_{\boldsymbol{p}}\sqrt{\frac{n_{\boldsymbol{p}}}{1+n_{\boldsymbol{p}}}}\,,&\quad\text{if
}\boldsymbol{p}\neq\pm\boldsymbol{q}\,,\\\
\sigma_{\boldsymbol{p}}\sqrt{\frac{2n_{\boldsymbol{p}}+\omega-1}{2n_{\boldsymbol{p}}+3-\omega}}\,,&\quad\text{if
}\boldsymbol{p}=\pm\boldsymbol{q}\,.\end{cases}$ (4.95)

Using these results, the excited state universal functional follows from
Levy’s constrained search formalism as

$\begin{split}\mathcal{F}_{\boldsymbol{\omega}}(\boldsymbol{n})&=\min_{\hat{\Gamma}_{\boldsymbol{\omega}}\mapsto\hat{\gamma}}\mbox{Tr}\left[\hat{\Gamma}_{\boldsymbol{\omega}}\hat{W}_{\mathrm{B}}\right]\\\
\
&=\mathcal{F}_{0}(\boldsymbol{n})+\sum_{\pm\boldsymbol{q}}n|W_{\boldsymbol{q}}|\left(\sqrt{n_{\boldsymbol{q}}(1+n_{\boldsymbol{q}})}-\frac{1}{2}\sqrt{4n_{\boldsymbol{q}}(n_{\boldsymbol{q}}+1)+\omega(4-\omega)-3}\right)\\\
\
&=\mathcal{F}_{0}(\boldsymbol{n})+n|W_{\boldsymbol{q}}|\left(2\sqrt{n_{\boldsymbol{q}}(1+n_{\boldsymbol{q}})}-\sqrt{4n_{\boldsymbol{q}}(n_{\boldsymbol{q}}+1)+\omega(4-\omega)-3}\right)\,.\end{split}$
(4.96)

where $\mathcal{F}_{0}$ denotes the ground state functional

$\mathcal{F}_{0}(\boldsymbol{n})=n\sum_{\boldsymbol{p}\neq\boldsymbol{0}}W_{\boldsymbol{p}}\left(n_{\boldsymbol{p}}-\mathrm{sgn}(W_{\boldsymbol{p}})\sqrt{n_{\boldsymbol{p}}(n_{\boldsymbol{p}}+1)}\right)\,.$
(4.97)

Consequently, Eq. (4.96) reduces to $\mathcal{F}_{0}$ for $\omega=1$, as
required.

As a consistency check with the weighted sum of energies given in Eq. (4.91),
we minimize the energy functional

$\mathcal{E}_{\boldsymbol{\omega}}(\boldsymbol{n})=\sum_{\boldsymbol{p}\neq\boldsymbol{0}}\varepsilon(\boldsymbol{p})n_{\boldsymbol{p}}+\mathcal{F}_{\boldsymbol{\omega}}(\boldsymbol{n})$
(4.98)

with respect to the occupation number vector $\boldsymbol{n}$. The solution of
$\varepsilon(\boldsymbol{p})=-\partial\mathcal{F}_{\boldsymbol{\omega}}/\partial
n_{\boldsymbol{p}}$ for every $\boldsymbol{p}$ is denoted by
$\accentset{\rule{3.5pt}{0.5pt}}{n}_{\boldsymbol{p}}$, leading to

$\accentset{\rule{3.5pt}{0.5pt}}{n}_{\boldsymbol{p}}=\frac{1}{2}\left(\frac{\varepsilon(\boldsymbol{p})+nW_{\boldsymbol{p}}}{\sqrt{\varepsilon(\boldsymbol{p})(\varepsilon(\boldsymbol{p})+2nW_{\boldsymbol{p}})}}-1\right)$
(4.99)

and

$\accentset{\rule{3.5pt}{0.5pt}}{n}_{\boldsymbol{q}}=\frac{1}{2}\left(\frac{(\varepsilon(\boldsymbol{q})+nW_{\boldsymbol{q}})(2-\omega)}{\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}}-1\right)\,.$
(4.100)

Note that we treat $\boldsymbol{q}$ and $-\boldsymbol{q}$ separately at this
point which, as already explained above, does not make any difference, if done
consistently. After collecting the contributions from $\boldsymbol{q}$ and
$-\boldsymbol{q}$ we eventually obtain

$\begin{split}E_{\boldsymbol{\omega}}&=-\sum_{\boldsymbol{p}>\boldsymbol{0}}\left(nW_{\boldsymbol{p}}+\varepsilon(\boldsymbol{p})-\sqrt{\varepsilon(\boldsymbol{p})(\varepsilon(\boldsymbol{p})+2nW_{\boldsymbol{p}})}\right)+(1-\omega)\sqrt{\varepsilon(\boldsymbol{q})(\varepsilon(\boldsymbol{q})+2nW_{\boldsymbol{q}})}\,,\end{split}$
(4.101)

which is in agreement with Eq. (4.91) and thus verifies that
$\mathcal{F}_{\boldsymbol{\omega}}$ obtained in Eq. (4.96) is indeed correct.

##### Legendre-Fenchel transform

Since the energy $E_{\boldsymbol{\omega}}$ and and the universal functional
$\mathcal{F}_{\boldsymbol{\omega}}$ are related through the Legendre-Fenchel
transform by (see Sec. 2.5)