#  Orbital liquid in the $e_{g}$ orbital Hubbard model in $d=\infty$
dimensions

Louis Felix Feiner  Institute for Theoretical Physics, Utrecht University, PO
Box 80.000, 3508 TA Utrecht, Netherlands Andrzej M. Oleś  Max Planck Institute
for Solid State Research, Heisenbergstrasse 1, D-70569 Stuttgart, Germany
Institute of Theoretical Physics, Jagiellonian University, Prof. Stanisława
Łojasiewicza 11, PL-30348 Kraków, Poland

(31 October, 2022)

###### Abstract

We demonstrate that the three-dimensional $e_{g}$ orbital Hubbard model can be
generalized to arbitrary dimension $d$, and that the form of the result is
determined uniquely by the requirements that (i) the two-fold degeneracy of
the $e_{g}$ orbital be retained, and (ii) the cubic lattice be turned into a
hypercubic lattice. While the local Coulomb interaction $U$ is invariant for
each basis of orthogonal orbitals, the form of the kinetic energy depends on
the orbital basis and takes the most symmetric form for the so-called complex-
orbital basis. Characteristically, with respect to this basis, the model has
two hopping channels, one that is orbital-flavor conserving, and a second one
that is orbital-flavor non-conserving. We show that the noninteracting
electronic structure consists of two nondegenerate bands of plane-wave real-
orbital single-particle states for which the orbital depends on the wave
vector. Due to the latter feature each band is unpolarized at any filling, and
has a non-Gaussian density of states at $d=\infty$. The orbital liquid state
is obtained by filling these two bands up to the same Fermi energy. We
investigate the $e_{g}$ orbital Hubbard model in the limit $d\to\infty$,
treating the on-site Coulomb interaction $U$ within the Gutzwiller
approximation, thus determining the correlation energy of the orbital liquid
and the (disordered) para-orbital states. In perfect analogy with the case of
the spin Hubbard model, the Gutzwiller approximation is demonstrated to be
exact at $d=\infty$ for the orbital Hubbard model, because of the collapse of
electron correlations to a single site. At half-filling (one electron per site
on average, $n=1$) one finds a Brinkman-Rice type “metal-insulator” transition
in the orbital liquid, which is analogous to the transition for a paramagnetic
state in the spin model, but occurs at stronger Hubbard interaction $U$ due to
the enhanced kinetic energy provided by the non-conserving hopping channel. We
show that the orbital liquid is the ground state everywhere in the $(n,U)$
phase diagram except close to half-filling at sufficiently large $U$, where
ferro-orbital order with real orbitals occupied is favored. The latter feature
is shown to be specific for $d=\infty$, being of mathematical nature due to
the exponential tails in the density of states.

##  I Introduction

Our understanding of electron correlations in solids stems to a considerable
extent from the Hubbard model [1]. This describes electrons which propagate on
a lattice and interact by an on-site Coulomb interaction $U$ when two
electrons with opposite spins occupy the same site. At half-filling (one
electron per site) it provides a picture of a metal-insulator transition due
to electron localization and of magnetism caused by superexchange according to
ideas going back to Anderson [2]. Away from half-filling it represents the
competition between the kinetic energy of the electrons and the localizing
effect of the Coulomb interaction. An important development regarding the
Hubbard model was the discovery by Metzner and Vollhardt that in the limit of
infinite dimension ($d\\!=\\!\infty$) only the on-site correlations survive
[3, 4]. They thus proved that in this limit the Gutzwiller approximation to
the variational ground state wave function [5, 6, 7] becomes exact [8], and
that a similar property applies [9] to the resonating valence-bond (RVB) state
[10]. It was then shown that the ground state phase diagram of the doped
Hubbard model at $d=\infty$ includes ferromagnetic (FM) and antiferromagnetic
(AF) order stabilized at large Coulomb interaction [11].

The $e_{g}$ orbital Hubbard model [12, 13, 14, 15, 16], discussed in the
present paper, describes spinless fermions which propagate on a lattice of
atoms with twofold degenerate orbitals and which interact by an on-site
Coulomb interaction $U$ when two electrons occupy both orbitals on a site.
Importantly, the hopping strength $t$ is dependent on the orientation of the
orbitals with respect to the hopping direction. Obviously, with twofold
orbital degeneracy replacing twofold spin degeneracy, the orbital model
closely resembles the spin model, but there is a crucial difference: in the
spin model the local symmetry is SU(2) but in the orbital model it is the
threefold rotation group $C_{3}$ and moreover this is linked to the hopping.
As an immediate consequence the model contains, in addition to the familiar
orbital-conserving hopping channel (analogous to the spin-conserving hopping
channel in the (spin) Hubbard model), also an orbital-flipping hopping
channel. This leads in the 3D case to the existence of an orbital liquid (OL)
state [12, 13], in analogy to a spin liquid [17, 18], which is disordered in a
way different from paramagnetism. Our aim is to analyze the limit $d\to\infty$
along the lines set out by Metzner and Vollhardt and explore whether this
gives further insight into the properties of the $e_{g}$ orbital Hubbard model
and the OL state in particular.

The $e_{g}$ orbital Hubbard model comes from the field of spin-orbital
physics, which has developed over the last two decades into a very active and
challenging part of solid state physics which unifies frustrated magnetism and
the phenomena in strongly correlated electron systems. It arose from the
pioneering ideas of Kugel and Khomskii [19] who recognized that in transition
metal oxides with partly filled degenerate $3d$ orbitals, both the spin and
the orbital degrees of freedom are quantum variables and have to be treated on
equal footing. At large local Coulomb interaction $U$ most of these compounds
are Mott (or charge-transfer) insulators due to electron localization, and
superexchange arises from coupled spins and orbitals together [20, 21, 22].
More recently similar phenomena have been found in ultracold fermion systems
in optical lattices [23].

In the degenerate-orbital case [24], quantum fluctuations in the insulating
state are enhanced with respect to the nondegenerate case [25, 26, 27]. In
some cases these joint fluctuations could even trigger a new state of quantum
matter—a spin-orbital liquid in model [28] or real systems [29, 30, 31]. In
most cases, however, one finds spin-orbital order due to spin interactions of
a Heisenberg form with SU(2) symmetry, coupled to anisotropic orbital
superexchange [32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47,
48, 49, 50, 51], in agreement with the Goodenough-Kanamori rules [52]. These
interactions give in general entangled ground states [53, 54, 55, 56, 57, 58].
We remark that the physical properties of such systems are very rich and
depend on whether the orbital degrees of freedom are of $e_{g}$, $t_{2g}$, or
$p$ type. The main difference is the spatial character of the orbitals which
makes their interactions directional.

An issue of high interest is how such systems, characterized by the presence
of orbital and spin degrees of freedom, behave under doping, i.e., how their
properties compare with the more familiar behavior of doped spin systems, such
as the high-$T_{c}$ cuprates [59, 60, 61]. In the latter systems remarkable
progress has been achieved by employing dynamical mean-field theory (DMFT)
[62], which gave very valuable insights into the transport properties [63] and
spectral functions [64]. This theory, which is currently used widely for
electronic structure calculations of strongly correlated materials where the
one-electron description breaks down [65], is in fact based on the above-
mentioned discovery by Metzner and Vollhardt regarding the Hubbard model [3,
4].

A realistic description of transition metal oxides with degenerate $3d$
orbitals requires the treatment of spin-orbital physics in at least a
degenerate Hubbard model, involving both orbital and spin degrees of freedom
[22]. However, in particular situations the problem becomes somewhat
simplified. In spin-orbital systems with FM order and weak spin-orbit
coupling, the spins may be ignored because the spin flavor is conserved in the
hopping processes, and the spins disentangle from the orbitals [66, 67] and do
not contribute to the dynamics. The degenerate Hubbard model appropriate for
the spin-orbital system [68] then reduces to an orbital Hubbard model with
direction-dependent hopping.

When the active orbitals are of $t_{2g}$ type this orbital Hubbard model is
rather similar to the standard spin Hubbard model because the orbital flavor
is then conserved in the hopping processes (like for spin) but the hopping is
two-dimensional (2D) [29, 69]. This gives almost localized orbital polarons
for single holes [70, 71] or orbital stripes due to self-organization of the
doped Mott insulator [72, 73].

However, in doped FM systems with active $e_{g}$ orbitals [74, 75] the orbital
flavor is not conserved, in contrast to the $t_{2g}$ case above, and the
kinetic energy consists of all possible hopping processes including those
which change the orbital flavor. Obviously, this represents a qualitative
difference with the spin Hubbard model. Moreover, it suggests that disordered
phases are favored more in doped $e_{g}$ than in doped $t_{2g}$ systems.

In doped manganites the FM metallic phase occurs due to the kinetic energy
gain by the double exchange mechanism [76], because antiferromagnetic (AF)
bonds hinder electron hopping while hopping is unrenormalized along FM bonds,
and the $e_{g}$ electrons reorient the $t_{2g}$ spins coupled to them by
Hund’s rule exchange. Double exchange for strongly correlated $e_{g}$
electrons is indeed responsible for the spectacular properties of doped
perovskite manganites [77, 78, 79], including colossal magnetoresistance in
the FM metallic phase [80]. At low doping the $A$-type AF order in
La1-xSrxMnO3 favors $e_{g}$ electron transport within 2D $(a,b)$ planes, and
gradually changes the spin order to FM along the $c$ axis. Increasing doping
generates orbital polarons [81]. Already a single polaron triggers [82] an
insulator-to-metal transition that occurs here to a FM metallic phase.
Eventually, the orbitals decouple from the spins and a disordered orbital
liquid (OL) phase arises in a 3D system [12, 13, 14, 15]. This OL phase may
explain the cubic dispersion relation of magnons in the FM phase of some
manganites which indicates that magnetic interactions are isotropic [83, 84]
unless they couple strongly to orbital excitations [57]. Specifically this
compound, La1-xSrxMnO3, provides a prime example of the physical situation
represented by the orbital model studied in this paper, with the spins
integrated out and the dynamics determined solely by the orbitals.

So the $e_{g}$ orbital Hubbard model, while still directly relevant for a
class of doped FM insulators, is the simplest nontrivial multi-orbital model
representing interacting lattice fermions, and as such presents an opportunity
for studying fundamental issues in orbital dynamics. Moreover, because of
their mathematical similarity one can readily compare its behavior with that
of the (spin) Hubbard model, at the same time exploiting the wealth of results
available for the latter one.

As stated, our overall aim is to investigate whether the limit of infinite
dimension elucidates the characteristic features of the $e_{g}$ orbital
Hubbard model, like it does for the spin Hubbard model [11]. In particular we
want to establish which features found in the 3D case [12, 13, 14, 15] are
characteristic for the model as such, i.e., continue to hold for arbitrary
dimension up to $d=\infty$, and which are specific for $d=3$ per se or are
caused by the fact that at finite dimension the Gutzwiller approach is an
approximation.

The purpose of this paper is therefore (i) to derive the orbital $e_{g}$ model
at arbitrary dimension $d$ from the requirements that the twofold degeneracy
of the $e_{g}$ orbitals be retained and that the cubic lattice symmetry be
turned into a hypercubic symmetry; (ii) to highlight the fundamental
difference with the spin Hubbard model which is that the orbital flavor is not
conserved and the kinetic energy is therefore potentially larger, which should
hinder orbital order and favor states with disordered orbitals; (iii) to
demonstrate that the Gutzwiller approximation is again exact in the limit
$d\to\infty$ like it is for the (spin) Hubbard model, and next to investigate
(iv) in what way the OL state avoids orbital polarization, and the role played
therein by symmetry; (v) whether at half-filling electron localization induces
a metal-insulator transition in the OL state and, if so, at which critical
value of the on-site interaction $U$ (vi) what the nature is of the ground
state (long-range order of ferro type or alternating type, or rather
disordered of para-orbital type or orbital-liquid type), dependent on the
strength of the interaction $U$ and the particle density $n$.

The paper is organized as follows. In Sec. II.1 we introduce the orbital
Hubbard model for spin-polarized $e_{g}$ electrons for a cubic lattice
($d=3$), and discuss its symmetry properties. We show that the cubic symmetry
of the hopping may be better appreciated when a particular basis consisting of
two orbitals with complex coefficients is used. This serves to introduce the
orbital model at general dimension $d$ and at $d=\infty$ in Sec. II.2. Next we
use the translational invariance to present the noninteracting Hamiltonian in
momentum space in Sec. II.3. Band dispersion and densities of states for
complex and real single-particle states are presented in Secs. III.1 and
III.2, respectively. In Sec. III.3 we address the single-particle eigenstates
and show that they exhibit a generic splitting into a lower and an upper band.
The orbital model is compared at dimension $d=\infty$ with the spin model in
Sec. III.4 and we argue from the densities of states that its general feature
is an enhanced kinetic energy. In Sec. III.5 we analyze the kinetic energy in
dependence of filling and compare it for the simplest symmetry-broken states,
disordered (para-orbital) states, and the OL state built from the single-
particle eigenstates.

The electron correlations are described using the Gutzwiller approximation
within a generalization of the Metzner approach [8, 9], analyzed in Sec. IV.1,
and we evaluate the renormalized propagator at $d=\infty$ using the collapse
of diagrams to a single site in Sec. IV.2. The general formalism for uniform
and two-sublattice states is introduced in Sec. IV.3. In Sec. V we treat
specific trial variational states: (i) ordered states in Sec. V.1, (ii) the OL
phase in Sec. V.2, concluding the demonstration that the Metzner approach
remains valid in the orbital case. In the next Section VI we compare the OL in
the orbital Hubbard model with the paramagnet in the spin Hubbard model, and
show that a Brinkman-Rice transition takes place in the OL just like in the
paramagnet but for a considerably larger value of $U$. The results of the
numerical analysis and the phase diagram of the orbital Hubbard model at
$d=\infty$ are presented in Sec. VII. At the end we focus on some general
aspects of the orbital physics and suggest possible extensions of the present
study (Sec. VIII.1). The paper is concluded in Sec. VIII.2 by pointing out the
differences between the spin and orbital Hubbard model in infinite dimension.
In the Appendix we present a proof that the orbital liquid phase is
unpolarized.

##  II The $e_{g}$ orbital Hubbard model

###  II.1 The model at dimension $d=3$

The usual choice of basis for $e_{g}$ orbitals in a 3D cubic lattice is to
take

$\textstyle{|z\rangle\equiv\frac{1}{\sqrt{6}}(2z^{2}-x^{2}-y^{2}),\hskip
11.38092pt|\bar{z}\rangle\equiv\frac{1}{\sqrt{2}}(x^{2}-y^{2}),}$ (1)

called real orbitals. However, because this basis is the natural one only for
the bonds parallel to the $c$ axis but not for those in the $(a,b)$ plane, the
kinetic energy takes then a rather nonsymmetric form, having a very different
appearance depending on the bond direction. It is thus preferred to use
instead the basis of complex orbitals at each site [85],

$\textstyle{|+\rangle=\frac{1}{\sqrt{2}}\big{(}|z\rangle+i|\bar{z}\rangle\big{)},\hskip
14.22636pt|-\rangle=\frac{1}{\sqrt{2}}\big{(}|z\rangle-i|\bar{z}\rangle\big{)},}$
(2)

corresponding to “up” and “down” pseudospin flavors, with the local pseudospin
operators defined as

$\displaystyle{\hat{T}}_{i}^{+}$ $\displaystyle=$
$\displaystyle{\hat{c}}_{i+}^{\dagger}{\hat{c}}_{i-},\hskip
34.14322pt{\hat{T}}_{i}^{-}={\hat{c}}_{i-}^{\dagger}{\hat{c}}_{i+},$
$\displaystyle{\hat{T}}_{i}^{z}$ $\displaystyle=$
$\displaystyle\textstyle{\frac{1}{2}}({\hat{c}}_{i+}^{\dagger}{\hat{c}}_{i+}-{\hat{c}}_{i-}^{\dagger}{\hat{c}}_{i-})=\textstyle{\frac{1}{2}}({\hat{n}}_{i+}-{\hat{n}}_{i-}).$
(3)

For later reference it is convenient to introduce also electron creation
operators ${\hat{c}}_{i}^{\dagger}(\psi_{i},\theta_{i})$ which create $e_{g}$
electrons in orbital coherent states
$|\Omega_{i}\rangle\equiv|\Omega_{i}(\psi_{i},\theta_{i})\rangle$ at site $i$,
defined as

$\displaystyle|\Omega_{i}\rangle\equiv\operatorname{e}^{-i\theta_{i}/2}\cos\left(\frac{\psi_{i}}{2}\right)|i+\rangle+\operatorname{e}^{+i\theta_{i}/2}\sin\left(\frac{\psi_{i}}{2}\right)|i-\rangle,$
(4)

in analogy with the well-known spin coherent states [86]. The local pseudospin
operator in this state is

$\langle\Omega_{i}|{\bm{\hat{T}}}_{i}|\Omega_{i}\rangle=\textstyle{\frac{1}{2}}(\sin\psi_{i}\cos\theta_{i},\sin\psi_{i}\sin\theta_{i},\cos\psi_{i}),$
(5)

i.e., it behaves like a vector [87].

Figure 1:  Schematic representation of the real orbital wave functions in the
’equatorial’ plane at dimensions (a) $d=3$, (c) $d=7$, and (b) of the
pseudospin vector at any dimension. The solid lines in (a) correspond to
directional $(3z^{2}-r^{2})$-like orbitals (the azimuthal angle $\theta_{i}/2$
in Eq. (6) equals one of the $\chi_{\alpha}/2$), while the dashed lines
correspond to the $(x^{2}-y^{2})$-like orbitals orthogonal to them
($\theta_{i}/2$ equalling $\chi_{\alpha}/2+\pi/2$). In (c) the solid lines
correspond to the $|x_{n}\rangle_{d}$ orbitals (the azimuthal angle
$\theta_{i}/2$ equals one of the $\chi_{n}/2$), while only the dashed line
corresponding to the central orthogonal orbital $|\bar{x}_{0}\rangle_{d}$ is
shown.

The parameter space of the coherent orbital (4) is a sphere, with the “poles”
($\psi_{i}=0$ and $\psi_{i}=\pi$) corresponding to the complex orbitals
$|i+\rangle$ and $|i-\rangle$, and the “equator” ($\psi_{i}=\pi/2$)
corresponding to (linear combinations of) the real orbitals,

$\left|\Omega_{i}\left(\frac{\pi}{2},\theta_{i}\right)\right\rangle=\cos\left(\frac{\theta_{i}}{2}\right)\left|iz\right\rangle+\sin\left(\frac{\theta_{i}}{2}\right)\left|i\bar{z}\right\rangle.$
(6)

The three real bases, $\\{|ix\rangle,|i\bar{x}\rangle\\}$,
$\\{|iy\rangle,|i\bar{y}\rangle\\}$, and $\\{|iz\rangle,|i\bar{z}\rangle\\}$,
associated with the cubic axes $a$, $b$, and $c$, respectively, are obtained
by setting the in-plane angle $\theta_{i}/2$ equal to $\chi_{a}/2=+2\pi/3$,
$\chi_{b}/2=-2\pi/3$, and $\chi_{c}/2=0$, for $|ix\rangle,|iy\rangle$, and
$|iz\rangle,$ and to $\chi_{\alpha}+\pi/2$ ($\alpha=a,b,c$) for
$|i\bar{x}\rangle,|i\bar{y}\rangle,|i\bar{z}\rangle$ , as illustrated in Fig.
1(a).

In the complex-orbital representation the orbital Hubbard model for $e_{g}$
electrons in $d=3$ takes the form [13]

$\displaystyle{\cal H}_{3}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\;t\sum_{\langle
ij\rangle\parallel{\hat{\alpha}}}\Big{\\{}\Big{(}{\hat{c}}_{i+}^{\dagger}{\hat{c}}_{j+}+{\hat{c}}_{i-}^{\dagger}{\hat{c}}_{j-}\Big{)}$
(7) $\displaystyle\hskip
19.91684pt+\gamma\Big{(}\operatorname{e}^{-i\chi_{\alpha}}{\hat{c}}_{i+}^{\dagger}{\hat{c}}_{j-}+\operatorname{e}^{+i\chi_{\alpha}}{\hat{c}}_{i-}^{\dagger}{\hat{c}}_{j+}\Big{)}\Big{\\}}$
$\displaystyle+$ $\displaystyle U\sum_{i}{\hat{n}}_{i+}{\hat{n}}_{i-},$

where the parameter $\gamma$ actually takes the value $\gamma=1$. This
parameter is introduced here as a device by which one may interpolate between
the standard 2-flavor Hubbard model with hopping $t/2$ at $\gamma=0$ and the
orbital $e_{g}$ model at $\gamma=1$. The appearance of the phase factors
$\operatorname{e}^{\pm i\chi_{\alpha}}$ is characteristic for the orbital
problem — they occur because the orbitals have an actual shape in real space
so that each hopping process depends on the bond direction and on the orbitals
between which hopping occurs. Moreover, Eq. (7) exhibits the crucial feature
that (except at $\gamma=0$) orbital flavor is not conserved, or, equivalently,
that pseudospin is not conserved, compare Eq. (3).

The model (7) consists of the kinetic energy $H_{\rm kin}\propto t$ and the
interorbital Coulomb interaction $H_{\rm int}\propto U$. The interaction is
invariant under any local basis transformation to a pair of orthogonal
orbitals, i.e., it gives an energy $U$ when a double occupancy occurs in any
representation, i.e., either when both real orbitals are simultaneously
occupied or when both complex orbitals are occupied,

$H_{\rm
int}=U\sum_{i}{\hat{n}}_{iz}{\hat{n}}_{i\bar{z}}=U\sum_{i}{\hat{n}}_{i+}{\hat{n}}_{i-}.$
(8)

We emphasize that this simple Hubbard-like form of the Coulomb interaction
operator obtained here corresponds to the high-spin ($S=1$) charge
excitations, $d^{1}_{i}d^{1}_{j}\rightarrow d^{2}_{i}d^{0}_{j}$, — these are
the only excited states in a FM system. The parameter $U$ stands then for the
onsite interorbital Coulomb repulsion, and is in fact lower than the
intraorbital Coulomb element (Kanamori parameter) $U_{0}$ by $3J_{H}$ due to
Hund’s exchange in the triplet state, i.e., $U=U_{0}-3J_{H}$ [33, 16].

Importantly, as noted in [13], the equivalence of the cubic axes shows up in
the kinetic energy term in Eq. (7) in a transparent way, viz. as a formal
invariance under threefold rotations of the phase angles
$\chi_{\alpha}\mapsto\chi_{\alpha}\\!-\\!4\pi/3$ in conjunction with a phase
shift of the electron operators
${\hat{c}}_{i\pm}^{\dagger}\mapsto\operatorname{e}^{\mp
2\pi/3}{\hat{c}}_{i\pm}^{\dagger}$. So the relevant symmetry group is $C_{3}$,
with the $e_{g}$ doublets $\\{|iz\rangle,|i\bar{z}\rangle\\}$ transforming as
real $E$ representations, or equivalently the pairs
$\\{|i+\rangle,|i-\rangle\\}$ transforming as conjugate complex $A_{1}$
representations.

###  II.2 Generalization to large dimension $d$

This form, Eq. (7), therefore lends itself to a natural generalization of the
$e_{g}$ orbital Hubbard model with two orbital flavors
$\\{|+\rangle,|-\rangle\\}$ to large dimension $d$ as follows,

$\displaystyle{\cal H}_{d}$ $\displaystyle=$
$\displaystyle-\;\frac{1}{\sqrt{2d}}\,\frac{t}{2}\sum_{\langle
ij\rangle\parallel\hat{x}_{n}}\Big{\\{}\Big{(}{\hat{c}}_{i+}^{\dagger}{\hat{c}}_{j+}+{\hat{c}}_{i-}^{\dagger}{\hat{c}}_{j-}\Big{)}$
(9) $\displaystyle\hskip
42.67912pt+\gamma\Big{(}\operatorname{e}^{-i\chi_{n}}{\hat{c}}_{i+}^{\dagger}{\hat{c}}_{j-}+\operatorname{e}^{+i\chi_{n}}{\hat{c}}_{i-}^{\dagger}{\hat{c}}_{j+}\Big{)}\Big{\\}}$
$\displaystyle+U\sum_{i}{\hat{n}}_{i+}{\hat{n}}_{i-},$

with $\chi_{n}=4\pi n/d$, where $n=-m,-m+1,\cdots,m-1,m$, and $d\equiv 2m+1$
(i.e., taken odd for convenience but $d$ could also be even with only slightly
different derivations below), and we note here that
$\sum_{n}\operatorname{e}^{i\chi_{n}}=0$. So instead of three cubic axes we
now have $d$ hypercubic axes $\\{\hat{x}_{n}\\}$, labelled by the index $n$.
Importantly, the number of orbital flavors remains two (as for spin $S=1/2$).

The real basis $\\{|x_{n}\rangle_{d},|\bar{x}_{n}\rangle_{d}\\}$ associated
with axis $\hat{x}_{n}$ is given by Eq. (6) with $\theta_{i}/2$ set equal to
$\chi_{n}/2$ and to $\chi_{n}/2+\pi/2$, respectively. This is illustrated in
Fig. 1(c) for $d=7$ in comparison with the case for $d=3$ shown in Fig. 1(a).
Figure 1(b) shows the associated behavior of the pseudospin vector, valid both
for $d=3$ and for general $d$.

The hopping parameter $t$ is scaled in the standard fashion [3, 8, 11] by
$\sqrt{2d}$ in order that the average kinetic energy remains finite as $d$
becomes large. Below we will use the abbreviation $\tilde{t}=t/\sqrt{2d}$
whenever convenient. Again, the parameter $\gamma$ takes the value $\gamma=1$
for the orbital problem, in which the orbital flavor (pseudospin) is not
conserved, while $\gamma=0$ for the corresponding spin problem with hopping
$t/2$, in which the (spin) flavor is conserved. The interorbital Coulomb
interaction $\propto U$ is the same for any dimension $d$ as it stands for the
on-site Coulomb repulsion when both orbital (or spin) flavors are occupied.

Obviously, the equivalence of the hypercubic axes manifests itself in the
Hamiltonian ${\cal H}_{d}$ as a $d$-fold rotational symmetry, i.e., by being
invariant under phase shifts
${\hat{c}}_{i\pm}^{\dagger}\mapsto\operatorname{e}^{\pm
i2\pi/d}{\hat{c}}_{i\pm}^{\dagger}$ of the electron operators in conjunction
with shifts of the angles $\chi_{n}\mapsto\chi_{n}\\!-\\!4\pi/d$, in perfect
analogy with the 3-fold rotational symmetry of ${\cal H}_{3}$. So ${\cal
H}_{d}$ is invariant under $C_{d}$, with the $e_{g}$ doublets
$\\{|ix_{n}\rangle_{d},|i\bar{x}_{n}\rangle_{d}\\}$ still transforming as real
$E$ representations. In the limit $d\to\infty$ the symmetry group $C_{d}$
turns into $C_{\infty}\equiv\mathbb{Z}$, which for practical purposes
approaches continuous rotational symmetry U(1).

The generalization to dimension $d=2m+1$ for the explicit form of the
normalized real basis orbitals associated with the central axis $\hat{x}_{0}$
in terms of the coordinates $x_{n}$ is obtained straightforwardly,

$\displaystyle|x_{0}\rangle_{d}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{\sqrt{2d}}\left(\sum_{n=-m}^{m}\cos(\chi_{n}/2)\,x_{n}^{2}\right),$
(10) $\displaystyle|\bar{x}_{0}\rangle_{d}$ $\displaystyle\equiv$
$\displaystyle\frac{2}{\sqrt{2d}}\left(\sum_{n=-m}^{m}\sin(\chi_{n}/2)\,x_{n}^{2}\right),$
(11)

of which Eq. (1) is seen to be the special case for $d=3$. We observe that for
general dimension $d$, the “directional” orbital $|x_{0}\rangle_{d}$ has a
less pronounced directional shape than $|z\rangle=|x_{0}\rangle_{3}$ in the 3D
case, because of the many finite contributions from all the coordinates
$x_{n}$ with $n\neq 0$, although the square of the central coordinate $x_{0}$
still has the largest coefficient. The orthogonal orbital
$|\bar{x}_{0}\rangle_{d}$ is perpendicular to the $\hat{x}_{0}$ axis as in the
3D case (the coefficient of ${x_{0}}^{2}$ in $|\bar{x}_{0}\rangle_{d}$ is zero
as $\chi_{0}=0$). When $d$ gets large the shape of $|\bar{x}_{0}\rangle_{d}$
resembles that of $|x_{0}\rangle_{d}$ but with permuted coordinates (because
Eq. (11) is obtained by replacing $\chi_{n}/2$ in Eq. (10) by
$\chi_{n}/2\\!-\\!\pi/2$). Clearly, at each site all real-orbital basis sets
$\\{|x_{n}\rangle_{d},|\bar{x}_{n}\rangle_{d}\\}$ are equivalent. We have made
the arbitrary choice to use the set
$\\{|x_{0}\rangle_{d},|\bar{x}_{0}\rangle_{d}\\}$ as the reference basis for
real orbitals, and for convenience to drop the subscript $d$.

###  II.3 Momentum space representation

Using translational invariance, we introduce creation and annihilation
operators in momentum space,

$\hat{c}_{{\bf
k}\lambda}^{\dagger}=\frac{1}{\sqrt{N}}\sum_{j}\mathrm{e}^{i{\bf k}{\rm
R}_{j}}\,\hat{c}_{j\lambda}^{\dagger},$ (12)

where $\lambda$ labels the orbital flavors, e.g. $\\{x_{0},\bar{x}_{0}\\}$ or
$\\{+,-\\}$, and $N$ is the number of sites. The free-electron part of ${\cal
H}_{d}$, Eq. (9), describing the kinetic energy ${\cal H}_{\rm kin}$, may be
represented upon Fourier transformation in the form

${\cal H}_{\rm kin}=-\tilde{t}\,\sum_{\bf k}\left({\hat{c}}_{{\bf
k}+}^{\dagger}\;{\hat{c}}_{{\bf
k}-}^{\dagger}\right)\\!\left(\begin{array}[]{cc}A_{\bf k}&\gamma G_{\bf
k}^{\ast}\\\\[5.69046pt] \gamma G_{\bf k}&A_{\bf
k}\end{array}\right)\\!\left(\begin{array}[]{c}{\hat{c}}_{{\bf
k}+}\\\\[5.69046pt] {\hat{c}}_{{\bf k}-}\end{array}\right),$ (13)

where the dispersion is given by orbital-conserving and orbital-non-conserving
terms,

$\displaystyle A_{\bf k}$ $\displaystyle=$ $\displaystyle\sum_{n}\cos
k_{n}\;,$ (14) $\displaystyle G_{\bf k}$ $\displaystyle=$
$\displaystyle\sum_{n}e^{+i\chi_{n}}\cos k_{n}\equiv|G_{\bf
k}|\,\operatorname{e}^{+i\phi_{\bf k}}\;.$ (15)

It will turn out convenient to introduce the following definitions,

$\displaystyle B_{\bf k}\\!$ $\displaystyle\equiv$ $\displaystyle|G_{\bf
k}|\\!=\\!\\!\left(\sum_{n,n^{\prime}}\cos(\chi_{n}\\!-\chi_{n^{\prime}})\cos
k_{n}\\!\cos k_{n^{\prime}}\\!\\!\right)^{\\!\\!\\!{1/2}}\\!\\!,$ (16)
$\displaystyle C_{\bf k}\\!$ $\displaystyle\equiv$ $\displaystyle\Re(G_{\bf
k})=B_{\bf k}\cos\phi_{\bf k}=\sum_{n}{\cos\chi_{n}\cos k_{n}}\,,$ (17)
$\displaystyle D_{\bf k}\\!$ $\displaystyle\equiv$ $\displaystyle\Im(G_{\bf
k})=B_{\bf k}\sin\phi_{\bf k}=\sum_{n}{\sin\chi_{n}\cos k_{n}}\,.$ (18)

One may now absorb the phase factors from the offdiagonal elements of the
matrix in Eq. (13) into the operators, by defining what we will call “phased”
complex-orbital single-particle operators,

$\hat{\tilde{c}}_{{\bf k}\pm}^{\dagger}=\mathrm{e}^{\mp i{\phi_{\bf
k}}/2}\,\hat{c}_{{\bf k}\pm}^{\dagger},$ (19)

and rewrite ${\cal H}_{\rm kin}$ as

${\cal H}_{\rm kin}=-\tilde{t}\,\sum_{\bf k}\left({\hat{\tilde{c}}}_{{\bf
k}+}^{\dagger}\;{\hat{\tilde{c}}}_{{\bf
k}-}^{\dagger}\right)\\!\left(\begin{array}[]{cc}A_{\bf k}&\gamma B_{\bf
k}\\\\[5.69046pt] \gamma B_{\bf k}&A_{\bf
k}\end{array}\right)\\!\left(\begin{array}[]{c}{\hat{\tilde{c}}}_{{\bf
k}+}\\\\[5.69046pt] {\hat{\tilde{c}}}_{{\bf k}-}\end{array}\right).$ (20)

If desired, the kinetic energy may also be rewritten from Eq. (13) in terms of
the real-orbital operators,

$\displaystyle\hskip 2.84544pt{\cal H}_{\rm kin}$ $\displaystyle=$
$\displaystyle-\tilde{t}\,\sum_{\bf k}\left[A_{\bf k}\,{\hat{\cal A}}_{\bf
k}-\gamma\left(C_{\bf k}\,{\hat{\cal C}}_{\bf k}+D_{\bf k}\,{\hat{\cal
D}}_{\bf k}\,\right)\right]$ $\displaystyle=$
$\displaystyle-\tilde{t}\,\sum_{\bf k}\left[A_{\bf k}\,{\hat{\cal A}}_{\bf
k}-\gamma B_{\bf k}\left(\cos\phi_{\bf k}\,{\hat{\cal C}}_{\bf
k}+\sin\phi_{\bf k}\,{\hat{\cal D}}_{\bf k}\,\right)\right],$

where

$\displaystyle{\hat{\cal A}}_{\bf k}$ $\displaystyle=$
$\displaystyle\left({\hat{c}}_{{\bf k}x_{0}}^{\dagger}{\hat{c}}_{{\bf
k}x_{0}}+{\hat{c}}_{{\bf k}{\bar{x}}_{0}}^{\dagger}{\hat{c}}_{{\bf
k}{\bar{x}}_{0}}\right),$ (22) $\displaystyle{\hat{\cal C}}_{\bf k}$
$\displaystyle=$ $\displaystyle\left({\hat{c}}_{{\bf
k}{\bar{x}}_{0}}^{\dagger}{\hat{c}}_{{\bf k}{\bar{x}}_{0}}-{\hat{c}}_{{\bf
k}x_{0}}^{\dagger}{\hat{c}}_{{\bf k}x_{0}}\right),$ (23)
$\displaystyle{\hat{\cal D}}_{\bf k}$ $\displaystyle=$
$\displaystyle\left({\hat{c}}_{{\bf k}x_{0}}^{\dagger}{\hat{c}}_{{\bf
k}{\bar{x}}_{0}}+{\hat{c}}_{{\bf k}{\bar{x}}_{0}}^{\dagger}{\hat{c}}_{{\bf
k}x_{0}}\right).$ (24)

The invariance with respect to the symmetry group $C_{d}$ of this real-orbital
representation of ${\cal H}_{\rm kin}$ is warranted by $A_{\bf k}$ and
${\hat{\cal A}}_{\bf k}$ both transforming as $A_{1}$, and by $\\{C_{\bf
k},D_{\bf k}\\}$ and $\\{{\hat{\cal C}}_{\bf k},{\hat{\cal D}}_{\bf k}\\}$
both transforming as real $E$ representations.

We note that whereas the kinetic-energy Hamiltonian (13), (20), or (II.3), is
block-diagonal in ${\bf k}$-space, the orbital flavors $\\{+,-\\}$ or
$\\{x_{0},{\bar{x}}_{0}\\}$ are still coupled. The further analysis of the
kinetic energy depends on whether some form of symmetry breaking takes place,
or whether we deal with an unpolarized orbital state.

##  III Noninteracting electrons

Below, in Sec. V, we will calculate the electron correlations due to the on-
site Hubbard repulsion $U$ using the type of wave function introduced by
Gutzwiller [6, 7],

$|\Psi\rangle=g^{\hat{D}}\big{|}\Phi_{0}\big{\rangle}=\prod_{i}\left(1-\left(1-g\right)\hat{D}_{i}\right)\big{|}\Phi_{0}\big{\rangle}\,,$
(25)

where $|\Phi_{0}\rangle$ is an arbitrary uncorrelated trial state,

$\hat{D}=\sum_{i}\hat{D}_{i},\quad\quad\quad\hat{D}_{i}=\hat{n}_{i+}\hat{n}_{i-},$
(26)

and $g$ is a variational parameter, $0\leq g\leq 1$. Frequently
$|\Phi_{0}\rangle$ is chosen to be a Gutzwiller wave function proper, i.e., a
Fermi-sea state of noninteracting particles,

$|\Phi_{0}\rangle=\prod_{{\bf k}\in{\cal K}_{\lambda}}\hat{c}_{{\bf
k}\lambda}^{\dagger}\prod_{{\bf q}\in{\cal K}_{\bar{\lambda}}}\hat{c}_{{\bf
q}\bar{\lambda}}^{\dagger}\,|0\rangle,$ (27)

where ${\cal K}_{\lambda}$ and ${\cal K}_{\bar{\lambda}}$ are the parts of the
$d$-dimensional Brillouin zone occupied by particles of opposite orbital
flavors $\lambda$ and $\bar{\lambda}$, determined from the condition that the
respective single-particle kinetic energies be smaller than Fermi energies
$E_{\mathrm{F},\lambda}$ and $E_{\mathrm{F},\bar{\lambda}}$.

For the orbital model, because it has only $C_{\infty}$ symmetry, even for a
Fermi-sea-like $|\Phi_{0}\rangle$ the free-particle kinetic energy depends on
the particular set of orbitals that are occupied, as in the 3D case [13]. This
is fundamentally different from the situation for the spin Hubbard model,
where because of the SU(2) symmetry of the spins one has just two (spin)
flavors quantized along an arbitrary direction in spin space, and both kinetic
energy and renormalization depend only on the filling of the spin subbands.
Therefore, we analyze in the present Section the dispersions and the densities
of states (DOSs) for the single-particle states used later on in building
variational states $|\Phi_{0}\rangle$.

###  III.1 Complex-orbital single-particle states

Figure 2:  DOSs of the various single-particle states: (a) for the complex-
orbital single-particle states, $\rho_{\rm c}(E)$, $\rho_{\rm c,ao}(E)$, and
$\rho_{\tilde{\rm c},{\rm ao}}(E)$; (b) for the real-orbital single-particle
states, $\rho_{\rm r}(E)$, $\rho_{\rm r,ao}(E)$, and $2\rho_{\rm r}(E)$, and
for the single-particle eigenstates, $\rho_{\ell}(E)$, $\rho_{u}(E)$, and
$\rho_{\ell+u}(E)$.

Particularly simple is the dispersion for the single-particle states with
complex orbitals as defined by Eq. (12)) for $\lambda\\!=\\!+$ or
$\lambda\\!=\\!-$,

$\varepsilon_{+}({\bf k})=\varepsilon_{-}({\bf k})\equiv\varepsilon_{\rm
c}({\bf k})=-\tilde{t}\,A_{\bf k}\,.$ (28)

We also consider single-particle states adapted to G-type partitioning of the
lattice into A and B sublattices. With the sublattice operators given by

$\displaystyle\hat{c}_{\mathrm{A}{\bf k}\lambda}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\hat{c}_{{\bf
k}\lambda}^{\dagger}+\hat{c}_{{\bf k}+{\bf Q}\lambda}^{\dagger}\right),$ (29)
$\displaystyle\hat{c}_{\mathrm{B}{\bf k}\lambda}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\hat{c}_{{\bf
k}\lambda}^{\dagger}-\hat{c}_{{\bf k}+{\bf Q}\lambda}^{\dagger}\right),$ (30)

where ${\bf Q}=(\pi,\pi,\cdots,\pi)$ and $\lambda$ is equal to either $+$ or
$-$, we introduce single-particle operators with alternating orbitals of
opposite flavor $\lambda$ (on A sites) and $\bar{\lambda}$ (on B sites)
according to

$\displaystyle\hat{c}_{{\bf k}\lambda\bar{\lambda},\ell}^{\dagger}$
$\displaystyle=$ $\displaystyle\frac{1}{\sqrt{2}}\left(\hat{c}_{\mathrm{A}{\bf
k}\lambda}^{\dagger}-\hat{c}_{\mathrm{B}{\bf
k}\bar{\lambda}}^{\dagger}\right),$ (31) $\displaystyle\hat{c}_{{\bf
k}\lambda\bar{\lambda},u}^{\dagger}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\left(\hat{c}_{\mathrm{A}{\bf
k}\lambda}^{\dagger}+\hat{c}_{\mathrm{B}{\bf
k}\bar{\lambda}}^{\dagger}\right).$ (32)

From Eq. (13) and making use of the relations $A_{\bf k+Q}=-A_{\bf k}$,
$B_{\bf k+Q}=B_{\bf k}$, and $\phi_{\bf k+Q}=\phi_{\bf k}+\pi$, one readily
finds that they correspond to a lower ($\ell$) and an upper ($u$) band in the
(reduced) Brillouin zone, both doubly degenerate, with dispersions given by

$\displaystyle\varepsilon_{+-,\ell}({\bf k})$ $\displaystyle=$
$\displaystyle-\varepsilon_{+-,u}({\bf k})=\varepsilon_{-+,\ell}({\bf
k})=-\varepsilon_{-+,u}({\bf k})$ (33) $\displaystyle\equiv$
$\displaystyle\varepsilon_{\rm c,ao}({\bf k})=-\gamma\tilde{t}\,C_{\bf
k}=-\gamma\tilde{t}\,B_{\bf k}\cos\phi_{\bf k},$

where the first (second) subscript in the first line gives the occupation of
the A (B) sublattice. Following the same procedure we also obtain operators
for single-particle states with phased complex orbitals and alternating
flavor, and find their dispersions,

$\displaystyle{\tilde{\varepsilon}}_{+-,\ell}({\bf k})$ $\displaystyle=$
$\displaystyle-{\tilde{\varepsilon}}_{+-,u}({\bf
k})={\tilde{\varepsilon}}_{-+,\ell}({\bf
k})=-{\tilde{\varepsilon}}_{-+,u}({\bf k})$ (34) $\displaystyle\equiv$
$\displaystyle{\tilde{\varepsilon}}_{\rm c,ao}({\bf
k})=-\gamma\tilde{t}\,B_{\bf k}\,.$

It is noteworthy that the dispersions (33) and (34) of the two alternating
single-particle states are proportional to $\gamma$, i.e. they originate
entirely from the orbital-flavor non-conserving part of the kinetic energy.

We now proceed to the DOSs associated with the above single-particle states.
In the limit $d\to\infty$ an exact analytical expression can be derived for
all single-particle states considered, also for those discussed below. This
was well known to be the case for the rather simple dispersion, identical to
Eq. (28), of the spin Hubbard model [3], but it holds generally because, from
a mathematical point of view [88], DOSs are probability distributions and in
the limit $d\to\infty$ are governed by the law of large numbers. It follows
that every DOS is either of Gaussian type,

$\rho(E)=\frac{1}{\sqrt{2\pi}}\;\frac{1}{wt}\;\mathrm{e}^{-\frac{1}{2}\left(\frac{E}{wt}\right)^{2}}\,,$
(35)

fully characterized by its (dimensionless) width $w$, or is derivable fairly
simply from Gaussians. For the above dispersions Eqs. (28) and (33) the DOSs
are Gaussians with widths $w_{\rm c}=1/2$ and $w_{\rm c,ao}=\gamma/2\sqrt{2}$,
respectively, while for the dispersion (34) one finds that the DOS is a
symmetrized Rayleigh distribution [88],

$\rho_{\tilde{\rm c},{\rm ao}}(E)=\frac{2}{\gamma t}\;\left|\frac{2E}{\gamma
t}\right|\;\mathrm{e}^{-\left(\frac{2E}{\gamma t}\right)^{2}}\,.$ (36)

These three complex-orbital DOSs are shown for the full orbital case, i.e.
$\gamma=1$, in Fig. 2(a) by the blue, green, and purple line, respectively.

###  III.2 Real-orbital single-particle states

For the single-particle states with real orbitals as defined by Eq. (12) with
$\lambda\\!=\\!x_{0}$ or $\lambda\\!=\\!\bar{x}_{0}$, and for the single-
particle states describing alternating real flavors $x_{0}$ and $\bar{x}_{0}$
on the sublattices defined in the same way as above by Eqs. (29)–(32), one
finds the dispersions

$\displaystyle\varepsilon_{x_{0}}({\bf k})$ $\displaystyle=$
$\displaystyle-\tilde{t}\,\left(A_{\bf k}+\gamma C_{\bf
k}\right)\equiv\varepsilon_{r}({\bf k})\,,$ (37)
$\displaystyle\varepsilon_{\bar{x}_{0}}({\bf k})$ $\displaystyle=$
$\displaystyle-\tilde{t}\,\left(A_{\bf k}-\gamma C_{\bf
k}\right)\equiv\varepsilon_{r}({\bf-k})\,,$ (38)
$\displaystyle\varepsilon_{{x_{0}}{\bar{x}_{0}},\ell}({\bf k})$
$\displaystyle=$ $\displaystyle-\varepsilon_{{x_{0}}{\bar{x}_{0}},u}({\bf
k})=\varepsilon_{{\bar{x}_{0}}{x_{0}},\ell}({\bf
k})=-\varepsilon_{{\bar{x}_{0}}{x_{0}},u}({\bf k})$ (39) $\displaystyle\equiv$
$\displaystyle\varepsilon_{\rm r,ao}({\bf k})=\gamma\tilde{t}D_{\bf
k}=\gamma\tilde{t}\,B_{\bf k}\sin\phi_{\bf k}.$

One observes that, as already indicated at the very end of Sec. II.2, the sets
of real-orbital single-particle states $\\{|{\bf k}x_{0}\rangle\\}$ and
$\\{|{\bf k}\bar{x}_{0}\rangle\\}$ are equivalent in the limit $d\to\infty$,
in the sense that they have the same energies be it at opposite ${\bf k}$. As
above for the complex-orbital single-particle states the alternating real-
orbital single-particle states correspond to a lower ($\ell$) and an upper
($u$) band in the (reduced) Brillouin zone, both doubly degenerate, and their
dispersions are proportional to $\gamma$, stemming entirely from the orbital-
flavor non-conserving part of the kinetic energy.

The corresponding DOSs are again Gaussians, with widths $w_{\rm r}=\sigma/2$
and $w_{\rm r,ao}=\gamma/2\sqrt{2}$, respectively, where
$\sigma=\sqrt{1+{\gamma^{2}}/2}$). They are shown, for $\gamma=1$, in Fig.
2(b) by the blue and green line respectively. One may note that $w_{\rm
r,ao}=w_{\rm c,ao}$, implying that the DOS of the alternating real-orbital
single-particle states is identical to that of its complex-orbital counterpart
shown in the upper panel.

###  III.3 Single-particle eigenstates

Figure 3:  Evolution of the DOS for the single-particle eigenstates $\ell$ and
$u$ of the $e_{g}$ orbital model at $U=0$ with increasing dimension: $d=3$
(blue), $d=5$ (green), $d=11$ (red), and $d=\infty$ (black). The lines at
finite $d$ have been calculated by ${\bf k}$-space sampling, those at
$d=\infty$ from Eq. (III.3).

The creation operators for the single-particle eigen-states of the kinetic
energy are seen from Eq. (20) to be simply proportional to the sum and
difference of the phased single-particle operators,

$\displaystyle{\hat{d}}_{{\bf
k},{\ell}}^{\dagger}=\frac{1}{\sqrt{2}}\left({\hat{\tilde{c}}}_{{\bf
k}+}^{\dagger}+{\hat{\tilde{c}}}_{{\bf k}-}^{\dagger}\right),$ (40)
$\displaystyle{\hat{d}}_{{\bf
k},{u}}^{\dagger}=\frac{-i}{\sqrt{2}}\left({\hat{\tilde{c}}}_{{\bf
k}+}^{\dagger}-{\hat{\tilde{c}}}_{{\bf k}-}^{\dagger}\right),$ (41)

and the corresponding energies of the lower (${\ell}$) and upper ($u$)
subbands are

$\varepsilon_{{\ell},u}({\bf k})=-t\left(A_{\bf k}\pm\gamma B_{\bf
k}\right)\;,$ (42)

where the $+$ ($-$) sign corresponds to the lower (upper) subband. We
calculate the DOSs for the two subbands by evaluating the convolution of the
DOSs corresponding to $A_{\bf k}$ and to the two branches of $B_{\bf k}$.
Strictly this requires $A_{\bf k}$ and $B_{\bf k}$ to be uncorrelated [88],
which is obviously not the case because both depend on ${\bf k}$, but we
assume that this condition becomes irrelevant in the limit $d\to\infty$, and
check the validity of this assumption afterwards (see below). We then obtain
for the subband DOSs

$\displaystyle\rho_{\ell,u}(E)$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2\pi}}\frac{2}{\sigma^{2}t}\;\mathrm{e}^{-\frac{1}{2}(\frac{2E}{t})^{2}}$
$\displaystyle\times$ $\displaystyle\left\\{1+\sqrt{\pi}\frac{\gamma E}{\sigma
t}\,\mathrm{e}^{\left(\frac{\gamma E}{\sigma
t}\right)^{2}}\left[\mathrm{erf}{\left(\frac{\gamma E}{\sigma t}\right)}\mp
1\right]\right\\}\,,$

where the upper (lower) sign corresponds to the $\ell$ ($u$) subband, and for
their sum, which we will need later,

$\displaystyle\rho_{\ell+u}(E)$ $\displaystyle=$
$\displaystyle\frac{2}{\sqrt{2\pi}}\frac{2}{\sigma^{2}t}\;\mathrm{e}^{-\frac{1}{2}(\frac{2E}{t})^{2}}$
(44) $\displaystyle\times$ $\displaystyle\left\\{1+\sqrt{\pi}\,\frac{\gamma
E}{\sigma t}\,\mathrm{e}^{(\frac{\gamma E}{\sigma
t})^{2}}\mathrm{erf}{\left(\frac{\gamma E}{\sigma t}\right)}\right\\}\,.$

These DOSs are also shown, for $\gamma=1$, in Fig. 2(b) by the dash-dotted red
and orange lines and by the dotted red line, respectively.

We emphasize that the subbands $\\{\varepsilon_{\ell}({\bf
k}),\varepsilon_{u}({\bf k})\\}$ are non-degenerate except for the ${\bf
k}$-points for which $B_{\bf k}=0$. This feature found already at dimension
$d=3$ is generic and persists at increasing dimension up to $d=\infty$,
contrary to the suggestion made in Ref. [89]. Figures 3 and 4 demonstrate this
explicitly. In Fig. 3 the DOSs of the two subbands are seen to become smoother
with increasing dimension but to remain distinct although they partly overlap.
In Fig. 4 the evolution of the subband DOSs with $\gamma$ is illustrated: in
the spin Hubbard model, i.e., at $\gamma=0$ (not shown) they necessarily
coincide, compare Eq. (42), while they are seen to get separated with their
maxima moving away from each other with increasing $\gamma$ and their maxima
finally reaching a maximum distance of $\simeq 0.84$ at $\gamma=1$ (not
shown), i.e., in the full orbital model. This figure also demonstrates the
validity of our assumption made in deriving Eq. (III.3): the colored lines,
calculated by ${\bf k}$-space sampling of the dispersion Eq. (42) at $10^{6}$
${\bf k}$-values, are indistinguishable from the black lines, calculated from
the analytical $d=\infty$ expression Eq. (III.3).

Figure 4:  Densities of states $\rho_{\ell,u}(E)$ for the single-particle
eigenstates of the $e_{g}$ orbital model at $U=0$, for three values of the
parameter $\gamma$ ($\gamma=0.25$, $0.50$, and $0.75$, with increasing
distance between the maxima), as given by the analytical expression (III.3)
for $d\to\infty$ (thin black lines), and as calculated numerically at
dimension $d=21$, where green lines are used for $|{\bf k},\ell\rangle$,
orange lines are used for $|{\bf k},u\rangle$).

We point out that the single-particle eigenstates given by Eqs. (40)-(41) are
in fact real and could be considered as the “phased” real-orbital single-
particle operators which are the counterpart to the phased complex-orbital
single-particle operators defined by Eq. (19) above,

$\displaystyle{\hat{d}}_{{\bf k},{\ell}}^{\dagger}$ $\displaystyle=$
$\displaystyle\>\>\>\,\cos(\phi_{\bf k}/2)\,{\hat{c}}_{{\bf
k}x_{0}}^{\dagger}+\sin(\phi_{\bf k}/2)\,{\hat{c}}_{{\bf
k}\bar{x}_{0}}^{\dagger},$ (45) $\displaystyle{\hat{d}}_{{\bf
k},{u}}^{\dagger}$ $\displaystyle=$ $\displaystyle-\sin(\phi_{\bf
k}/2)\,{\hat{c}}_{{\bf k}x_{0}}^{\dagger}+\cos(\phi_{\bf
k}/2)\,{\hat{c}}_{{\bf k}\bar{x}_{0}}^{\dagger},$ (46)

as one readily verifies from Eqs. (19) and (1), or alternatively by
recognizing that Eqs. (40) and (41) have the form of coherent orbital states
like Eq. (4) with $\psi=0$ and thus are given by Eq. (6). So they are
represented in the “equatorial” plane, compare Fig. 1(c), by an angle $\theta$
equal to the phase angle $\phi_{\bf k}$, which is entirely determined by the
wavevector and therefore is in general not equal to any of the $\chi_{n}$. An
important consequence is that each of these single-particle states carries a
finite orbital polarization, i.e., contributes at each lattice site to the
value of the pseudospin in the $xy$-plane, see Fig. 1(b),

$\displaystyle\langle{\bf k},{\ell}|\hat{T}_{i}^{x}|{\bf k},{\ell}\rangle$
$\displaystyle=$ $\displaystyle-\langle{\bf k},u|\hat{T}_{i}^{x}|{\bf
k},u\rangle=\frac{1}{2N}\cos\phi_{\bf k},$ (47) $\displaystyle\langle{\bf
k},{\ell}|\hat{T}_{i}^{y}|{\bf k},{\ell}\rangle$ $\displaystyle=$
$\displaystyle-\langle{\bf k},u|\hat{T}_{i}^{y}|{\bf
k},u\rangle=\frac{1}{2N}\sin\phi_{\bf k},$ (48) $\displaystyle\langle{\bf
k},{\ell}|\hat{T}_{i}^{z}|{\bf k},{\ell}\rangle$ $\displaystyle=$
$\displaystyle\langle{\bf k},u|\hat{T}_{i}^{z}|{\bf k},u\rangle=0.$ (49)

It is noteworthy that while “$\ell$” and “$u$” may properly be called
(sub)bands, the polarization direction of the states in each band varies with
${\bf k}$, in contrast to the familiar feature of bands of particles carrying
spin where the spin direction is the same for all states in a band.

Figure 5:  Total DOS $\rho_{\ell+u}(E)$ (44) (solid lines) and the electron
density per orbital/spin flavor $n(E)$ (dashed lines) in the unpolarized
state, as obtained for the orbital $e_{g}$ Hubbard model ($\gamma=1$, red
lines) and the spin Hubbard model ($\gamma=0$, blue lines) at $d=\infty$.

###  III.4 Comparison with the spin Hubbard model

The expression (44) derived for the combined DOS of the two subbands
$\rho_{\ell+u}(E)$ gives us the opportunity to compare the distribution of the
single-particle eigenstates in the $e_{g}$ orbital Hubbard with the one in the
spin Hubbard model, see Fig. 5. In the orbital model the total DOS is broader
and the total electron density $n(E)$ grows slower with increasing Fermi
energy $E$. The plot already indicates that the kinetic energy in the orbital
model at a particular density $n$, obtained by filling the two subbands up to
the same Fermi energy, is lower than the kinetic energy obtained at the same
density in the spin model. This difference is a manifestation of the flavor
non-conserving hopping which opens an extra channel for the kinetic energy and
thus enhances it in the orbital model.

###  III.5 Kinetic energies

Before taking electron correlations into account, we examine the kinetic
energy of the various trial states $|\Phi_{0}\rangle$ built from the single-
particle states discussed above, and identify the most favorable states in the
absence of Coulomb interactions. We consider on the one hand uniformly
polarized ordered states, obtained by filling a single band (or, in the cases
of alternating order, subsequently the lower and upper band corresponding to
the same sublattice occupation), and on the other hand unpolarized “para-
orbital” states, obtained by filling two degenerate bands up to the same Fermi
energy, compare Eq. (27). For each such trial state the particle density $n$
and the kinetic energy $E_{\rm kin}$ as a function of the Fermi energy $E_{\rm
F}$ are given by

$\displaystyle n(E_{\rm F})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{E_{\rm F}}\rho(E)\,\mathrm{d}E\,,$ (50)
$\displaystyle E_{\mathrm{kin}}(E_{\rm F})$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{E_{\rm F}}E\,\rho(E)\,\mathrm{d}E\,,$ (51)

where $\rho(E)$ is the DOS of the corresponding band(s) of single-particle
states. Only the density range $0\leq n\leq 1$ will be considered since $n$
and $2-n$ are equivalent because of particle-hole symmetry.

Figure 6:  Kinetic energies $E_{\rm kin}$ in the limit $d\to\infty$ of the
orbital ordered and the para-orbital states of the $e_{g}$ Hubbard model (9)
(i.e. $\gamma=1$) at $U=0$ for increasing electron density $n$, compared with
that of the orbital liquid state: (a) phases with complex orbitals: ferro-
orbital FOc (solid blue line), alternating-orbital AOc (solid green line),
alternating phased orbital AO$\rm\tilde{c}$ (solid cyan line), para-orbital
POc (dashed blue line), para alternating orbital PAOc (dashed green line),
para alternating phased orbital PAO$\rm\tilde{c}$ (dashed cyan line), and
orbital liquid OL (dotted red line); (b) phases with real orbitals: ferro-
orbital FOr (solid blue line), alternating-orbital AOr (solid green line),
para-orbital POr (dashed blue line), para alternating orbital PAOr (dashed
green line), and orbital liquid OL [dotted red line, the same as in (a)].

We begin with the ordered orbital phases built from complex-orbital single-
particle states. The simplest one is the complex ferro-orbital state
$\big{|}\Phi_{0}^{\rm FOc}\big{\rangle}$ obtained by filling the $|{\bf
k}+\rangle$-band (or the $|{\bf k}-\rangle$-band). The complex alternating-
orbital state $\big{|}\Phi_{0}^{\rm AOc}\big{\rangle}$ is, compare Eqs. (31)
and (32), obtained by filling consecutively the $|{\bf k},+-,\ell\rangle$-band
and the $|{\bf k^{\prime}},+-,u\rangle$-band, and the phased complex
alternating-orbital state $\big{|}\Phi_{0}^{\rm AO\tilde{\rm c}}\big{\rangle}$
is built similarly from the $\hat{\tilde{c}}_{{\bf k}+-,\ell}^{\dagger}$ and
$\hat{\tilde{c}}_{{\bf k^{\prime}}+-,u}^{\dagger}$ single-particle states.

Moving on to the unpolarized states we have the complex para-orbital state
POc, obtained by filling simultaneously the $|{\bf k}+\rangle$\- and the
$|{\bf k}-\rangle$-band, the complex alternating para-orbital state PAOc, by
filling simultaneously the $|{\bf k},+-,\ell\rangle$-band and the $|{\bf
k^{\prime}},-+,\ell\rangle$-band, compare Eq. (31), i.e., from the lower band
only, and finally the similarly defined complex phased alternating para-
orbital state PAO$\tilde{\rm c}$. Figure 6(a) shows the kinetic energies of
these six states plotted versus the particle density.

Next we consider the trial states built from real-orbital single-particle
states, again beginning with the ordered phases. Analogous to the complex
trial states above we have the real ferro-orbital state $\big{|}\Phi_{0}^{\rm
FOr}\big{\rangle}$ obtained by filling the $|{\bf k}x_{0}\rangle$-band (or
alternatively the $|{\bf k}{\bar{x}}_{0}\rangle$-band) and the real
alternating-orbital state $\big{|}\Phi_{0}^{\rm AOr}\big{\rangle}$ obtained by
filling consecutively the $|{\bf k},x_{0}{\bar{x}}_{0},\ell\rangle$-band and
the $|{\bf k^{\prime}},x_{0}{\bar{x}}_{0},u\rangle$-band, compare Eqs. (31)
and (32). Here the unpolarized states are the real para-orbital state POr,
with equally filled $|{\bf k}x_{0}\rangle$-band and $|{\bf
k}{\bar{x}}_{0}\rangle$-band, and the real alternating para-orbital state
PAOr, with the $|{\bf k},x_{0}{\bar{x}}_{0},\ell\rangle$-band and the $|{\bf
k^{\prime}},{\bar{x}}_{0}x_{o},\ell\rangle$-band filled equally. The kinetic
energies of these four states plotted versus particle density are shown in
Fig. 6(b).

The $n$-dependence of the various kinetic energies can be understood as
follows. For all ordered states for which the single-particle DOS is Gaussian,
which includes the FOc, AOc, FOr, and AOr states, the particle density and
kinetic energy take the form

$\displaystyle n(E_{\rm F})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1+\operatorname{erf}\\!\left(\frac{E_{\rm
F}}{wt\sqrt{2}}\right)\right]\,,$ (52) $\displaystyle E_{\mathrm{kin}}(E_{\rm
F})$ $\displaystyle=$
$\displaystyle-\frac{wt}{\sqrt{2\pi}}\;\mathrm{exp}\\!\left(\\!-\frac{1}{2}\\!\left(\frac{E_{\rm
F}}{wt}\right)^{\\!\\!2}\right)\,,$ (53)

where $\operatorname{erf}$ is the error function. Formally, inversion of Eq.
(52) gives

$E_{\rm F}(n)=wt\sqrt{2}\>\mathrm{inverf}(2n-1)\,,$ (54)

where $\mathrm{inverf}$ is the inverse error function, and inserting this into
Eq. (53) one obtains

$E_{\mathrm{kin}}(n)=-\frac{wt}{\sqrt{2\pi}}\;\mathrm{exp}\\!\left(-\left(\mathrm{inverf}(2n-1)\right)^{2}\right)\,.$
(55)

This demonstrates that the shape of $E_{\mathrm{kin}}(n)$ is the same for the
states mentioned above, as is evident in Fig. 6, and that the only difference
is in the width prefactor $w$, which also determines the value of the minimum,
$-wt/\sqrt{2\pi}$, at $n=1/2$. The only exception is the AO$\tilde{\rm c}$
state for which the Rayleigh distribution type DOS Eq. (36) leads to

$\displaystyle n^{{\rm AO}\tilde{\rm c}}(E_{\rm F})$ $\displaystyle=$
$\displaystyle\frac{1}{2}\\!\left[1+\mathrm{sgn}(E_{\rm
F})\\!\left(\\!1-\mathrm{e}^{-\left(\\!\frac{2E_{\rm F}}{\gamma
t}\\!\right)^{2}}\right)\\!\right]\\!\\!,$ (56) $\displaystyle E^{{\rm
AO}\tilde{\rm c}}_{\mathrm{kin}}(E_{\rm F})$ $\displaystyle=$
$\displaystyle-\frac{t}{\sqrt{2\pi}}\frac{\gamma\pi}{4\sqrt{2}}\;\left[1-\operatorname{erf}\\!\left(\\!\frac{2|E_{\rm
F}|}{\gamma t}\\!\right)\right.$ (57) $\displaystyle\hskip
36.98866pt+\left.\frac{2}{\sqrt{\pi}}\frac{2|E_{\rm F}|}{\gamma
t}\mathrm{e}^{-\left(\\!\frac{2E_{\rm F}}{\gamma t}\\!\right)^{2}}\right]\\!.$

As Fig. 6(a) shows, this makes the AO$\tilde{\rm c}$ state have lower kinetic
energy than the FOc state close to quarter-filling, i.e., in the interval
$0.25<n<0.75$. This can be understood from the fact that the DOS of the phased
complex-orbital single-particle states has more weight at more negative
energies, see Fig. 2(a). However, as seen in Fig. 6(b), it is the FOr state
that clearly has the lowest kinetic energy of all ordered states due to its
large width, with minimum energy $-0.2443~{}t$, whereas the minimum for the
AO$\tilde{\rm c}$ state is only $-0.2216~{}t$, see Eq. (57). The basic reason
for the difference is that the FOr state benefits fully from the orbital-
flavor conserving hopping channel and also partially from the non-conserving
channel, compare Eq. (37), while the AO$\tilde{\rm c}$ state benefits fully
from the orbital-flavor non-conserving channel but not at all from the
conserving channel, cf. Eq. (34).

As regards the unpolarized states one may observe that for every para-orbital
state PX its kinetic energy is related to that of the associated ordered state
X by

$E_{\mathrm{kin}}^{\rm PX}(n)=2\,E_{\mathrm{kin}}^{\rm X}(n/2)\,,$ (58)

which follows directly from the defining equations (50) and (51). Equation
(58) implies that the relative magnitude of the kinetic energies of the
unpolarized states is the same as that of the ordered states, as indeed seen
in Fig. 6. In particular, the POr state has the lowest kinetic energy of all
para-orbital disordered phases, for the same reason as the FOr state is the
lowest-energy ordered state. Also, the para-orbital states have considerably
lower kinetic energy than their ordered counterparts, but they will get
renormalized at finite $U$ whereas the ordered states will not.

Finally we consider the trial state built from the single-particle eigenstates
by filling the lower and the upper subband up to the same Fermi energy,
written down here explicitly, compare Eqs. (40) – (42),

$|\Phi_{0}^{\rm OL}\rangle=\prod_{{\bf k};\varepsilon_{\ell}({\bf k})<E_{\rm
F}}\,\prod_{{\bf k^{\prime}};\varepsilon_{u}({\bf k^{\prime}})<E_{\rm
F}}\,\hat{d}_{{\bf k},{\ell}}^{\dagger}\,\hat{d}_{{\bf
k^{\prime}},u}^{\dagger}|0\rangle\,.$ (59)

This state we call the OL state because it is unpolarized in the following
remarkable way.

Whereas all single-particle states in $|\Phi_{0}^{\rm OL}\rangle$ are
polarized individually, compare Eqs. (47) and (48), the variation with ${\bf
k}$ makes their contributions to the pseudospin add up to zero at each lattice
site, at any filling $n(E_{\rm F})$ and for each subband $\ell$ and $u$
separately. This follows from the following considerations (for a more
detailed proof, see the Appendix): (i) for every single-particle state $|{\bf
k},\ell\rangle$ from the lowest subband which is occupied in $|\Phi_{0}^{\rm
OL}\rangle$ because it satisfies $\varepsilon_{\ell}({\bf k})<E_{\rm F}$, all
single-particle states $|{\bf k}^{\prime}_{n},\ell\rangle$ with wave vectors
${\bf k}^{\prime}_{n}$ generated from ${\bf k}$ by successive cyclic
permutations of its components are also occupied; this holds because both
$A_{\bf k}$ and $B_{\bf k}$ are invariant under such permutations, compare
Eqs. (14) and (16), implying that $\varepsilon_{\ell}({\bf
k}^{\prime}_{n})=\varepsilon_{\ell}({\bf k})$, and it follows that
$\varepsilon_{\ell}({\bf k}^{\prime}_{n})<E_{\rm F}$; the same reasoning
applies for the upper subband $u$; (ii) under each such transformation ${\bf
k}\mapsto{{\bf k}^{\prime}_{n}}$ the corresponding phase angle is transformed
as $\phi_{\bf k}\mapsto\phi_{{\bf k}}+4n\pi/d$, so summing $\exp({i\phi_{{\bf
k}^{\prime}_{n}})}$ over all ${\bf k}^{\prime}_{n}$ yields $0$, and the
compensation of the contributions to the pseudospin follows, compare again
Eqs. (47) and (48). Therefore, even though the densities in the two partially
filled subbands are different, $n_{\ell}(E_{\rm F})>n_{u}(E_{\rm F})$, as
recognized from the DOSs shown in Fig. 2(b), the OL state is unpolarized,
because each subband is unpolarized by itself at any filling. This mechanism
is entirely different from the compensation of two oppositely polarized bands
at equal filling as occurs for para-orbital (or paramagnetic) states.

We emphasize that the reasoning above is not specific for $d=\infty$ but holds
in any dimension, and so the resulting absence of polarization in both bands
is a characteristic feature of the OL per se. It is associated with the fact
that the OL transforms as $A_{1}$ under the symmetry group $C_{d}$ (is
invariant under U(1) in the limit $d\to\infty$) and implies that the OL phase
protects the (hyper)cubic symmetry. — We further point out that there is
necessarily double occupancy in $|\Phi_{0}^{\rm OL}\rangle$. This is obvious
from the fact that the single-particle states in $|\Phi_{0}^{\rm OL}\rangle$
do not correspond to a single common direction in the ‘equatorial’ plane of
Fig. 1(c) because of the variation of $\phi_{\bf k}$ with the wave vector
${\bf k}$.

The relevant DOS from which to obtain the particle density and the kinetic
energy in the OL state is $\rho_{\ell+u}(E)$, Eq. (44). One finds

$\displaystyle n_{\rm OL}(E_{\rm F})$ $\displaystyle=$ $\displaystyle
1+\mathrm{erf}\left(\frac{\sqrt{2}E_{\rm
F}}{t}\right)-\frac{\gamma}{\sigma\sqrt{2}}\mathrm{e}^{-2\left(\frac{E_{\rm
F}}{\sigma t}\right)^{2}}\mathrm{erf}\left(\frac{\gamma E_{\rm F}}{\sigma
t}\right)\,,$ (60) $\displaystyle E_{{\rm kin};\rm OL}(E_{\rm F})$
$\displaystyle=$
$\displaystyle-\frac{t}{\sqrt{2\pi}}\left\\{\,\mathrm{e}^{-2\left(\frac{E_{\rm
F}}{t}\right)^{2}}+\gamma\sqrt{\pi}\left[\frac{E_{\rm F}}{\sigma
t}\,\mathrm{e}^{-2\left(\frac{E_{\rm F}}{\sigma
t}\right)^{2}}\textrm{erf}\,\left(\frac{\gamma E_{\rm F}}{\sigma
t}\right)+\mathrm{F}\left(\gamma;\frac{|E_{\rm F}|}{\sigma
t}\right)\right]\right\\},$ (61)

where

$\mathrm{F}(\gamma;x)=\int_{x}^{\infty}\mathrm{e}^{-2y^{2}}\mathrm{erf}(\gamma
y)\,\mathrm{d}y.$ (62)

The function $\mathrm{F}(\gamma;x)$ has to be calculated numerically. Since
$\mathrm{F}(1;0)=\frac{1}{\sqrt{2\pi}}\left[\frac{\pi}{2}-\tan^{-1}(\sqrt{2})\right]=0.245541$
[90], we have rewritten Eq. (62) at $\gamma=1$ as

$\mathrm{F}(1;x)=0.245541-\int_{0}^{x}\mathrm{e}^{-2y^{2}}\mathrm{erf}(y)\,\mathrm{d}y.$
(63)

to perform the calculation for the full orbital case.

The kinetic energy of the OL state is included in both Fig. 6(a) and Fig.
6(b), and is seen to be lower than that of any other trial state, remarkably
including the lowest-energy complex-orbital and real-orbital para states. The
explanation is similar to that above for the FOr and the POr states. The OL
state fully captures all kinetic energy available in the orbital Hubbard model
Eq. (9), both from the orbital-flavor conserving hopping channel and from the
non-conserving hopping channel, cf. Eq. (42). This is a strong hint that the
OL phase could be favored also in the presence of electron interactions. To
establish the phase with the lowest overall energy we need to consider the
energy renormalization of the para-orbital states and the OL state in the
presence of the Coulomb interaction, i.e., at finite $U$.

##  IV Correlations by the Gutzwiller method

The original analysis of the Gutzwiller approach in the limit $d\to\infty$,
including the proof that the Gutzwiller approximation becomes exact in that
limit, was carried out specifically for the spin Hubbard model, implying the
implicit assumption of SU(2) symmetry [3, 8, 9]. We therefore need to
investigate which modifications are required in the orbital case where we have
only $C_{d}$ symmetry, and check if the proof of exactness still holds.

###  IV.1 General derivation at dimension $d$

Electron correlations due to the on-site Hubbard repulsion $U$ will be
implemented using the wave function introduced by Gutzwiller,

$|\Psi(g)\rangle=g^{\hat{D}}\left|\Phi_{0}\right\rangle=\prod_{i}\left(1-(1-g)\hat{D}_{i}\right)\left|\Phi_{0}\right\rangle\,,$
(64)

with

$\hat{D}=\sum_{i}\hat{D}_{i},\quad\quad\hat{D}_{i}=\hat{n}_{i+}\hat{n}_{i-},$
(65)

where the variational parameter $g$ is determined by minimizing the energy [6,
7, 8, 9],

$E(g)=\frac{\langle\Psi(g)|{\cal
H}_{d}|\Psi(g)\rangle}{\langle\Psi(g)|\Psi(g)\rangle}\,.$ (66)

First, we consider here the general case and derive the propagator and self-
energy. In the present case the fermion (orbital) flavor is not conserved and
we introduce a propagator,

${\cal
P}_{ij;\alpha\beta}=\frac{\big{\langle}\Psi\big{|}\hat{c}_{i\alpha}^{\dagger}\hat{c}_{j\beta}\big{|}\Psi\big{\rangle}}{\langle\Psi(g)|\Psi(g)\rangle}\,.$
(67)

which we allow to have, at least in principle, nonzero offdiagonal matrix
elements with respect to flavor (i.e., $\alpha\neq\beta$). Following Metzner
[9], one finds the following expression for the matrix elements,

$\displaystyle{\cal P}_{ij;\alpha\beta}$ $\displaystyle=$
$\displaystyle\big{\\{}c_{i\alpha}^{\dagger}c_{j\beta}[1-(1-g)(n_{i\bar{\alpha}}+n_{j\bar{\beta}})$
(68) $\displaystyle+$
$\displaystyle(1-g)^{2}(n_{i\bar{\alpha}}n_{j\bar{\beta}}+\delta_{ij}\delta_{\alpha\beta}n_{i\bar{\alpha}})X\big{\\}}_{0}^{c}\,,$

where

$X\equiv 1+\sum_{m=1}^{\infty}\frac{(g^{2}-1)^{m}}{m!}\sum_{g_{1}\dots
g_{m}}{D}_{g_{1}},\dots,{D}_{g_{m}}\,.$ (69)

Similarly, the bare propagator is defined by,

${\cal
P}_{ij;\alpha\beta}^{0}=\big{\langle}\Phi_{0}\big{|}\hat{c}_{i\alpha}^{\dagger}\hat{c}_{j\beta}\big{|}\Phi_{0}\big{\rangle}=\left\\{c_{i\alpha}^{\dagger}c_{j\beta}\right\\}^{c}_{0}\,.$
(70)

Here ${\cal P}^{0}$ is, like ${\cal P}$, an $N\times N\times 2\times 2$ matrix
with off-diagonal elements with respect to orbital flavor generated by the
off-diagonal hopping in orbital space. We emphasize that the curly brackets
$\\{\dots\\}_{0}^{c}$ in Eqs. (68) and (70) denote the sum over all connected
products of anticommuting contractions, evaluated in the trial state
$|\Phi_{0}\rangle$, and thus the ($c_{i\alpha}$, etc.) symbols are Grassman
variables and not the fermion operators ($\hat{c}_{i\alpha}$, etc.).

Next we define the self-energy ${\cal S}$. Here ${\cal S}_{ij}$ is a $2\times
2$ matrix with off-diagonal elements in orbital space. The elements of the
self-energy ${\cal S}$ are therefore different from those for the RVB wave
function [9], and include the specific last term below,

$\displaystyle{\cal S}_{ij;\alpha\beta}\\!$ $\displaystyle=$
$\displaystyle\left\\{\left[(1-g)^{2}c_{i\alpha}^{\dagger}n_{i\bar{\alpha}}c_{j\beta}n_{j\bar{\beta}}\right.\right.$
(71) $\displaystyle-$
$\displaystyle\\!\left.\left.\delta_{ij}(g^{2}-1)\\!\left(\delta_{\alpha\beta}n_{i\bar{\alpha}}\\!-\delta_{\alpha\bar{\beta}}c_{i\alpha}^{\dagger}c_{j\beta}\right)\right]\\!X\right\\}_{0}^{c}\\!.$

The self-energy gives the double occupancy,

$\overline{d}_{i}=\frac{1}{2}\,\frac{g^{2}}{1-g^{2}}{\rm
Tr}_{\alpha}\left({\cal S}{\cal P}^{0}\right)_{ii}\,.$ (72)

The further derivation, following closely the case of the RVB state considered
by Metzner [9], leads after some algebraic manipulations to the propagator in
real space between sites $i$ and $j$,

$\displaystyle{\cal P}_{ij}$ $\displaystyle=$ $\displaystyle{\cal
P}_{ij}^{0}+\delta_{ij}\left[\frac{1-g}{1+g}\left({\cal S}{\cal
P}^{0}\right)_{ii}\delta_{\alpha\beta}-\frac{{\cal S}_{ii}}{(1+g)^{2}}\right]$
(73) $\displaystyle+$ $\displaystyle\left[\left({\cal
P}^{0}-\frac{1}{1+g}\right){\cal S}\left({\cal
P}^{0}-\frac{1}{1+g}\right)\right]_{ij}\,.$

Here in the first local term at site $i$ only the diagonal elements of the
matrix ${\cal S}{\cal P}^{0}$ contribute, as indicated by
$\delta_{\alpha\beta}$.

The first term in Eq. (68) is related to the self-energy ${\cal S}$ by

$\big{\\{}c_{i\alpha}^{\dagger}c_{j\beta}X\big{\\}}_{0}^{c}\equiv\left({\cal
P}^{0}+{\cal P}^{0}{\cal S}{\cal P}^{0}\right)_{ij;\alpha\beta}\,.$ (74)

Finally, we define a “proper self-energy” ${\cal S}^{*}$ as the sum over all
one-particle irreducible diagrams contributing to ${\cal S}$—the two self-
energies are then related by a Dyson equation,

${\cal S}={\cal S}^{*}+{\cal S}^{*}\circ{\cal P}^{0}\circ{\cal S}\,.$ (75)

After introducing a renormalized propagator

$\overline{\cal P}={\cal P}^{0}+{\cal P}^{0}\circ{\cal S}\circ{\cal P}^{0}\,,$
(76)

one finds that, owing to Eq. (75), it satisfies

$\overline{\cal P}={\cal P}^{0}+{\cal P}^{0}\circ{\cal
S}^{*}\circ\overline{\cal P}\,.$ (77)

The site-diagonal matrix elements of ${\cal S}$ are related to those of
$\overline{\cal P}$ as follows,

$\displaystyle{\cal S}_{ii;\alpha\alpha}$ $\displaystyle=$
$\displaystyle\;\;(1-g^{2}){\overline{\cal
P}}_{ii;\bar{\alpha}\bar{\alpha}}\,,$ (78) $\displaystyle{\cal
S}_{ii;\alpha\bar{\alpha}}$ $\displaystyle=$
$\displaystyle-(1-g^{2})\overline{\cal P}_{ii;\alpha\bar{\alpha}}\,.$ (79)

###  IV.2 Single-site collapse at $d=\infty$

Similar to the spin Hubbard model at $d\to\infty$, the diagrams describing the
perturbation expansion of ${\cal S}^{*}$ with respect to the local interaction
$U$ collapse. Then the lattice sums reduce to summations of terms in which all
indices coincide on the same lattice site, and the proper self-energy becomes
site-diagonal, i.e., ${\cal S}^{*}_{ij}=\delta_{ij}{\cal S}^{*}_{ii}$, The
collapse permits the matrix elements of ${\cal S}^{*}_{ii}$ (with respect to
the orbital label) to be expressed in terms of those of the on-site matrix
elements of $\overline{\cal P}_{ii}$ defined above. One finds in matrix form,

${\cal S}^{*}_{ii}=f\left(\begin{array}[]{cc}-\overline{\cal
P}_{ii;--}&+\overline{\cal P}_{ii;+-}\\\ +\overline{\cal
P}_{ii;-+}&-\overline{\cal P}_{ii;++}\end{array}\right)\,.$ (80)

We emphasize that ${\cal S}^{*}_{ii}$ is not yet necessarily independent of
$i$. Here the factor $f$ is

$f=\frac{1-\left[1+4(1-g^{2})\det{\overline{\cal
P}_{ii}}\right]^{1/2}}{2\det\\{\overline{\cal P}_{ii}\\}}\,,$ (81)

and upon inversion,

$\overline{\cal P}_{ii}=\frac{1}{(1-g^{2})-\det{{\cal
S}^{*}_{ii}}}\left(\begin{array}[]{cc}+{\cal S}^{*}_{ii;--}&-{\cal
S}^{*}_{ii;+-}\\\ -{\cal S}^{*}_{ii;-+}&+{\cal
S}^{*}_{ii;++}\end{array}\right)\,.$ (82)

Equating this to the site-diagonal part of $\overline{\cal P}$ as obtained
from Eq. (77) generates two equations which together determine ${\cal S^{*}}$,
for any given ${\cal P}^{0}$ and $g$. Next, Eq. (75) gives the self-energy
${\cal S}$, which finally determines the propagator via Eq. (73) and the
double-occupancy using Eq. (72), and thus the total energy at any given
electron filling $n$.

###  IV.3 Momentum space representation

As we shall be mainly interested in translation invariant states, it is
expedient to Fourier transform all variables of interest upon which matrix
products with respect to site indices turn into ordinary products in ${\bf
k}$-space, e.g. Eq. (73) becomes

$\displaystyle{\cal P}({\bf k})$ $\displaystyle=$ $\displaystyle{\cal
P}^{0}({\bf k})$ $\displaystyle+$ $\displaystyle\sum_{\bf
k^{\prime}}\left[\frac{1-g}{1+g}\left({\cal S}({\bf k^{\prime}}){\cal
P}^{0}({\bf k^{\prime}})\right)_{\alpha\alpha}\delta_{\alpha\beta}-\frac{{\cal
S}({\bf k^{\prime}})}{(1+g)^{2}}\right]$ $\displaystyle+$
$\displaystyle\left[\left({\cal P}^{0}({\bf k})-\frac{1}{1+g}\right){\cal
S}({\bf k})\left({\cal P}^{0}({\bf k})-\frac{1}{1+g}\right)\right].$

Consider now a homogeneous state where ${\cal S}^{*}_{ii}\equiv S^{*}$ is
site-independent, and we may use a shorthand notation for its matrix elements,

${S}^{*}=\left(\begin{array}[]{cc}S^{*}_{+}&R^{*}\\\
(R^{*})^{+}&S^{*}_{-}\end{array}\right)=\left(\begin{array}[]{cc}S^{*}_{+}&|R^{*}|\mathrm{e}^{i\phi}\\\
|R^{*}|\mathrm{e}^{-i\phi}&S^{*}_{-}\end{array}\right)\,.$ (84)

Using the Dyson equation for $\overline{\cal P}$ in Eq. (77), $\overline{\cal
P}({\bf k})$ can be expressed (in ${\bf k}$-space) explicitly in terms of
${\cal P}^{0}({\bf k})$ and ${\cal S}^{*}({\bf k})$,

$\overline{\cal P}({\bf k})=\frac{1}{N_{\bf k}}\left[{\cal P}^{0}({\bf
k})-\det{{\cal P}^{0}({\bf k})}\left(\begin{array}[]{cc}S^{*}_{-}&-R^{*}\\\
-(R^{*})^{+}&S^{*}_{+}\end{array}\right)\right]\,,$ (85)

where the normalizing prefactor $N_{\bf k}$ is

$N_{{\bf k}}=1-\mathrm{Tr}_{\alpha}\left\\{S^{*}{\cal P}^{0}({\bf
k})\right\\}+(\det{S^{*}})\det{{\cal P}^{0}({\bf k})}\,.$ (86)

It depends on the occupancy of the single-particle states at wave vector ${\bf
k}$ in the many-particle trial state $|\Phi\rangle$. Upon solving for
$S^{*}_{+}$, $S^{*}_{-}$, $R^{*}$, and $(R^{*})^{+}$, one obtains the self-
energy from

${\cal S}({\bf k})=\frac{1}{N_{\bf
k}}\left[S^{*}-\det{S^{*}}\left(\begin{array}[]{cc}{\cal P}^{0}({\bf
k})_{--}&-\,{\cal P}^{0}({\bf k}_{+-}\\\ -\,{\cal P}^{0}({\bf k})_{-+}&{\cal
P}^{0}({\bf k})_{++}\end{array}\right)\right]\,,\\\ $ (87)

upon which it is straightforward to write down explicitly the expression for
the propagator ${\cal P}({\bf k})$ and calculate all quantities of physical
interest.

##  V Gutzwiller renormalization

###  V.1 Ordered states

As an example, we consider an ordered, i.e. polarized, state with complex
orbitals, Oc,

$\left|\Phi_{0}^{\rm Oc}\right\rangle=\prod_{{\bf k}\in{\cal
K}_{+}}{\hat{c}}_{{\bf k}+}^{\dagger}\prod_{{\bf q}\in{\cal
K}_{-}}{\hat{c}}_{{\bf q}-}^{\dagger}\,|0\rangle\,.$ (88)

It contains fixed numbers of electrons $n_{+}^{0}$ and $n_{-}^{0}$ in
$|i+\rangle$ and $|i-\rangle$ states ($n_{+}^{0}\geq n_{-}^{0}$),
respectively,

$\displaystyle n_{+}^{0}$ $\displaystyle=$ $\displaystyle\sum_{{\bf k}\in{\cal
K}_{+}}n_{{\bf k},+}^{0},$ $\displaystyle n_{-}^{0}$ $\displaystyle=$
$\displaystyle\sum_{{\bf q}\in{\cal K}_{-}}n_{{\bf q},-}^{0},$ (89)

which determine the total electron number $n_{0}$ and orbital polarization
$m_{0}$,

$\displaystyle n_{0}$ $\displaystyle=$ $\displaystyle n_{+}^{0}+n_{-}^{0},$
$\displaystyle m_{0}$ $\displaystyle=$ $\displaystyle n_{+}^{0}-n_{-}^{0}.$
(90)

The state (88) includes both full polarization, $n_{-}^{0}=0$ and
$n_{0}=m_{0}=n_{+}^{0}$, i.e., the complex ferro-orbital state FOc, and zero
polarization, $n_{0\pm}=\frac{1}{2}n_{0}$, $m_{0}=0$, i.e., the complex para-
orbital state POc.

The bare propagator is a diagonal $2\times 2$ matrix in the orbital
$\\{|+\rangle,|-\rangle\\}$ basis, i.e., ${\cal P}_{\alpha\beta}^{0}({\bf
k})=\delta_{\alpha\beta}n_{{\bf k}\alpha}^{0}$, and thus $\det\\{{\cal
P}^{0}({\bf k})\\}=n_{{\bf k}+}^{0}n_{{\bf k}-}^{0}$. Consequently, the
explicit form of the matrix elements of $\overline{\cal P}({\bf k})$ and
$\overline{\cal P}_{ii}\equiv\overline{P}$ simplify, and the equations
following from Eqs. (85) and (82) reduce to

$\displaystyle\overline{P}_{++}\\!$ $\displaystyle=$
$\displaystyle\frac{1}{N}\\!\sum_{\bf k}\overline{\cal P}({\bf k})_{++}$ (91)
$\displaystyle=$ $\displaystyle\frac{1}{N}\\!\sum_{\bf
k}\\!\left[\frac{n_{{\bf k}+}^{0}(1\\!-\\!n_{{\bf
k}-}^{0})}{1-S_{+}^{*}}\\!+\\!\frac{n_{{\bf k}+}^{0}n_{{\bf
k}-}^{0}(1-S_{-}^{*})}{(1\\!-\\!S_{+}^{*})(1\\!-\\!S_{-}^{*})-|R^{*}|^{2}}\right]$
$\displaystyle=$
$\displaystyle\frac{S_{-}^{*}}{1-g^{2}-(S_{+}^{*}S_{-}^{*}-|R^{*}|^{2})},$
$\displaystyle\overline{P}_{+-}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}\overline{\cal P}({\bf k})_{+-}$ (92)
$\displaystyle=$ $\displaystyle\frac{1}{N}\sum_{\bf k}\frac{n_{{\bf
k}+}^{0}n_{{\bf k}-}^{0}R^{*}}{(1-S_{+}^{*})(1-S_{-}^{*})-|R^{*}|^{2}}$
$\displaystyle=$
$\displaystyle-\frac{R*}{1-g^{2}-(S_{+}^{*}S_{-}^{*}-|R^{*}|^{2})},$
$\displaystyle\overline{P}_{--}\\!$ $\displaystyle=$
$\displaystyle\frac{1}{N}\\!\sum_{\bf k}\overline{\cal P}({\bf k})_{--}$ (93)
$\displaystyle=$ $\displaystyle\frac{1}{N}\\!\sum_{\bf
k}\\!\left[\frac{n_{{\bf k}-}^{0}(1\\!-\\!n_{{\bf
k}+}^{0})}{1-S_{-}^{*}}\\!+\\!\frac{n_{{\bf k}+}^{0}n_{{\bf
k}-}^{0}(1-S_{+}^{*})}{(1\\!-\\!S_{+}^{*})(1\\!-\\!S_{-}^{*})-|R^{*}|^{2}}\right]$
$\displaystyle=$
$\displaystyle\frac{S_{+}^{*}}{1-g^{2}-(S_{+}^{*}S_{-}^{*}-|R^{*}|^{2})},$

From Eq. (92) it follows that $R^{*}=0$, and Eqs. (91) and (93) are simplified
to

$\displaystyle\frac{n_{+}^{0s}}{1-S_{+}^{*}}+\frac{d^{0}}{1-S_{+}^{*}}$
$\displaystyle=$ $\displaystyle\frac{S_{-}^{*}}{1-g^{2}-S_{+}^{*}S_{-}^{*}},$
(94) $\displaystyle\frac{n_{-}^{0s}}{1-S_{-}^{*}}+\frac{d^{0}}{1-S_{-}^{*}}$
$\displaystyle=$ $\displaystyle\frac{S_{+}^{*}}{1-g^{2}-S_{+}^{*}S_{-}^{*}},$
(95)

where $n_{\alpha}^{0s}$ is the number of singly occupied ${\bf k}$-points of
flavor $\alpha\in\\{+,-\\}$ and $d^{0}$ is the number of doubly occupied ${\bf
k}$-points, hence

$n_{\alpha}^{0}=n_{\alpha}^{0s}+d^{0}.$ (96)

After some algebraic manipulation one finds

$n^{0}_{+}\,\frac{S_{+}^{*}}{1-S_{+}^{*}}=n^{0}_{-}\,\frac{S_{-}^{*}}{1-S_{-}^{*}},$
(97)

and the solution for $S_{+}^{*}$ and $S_{-}^{*}$ can be found.

Using the simplified (because of $R^{*}=0$) expressions for ${\cal S}({\bf
k})$ and $\overline{\cal P}({\bf k})$ one finds,

${\cal P}({\bf k})_{\alpha\beta}=\delta_{\alpha\beta}n_{{\bf
k}\alpha}=\delta_{\alpha\beta}\left[q_{\alpha}(g)n_{{\bf
k}\alpha}^{0}+b_{\alpha}\right],$ (98)

with

$\displaystyle q_{\alpha}(g)$ $\displaystyle\equiv$ $\displaystyle
1-\frac{S_{\alpha}^{*}}{1+g^{2}}\left(1-\frac{g^{2}}{1-S_{\alpha}^{*}}\right),$
(99) $\displaystyle b_{\alpha}$ $\displaystyle\equiv$
$\displaystyle\left[1-q_{\alpha}(g)\right]n_{\alpha}^{0}.$ (100)

This is the Gutzwiller approximation for the ordered state Eq. (88). Note that

$n_{\alpha}=\frac{1}{N}\sum_{\bf k}n_{{\bf
k}\alpha}=q_{\alpha}n_{\alpha}^{0}+b_{\alpha}=n_{\alpha}^{0},$ (101)

i.e., Luttinger theorem is satisfied. The average double occupancy is given by

$\overline{d}=\frac{g^{2}}{1-g^{2}}\frac{S_{+}^{*}}{1-S_{+}^{*}}\,n_{0+}=\frac{g^{2}}{1-g^{2}}\frac{S_{-}^{*}}{1-S_{-}^{*}}\,n_{0-}\,.$
(102)

One recognizes that the above results are identical to those obtained in the
spin case, i.e., for SU(2) symmetry.

For any other trial state built from a single pair of orthogonal orbitals, cf.
Eq. (4),

$\left|\Phi_{0}^{\psi\theta}\right\rangle=\prod_{{\bf k}\in{\cal
K}}{\hat{c}}_{\bf k}^{\dagger}(\psi,\theta)\prod_{{\bf q}\in{\cal
Q}}{\hat{c}}_{\bf k}^{\dagger}(\psi+\pi,\theta+\pi)|0\rangle,$ (103)

one obtains exactly analogous results when similarly writing all expressions
in terms of the corresponding orbital basis. So, e.g. for the real Or states
built from $\\{|{\bf k}x_{0}\rangle\\}$ and $\\{|{\bf q}\bar{x}_{0}\rangle\\}$
single-particle states, compare Eq. (88), one has

$n_{{\bf k}x_{0}}=q_{x_{0}}(g)n_{{\bf k}x_{0}}^{0}+b_{x_{0}}\,,$ (104)

compare Eq. (98). The renormalization $q_{x_{0}}(g)$ and the background
constant density $b_{x_{0}}$ are straightforwardly obtained from
$q_{\alpha}(g)$ and $b_{\alpha}$ by the corresponding orbital transformation.

###  V.2 Orbital liquid state

Next we consider the OL state introduced above, see Eq. (59),

$|\Phi_{0}^{\rm OL}\rangle=\prod_{{\bf k};\varepsilon_{\ell}({\bf k})<E_{\rm
F}}\,\prod_{{\bf k^{\prime}};\varepsilon_{u}({\bf k^{\prime}})<E_{\rm
F}}\,\hat{d}_{{\bf k},{\ell}}^{\dagger}\,\hat{d}_{{\bf
k^{\prime}},{u}}^{\dagger}|0\rangle\,,$ (105)

where ${\bf k}$ and ${\bf k^{\prime}}$ label the occupied states from the
lower (${\ell}$) and upper ($u$) subband, and both subbands are occupied up to
the same Fermi energy $E_{\rm F}$.

The free propagator is, in the $(+,-)$ basis,

${\cal P}^{0}_{\bf k}=\frac{1}{2}\left(\begin{array}[]{cc}n_{{\bf
k},{\ell}}^{0}+n_{{\bf k},u}^{0}&\left(n_{{\bf k},{\ell}}^{0}-n_{{\bf
k},u}^{0}\right)\,e^{i\phi_{\bf k}}\\\ \left(n_{{\bf k},{\ell}}^{0}-n_{{\bf
k},u}^{0}\right)\,e^{-i\phi_{\bf k}}&n_{{\bf k},{\ell}}^{0}+n_{{\bf
k},u}^{0}\end{array}\right),$ (106)

This leads to the following expressions for the on-site matrix of
$\overline{\cal P}$ in real space, $\overline{\cal P}_{ii}\equiv\overline{P}$,

$\displaystyle\overline{P}_{++}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}\left[\frac{n_{{\bf
k},\ell}^{0}(1-n_{{\bf k},u}^{0})}{2N_{{\bf k}s}}+\frac{n_{{\bf
k},\ell}^{0}n_{{\bf k},u}^{0}(1-S^{*}_{-})}{N_{{\bf k}d}}\right],$
$\displaystyle\overline{P}_{+-}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}\left[\frac{n_{{\bf
k},\ell}^{0}(1-n_{{\bf k},u}^{0})\mathrm{e}^{i\phi_{\bf k}}}{2N_{{\bf
k}s}}+\frac{n_{{\bf k},\ell}^{0}n_{{\bf
k},u}^{0}|R^{*}|\mathrm{e}^{i\phi}}{N_{{\bf k}d}}\right],$
$\displaystyle\overline{P}_{--}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}\left[\frac{n_{{\bf
k},\ell}^{0}(1-n_{{\bf k},u}^{0})}{2N_{{\bf k}s}}+\frac{n_{{\bf
k},\ell}^{0}n_{{\bf k},u}^{0}(1-S^{*}_{+})}{N_{{\bf k}d}}\right]\,.$

where the normalizing prefactor $N_{\bf k}$ depends on the occupancy of the
subbands at wave vector ${\bf k}$: it equals $N_{{\bf k}s}$ if only a single
(viz. the lower) subband is occupied and equals $N_{{\bf k}d}$ if both $|{\bf
k},\ell\rangle$ and $|{\bf k},{\rm u}\rangle$ are occupied, with

$\displaystyle N_{{\bf k}s}$ $\displaystyle=$ $\displaystyle
1-\frac{1}{2}(S^{*}_{+}+S^{*}_{-})-|R^{*}|\cos(\phi-\phi_{\bf k})$ (110)
$\displaystyle N_{{\bf k}d}$ $\displaystyle=$ $\displaystyle
1-S^{*}_{+}-S^{*}_{-}+\det{S^{*}}$ (111) $\displaystyle=$
$\displaystyle(1-S^{*}_{+})(1-S^{*}_{-})-|R^{*}|^{2}\,.$

Setting the expressions above equal to those obtained in Eq. (82) one obtains
a set of four equations from which to determine $S^{*}_{+}$, $S^{*}_{-}$,
$|R^{*}|$, and $\phi$.

Because of the occurrence of the phases $\phi-\phi_{\bf k}$ in Eqs. (V.2) -
(V.2) these equations cannot be solved in full generality and a simplifying
assumption is necessary. Thus we make here the ansatz that the off-diagonal
elements of ${\cal S}^{*}$ vanish, i.e., $R^{*}=0$. Then the terms that
involve summation over the phase factors $\\{\phi_{\bf k}\\}$ drop out and the
structure of the equations needed to determine ${\cal S}^{*}$ simplifies.
Because the OL state is not polarized the two flavors are equivalent and the
diagonal elements of ${\cal S}^{*}_{ii}$ must be equal,
$S^{*}_{+}=S^{*}_{-}\equiv S^{*}$. Defining single and double occupancy in
${\bf k}$-space per ${\bf k}$-point,

$\displaystyle\overline{S}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}(1-n_{{\bf k},u})\,n_{{\bf k},{\ell}}\,,$
(112) $\displaystyle\overline{D}$ $\displaystyle=$
$\displaystyle\frac{1}{N}\sum_{\bf k}n_{{\bf k},u}\,n_{{\bf k},{\ell}}\,,$
(113)

one obtains a single equation for $S^{*}$ in terms of $\overline{S}$ and
$\overline{D}$,

$\frac{1}{2}\frac{\overline{S}}{1-S^{*}}+\frac{\overline{D}}{1-S^{*}}=\frac{S^{*}}{1-(S^{*})^{2}-g^{2}}\,.\\\
$ (114)

Using the number of electrons per site $n=\overline{S}+2\overline{D}$, one
finds finally for the orbital liquid state

$S^{*}=\frac{1}{2-n}\left(1-\sqrt{1-n(2-n)(1-g^{2})}\right)\,.$ (115)

Inserting this result into the expressions for the renormalization factor and
double occupancy in terms of $S^{*}$,

$\displaystyle q$ $\displaystyle=$ $\displaystyle
1-\frac{S^{*}}{(1+g)^{2}}\,\left(1-\frac{g^{2}}{1-S^{*}}\right)\,,$ (116)
$\displaystyle\overline{d}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\frac{g^{2}n}{1-g^{2}}\,\frac{S^{*}}{1-S^{*}}\,,$
(117)

one finally obtains

$\displaystyle q(g)$ $\displaystyle=$ $\displaystyle
1-\frac{1}{(1+g)^{2}}\,\frac{\left(1-\sqrt{1-n(2-n)(1-g^{2})}\right)^{2}}{n(2-n)}\,,$
$\displaystyle\overline{d}(g)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{1-n(2-n)(1\\!-\\!g^{2})}+n(1\\!-\\!g^{2})-1}{2(1-g^{2})}.$
(119)

This expression for the kinetic energy renormalization (V.2) is identical to
that which follows from Eq. (99) when the ordered state is unpolarized, i.e.,
for $q_{+}(g)=q_{-}(g)\equiv q(g)$. At $U=\infty$ one finds $g=0$ and

$q(0)=\frac{1-n}{1-n/2},$ (120)

which reproduces the kinetic energy renormalization in the spin model.

##  VI Brinkman-Rice transition

In the following we focus on the OL state, and in particular compare its role
in the orbital Hubbard model with that of the paramagnetic phase in the spin
Hubbard model. For the latter purpose it is convenient to introduce
$t^{\prime}\equiv t/2$ – then the results obtained here at $d=\infty$ for
orbital phases are directly comparable with those for magnetic phases [11].

Figure 7:  Total energy $E_{\rm tot}$ (121) versus electron density $n$ for
the OL phase in the orbital $e_{g}$ model with $t^{\prime}\equiv t/2$. The
different curves stand for $U=0$ (black dashed line), metallic regime with
$U=2t^{\prime}$, $4t^{\prime}$, $6t^{\prime}$ (blue solid lines) and
insulating regime $U=10t^{\prime}$ (red heavy line).

Altogether the energy of the OL is determined by

$E_{\rm tot}^{\rm OL}(n;g)=q(n;g)E_{\rm kin}^{\rm OL}(n)+U\overline{d}(n;g),$
(121)

and one has to use Eqs. (V.2) and (119) with the optimal value of the
variational parameter $g$. The result for the density dependence of $E_{\rm
tot}$ for the OL phase is shown in Fig. 7 for a number of values of $U$. One
observes that the renormalization of the kinetic energy due to electron
correlations induced by $U$ is largest at half-filling, i.e., at $n=1$.
Indeed, here one might expect an insulating state to appear eventually at
sufficiently large $U$, with completely suppressed electron dynamics, see
below. Further, with increasing $U$ the total energy (121) gradually increases
at any density $n$ and its minimum moves towards $n=0.5$. The pattern seen in
Fig. 7 is entirely similar to that shown by the paramagnetic phase in the spin
Hubbard model when treated likewise by the Gutzwiller method, because it is
dominated by the behavior of $q(n;g)$ and $\overline{d}(n;g)$ and hardly
affected by the slight difference in the shape of $E_{\rm kin}(n)$.

Figure 8:  Dependence on $U/t^{\prime}$, at electron density $n=1$, of (a) the
variational parameter $g$, (b) the double occupancy $\overline{d}$, and (c)
the total energy $E_{\rm tot}/t^{\prime}$, for the spin Hubbard model (blue
dashed lines) and for the orbital $e_{g}$ Hubbard model (red full lines), at
$d=\infty$. For easier comparison with the spin case, we use here the hopping
parameter $t^{\prime}\equiv t/2$ for both spin and orbital Hubbard model.

Figure 8 shows in some more detail what happens at $n=1$, demonstrating that
for the OL in the orbital model the dependence of the variational parameter
$g$, the double occupancy $\overline{d}$, and the total energy $E_{\rm tot}$
on $U$ is qualitatively the same as for the paramagnetic state in the spin
model. In particular double occupancy gets fully suppressed ($\overline{d}=0$)
at and above a critical value $U_{\rm BR}$, while simultaneously $g=0$
implying $q=0$, i.e., the electrons are entirely localized and itinerancy
disappears. This has been recognized in the spin case by Brinkman and Rice
[91] as a qualitative criterion for the onset of the transition from a
metallic to an insulating state [92], which in fact becomes exact at
$d=\infty$ [93].

Figure 9:  Evolution of the OL phase from free electrons to a strongly
correlated metal at $n=0.9$ for increasing $U/t^{\prime}$ as obtained for the
orbital $e_{g}$ Hubbard model (solid lines) at $d=\infty$: (a) the Gutzwiller
renormalization factor $q$, (b) the double occupancy $\overline{d}$, and (c)
the total energy $E_{\rm tot}$. For easier comparison with the spin case, we
use here the hopping parameter $t^{\prime}\equiv t/2$ for the orbital Hubbard
model (9).

However, quantitatively there is a considerable difference: for the OL in the
orbital $e_{g}$ Hubbard model this transition occurs at $U_{\rm
BR}=9.161\,t^{\prime}$, enhanced from $U_{\rm BR}=6.383\,t^{\prime}$ for the
paramagnetic state in the spin Hubbard model (equivalent to the ($\gamma=0$)
POc state in the orbital model), because $U_{\rm BR}$ is generally given by

$U_{\rm BR}=8\,|E_{\rm kin}(E_{\rm F}=0)|\,,$ (122)

which is larger by the factor $1+\sqrt{\pi}F(1;0)=1.435$ for the OL, compare
Eqs. (61) and (53). Clearly this is entirely due to the extra contribution to
the kinetic energy coming from the orbital-flavor non-conserving hopping at
$\gamma=1$. We note in passing that the real para-orbital phase POr shows a
Brinkman-Rice transition at the intermediate value $U_{\rm
BR}=7.818\,t^{\prime}$, enhanced only by a factor $w_{\rm r}/w_{\rm
c}=\sigma(\gamma=1)=\sqrt{3/2}=1.225$, because it profits only partially from
the additional kinetic energy available in the orbital model.

At $n<1$ the renormalization factor $q(g)$, Eq. (V.2), remains finite at large
$U$, and thus the electrons remain itinerant, be it with reduced kinetic
energy, see Fig. 9. Increasing $U/t^{\prime}$ defines here three distinct
regimes: (i) at very small $U$ ($U\lesssim t^{\prime}$) correlations are
virtually absent, $q\approx 1$ and $\overline{d}\approx(n/2)^{2}$, (ii) at
small $U$ ($U<U_{\rm BR}$) the correlations gradually develop, $q(g)$ is
reduced from 1.0 down to $\approx 0.2$, double occupancy is strongly reduced
($\overline{d}\to 0$), and the total energy increases, and (iii) at large $U$
($U>U_{\rm BR}$) the correlations dominate, double occupancy becomes excluded
and the electrons move in the strongly correlated OL without creating double
occupancies.

##  VII Phase diagram

Figure 10:  Total energy $E_{\rm tot}$ as a function of density $n$ at $U=10t$
for the OL state (red line), the real-orbital POr state (blue line), the real-
orbital FOr state (black line), the complex-orbital FOc state (green line),
and the phased complex-orbital AO$\tilde{{\rm c}}$ state (aubergine line): (a)
for OL, POr, FOr, FOc, and AO$\tilde{{\rm c}}$ in the range $0\leq n\leq 1$;
(b) zoom-in on the range close to $n=1$ for OL, POr, and FOr only.

From the unrenormalized kinetic energies shown in Fig. 6 we concluded that the
uniformly polarized FOr phase and possibly the disordered POr phase are the
only candidates competing with the OL phase for being the ground state after
the energies of the POr and OL phase have been renormalized by correlations,
for the following reasons: First, as double occupancies do not occur in the
FOr phase or the other uniform phases with broken symmetry, these phases are
not affected by renormalization and so the sequence of their energies is
fixed. Second, although the energies of the disordered phases are changed by
accounting for correlations, their sequential order also cannot change if all
of them are subjected to the same renormalization procedure (at least if the
procedure is only sensitive to energies as is the case for the Gutzwiller
method). This holds because for any phase its total renormalized energy is
given by an expression like Eq. (121) with the variational parameter(s)
optimized for the lowest outcome; if the kinetic energy in this equation is
replaced by a lower value associated with a different phase, the total energy
is lowered even if the variational parameter $g$ is kept fixed; subsequently
allowing $g$ to optimize can only lead to a further lowering of the total
energy. This argument should not only hold for the para-orbital phases but
also be valid between the POr and OL phases, implying that the only
competition of interest should be between the FOr phase and the OL phase.

Figure 11:  Total energy $E_{\rm tot}$ as a function of the Hubbard
interaction parameter $U/t$ for the OL phase (red line), real-orbital FOr
state (black horizontal line, independent of $U$), and real-orbital POr state
(blue line) for fixed $n=0.95$.

Numerical results are shown in Fig. 10(a) for the dependence upon density of
the total energies of various phases at a particular value of $U$, viz.
$U=10t$. The states other than the FOr and the OL states are included here
just to demonstrate that they do not play a role, as argued above. Yet indeed,
a competition between the FOr phase and the OL does take place and the total
energies of these phases are seen to be quite close to one another. In fact
they cross and, whereas the energy of the OL is well below that of the FOr
state over most of the density range even for this rather large value of $U$,
the OL phase is unstable against the FOr phase close to $n=1$, as shown in
more detail in Fig. 10(b). We point out that the latter result is different
from that obtained for the 3D orbital Hubbard model (at $d=3$) where the OL is
stable in the entire regime of $n$ also at large $U$ [13].

The transition between the OL phase and the FOr phase can be determined more
accurately by calculating at fixed $n$ the total energy of the OL as a
function of $U$ and comparing this with the $U$-independent total energy of
the FOr state, as shown in Fig. 11. One thus obtains $U_{c}(n)$ where the two
phases have equal energy. This defines the OL-FOr phase boundary because the
transition is first-order, since the states $|\Psi^{\rm OL}\rangle$ and
$|\Psi^{\rm FOr}\rangle$ cannot transform continuously into one another
because of the globally different structure of $|\Phi_{0}^{\rm OL}\rangle$ and
$|\Phi_{0}^{\rm FOr}\rangle$. This stands in clear contrast to the
paramagnetic-ferromagnetic transition in the spin Hubbard model, which is
second order because the magnetization can evolve gradually from the
paramagnetic state.

The resulting phase diagram of the orbital $e_{g}$ Hubbard model in the ($n$,
$U$) plane as obtained in the Gutzwiller approximation is shown in Fig. 12.
The FOr phase is more stable than the OL phase for $U>U_{c}(n)$ if $n>n_{c}$.
Here $n_{c}=0.8746$ is the critical value of the particle density at which the
energies of the phases are equal at $U=\infty$, and below which the OL is
therefore always stable, see Fig. 12. It is obtained from a numerical
comparison using the analytical expressions for the kinetic energies of the
two phases, Eq. (55) and Eq. (61), and for the renormalization factor at
$U=\infty$, Eq. (120).

Figure 12: Phase diagram of the OL state versus the FOr state at $d=\infty$.
The FOr phase is more stable than the OL phase for $U>U_{c}(n)$ if $n>n_{c}$.
The vertical red dashed line indicates the phase boundary for $U=\infty$ at
$n_{c}$; $U_{c}(1)=4.5805\,t$ is equal to the critical value of the
interaction for the Brinkman-Rice transition [91] in the OL state.

Next we consider the intersection of the phase boundary with the $U/t$-axis at
$n=1$. With our simple representation of real-orbital ferro order by the state
$\big{|}\Psi^{\rm FOr}\big{\rangle}=\big{|}\Phi_{0}^{\rm FOr}\big{\rangle}$,
which has fully saturated orbital order such that $E_{\rm tot}^{\rm
FOr}(1)=E_{\rm kin}^{\rm FOr}(1)=0$ independent of $U$, this intersection
necessarily coincides with the Brinkman-Rice transition point where the energy
of the OL attains $0$, so $U_{c}(1)=U_{\rm BR}=4.5805\,t$. However, as pointed
out by Fazekas et al. for the spin case [11], precisely close to $n=1$ a small
reduction of the polarization must occur involving the exponential tails of
the DOSs of the $x_{0}$\- and ${\bar{x}}_{0}$-bands. If this would be taken
into account it would produce some lowering of the energy of the FOr phase
thus leading to a slight decrease of $U_{c}(1)$.

Allowing some polarization of the FOr state at general density, at the price
of creating double occupancy, would similarly give some energy lowering.
However, this is severely limited by the restriction that upon transferring
one particle from the majority to the minority band the gain in kinetic energy
must exceed $U$. This condition is only met as long as the Fermi energy of the
minority band is in the exponential tail of the DOS, i.e. only for a very
small fraction of the particles, and also $U\lesssim 4t$, compare Fig. 2(b).
Thus our use of the fully ordered state $\big{|}\Psi^{\rm FOr}\big{\rangle}$
neglects at most a minor expansion of the FOr regime in the phase diagram.

We have not attempted to construct a more sophisticated trial state for any of
the AO phases, such as has been shown to be important to obtain a good
description of the AF phase in the spin Hubbard model [8, 11]. So, strictly we
cannot exclude that alternating orbital order would show up in the orbital
model at and/or close to half-filling by outcompeting the ferro orbital order,
as occurs between antiferromagnetism and ferromagnetism in the spin case.
However, this seems unlikely in the orbital case because of the large
difference in kinetic energy between the fully ordered states, see Fig. 6,
viz. by a factor $w_{\rm r}/w_{\rm r,ao}=\sigma\sqrt{2}=\sqrt{3}=1.73$ between
the FOr and AOr states, and (because the kinetic energy of the AO${\tilde{c}}$
state is very close to that of the FOc state) by a factor $\simeq w_{\rm
r}/w_{\rm c}=\sigma=\sqrt{3/2}=1.22$ between the FOr and AO${\tilde{c}}$
states. As shown above, this difference is due to the more efficient use of
both pseudospin-conserving and non-conserving hopping channels by real-orbital
ferro order than by alternating orbital order. What we do have established is
that the OL state is the ground state over most of the ($n,U$) phase diagram,
as illustrated by Fig. 12.

Finally, we point out that the phase diagram obtained here at $d=\infty$ is
qualitatively different from that found before at $d=3$ [13]. In the latter
case the lowest-energy ordered FOx phase was completely eliminated from the
phase diagram by the non-conserving hopping, and the OL phase was the only
phase present. We propose that this is caused by the different energy ranges:
$\left[-3t,3t\right]$ at $d=3$ versus $\left[-\infty,\infty\right]$ at
$d=\infty$. The argumentation is as follows:

Let us consider, at $U=\infty$, the derivatives with respect to $n$ of the
total energy of the candidate phases (viz. FOr and OL at $d=\infty$, FOx and
OL at $d=3$) at $n=1$, as these determine which phase has the lowest energy in
the immediate vicinity of half-filling, compare Fig. 10. Since
$\overline{d}=0$ at $U=\infty$, we have $E_{\rm tot}^{\rm FO}(n)=E_{\rm
kin}^{\rm FO}(n)$ and $E_{\rm tot}^{\rm OL}(n)=q(n)E_{\rm kin}^{\rm OL}(n)$,
compare Eq. (121), and therefore

$\displaystyle\frac{dE_{\rm tot}^{\rm FO}}{dn}(n)$ $\displaystyle=$
$\displaystyle\frac{dE_{\rm kin}^{\rm FO}}{dn}(n)=E_{\rm F}^{\rm FO}(n)\,,$
(123) $\displaystyle\frac{dE_{\rm tot}^{\rm OL}}{dn}(n)$ $\displaystyle=$
$\displaystyle\frac{dq}{dn}(n)E_{\rm kin}^{\rm OL}(n)+q(n)\frac{dE_{\rm
kin}^{\rm OL}}{dn}(n)$ (124) $\displaystyle=$
$\displaystyle\frac{dq}{dn}(n)E_{\rm kin}^{\rm OL}(n)+q(n)E_{\rm F}^{\rm
OL}(n)\,.$

Both expressions can be evaluated, (i) at $d=3$: by consulting [13], for the
FOx phase using Eq. (3.23) therein (from which one obtains $E_{\rm kin}^{\rm
FOx}(1-x)\eqsim 3t(1-x+\cdots)$ and thus $E_{\rm F}^{\rm FOx}(1)=3t$), for the
OL phase making use of the data plotted in Fig. 5 therein [from which one
obtains $E_{\rm kin}^{\rm OL}(1)\eqsim-1.52t$ and $E_{\rm F}^{\rm OL}(1)=0$];
(ii) at $d=\infty$: for the FOr phase using Eqs. (54) in the present paper
with $w_{\rm r}=\sqrt{3/2}/2$, for the OL phase using Eq. (61) and so obtain
$E_{\rm kin}^{\rm OL}(1)\eqsim-0.5726t$ ; (iii) inserting $q(1)=0$ and
$dq/dn(1)=-2$.

The results are as follows:

$\displaystyle d$ $\displaystyle=$ $\displaystyle 3:\frac{dE_{\rm tot}^{\rm
FOx}}{dn}(1)=3.00\,t,$ (125) $\displaystyle d$ $\displaystyle=$ $\displaystyle
3:\frac{dE_{\rm tot}^{\rm OL}}{dn}(1)=3.04\,t,$ (126) $\displaystyle d$
$\displaystyle=$ $\displaystyle\infty\\!:\frac{dE_{\rm tot}^{\rm
FOr}}{dn}(1-x)\\!=\\!0.8660\>\mathrm{inverf}(1-2x)\,t,$ (127) $\displaystyle
d$ $\displaystyle=$ $\displaystyle\infty\\!:\frac{dE_{\rm tot}^{\rm
OL}}{dn}(1)=1.1451\,t\,.$ (128)

So, for $d=3$ the derivative is slightly larger for the OL, implying that
close to $n=1$ this is the lowest state and then will remain so in the entire
phase diagram, as we found before [13]. However, for $d=\infty$ the slope for
the FOr phase initially diverges (the function $\mathrm{inverf}(y)\to\infty$
at $y\to 1$), so the FOr state is the ground state close to $n=1$, as we found
in the present paper. This divergence is not seen in the plots because it is
very steep. Thus already for $x=1-n$ larger than $0.03$ the derivative is of
order 1 and with further increasing $x$ the FOr phase begins to be overtaken
by the OL. So the reason for the FOr phase being stable at $d=\infty$
apparently lies in the exponential tails in its DOS, particle density and
kinetic energy. As this is a mathematical rather than a physical phenomenon,
the above result suggests that stability of ordered phases at $d=\infty$ is
mathematically correct but physically spurious, and that the earlier
theoretical result at $d=3$ [13] describes the generic physics.

We thus conclude that regarding the phase diagram the $d=\infty$ case is not
representative for all dimensions: the $d=3$ case is qualitatively different.
This actually demonstrates that the effect of the non-conserving hopping terms
themselves is qualitatively different for different dimensions, which is not
the case for the conserving terms, since FM and AF phases appear in the phase
diagram of the spin Hubbard model in all dimensions.

##  VIII Discussion and conclusions

###  VIII.1 General aspects of the $d=\infty$ orbital model

A few aspects of the $e_{g}$ orbital Hubbard model at $d=\infty$ deserve some
further comments. They mainly concern the limitations of the Gutzwiller
approach in comparison with more powerful recent methods, in particular
selfconsistent DMFT. As mentioned briefly in Sec. I the development of the
DMFT approach [62] was made possible by the discovery by Metzner and Vollhardt
that in the limit $d=\infty$ only on-site correlations survive [8]. Nowadays
DMFT is recognized as a standard method to study electron correlation effects
in the electronic structure [94, 95, 96, 97]. This method has been
successfully applied, inter alia, to the Falicov-Kimball model [98] and to
nonequilibrium dynamics [99]. The ideas employed using the Gutzwiller wave
function helped to formulate the DMFT method in spin systems [95].

First of all, one may wonder to what extent the description of the metal-
insulator transition, as presented in Sec. VI, i.e., as a Brinkman-Rice
transition, is realistic. An additional aspect for the orbital model are
possible crystal-field splittings but it has been established that orbital-
selective Mott phases occur then for large enough interaction $U$ [100, 101].
However, when the orbitals remain equivalent as in the OL state considered in
the present paper, one may expect that the transition occurs in a similar way
as in the paramagnetic state in the spin model. So our treatment of the metal-
insulator transition in the orbital model at $d=\infty$ by the Gutzwiller
approach, leading to the Brinkman-Rice transition, then has the same status as
the treatment of the metal-insulator transition in the spin model at
$d=\infty$ by that same approach, and gives valuable insight how correlations
induce the system to approach the localized state [95].

Second, past calculations using the Gutzwiller wave function have led to a
better understanding of the strongly correlated regime of the Hubbard model.
In one dimension the metal-insulator transition is absent but several
quantities have been calculated exactly [102]: the double occupation, the
momentum distribution, as well as its discontinuity at the Fermi surface.
These quantities determine the expectation value of the 1D Hubbard Hamiltonian
for any symmetric and monotonically increasing dispersion. The Gutzwiller wave
function has also been found to predict ferromagnetic behavior for
sufficiently large interaction $U$ [103, 104]. This agrees with the present
result that the orbital Hubbard model also has a range of ferro-orbital states
close to $n=1$.

However, it has been known for some time that the Brinkman-Rice picture of the
metal-insulator transition, as given by the Gutzwiller approach both at $d=3$
and at $d=\infty$, is oversimplified, in particular yielding a poor
description of the dynamics on the insulating side [105]. Recent studies of
the metal-insulator transition in the spin Hubbard model at $d=\infty$ by DMFT
have demonstrated the existence of two different transitions, a metal-to-
insulator transition at $U_{c2}$ where the quasiparticle weight becomes zero
and an insulator-to-metal transition at $U_{c1}$ where the gap closes, both
first-order, and the occurrence of hysteresis between the two critical
interaction strengths [106, 107, 108]. This holds from $T=0$ up to a critical
temperature $T_{c}$ where $U_{c1}$ and $U_{c2}$ coalesce, while for $T>T_{c}$
there is a smooth crossover between metallic and insulating behavior. Notably,
the thermodynamics of the hysteresis region was recently explained at the two-
particle level [109]. Obviously, the above more-detailed understanding of the
Mott-Hubbard transition implies that this transition is considerably more
complex than the Brinkman-Rice picture suggests.

These differences between the Mott-Hubbard transition in the spin Hubbard
model as calculated by DMFT and the Brinkman-Rice picture are clearly due to
the fact that DMFT is coping more effectively with correlations than the
Gutzwiller approach. One would therefore expect that a DMFT treatment of the
$e_{g}$ orbital Hubbard model would produce similar changes to the metal-
insulator transition as in the spin Hubbard model, because this transition is
typically a correlation-induced phenomenon. By contrast, one would expect the
more prominent role of the disordered OL in the phase diagram (as compared to
the paramagnetic phase in the spin case) to remain largely unaffected, because
it is not due to correlations but to the additional kinetic energy provided by
the non-conserving hopping channel.

However, these are just speculations. We emphasize that orbital models such as
that in Eq. (9) [or Eq. (7)] have as yet not been studied in the context of
self-consistent DMFT calculations. Therefore, even the corresponding effective
impurity model for such calculations is not known, and its construction
appears not to be entirely trivial. A prime candidate is a two-orbital
Anderson impurity model with on-site $U$ and selfconsistently generated
orbital flipping terms. In view of the overriding importance of the orbital
non-conserving hopping channel and the associated $C_{d}$ symmetry
demonstrated above, these features should apparently be incorporated somehow.
The only way seems to be by imposing a specific structure on the bath modes
and their coupling to the impurity. It would be desirable if a DMFT treatment
along these lines were performed, both to see if the above speculations on the
$e_{g}$ orbital Hubbard model are born out and to further explore the
application of the DMFT method to a wider range of models.

The early ideas of Kugel and Khomskii [19] culminate in systems where spins
and orbitals become almost equivalent. This idealization happens for face-
sharing octahedra and one finds indeed a highly symmetric SU(4) model [110].
So far, studies of such systems have only been performed theoretically and
these have identified an interesting evolution of spin-orbital entanglement
with increasing spin-orbit coupling [66]. In spite of the absence of long-
range spin-orbital order such systems exhibit features typical of those
manifesting themselves at phase transitions [111]. The experimental search for
quantum spin-orbital liquids identified FeSc2S4 as a quantum material with
long-range entanglement [112]. Such states are primarily realized in the
Kitaev materials where spins and orbitals lose their identity and are strongly
entangled by large spin-orbit coupling [113]. The challenge in the theory is
to find a quantum spin-orbital liquid at infinite dimension.

###  VIII.2 Summary

Summarizing, we have presented a generalization of the $e_{g}$ orbital model
to infinite dimension $d=\infty$, preserving the two-fold degeneracy of the
orbitals and turning the lattice symmetry from cubic into hypercubic. It is
quite remarkable that the two-flavor orbital model is manifestly different
from the corresponding spin model at $d=\infty$ by the presence of orbital-
flipping hopping terms. At the same time it is somewhat surprising that the
two orbital flavors are equivalent in the limit of $d=\infty$ when they are so
different in the 3D model, including the physics of the manganites [80]. We
have further shown that the Gutzwiller approximation becomes exact in the
limit $d\to\infty$ for the orbital Hubbard model in perfect analogy with the
spin Hubbard model.

We conclude that the peculiar features of the orbital Hubbard model (9) are
due to the fact that the extra hopping terms $\propto\gamma t$ do not conserve
the pseudospin as they mix the two orbital flavors. They thus induce the
following distinctive features: (i) there are two classes of single-particle
plane-wave states: those with complex orbitals and those with real orbitals,
which behave differently; the same holds for the Fermi-sea-like multiparticle
states built from them; (ii) the single-particle eigenstates of the kinetic
energy form two distinct bands with different dispersions; (iii) each such
eigenstate carries a nonzero polarization, i.e., contributes to the pseudospin
at all sites; yet when summed over all eigenstates in an energy shell, these
contributions add up to zero; by this mechanism the orbital liquid state (OL),
which is obtained by filling the two bands of single-particle eigenstates to
the same Fermi energy, is unpolarized at any filling. The above qualitative
features are generic for the $e_{g}$ orbital Hubbard model, i.e., they hold at
any dimension, from $d=3$ up to $d=\infty$.

There are also a number of features induced by the non-conserving hopping
channel that are to some extent quantitative and are most pronounced at
$d=\infty$. (iv) The extra hopping channel lowers the kinetic energy with
respect to the spin case and thus makes the Brinkman-Rice transition to an
insulator occur at a larger value of $U$; (v.a) phases with alternating
orbitals are much less favored than antiferromagnetism for fermions on a
bipartite lattice like the hypercubic lattices considered here, because they
benefit relatively less from the pseudospin-conserving channel, whereas (v.b)
unlike ferromagnetic states in the spin Hubbard model, ferro-orbital states
are not eigenstates of the orbital Hubbard model (9) and benefit energetically
more from the orbital-mixing term; (v.c) together this makes the FOr state the
main competitor of the orbital liquid phase, and the $d=\infty$ phase diagram
confirms that FOr order indeed occurs at large $U$ close to half-filling
[114], while the orbital liquid phase is the ground state elsewhere; (v.d) the
stability of the FOr is apparently due to the exponential tails in the DOS,
particle density and kinetic energy, and so is specific for $d=\infty$,
whereas at $d=3$ the OL phase is the lowest-energy state in the entire phase
diagram. Perhaps feature (v.a) is the most remarkable quantitative consequence
of orbital physics. Having richer hopping processes both with and without the
restriction that the orbital flavor is conserved, one has to accept that
alternating orbital order is more difficult to realize than in the spin model
[11]. Feature (v.b) makes it more difficult to realize fully polarized ferro-
orbital states so it may sound surprising that nevertheless the FOr phase is
more stable than the orbital liquid in a range of electron densities $n\geq
n_{c}$. Note that the theorem formulated by Nagaoka [115] at $U=\infty$ for
the spin case does not apply to the orbital Hubbard model.

One of the attractive ideas in this field has been the possible existence of
polarized ferro-orbital states with partly filled complex orbitals [116]. We
have established that unfortunately such states cannot be realized. We have
shown that they are unstable and one has to consider instead FOr states.

###### Acknowledgements.

We thank Walter Metzner for his interest and for very insightful discussions.
L.F.F. thanks the Max Planck Institute for Solid State Research in Stuttgart
and the Jagiellonian University in Kraków for their kind hospitality. A.M.O.
kindly acknowledges Narodowe Centrum Nauki (NCN, Poland) Project No.
2016/23/B/ST3/00839 and is grateful for support via the Alexander von Humboldt
Foundation Fellowship [117] (Humboldt-Forschungspreis). *

##  Appendix A  Proof that the orbital liquid phase is unpolarized

Consider the cyclic permutation $\mathrm{C}$ defined by

$\displaystyle\mathrm{C}{\bf k}$ $\displaystyle=$
$\displaystyle\mathrm{C}(k_{-m},k_{-m+1},....,k_{0},....,k_{m})$ (129)
$\displaystyle\equiv(k_{m},k_{-m},....,k_{-1},....,k_{m-1}).$

Then we have, compare Eqs. (15)–(16),

$\displaystyle G_{{\mathrm{C}\bf k}}$ $\displaystyle=$ $\displaystyle
B_{\mathrm{C}{\bf k}}\,\mathrm{exp}(i\phi_{\mathrm{C}{\bf
k}})=\sum_{n}\cos(\mathrm{C}{\bf k})_{n}\,\mathrm{e}^{i\chi_{n}}$ (130)
$\displaystyle=$ $\displaystyle\sum_{n}\cos
k_{n-1}\,\mathrm{e}^{i\chi_{n}}=\sum_{n}\cos k_{n}\,\mathrm{e}^{i\chi_{n+1}}$
$\displaystyle=$ $\displaystyle\mathrm{e}^{i4\pi/d}\,G_{\bf k}=B_{\bf
k}\,\mathrm{e}^{i4\pi/d}\,\mathrm{exp}(i\phi_{\bf k}),$

and it follows that $B_{\mathrm{C}{\bf k}}=B_{\bf k}$ and
$\phi_{\mathrm{C}{\bf k}}=\phi_{\bf k}+4\pi/d$. The set of ${\bf k}$-vectors
generated from ${\bf k}$ by successive permutations is $\Pi({\bf
k})=\\{\mathrm{C}^{n}{\bf k}\,|-m\leq n\leq m\\}$, so

$\displaystyle\sum_{{\bf k}^{\prime}\in\Pi({\bf
k})}\\!\\!\\!\mathrm{e}^{i\phi_{{\bf k}^{\prime}}}$ $\displaystyle=$
$\displaystyle\sum_{n}\mathrm{exp}(i\phi_{\mathrm{C}^{n}{\bf
k}})=\sum_{n}\mathrm{exp}(i(\phi_{\bf k}+4n\pi/d))$ (131) $\displaystyle=$
$\displaystyle\mathrm{e}^{i\phi_{\bf
k}}\,\sum_{n}\mathrm{e}^{i4n\pi/d}=\mathrm{e}^{i\phi_{\bf
k}}\,\sum_{n}\mathrm{e}^{i\chi_{n}}=0\,.$

## References

  * [1] J. Hubbard, Electron Correlations in Narrow Energy Bands, Proceedings of the Royal Society of London A 276, 238 (1963).
  * [2] P. W. Anderson, New Approach to the Theory of Superexchange Interactions, Phys. Rev. 115, 2 (1959).
  * [3] W. Metzner and D. Vollhardt, Correlated Lattice Fermions in $d=\infty$ Dimension, Phys. Rev. Lett. 62, 324 (1989).
  * [4] W. Metzner and D. Vollhardt, Analytic calculation of ground-state properties of correlated fermions with the Gutzwiller wave-function, Phys. Rev. B 37, 7382 (1988).
  * [5] M. C. Gutzwiller, Effect of correlation on the ferromagnetism of transition metals, Phys. Rev. Lett. 10, 159 (1963).
  * [6] M. C. Gutzwiller, Correlation of Electrons in a Narrow $s$ Band, Phys. Rev. 137, A1726 (1965).
  * [7] D. Vollhardt, Normal He3 — An almost localized Fermi liquid, Rev. Mod. Phys. 56, 99 (1984).
  * [8] W. Metzner, Variational theory for correlated lattice fermions in high dimension, Z. Phys. B 77, 253 (1989).
  * [9] W. Metzner, Analytic evaluation of resonating valence bond states, Z. Phys. B 82, 183 (1991).
  * [10] P. W. Anderson, Resonating valence bonds: A new kind of insulator?, Mater. Res. Bull. 8, 153-160 (1973).
  * [11] P. Fazekas, B. Menge, and E. Müller-Hartmann, Ground state phase diagram of the infinite dimensional Hubbard model: A variational study, Z. Phys. B 78, 69 (1990).
  * [12] S. Ishihara, M. Yamanaka, and N. Nagaosa, Orbital liquid in perovskite transition-metal oxides, Phys. Rev. B 56, 686 (1997).
  * [13] L. F. Feiner and A. M. Oleś, Orbital liquid in ferromagnetic manganites: The orbital Hubbard model for $e_{g}$ electrons, Phys. Rev. B 71, 144422 (2005).
  * [14] T. Tanaka, M. Matsumoto, and S. Ishihara, Randomly diluted $e_{g}$ orbital-ordered systems, Phys. Rev. Lett. 95, 267204 (2005).
  * [15] T. Tanaka and S. Ishihara, Dilution effects in two-dimensional quantum orbital systems, Phys. Rev. Lett. 98, 256402 (2007).
  * [16] P. Czarnik, J. Dziarmaga, and A. M. Oleś, Overcoming the sign problem at finite temperature: Quantum tensor network for the orbital $e_{g}$ model on an infinite square lattice, Phys. Rev. B 96, 014420 (2017).
  * [17] L. Balents, Spin liquids in frustrated magnets, Nature (London) 464, 199-208 (2010).
  * [18] L. Savary and L. Balents, Quantum spin liquids: A review, Rep. Prog. Phys. 80, 016502 (2017).
  * [19] K. I. Kugel and D. I. Khomskii, The Jahn-Teller effect and magnetism: Transition metal compounds, Usp. Fiz. Nauk 136, 621 (1982) [Sov. Phys. Usp. 25, 231 (1982)].
  * [20] D. I. Khomskii, Transition Metal Oxides (Cambridge University Press, Cambridge, 2014).
  * [21] M. Imada, A. Fujimori, and Y. Tokura, Metal-insulator transitions, Rev. Mod. Phys. 70, 1039 (1998).
  * [22] Y. Tokura and N. Nagaosa, Orbital Physics in Transition-Metal Oxides, Science 288, 462-468 (2000).
  * [23] A. V. Gorshkov, M. Hermele, V. Gurarie, C. Xu, P. S. Julienne, J. Ye, P. Zoller, E. Demler, M. D. Lukin, and A. M. Rey, Two-orbital SU(N) magnetism with ultracold alkaline-earth atoms, Nature Phys. 6, 289 (2010).
  * [24] A. Koga, N. Kawakami, T. M. Rice, and M. Sigrist, Orbital-Selective Mott Transitions in the Degenerate Hubbard Model, Phys. Rev. Lett. 92, 216402 (2004).
  * [25] L. F. Feiner, A. M. Oleś, and J. Zaanen, Quantum Melting of Magnetic Order due to Orbital Fluctuations, Phys. Rev. Lett. 78, 2799 (1997).
  * [26] L. F. Feiner, A. M. Oleś, and J. Zaanen, Quantum disorder versus order-out-of-disorder in the Kugel-Khomskii model, J. Phys.: Condens. Matter 10, L555 (1998).
  * [27] G. Khaliullin and V. Oudovenko, Spin and orbital excitation spectrum in the Kugel-Khomskii model, Phys. Rev. B 56, R14243 (1997).
  * [28] P. Corboz, M. Lajkó, A. M. Laüchli, K. Penc, and F. Mila, Spin-Orbital Quantum Liquid on the Honeycomb Lattice, Phys. Rev. X 2, 041013 (2012).
  * [29] G. Khaliullin and S. Maekawa, Orbital liquid in three-dimensional Mott insulator: LaTiO3, Phys. Rev. Lett. 85, 3950 (2000).
  * [30] G. Khaliullin, Orbital order and fluctuations in Mott insulators, Prog. Theor. Phys. Suppl. 160, 155 (2005).
  * [31] K. Kitagawa, T. Takayama, Y. Matsumoto, A. Kato, R. Takano, Y. Kishimoto, S. Bette, R. Dinnebier, G. Jackeli, and H. Takagi, A spin-orbital-entangled quantum liquid on a honeycomb lattice, Nature (London) 554, 341 (2018).
  * [32] L. F. Feiner and A. M. Oleś, Electronic origin of magnetic and orbital ordering in insulating LaMnO3, Phys. Rev. B 59, 3295 (1999).
  * [33] J. van den Brink, P. Horsch, F. Mack, and A. M. Oleś, Orbital dynamics in ferromagnetic transition-metal oxides, Phys. Rev. B 59, 6795 (1999).
  * [34] A. M. Oleś, G. Khaliullin, P. Horsch, and L. F. Feiner, Fingerprints of spin-orbital physics in cubic Mott insulators: Magnetic exchange interactions and optical spectral weights, Phys. Rev. B 72, 214431 (2005).
  * [35] A. Reitsma, L. F. Feiner, and A. M. Oleś, Orbital and spin physics in LiNiO2 and NaNiO2, New J. Phys. 7, 121 (2005).
  * [36] J. Chaloupka and G. Khaliullin, Orbital order and possible superconductivity in LaNiO3/LaMO3 superlattices, Phys. Rev. Lett. 100, 016404 (2008).
  * [37] V. I. Solovyev, Spin-orbital superexchange physics emerging from interacting oxygen molecules in KO2, New J. Phys. 10, 013035 (2008).
  * [38] B. Normand and A. M. Oleś, Frustration and entanglement in the $t_{2g}$ spin-orbital model on a triangular lattice: valence-bond and generalized liquid states, Phys. Rev. B 78, 094427 (2008).
  * [39] F. Krüger, S. Kumar, J. Zaanen, and J. van den Brink, Spin-orbital frustrations and anomalous metallic state in iron-pnictide superconductors, Phys. Rev. B 79, 054504 (2009).
  * [40] W. Lv, F. Krüger, and P. Phillips, Orbital ordering and unfrustrated $(\pi,0)$ magnetism from degenerate double exchange in the iron pnictides, Phys. Rev. B 82, 045125 (2010).
  * [41] B. Normand, Multicolored quantum dimer models, resonating valence-bond states, color visons, and the triangular-lattice $t_{2g}$ spin-orbital system, Phys. Rev. B 83, 094427 (2011).
  * [42] J. Chaloupka and A. M. Oleś, Spin-orbital resonating valence-bond liquid on a triangular lattice: Evidence from finite cluster diagonalization, Phys. Rev. B 83, 094427 (2011).
  * [43] K. Wohlfeld, M. Daghofer, and A. M. Oleś, Spin-orbital physics for $p$ orbitals in alkali hyperoxides—Generalization of the Goodenough-Kanamori rules, EPL (Europhysics Letters) 96, 27001 (2011).
  * [44] K. Wohlfeld, M. Daghofer, S. Nishimoto, G. Khaliullin, and J. van den Brink, Intrinsic Coupling of Orbital Excitations to Spin Fluctuations in Mott Insulators, Phys. Rev. Lett. 107, 147201 (2011).
  * [45] J. Schlappa, K. Wohlfeld, K. J. Zhou, M. Mourigal, M. W. Haverkort, V.N. Strocov, L. Hozoi, C. Monney, S. Nishimoto, S. Singh, A. Revcolevschi, J.-S. Caux, L. Patthey, H. M. Ronnow, J. van den Brink, and T. Schmitt, Spin-orbital separation in the quasi-one-dimensional Mott insulator Sr2CuO3, Nature (London) 485, 82-U108 (2012).
  * [46] K. Wohlfeld, S. Nishimoto, M. W. Haverkort, and J. van den Brink, Microscopic origin of spin-orbital separation in Sr2CuO3, Phys. Rev. B 88, 195138 (2013).
  * [47] V. Bisogni, K. Wohlfeld, S. Nishimoto, C. Monney, J. Trinckauf, K. Zhou, R. Kraus, K. Koepernik, C. Sekar, V. Strocov, B. Büchner, T. Schmitt, J. van den Brink, and J. Geck, Orbital Control of Effective Dimensionality: From Spin-Orbital Fractionalization to Confinement in the Anisotropic Ladder System CaCu2O3, Phys. Rev. Lett. 114, 096402 (2015).
  * [48] C.-C. Chen, M. van Veenendaal, T. P. Devereaux, and K. Wohlfeld, Fractionalization, entanglement, and separation: Understanding the collective excitations in a spin-orbital chain, Phys. Rev. B 91, 165102 (2015).
  * [49] W. Brzezicki, A. M. Oleś, and M. Cuoco, Spin-Orbital Order Modified by Orbital Dilution in Transition-Metal Oxides: From Spin Defects to Frustrated Spins Polarizing Host Orbitals, Phys. Rev. X 5, 011037 (2015).
  * [50] W. Brzezicki, M. Cuoco, and A. M. Oleś, Novel Spin-Orbital Phases induced by Orbital Dilution, J. Supercond. Novel Magn. 29, 536 (2016).
  * [51] W. Brzezicki, M. Cuoco, and A. M. Oleś, Exotic Spin-Orbital Physics in Hybrid Oxides, J. Supercond. Novel Magn. 30, 129 (2017).
  * [52] J. B. Goodenough, Magnetism and the Chemical Bond (Interscience, New York, 1963).
  * [53] A. M. Oleś, P. Horsch, L. F. Feiner, and G. Khaliullin, Spin-Orbital Entanglement and Violation of the Goodenough-Kanamori Rules, Phys. Rev. Lett. 96, 147205 (2006).
  * [54] A. M. Oleś, Fingerprints of spin-orbital entanglement in transition metal oxides, J. Phys.: Condensed Matter 24, 313201 (2012).
  * [55] W. Brzezicki, J. Dziarmaga, and A. M. Oleś, Noncollinear Magnetic Order Stabilized by Entangled Spin-Orbital Fluctuations, Phys. Rev. Lett. 109, 237201 (2012).
  * [56] W. Brzezicki, J. Dziarmaga, and A. M. Oleś, Exotic spin orders driven by orbital fluctuations in the Kugel-Khomskii model, Phys. Rev. B 87, 064407 (2013).
  * [57] M. Snamina and A. M. Oleś, Magnon dressing by orbital excitations in ferromagnetic planes of K2CuF4 and LaMnO3, New J. Phys. 21, 023018 (2019).
  * [58] W. Brzezicki, Spin, orbital and topological order in models of strongly correlated electrons, J. Phys.: Condensed Matter 32, 023001 (2020).
  * [59] P. A. Lee, N. Nagaosa, and X.-G. Wen, Doping a Mott insulator: Physics of high-temperature superconductivity, Rev. Mod. Phys. 78, 17 (2006).
  * [60] M. Ogata and H. Fukuyama, The $t$-$J$ model for the oxide high-$T_{c}$ superconductors, Rep. Prog. Phys. 71, 036501 (2008).
  * [61] B. Keimer, S. A. Kivelson, M. R. Norman, S. Uchida, and J. Zaanen, High Temperature Superconductivity in the Cuprates, Nature 518, 179 (2015).
  * [62] A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg, Dynamical mean-field theory of strongly correlated fermion systems and the limit of infinite dimensions, Rev. Mod. Phys. 68, 13 (1996).
  * [63] J. Merino and R. H. McKenzie, Transport properties of strongly correlated metals: A dynamical mean-field approach, Phys. Rev. B 61, 7996 (2000).
  * [64] B. Kyung, S. S. Kancharla, D. Sénéchal, A.-M. S. Tremblay, M. Civelli, and G. Kotliar, Pseudogap induced by short-range spin correlations in a doped Mott insulator, Phys. Rev. B 73, 165114 (2006).
  * [65] G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko, O. Parcollet, and C. A. Marianetti, Electronic structure calculations with dynamical mean-field theory, Rev. Mod. Phys. 78, 865 (2006).
  * [66] D. Gotfryd, E. Pärschke, J. Chaloupka, A. M. Oleś, and K. Wohlfeld, How Spin-Orbital Entanglement depends on the Spin-Orbit Coupling in a Mott Insulator, Phys. Rev. Research 2, 013353 (2020).
  * [67] D. Gotfryd, E. Pärschke, K. Wohlfeld, and A. M. Oleś, Evolution of Spin-Orbital Entanglement with Increasing Ising Spin-Orbit Coupling, Condensed Matter 5, 53 (2020).
  * [68] A. M. Oleś, Antiferromagnetism and correlation of electrons in transition metals, Phys. Rev. B 28, 327 (1983).
  * [69] A. B. Harris, T. Yildirim, A. Aharony, Ora Entin-Wohlman, and I. Y. Korenblit, Unusual Symmetries in the Kugel-Khomskii Hamiltonian, Phys. Rev. Lett. 91, 087206 (2003).
  * [70] M. Daghofer, K. Wohlfeld, A. M. Oleś, E. Arrigoni, and P. Horsch, Absence of Hole Confinement in Transition-Metal Oxides with Orbital Degeneracy, Phys. Rev. Lett. 100, 066403 (2008).
  * [71] K. Wohlfeld, M. Daghofer, A. M. Oleś, and P. Horsch, Spectral properties of orbital polarons in Mott insulators, Phys. Rev. B 78, 214423 (2008).
  * [72] P. Wróbel and A. M. Oleś, Ferro-Orbitally Ordered Stripes in Systems with Alternating Orbital Order, Phys. Rev. Lett. 104, 206401 (2010).
  * [73] P. Wróbel, R. Eder, and A. M. Oleś, Optical conductivity due to orbital polarons in systems with orbital degeneracy, Phys. Rev. B 86, 064415 (2012).
  * [74] P. Horsch, J. Jaklič, and F. Mack, Optical conductivity of colossal-magnetoresistance compounds: Role of orbital degeneracy in the ferromagnetic phase, Phys. Rev. B 59, 6217 (1999).
  * [75] P. Horsch and F. Mack, Optical conductivity in doped manganites with planar $x^{2}-y^{2}$ orbital order, Phys. Rev. Lett. 82, 3160 (1999).
  * [76] C. Zener, Interaction Between the $d$ Shells in the Transition Metals, Phys. Rev. 81, 440 (1951).
  * [77] E. Dagotto, T. Hotta, and A. Moreo, Colossal magnetoresistant materials: The key role of phase separation, Phys. Rep. 344, 1 (2001).
  * [78] E. Dagotto, Open questions in CMR manganites, relevance of clustered states and analogies with other compounds including the cuprates, New J. Phys. 7, 67 (2005).
  * [79] T. V. Ramakrishnan, Modelling colossal magnetoresistance manganites, J. Phys.: Condens. Matter 19, 125211 (2007).
  * [80] Y. Tokura, Critical features of colossal magnetoresistive manganites, Rep. Prog. Phys. 69, 797 (2006).
  * [81] R. Kilian and G. Khaliullin, Orbital polarons in the metal-insulator transition of manganites, Phys. Rev. B 60, 13458 (1999).
  * [82] M. Daghofer, A. M. Oleś, and W. von der Linden, Orbital polarons versus itinerant $e_{g}$ electrons in doped manganites, Phys. Rev. B 70, 184430 (2004).
  * [83] G. Khaliullin and R. Kilian, Theory of anomalous magnon softening in ferromagnetic manganites, Phys. Rev. B 61, 3494 (2000).
  * [84] A. M. Oleś and L. F. Feiner, Why spin excitations in metallic ferromagnetic manganites are isotropic, Phys. Rev. B 65, 052414 (2002).
  * [85] The sign convention taken here is natural for the purely electronic problem, and is opposite to that in [13] which is commonly used in the Jahn-Teller problem for $e_{g}$ doublets, cf. R. Englman, The Jahn-Teller Effect in Molecules and Crystals (John Wiley, London, 1972).
  * [86] J. R. Klauder, Path integrals and stationary-phase approximations, Phys. Rev. D 19, 2349 (1979).
  * [87] Note that $\theta_{i}$ is the azimuthal, not the polar angle.
  * [88] N. C. Petroni, Probability and Stochastic Processes for Physicists (Springer, Berlin, 2021).
  * [89] J. Bünemann, W. Weber, and F. Gebhard, Multiband Gutzwiller wave functions for general on-site interactions, Phys. Rev. B 57, 6896 (1998).
  * [90] E. W. Ng and M. Geller, A Table of Integrals of the Error Functions, J. Res. NBS B Math. Sci. 73B, 1 (1969).
  * [91] W. F. Brinkman and T. M. Rice, Application of Gutzwiller’s variational method to the metal-insulator transition, Phys. Rev. B 2, 4302 (1970).
  * [92] This is an idealization and at finite dimension the condition $d=0$ is satisfied only at $U=\infty$, because the Brinkman-Rice criterion neglects virtual excitations, which generate antiferromagnetism at finite $U$ in the spin Hubbard model.
  * [93] P. G. J. van Dongen, F. Gebhard, and D. Vollhardt, Variational evaluation of correlation functions for lattice electrons in high dimensions, Z. Phys. B 76, 199-210 (1989).
  * [94] D. Vollhardt, Dynamical Mean-Field Approach for Strongly Correlated Materials, in: Lecture Notes of the Autumn School on Correlated Electrons 2011, Vol. 1, Eds. E. Pavarini, E. Koch, D. Vollhardt, and A. Lichtenstein (Verlag des Forschungszentrum Jülich, 2011).