Rotationally-invariant slave-boson formalism and momentum dependence of
  the quasiparticle weight
arXiv:0704.1434v2  [cond-mat.str-el]  2 Aug 2007
Rotationally-invariant slave-boson formalism
and momentum dependence of the quasiparticle weight
Frank Lechermann,1, 2, ∗ Antoine Georges,2 Gabriel Kotliar,3, 2 and Olivier Parcollet4
1I. Institut für Theoretische Physik, Universität Hamburg,Jungiusstrasse 9,20355 Hamburg, Germany
2Centre de Physique Théorique, École Polytechnique, 91128 Palaiseau Cedex, France
3Serin Physics Laboratories, Rutgers University, Piscataway, NJ, USA
4Service de Physique Théorique, CEA/DSM/SPhT-CNRS/SPM/URA 2306 CEA Saclay, F-91191 Gif-Sur-Yvette, France
We generalize the rotationally-invariant formulation of the slave-boson formalism to multiorbital
models, with arbitrary interactions, crystal fields, and multiplet structure. This allows for the study
of multiplet effects on the nature of low-energy quasiparticles. Non-diagonal components of the
matrix of quasiparticle weights can be calculated within this framework. When combined with
cluster extensions of dynamical mean-field theory, this method allows us to address the effects of
spatial correlations, such as the generation of the superexchange and the momentum dependence
of the quasiparticle weight. We illustrate the method on a two-band Hubbard model, a Hubbard
model made of two coupled layers, and a two-dimensional single-band Hubbard model (within a
two-site cellular dynamical mean-field approximation).
PACS numbers: 71.10.-w,71.10.Fd,71.30.+h,74.25.Jb
I. INTRODUCTION AND MOTIVATIONS
A. General motivations
The method of introducing auxiliary bosons in order
to facilitate the description of interacting fermionic sys-
tems is an important technique in theoretical many-body
physics. In this regard, the so-called slave boson (SB)
approach is a very useful tool in dealing with models of
strongly correlated electrons. Slave boson mean-field the-
ory (SBMFT), i.e., at the saddle-point level, is the sim-
plest possible realization of a Landau Fermi liquid (for
a review, see e.g. Ref.1). Within SBMFT, a simplified
description of the low-energy quasiparticles is obtained,
while high-energy (incoherent) excitations are associated
with fluctuations around the saddle point. In particu-
lar, two essential features are captured by SBMFT: (i)
the Fermi surface (FS) of the interacting system (sat-
isfying Luttinger’s theorem) is determined by the zero-
frequency self-energy, which is in turn determined by the
Lagrange multipliers associated with the constraints and
(ii) the quasiparticle (QP) weight Z is determined by the
saddle-point values of the slave bosons. Hence SBMFT is
a well-tailored technique when attempting to understand
the low-energy physics emerging from more sophisticated
theoretical tools, such as for example dynamical mean-
field theory (DMFT), which deals with the full frequency
dependence of the self-energy.
In this paper, we are concerned with the construction
of a slave-boson formalism which is able to deal with the
two following problems.
1) In multi-orbital models, handle an arbitrary form
of the interaction hamiltonian, not restricted to density-
density terms, and possibly including interorbital hop-
pings or hybridizations. We aim in particular at describ-
ing correctly the multiplets (eigenstates of the atomic
hamiltonian), but we also want to be able to work in
an arbitrary basis set, not necessarily that of the atomic
multiplets (and of course, to obtain identical results, in-
dependent of the choice of basis).
2) Describe situations in which the QP weight is not
uniform along the Fermi surface, but instead varies as a
function of the momentum, i.e., Z=Z(k).
There are clear physical motivations for addressing
each of these issues. The first one is encountered when-
ever one wants to deal with a specific correlated material
in a realistic setting (see e.g. Ref. 2). Usually, more than
one band is relevant to the physics (e.g. a t2g triplet
or eg doublet for transition metal-oxides, or the full 7-
fold set of f -orbitals in rare-earth, actinides and their
compounds). The second issue is an outstanding one
in connection with cuprate superconductors. In those
materials, a strong differentiation in momentum space
is observed in the “normal” (i.e., non-superconducting)
state, especially in the underdoped regime (for a review,
see e.g. Ref. 3). For momenta close to the nodal re-
gions, i.e., close to the regions where the superconduct-
ing gap vanishes, reasonably long-lived QPs are found. In
contrast, in the antinodal directions, the angle-resolved
photoemission (ARPES) spectra reveal only a broad line-
shape with no well-defined QPs. The nature of the in-
cipient normal state in the underdoped regime (i.e., the
state achieved by suppressing the intervening supercon-
ductivity) has been a subject of debate. One possibil-
ity is that QPs would eventually emerge at low-enough
temperature in the antinodal region as well, but with a
much smaller QP weight ZAN≪ZN. Another possibility
is that coherent QPs simply do not emerge in the antin-
odal region. Anyhow, there is evidence from ARPES
and other experiments3 that the QP weight (whenever
it can be defined) has significant variation along the FS
and is larger at the nodes. Since the QP weight sets the
scale for the coherence temperature below which long-
lived QPs form, a smaller Z means a smaller coher-
ence temperature. Hence, if the temperature is higher
http://arxiv.org/abs/0704.1434v2
than the coherence scale associated with momenta close
to the nodes, and larger than the one associated with
the antinodes, QPs will be visible only in the nodal re-
gions. At this temperature, the FS will thus appear as
being formed of “Fermi arcs”, as indeed observed experi-
mentally4. Important differences between the nodal and
antinodal region in the superconducting state have also
been unraveled by recent experiments, in particular from
Raman scattering which revealed two different energy
scales with different doping dependence, associated with
each of these regions5. Momentum-space differentiation
of QP properties is therefore a key feature of cuprate su-
perconductors, but it is also an issue which is particularly
difficult to handle theoretically.
As we now explain, these two issues are actually closely
related one to the other. In a general multi-orbital model,
the self-energy is a matrix Σαβ (α, β are orbital indices).
Except when a particular symmetry dictates otherwise,
this matrix has in general off-diagonal (interorbital) com-
ponents and these off-diagonal components may have
a non-zero linear term in the low-frequency expansion,
hence yielding non-diagonal components of the matrix of
QP weights defined as:
. (1)
On the other hand, a momentum dependent QP weight
Z(k) means that, in real-space, Zij=Z(Ri−Rj) depends
on the separation between lattice sites (a momentum-
independent Z means that Zij=Z δij is purely local).
Hence, in both cases, one has to handle a QP weight
which is a matrix in either the orbital or the site in-
dices. The connection becomes very direct in the frame-
work of cluster extensions of DMFT (for reviews, see e.g.
Refs [6–9]). There, a lattice problem is mapped onto a
finite-size cluster which is self-consistently coupled to an
environment. This finite-size cluster can be viewed as a
multi-orbital (or molecular) quantum impurity problem,
in which each site plays the role of an atomic orbital. Re-
cently, numerical solutions of various forms of cluster ex-
tensions to the DMFT equations for the two-dimensional
Hubbard model have clearly revealed the phenomenon
of momentum-space differentiation10–12. Developing low-
energy analytical tools to interpret, understand, and gen-
eralize the results of these calculations is clearly an im-
portant and timely issue. The slave-boson methods de-
veloped in the present work are a step in this direction.
Obviously, the existence of off-diagonal components of
the Ẑ-matrix is a basis-set dependent issue. A proper
choice of orbital basis can be made, which diagonalizes
this matrix. In certain cases, this basis is dictated by
symmetry considerations, while in the absence of sym-
metries, the basis set in which Z is diagonal cannot be
guessed a priori. For instance, in a two-site cluster or
two-orbital model in which the two sites play equivalent
roles, even and odd combinations diagonalize not only
the Ẑ-matrix, but in fact the self-energy matrix itself
for all frequencies (see Sec. III). In such cases, it may
be favorable to work in this orbital basis set, and deal
only with diagonal QP weights. However, performing the
rotation into this orbital basis set will in general trans-
form the interacting hamiltonian into a more complicated
form. For example, starting from a density-density inter-
action, it may induce interaction terms which are not of
the density-density type (i.e., involve exchange, pair hop-
ping, etc.). For these reasons, it is essential to consider
slave-boson formalisms which can handle both arbitrary
interaction terms, and non-diagonal components of the
QP weight matrix: these two issues are indeed connected.
The formalism presented in this article builds on earlier
ideas of Li, Wölfle, Hirschfeld and Frésard13,14 (see Ap-
pendix. A), in which the SB formalism is formulated in a
fully rotationally-invariant manner (see also Refs [15,16]
in the framework of the Gutzwiller approximation), so
that the orbital basis set needs not be specified from
the beginning, and the final results are guaranteed to
be equivalent irrespectively of the chosen basis set.
B. Some notations
In this paper, we shall consider multi-orbital models of
correlated electrons with hamiltonians of the form:
H = Hkin +
Hloc[i] , (2)
with Hkin =
εαβ(k) d
kαdkβ . (3)
In these expressions, α, β label electronic species and run
from 1 to M (i.e., M is twice the number of atomic or-
bitals in the context of a multi-orbital model of electrons
with spin: α=(m,σ),σ=↑, ↓). The k vector runs over the
Brillouin zone of the lattice, whose sites are labeled by
i (in the context of cluster-DMFT, i will label clusters
and runs over the superlattice sites, thus k runs over the
reduced Brillouin zone of the superlattice, see Sec. III C).
The first term in (2) is the kinetic energy: εαβ(k) is the
Fourier transform of the (possibly off-diagonal) hoppings
and does not contain any local terms (i.e.,
εαβ(k)=0).
Hloc contains both the one-body local terms and the in-
teractions, assumed to be local. A general form for Hloc
is55:
Hloc =
ε0αβd
αdβ +
Uαβγδ d
βdδdγ . (4)
Fock states form a convenient basis set of the local
Hilbert space on each site. They are specified by se-
quences n=(n1, . . . , nM ), with nα=0, 1 (we consider a
single site and drop the site index):
|n〉 =
· · ·
|vac〉 . (5)
In the following, {|A〉} will denote an arbitrary basis set
of the local Hilbert space, specified by its components on
the Fock states:
|A〉 =
〈n|A〉 |n〉 , (6)
while |Γ〉 will denote the eigenstates of the local hamil-
tonian, i.e., the ‘atomic’ multiplets such that:
Hloc|Γ〉 = EΓ |Γ〉 . (7)
C. Slave bosons for density-density interactions: a
reminder
When the orbital densities nα are good quantum num-
bers for the local hamiltonian Hloc, i.e., when the eigen-
states of the latter are labeled by nα, a very simple slave
boson formalism can be constructed which is a direct
multi-orbital generalization of the 4-boson scheme intro-
duced by Kotliar and Ruckenstein17,18. While this is
standard material, we feel appropriate to briefly remind
the reader of how this scheme works, in order to con-
sider generalizations later on. We thus specialize in this
subsection to a local hamiltonian of the form:
Hloc =
ε0αn̂α +
Uαβ n̂α n̂β , (8)
so that the multiplets are the Fock states |n〉 themselves,
with eigenenergies:
ε0αnα +
Uαβ nα nβ . (9)
To each Fock state, one associates a boson creation op-
erator φ†n. Furthermore, auxiliary fermions f
α are intro-
duced which correspond to quasiparticle degrees of free-
dom. The (local) enlarged Hilbert space thus consists
of states which are built from tensor products of a QP
Fock state, times an arbitrary number of bosons. In con-
trast, the physical Hilbert space is generated by the basis
set consisting of the 2M states which contain exactly one
boson, and in which this boson matches the QP Fock
state. Thus, the states representing the original physical
states (5) in the enlarged Hilbert space, in a one-to-one
manner, are the following (“physical”) states:
|n 〉 ≡ φ†n|vac〉 ⊗ |n〉f . (10)
The underlining in |n 〉 allows one to distinguish between
the original Fock state of the physical electrons |n〉, and
its representative state in the enlarged Hilbert space. In
this expression, |n〉f stands for the QP Fock state:
|n〉f ≡
· · ·
|vac〉 . (11)
It is easily checked that a simple set of constraints
uniquely specifies the physical states among all the states
of the enlarged Hilbert space, namely:
φ†nφn = 1 (12)
nφn = f
αfα , ∀α . (13)
The first constraint imposes that only states with a single
boson are retained, while the second one insures that the
fermionic (QP) and bosonic contents match. Obviously,
the saddle-point values of the slave bosons will have a
simple interpretation, |φn|2 being the probability associ-
ated with the Fock space configuration n.
The operator:
d†α =
〈n|f †α|m〉φ†nφm f †α (14)
is a faithful representation of the physical electron cre-
ation operator on the representatives (10), namely:
d†α|n 〉 =
〈n′|d †a |n〉 |n′ 〉 , (15)
in which, in fact, the r.h.s is either zero (if nα= 1) or
composed of just a single state (with n′α= 1 and otherwise
β for β 6=α). This expression of the physical electron
operators is not unique however: obviously, one can for
example multiply this with any operator acting as the
identity on the physical states. This is true as long as the
constraint is treated exactly. When treated in the mean-
field approximation however, (i.e., at saddle point), these
equivalent expressions will not lead to the same results.
In fact, (14) suffers from a serious drawback namely it
does not yield the exact non-interacting (Uαβ=0) limit
at saddle-point. Instead, the expression:
〈n|f †α|m〉 [∆̂α]−
2φ†nφm [1− ∆̂α]−
2 f †α (16)
with ∆̂α[φ] ≡
nφn , (17)
turns out to satisfy this requirement, while having ex-
actly the same action as (14) when acting on physical
states. This choice of normalization is actually very natu-
ral given the probabilistic interpretation of |φn|2: the ex-
pression [∆̂α]
−1/2φ†n (resp. φm [1−∆̂α]−1/2) is actually a
probability amplitude, normalized over the restricted set
of physical states such that nα=1 (resp. nα=0). Hence
the combination of boson fields in (16) is a transition
probability between the statem withmα=0 and the state
n with nα=1.
Anyhow, whether the simplest expression (14) or the
normalized expression (16) is chosen for the physical op-
erator, the relation between the physical and QP single-
particle operators is of the form:
dα = r̂α[φ] fα (18)
It is important to note that the orbital index carried by
the physical operator is identical to that of the QP oper-
ator. An immediate consequence is that the self-energy
at the saddle-point level is a diagonal matrix in orbital
space Σαβ=δαβΣα, which reads:
Σα(ω) = Σα(0) + ω
, (19)
with Zα = |rα|2 (20)
Σα(0) = λα/|rα|2 − ε0α . (21)
In these expressions, rα is evaluated at saddle-point level,
and λα is the saddle-point value of the Lagrange multi-
pliers enforcing the constraint (13).
The expression (20) of the QP weight is an immediate
consequence of (18): at saddle-point level, rα becomes
a c-number and (18) implies that the physical electron
carries a spectral weight |rα|2. Hence, in order to describe
within SBMFT situations in which the QP weight is a
non-diagonal matrix, one must disentangle the orbital
indices carried by the physical electron and those carried
by the QP degrees of freedom. These operators will then
be related by a non-diagonal matrix:
dα = R̂αβ [φ] fβ . (22)
This is precisely what the formalism exposed in this ar-
ticle achieves. The physical significance of such a non-
diagonal relation is that creating a physical electron in
a given orbital may induce the creation of QPs in any
other orbital. Thinking of orbital as real-space indices
(within e.g. cluster-DMFT), this means that the cre-
ation of a physical electron on a given site induces QPs
on other sites in a non-local manner, corresponding to a
momentum-dependent Z(k).
D. Difficulties with naive generalizations to the
multi-orbital case
Let us come back to the general multi-orbital inter-
action (4). In order to motivate the fully rotationally-
invariant formalism exposed in the next section, let us
point out some difficulties arising when attempting to
generalize the simple SB formalism of the previous sec-
tion.
The central difference between the general interaction
(4) and the density-density form (8) is that the atomic
multiplets |Γ〉 are no longer Fock states. Thus, it would
seem natural to associate a slave boson φΓ to each of the
atomic multiplets. Indeed, Bünemann et al.19 (see also
[15]) have proposed generalized Gutzwiller wave func-
tions in which a variational parameter (a.k.a a probabil-
ity |φΓ|2) is associated with each atomic multiplet (see
also Ref. [20] and the recent work of Dai et al.21 in the
SB context). A slave-boson formulation requires a clear
identification of the physical states within the enlarged
Hilbert space. A natural idea is to define those in one-
to-one correspondence with the atomic multiplets, as:
|Γ〉 ?= φ†Γ|vac〉 ⊗
〈n|Γ〉 |n〉f . (23)
The local part of the hamiltonian has a simple repre-
sentation on these physical states Hloc=
However, a major difficulty is that there is no simple
constraint implementing the restriction to these physi-
cal states, and such that it is quadratic in the fermionic
(QP) degrees of freedom (which is essential in order to
yield a manageable saddle point). In particular, it is eas-
ily checked that the apparently natural constraint21:
f †αfα
〈Γ|n̂α|Γ〉φ†ΓφΓ (24)
is actually not satisfied by the states (23) as an operator
identity56. Further difficulties also arise when attempt-
ing to derive an expression for the physical creation op-
erators. These difficulties stem from the fact that two
atomic multiplets having particle numbers differing by
one unit cannot in general be related by the action of a
single-fermion creation.
One might also think of defining the physical states in
correspondence to the Fock states, as:
|n〉 ?= |n〉f ⊗
〈Γ|n〉φ†Γ|vac〉 (25)
which do satisfy the following quadratic constraint:
f †αfα =
〈Γ|n̂α|Γ′〉φ†ΓφΓ′ . (26)
However, another difficulty then arises. Namely, it is
not possible to write the local interaction hamiltonian
purely in terms of bosonic degrees of freedom, which is
the whole purpose of SB representations. In particular,
the obvious expression Hloc=
ΓφΓ which had the
correct action on states (23) no longer works for states
(25) since it leaves unchanged the fermionic content of
them.
After some thinking, one actually realizes that these
naive generalizations are all faced with the same problem,
namely that they do not embody the crucial conceptual
distinction between physical and QP degrees of freedom.
Both (23) and (25) assume a priori a definite relation be-
tween the physical and QP content of a state. The key
to a successful SB formalism is therefore to disentangle
physical and quasiparticle degrees of freedom, and letting
the variational principle at saddle point decide which re-
lationship actually exists between those.
We shall see however in Sec. III A that, provided the lo-
cal hamiltonian has enough symmetries, the rotationally-
invariant formalism of the present article does correspond
to assigning at saddle point a probability to each atomic
configuration (multiplet) |Γ〉, hence establishing contact
with the previous works of Refs. [19,21]. Yet for less sym-
metric hamiltonians, the general formalism of the present
article is requested.
II. ROTATIONALLY-INVARIANT
SLAVE-BOSON FORMALISM
A. Physical Hilbert space and constraints
In order to construct a SB formalism in which physi-
cal and QP states are disentangled, we shall associate a
slave boson φΓn to each pair of atomic multiplet |Γ〉 and
QP Fock state |n〉f . More generally, we can work in an
arbitrary basis set |A〉 of the local Hilbert space, not nec-
essarily that of the atomic multiplets, and consider slave
bosons φAn. As we shall see, the formalism introduced
in this article is such that two different choices of basis
sets are related by a unitary transformation and there-
fore lead to identical results. In particular, one could also
choose the physical Fock states |m〉d as the basis set A,
and work with slave bosons φmn which form the com-
ponents of a density matrix connecting the physical and
QP spaces. It is crucial however to keep in mind that the
first index (A) refers to physical-electron states, while the
second one (n) refers to quasiparticles.
A priori, a slave boson φAn can be introduced for any
pair (A, n). However, in this paper, we shall restrict our-
selves to phases which do not display an off-diagonal su-
perconducting long-range order, and hence one can re-
strict the φAn’s to pairs of states which have the same
total particle number on a given site (the local hamilto-
nian Hloc commutes with
αdα). The formalism is
easily extended to superconducting states14,15,22 by lift-
ing this assumption and modifying appropriately the ex-
pressions derived in this section. In the following, we
consider basis states A which are eigenstates of the local
particle number (denoted by NA), and hence a φAn is
introduced provided
α nα=NA.
The representation of such a basis state in the enlarged
Hilbert space is defined as:
|A〉 ≡ 1√
An|vac〉 ⊗ |n〉f . (27)
In this expression, DA denotes the dimension of the sub-
space of the Hilbert space with particle number identical
to that of A, i.e., DA≡ D(NA)=
. This insures a
proper normalization of the state. As before, the “under-
line” in |A〉 allows to distinguish this state, which lives in
the tensor product Hilbert space of QP and boson states,
from the physical electron state |A〉.
Having decided on the physical states, we need to iden-
tify a set of constraints which select these physical states
out of the enlarged Hilbert space in a necessary and suf-
ficient manner. It turns out that the following (M2 + 1)
constraints achieve this goal:
AnφAn = 1 (28)
An′φAn 〈n|f
αfα′ |n′〉 = f †α fα′ , ∀α .(29)
The first constraint is obvious and requires that the
physical states are single-boson states. It is easy to check
that the physical states satisfy the second set of con-
straints (29), but a little more subtle to actually prove
that this set of constraints is sufficient to uniquely select
the physical states (27) in the enlarged Hilbert space.
The detailed proof is given in Appendix B. Let us em-
phasize that the order of primed and unprimed indices
in (29) is of central importance.
B. Representation of the physical electron
operators
We now turn to the representation of the physical elec-
tron creation operator on the representatives (27) of the
physical states in the enlarged Hilbert space. We need to
find an operator which acts on these representatives ex-
actly as d †α acts on the physical basis |A〉. Namely, given
the matrix elements 〈A|d †α|B〉 such that
d †α |B〉 =
〈A|d †α|B〉 |A〉 , (30)
we want to find an operator d†α (in terms of the boson
and QP operators) such that
d†α |B〉 =
〈A|d †α|B〉 |A〉 . (31)
1. Proximate expression
As in the case of the density-density interactions dis-
cussed above (Sec. IC), the answer is not unique. We
first construct the generalization of expression (14) to
the present formalism (i.e., ignore at first the question of
the proper operators to be inserted in order to recover the
correct non-interacting limit). The following expression
is shown in Appendix C to satisfy (31):
β,AB,nm
〈A|d †α|B〉〈n|f
β |m〉√
NA(M −NB)
AnφBm f
β . (32)
We note that NA=NB+1 in this expression can take the
values 1, . . . ,M .
Hence, we see that within this formalism, the physical
and QP operators are indeed related by a non-diagonal
transformation (22):
dα = R̂[φ]αβ fβ (33)
with the R̂-matrix corresponding to (32) given by (R̂∗αβ
denotes the complex conjugate of R̂αβ):
R̂[φ]∗αβ =
AB,nm
〈A|d †α|B〉〈n|f
β |m〉√
NA(M −NB)
AnφBm . (34)
The action of (32) on physical states, and the proof that
it satisfies (31) are detailed in Appendix C.
2. Improved expression
The simple expression (32), although having the cor-
rect action on the physical states, suffers from the same
drawback than (14) in the case of density-density inter-
actions. Namely, at saddle-point level (i.e., with the con-
straint satisfied on average instead of exactly), the non-
interacting limit is not appropriately recovered. Thus,
one needs to generalize the improved expression (16) to
the present rotationally-invariant formalism. However,
care must be taken to do so in a way which respects
gauge invariance (i.e., the possibility of making an ar-
bitrary unitary rotation on the QP orbital indices, see
Sec. II C).
We consider the following operators, bilinear in the
bosonic fields:
AnφAm〈m|f
αfβ |n〉 (35)
AnφAm〈m|fβf
α|n〉 , (36)
which can be interpreted as particle- and hole- like QP
density matrices (note that when the constraint is satis-
fied exactly: ∆̂
αβ=δαβ−∆̂
αβ). We then choose to modify
the R-matrix in the following manner (see Appendix C):
R̂[φ]∗αβ =
AB,nm,γ
〈A|d †α|B〉〈n|f †γ |m〉φ
AnφBmMγβ , (37)
with Mγβ ≡
∣∣∣∣∣
∆̂(p)∆̂(h) + ∆̂(h)∆̂(p)
)]− 12
∣∣∣∣∣
.(38)
We chose to let the QP density matrices enter the M -
matrix in a symmetrized way in order to respect equiva-
lent treatment of particles and holes. Expression (37) can
be shown to be gauge-invariant, and turns out to yield
the correct non-interacting limit at saddle point. How-
ever, although it yields a saddle point satisfying all the
appropriate physical requirements, it is not fully justified
as an operator identity.
C. Gauge invariance
As usual in formalisms using slave particles, a gauge
symmetry is present which allows one to freely rotate the
QP orbital indices, independently on each lattice site.
Physical observables are of course gauge-invariant. Let
us consider an arbitrary SU(M) rotation of the QP op-
erators:
f †α =
Uαβ f̃
β (39)
This rotation induces a corresponding unitary transfor-
mation U(U) of the QP Fock states |n〉f . This unitary
transformation is characterized by the fact that the ex-
pectation value of fα in its Fock basis is an invariant
tensor: it is the same in every basis. Therefore (summa-
tion over repeated indices is implicit everywhere in the
following):
〈n|f †β|m〉 = Uββ′U(U)
nn′〈n′|f
β′ |m
′〉 U(U)mm′ (40)
n|f †αfβ|m
= Uαα′U
ββ′U∗(U)nn′
n′|f †α′fβ′ |m
U(U)mm′
(the second expression can be deduced from the first us-
ing closure relations). We can now check that if the slave
bosons transforms like
φAn = U(U)nn′ φ̃An′ (42)
then the constraints and the expressions of the phys-
ical electron operator (either (32) or (37)) are gauge-
invariant. Namely, the R-matrix obeys the following
transformation law:
R̂[φ]αβ = R̂[φ̃]αβ′ Uββ′ (43)
and therefore the physical electron operator is invariant:
dα = R̂[φ̃]αβ f̃β = R̂[φ]αβfβ (44)
D. Change of physical and quasiparticle basis sets
It is clear that the basis |A〉 of the local Hilbert space
(i.e., the physical basis states) can be chosen arbitrarily
in this formalism. Indeed, making a basis change from
|A〉 to |Ã〉, all the expressions above keep an identical
form provided the bosons corresponding to the new basis
are defined as:
〈A|Ã〉φ†An (45)
As mentioned above, it is often convenient to use the
eigenstates |Γ〉 of Hloc as a basis set.
Changing the basis states associated with quasiparti-
cles is a somewhat trickier issue. Up to now, we have
worked with Fock states |n〉f . A different basis set |Q〉f
can be used, provided however the unitary matrix 〈Q|n〉
is real, i.e., 〈Q|n〉=〈n|Q〉. Indeed, the matrix element
〈Q|n〉 appears in the transformation of the physical states
and of the constraint, while 〈n|Q〉 appears in the trans-
formation of the physical electron operator. When this
matrix elements are real, new bosons can be defined in
the transformed QP basis according to:
〈Q|n〉φ†An , (〈Q|n〉 = 〈n|Q〉) (46)
In particular, when the local hamiltonian is a real sym-
metric matrix, the same linear combinations of Fock
states which define the atomic multiplets |Γ〉 can be used
for QPs, and bosons φΓΓ′ can be considered. This is
sometimes a useful way of interpreting the formalism and
the results at saddle point (see Sec. III A).
E. Expression of the hamiltonian, free energy and
Green’s function
In this section, we derive the expression of the hamil-
tonian in terms of the slave boson and QP fermionic vari-
ables. We then construct the free-energy functional to be
minimized within a mean-field treatment, and express the
Green’s function and self-energy at saddle point.
We recall that the full hamiltonian (2) reads, in terms
of the physical electron variables: H=Hkin +
iHloc[i]
with Hkin=
αβ εαβ(k) d
kαdkβ the intersite kinetic
energy and Hloc the local part of the hamiltonian on a
given site i, with general form (4).
It is easily checked that the following bosonic operator
is a faithful representation of Hloc on the representatives
of the physical states in the enlarged Hilbert space:
H loc =
〈A|Hloc|B〉
AnφBn (47)
If the basis |Γ〉 of atomic multiplets is used, this simplifies
down to:
H loc =
ΓnφΓn (48)
Using the bosonic R-operators relating the physical elec-
tron to the QP operators, yields the following expression
of the kinetic energy:
H kin =
αα′ββ′
[R̂†]αα′εα′β′(k)R̂β′β f
kαfkβ (49)
A mean-field theory is obtained by condensing the slave
bosons into c-numbers 〈φAn〉 ≡ ϕAn. The constraints
are implemented by introducing Lagrange multipliers: λ0
associated with (28) and λαα′ ≡ [Λ]αα′ associated with
(29). The saddle point is obtained by extremalizing, over
the ϕAn’s and the Lagrange multipliers, the following
free-energy functional:
Ω[{ϕAn}; Λ, λ0] = (50)
= − 1
tr ln
1 + e−β(R
†(ϕ)ε(k)R(ϕ)+Λ)
ABnn′
ϕ∗An′
δnn′δAB λ0 + δnn′〈A|Hloc|B〉
Λαβ〈n|f †αfβ |n′〉
ϕBn (51)
The saddle-point equations, as well as technical aspects
of their numerical solution, are detailed in Appendix D.
Finally, we derive the expressions of the Green’s func-
tions Ĝ, the self-energy Σ̂ and the QP weight Ẑ at sad-
dle point. For the QPs, the one-particle Green’s func-
tion Gf,αβ(k, τ − τ ′)≡−〈f †kα(τ)fkβ(τ ′)〉 reads (in matrix
form):
f (k, ω) = ω −R
†(ϕ) ε(k)R(ϕ) −Λ (52)
and hence the physical electron Green’s function reads
(we drop the ϕ dependence for convenience):
d (k, ω) = [R
†]−1G−1f R
= ω (RR†)−1 − [R†]−1ΛR−1 − ε(k) , (53)
while the non-interacting Green’s function is (including
the one-body term present in Hloc):
d0 (k, ω) = ω11− ε
0 − ε(k) . (54)
The physical self-energy is thus:
Σd(ω) ≡ G−1d0 −G
1− [RR†]−1
+ [R†]−1ΛR−1 − ε0 . (55)
So that the matrix of QP weights is obtained in terms of
the R̂-matrix at saddle point as:
Z = RR† . (56)
This generalizes (20) to non-diagonal cases. It is easily
checked that these expressions of the physical quantities
Gd,Σd and Z are indeed gauge-invariant.
III. ILLUSTRATIVE RESULTS
In the following, we apply the above formalism to three
different model problems in strongly correlated physics.
First, we consider two popular models, namely the two-
band Hubbard model on a three-dimensional (3D) cubic
lattice, and a “bi-layer” model, coupling two Hubbard 3D
cubic lattices . Finally a two-site cluster (cluster-DMFT)
approximation to the single-band Hubbard model on
a two-dimensional (2D) square lattice is investigated.
Hence these models have in common that they all involve
two coupled orbitals (associated, in the cluster-DMFT
(CDMFT) framework, to the dimer made of two lattice
sites). The present formalism is of course not restricted
to two-orbital problems, however such models provide the
simplest examples where the power of the method may
be demonstrated.
A. Two-band Hubbard model
The Hubbard model involving two correlated bands,
without further onsite hybridization or crystal-field split-
ting, serves as one of the standard problems in condensed
matter theory. In contrast to the traditional single-band
model, the formal interaction term in Eq. (4) now gener-
ates in the most general fully SU(2) symmetric case four
energy parameters, i.e., the intraorbital Hubbard U and
the interorbital Hubbard U ′ as well as the two exchange
couplings J ,JC . Thus the present atomic hamiltonian
TABLE I: Eigenstates |Γ〉 of the SU(2) rotationally-invariant two-band Hubbard model. Spin values and energies are given for
the eigenstates. The last column shows the slave bosons for the description of the eigenstates in the SBMFT formalism.
No. Eigenstate |Γ〉 SΓ SzΓ EΓ φΓn
1 |00, 00〉 0 0 0 φ1,|00,00〉
2 | ↑ 0, 00〉 1
0 φ2,|↑0,00〉
3 |0 ↓, 00〉 1
0 φ3,|0↓,00〉
4 |00, ↑ 0〉 1
0 φ4,|00,↑0〉
5 |00, 0 ↓〉 1
0 φ5,|00,0↓〉
6 | ↑ 0, ↑ 0〉 1 1 U ′ − J φ6,|↑0,↑0〉
(| ↑ 0, 0 ↓〉+ |0 ↓, ↑ 0〉) 1 0 U ′ − J
φ7,|↑0,0↓〉, φ7,|0↓,↑0〉
8 |0 ↓, 0 ↓〉 1 -1 U ′ − J φ8,|0↓,0↓〉
(| ↑ 0, 0 ↓〉 − |0 ↓, ↑ 0〉) 0 0 U ′ + J
φ9,|↑0,0↓〉, φ9,|0↓,↑0〉
10 1√
(| ↑↓, 00〉 − |00, ↑↓〉) 0 0 U − JC
φ10,|↑↓,00〉, φ10,|00,↑↓〉
11 1√
(| ↑↓, 00〉+ |00, ↑↓〉) 0 0 U + JC
11,|↑↓,00〉, φ11,|00,↑↓〉
12 | ↑↓, ↑ 0〉 1
U + 2U ′ − J φ12,|↑↓,↑0〉
13 | ↑↓, 0 ↓〉 1
U + 2U ′ − J φ13,|↑↓,0↓〉
14 | ↑ 0, ↑↓〉 1
U + 2U ′ − J φ
14,|↑0,↑↓〉
15 |0 ↓, ↑↓〉 1
U + 2U ′ − J φ
15,|0↓,↑↓〉
16 | ↑↓, ↑↓〉 0 0 2U + 4U ′ − 2J φ16,|↑↓,↑↓〉
reads
Hloc = U
nα↑nα↓ + U
n1σn2σ′
n1σn2σ + J
2σ̄d1σ̄d2σ
1↓d2↓d2↑ + d
2↓d1↓d1↑
. (57)
The kinetic energy shall contain only intraband terms
for a basic tight-binding (TB) model for s-bands on a
3D simple cubic lattice with lattice constant a. Thus the
corresponding hamiltonian is written as
Hkin = −
α=1,2
iασdjασ , (58)
with the eigenvalues
εα(k) = −
µ=xyz
cos(kµa) , (59)
where tα denotes the hopping parameter for orbital
α=1,2. For convenience, we set a=1. The factor 1/3 in
eq. (58) is to normalize the total bandwidth to Wα=4tα.
Because of the cubic symmetry, U ′=U−2J may be used,
and furthermore we set J=JC . This model is similar to
the one considered by Bünemann et al.19 using a gener-
alized Gutzwiller approximation (see also [15]). Our sim-
pler TB description exhibits in principle perfect nesting,
however this issue is not relevant at the present level. The
3D two-band Hubbard model is studied to make contact
with the named previous work and in order to establish
the connection between the Gutzwiller and slave-boson
points of view.
When working in the SU(2) rotationally-invariant
case, the 24=16 atomic eigenstates |Γ〉 of the local hamil-
tonian (57) serve as the appropriate atomic basis (see Ta-
ble I), however also the simpler Fock basis (or any other)
may be used. Of course, in the Fock basis, a more com-
plicated energy matrix must be used in the saddle-point
equations (see Appendix D). It should be clear from
Tab. I that there are 20 nonzero slave-boson amplitudes
φΓn for the current problem. The S
z=0 triplet as well
as the three singlets are described with two φΓn, respec-
tively. In principle, even more φΓn may be introduced in
the beginning of the iteration cycle to minimize Ω, but at
convergence those one will come out to be strictly zero.
Of course, in high-symmetry situations there is still some
redundancy within the set of the 20 SBs. For instance,
for equal bandwidth at half filling (see Fig. 1), all the
one- and three-particle SBs are identical, as well as the
zero- and four particle SB. Moreover the Sz = ±1 triplet
SBs are equal because of the degeneracy. The two SBs
describing the Sz=0 triplet are also identical, with a mag-
nitude φ
(t,0)n
(t,±1)n/
2. Also the bosons describing
one specific singlet have the same absolute value, however
they carry the multiplet phase information, i.e., have plus
or minus sign. In conclusion, in the orbitally degenerate
case, the SB amplitudes at saddle point are of the form:
ϕΓn = 〈n|Γ〉 yΓ , (60)
in which the matrix element 〈Γ|n〉 is entirely determined
by Hloc and yΓ is a (coupling-dependent) amplitude, de-
pending only on the eigenstate Γ. This is more clearly in-
0 0.5 1.0 1.5 2.0
U/(2t)
φ1, |00,00>
φ2, |10,00>
φ10, |11,00>
φ6, |10,10>
φ11, |11,00>
FIG. 1: (color online) Inequivalent Slave-boson probabili-
ties |φΓn|2 for the two-band Hubbard model at half-filling
for equal bandwidth and J/U=0.2. Note that φ10,|↑↓,00〉 and
φ11,|↑↓,00〉 describe part of the singlet states, hence their over-
all amplitude is scaled by 1/
terpreted when atomic states are also used as basis states
for QPs (Sec. II D). Indeed, Eq. (60) means that:
ϕΓΓ′ = δΓΓ′ yΓ (61)
Hence, in this highly symmetric case, the saddle point is
indeed of the diagonal form considered in Refs. [19,21].
Once the symmetry is lowered, more SBs become inequiv-
alent and this relation does not hold anymore: there are
off-diagonal components even when the basis of atomic
states is used for both physical and QP states. In this
context, the present formalism becomes essential. Dif-
ferent bandwidths for each orbital, together with a finite
doping away from half-filling lead for instance to two dif-
ferent absolute values for the two SBs associated with the
singlets formed by the two doubly-occupied Fock states
(as seen at the end of this paragraph in Fig.6).
Since no interorbital hybridization is applied in this
section, the Ẑ-matrix is diagonal. We consider first the
simple case of equal bandwidths t1=t2=0.5 (note that
in all our applications, t sets the unit of energy), thus
Z11=Z22=Z. Figure 2 shows the variation of Z for
different ratios J/U in the half-filled case (n=2). The
critical coupling Uc for the Mott transition with J=0 ob-
tained from this slave-boson calculation is in accordance
with the result of the analytical formula given by Frésard
and Kotliar18. It is seen that an increased J lowers the
critical U and moreover changes the transition from sec-
ond to first order. Note that in this regard, Fig. 2 depicts
Z up to the spinodal boundary, i.e., the true transition
(following from an energy comparison) is expected to be
at slightly lower Uc. One can also observe the nonmono-
tonic character for the evolution of the critical Z at this
boundary when increasing J/U . We plot in Fig. 2 addi-
tionally the results when restricting the atomic hamilto-
nian to density-density terms only, in order to check for
the importance of the then neglected spin-flip and pair-
hopping terms. For larger J/U the critical Z from the
0 1 2 3 4
U/(2t)
rotationally-invariant
density-density
FIG. 2: (color online) Influence of J on the Mott
transition in the two-band Hubbard model at half fill-
ing (n=2) for equal bandwidth. From right to left:
J/U=0, 0.01, 0.02, 0.05, 0.10, 0.20, 0.45.
latter description is larger compared to the rotationally-
invariant one and moreover it is monotonically grow-
ing. The latter feature strengthens the first-order char-
acter in the density-density formulation for growing J/U ,
whereas for rotationally-invariant interactions this char-
acter is strongly weakened in that regime. Although for
J/U=0.45 the jump of Z is quite small, the transition is
however still first order in the present calculation. Fur-
thermore, there appears to be a crossover between the
two approaches concerning the reachable metallic spin-
odal boundary when increasing J/U .
At quarter filling (n=1) a continuous transition is ob-
tained for all the previous interactions (see Fig. 3). Com-
pared to the half-filled case, the density-density approx-
0 2 4 6 8
U/(2t)
rotationally-invariant
density-density
J/U=0.45
J/U=0.20
FIG. 3: (color online) Influence of J on the Mott transi-
tion in the two-band Hubbard model at quarter filling (n=1)
for equal bandwidth. The first three combined curves for
the two types of interactions (from left to right) belong to:
J/U=0, 0.05, 0.10. Arrows indicate the labelling for the two
larger J/U ratios.
0 0.2 0.4 0.6 0.8 1.0
J/(2t)
U = 1
U = 2
U = 3
FIG. 4: (color online) Influence of J for fixed U at n=2 (solid
lines) and n=1 (dotted-dashed lines) for equal bandwidth and
full SU(2) symmetry. The vertical dotted lines mark the limit
we set for J , respectively.
imation appears to be less severe for small J/U , but
leads to some differences compared to the rotationally-
invariant form for large J/U . Note that for J/U=0.45,
U ′−J in the local hamiltonian (57) becomes negative.
Thus a corresponding change of the ground state may
lead to the resulting nonmonotonic behavior for Uc then
observed in Fig. 3. The critical U for J=0 is smaller than
at half filling and with increasing J/U the transition is
shifted to larger Uc (with the above named exception for
J/U large). Hence J has a rather different influence on
the degree of correlation for the two fillings. While for
n=2 the Hund’s rule coupling substantially enhances the
correlations, seen by the decrease in Z, for n=1 the oppo-
site effect may be observed. This is also demonstrated in
Fig. 4 which displays the influence of J for fixed values of
U comparing half filling with quarter filling. The strong
decrease in Z upon increasing J was recently shown to be
0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0
filling
J/U = 0
J/U = 0.1
J/U = 0.2
FIG. 5: Filling dependence of Z for selected values of J/U
within the equal-bandwidth two-band model with full SU(2)
symmetry (U=1.75).
important for the physical properties of actinides, in par-
ticular regarding the distinct properties of δ-Plutonium
and Curium23. For each U shown in Fig. 4, the density-
density limiting value U/3 was used as an upper bound
for J . However, for U=2 and U=3 the system shows
already a first-order transition at half-filling below the
latter limit.
The QP residue Z is shown as a function of filling n
in Fig. 5 for U=1.75 and three ratios J/U . For J=0 it
is observed that Z(n) exhibits two minima, both located
at integer filling. The minimum at n=1 is deeper, cor-
responding to a lower value for Uc in the quarter-filled
case. Because of the filling-dependent effect of J seen in
Fig. 4, the nonmonotonic character of Z(n) is lifted for
growing J/U .
The two-band Hubbard model was already extensively
studied in the more elaborate DMFT framework in in-
finite dimensions. Such investigations reveal the same
qualitative change of the critical U for different integer
fillings24,25, of course with some minor quantitative dif-
ferences. Also the reduction24,26–29 of Uc and the onset
of a first-order Mott transition24,26,27,29 for finite J/U at
half-filling is in accordance. Concerning the latter effect,
the trend of weakening the first-order tendency for large
J/U is also reproduced and there is some discussion26,27
about the possibility of even changing back to a contin-
uous Mott transition in that regime. The increasing Uc
with growing J at quarter filling was also found by Song
and Zou29.
Finally, in Fig. 6 a comparison between the equal-
bandwidth and the different-bandwidth cases at noninte-
ger filling n=1.5 is displayed (J/U=0.2). ForW1=W2 the
model does not show a metal-insulator transition because
of the doping. Also the filling of both bands is iden-
0 1 2 3 4
U/(2t)
0 1 2 3 4
FIG. 6: (color online) Comparison of the two-band model
for W1=W2 (left) and W1=2W2 (right) at noninteger filling
n=1.5 and J/U=0.2. The ratio x is plotted for the singlet
states coupling the doubly occupied Fock states (see Tab I),
demonstrating that ϕΓΓ′ is no longer diagonal in this case.
tical and constant with increasing U (ns1=ns2=0.375),
and as stated earlier the SBs are still of the form given
by eq. (60). However, when breaking the symmetry be-
tween the two bands by considering different bandwidths,
the model behaves qualitatively rather differently. The
individual band fillings are not identical anymore, favor-
ing the larger-bandwidth band for U=0. With increas-
ing U the system manages to drive at least one band
insulating by transferring charge from the broader into
the narrower band, until the latter is filled with one elec-
tron30,31. Hence Z2 of the narrower band becomes zero at
an orbital-selective Mott transition (OSMT)30–35. This
asymmetric model has also a more sophisticated SB de-
scription, since for instance the SBs of the singlets built
out of the respective doubly-occupied Fock states have
now different amplitudes.
B. The Hubbard bilayer
Next, we consider a model consisting of two single-
band Hubbard models (two “layers”), coupled by an in-
terlayer hopping V . This rather popular model has al-
ready been subject of various studies36–39. For simplicity
and in order to make connection to the previous section,
each layer is described here by a 3D cubic lattice, with
an onsite repulsion U and an intralayer bandwidth Wα
(α = 1, 2), possibly different for the two layers. Hence
the local hamiltonian for this problem reads
Hloc = U
α=1,2
nα↑nα↓ + V
1σd2σ + d
2σd1σ
1σd1σ′d
2σ′d2σ , (62)
where the last term describes a possible spin-spin interac-
tion between the layers. However, for simplicity, we only
0 1 2 3
U/(2t)
V = 0
V = 0.25
FIG. 7: (color online) Half-filled bi-layer with equal band-
width for V=0 and V=0.25. For V=0 the filling per spin
within the two bands is identical (ns1=ns2=0.5), whereas
for finite V the the symmetry-adapted bonding/antibonding
states have different filling, denoted n+,n−.
0 0.5 1.0
U/(2t)
0 0.5 1.0 1.5 2.0
 = 0.75
 = 0.50
 = 0.25
FIG. 8: (color online) QP residues Zi and symmetry-adapted
fillings n+,n− for the half-filled bilayer. Left: W2/W1=0.5
and V=0. Right: various bandwidth ratios and V=0.1. In
the right part, the curves for smaller Z and n are associated
with the lower-bandwidth band.
present in this article results with J=0. Our choice of
kinetic energy is equivalent to the one in the last section,
i.e., given by Eqs. (58,59).
In the presence of V , an off-diagonal self-energy
Σ12(ω) is generated. Furthermore, away from half-filling
(n1+n2=2), this self-energy is expected to have a term
linear in ω at low frequency, and hence Z12 6=0. We note
that, when the bandwidths are equal (W1=W2), the bi-
layer model can be transformed into a two-orbital model
by a k-independent rotation to the bonding-antibonding
(or +,−) basis. In the latter basis, there is no hybridiza-
tion but instead a crystal-field splitting (=2V ) between
the two orbitals. The couplings of the two-orbital hamil-
tonian are given by (in the notation of the previous sec-
tion, and for J=0): Ueff=U
eff=Jeff=U/2. When the
bandwidths are different however, the interlayer hopping
cannot be eliminated without generating non-local inter-
dimer interactions.
Due to the reduced symmetry of the present model
in comparison to the two-band Hubbard model from the
previous section, the number of nonzero SBs φAn equals
now 36 (we use here the Fock basis for |A〉). We first con-
sider the simplest case of a half-filled system (n1=n2=1)
with equal bandwidths W1=W2 (and J=0). Results for
the intralayer QP weight and the orbital occupancies of
the bonding and antibonding bands are given in Fig. 7.
It is seen that the Mott transition is continuous for V=0
but becomes discontinuous in the presence of an inter-
layer hopping V 6=0. For V=0.25 the spinodal boundary
of the metallic regime is reached for U∼2.055. These re-
sults are consistent with findings in previous works36–40
within the DMFT framework.
Still focusing on the half-filled case, we display in Fig. 8
the QP weight as a function of U for different bandwidth
ratios W2/W1. When V=0, one has two independent
1.5 1.6 1.7 1.8 1.9 2.0
filling
-0.03
-0.02
-0.01
FIG. 9: (color online) Doped bilayer with equal bandwidth
and V=0.25 for U=2.054 (<Uc).
0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5
U/(2t)
-0.20
-0.15
-0.10
-0.05
FIG. 10: (color online) Bilayer at fixed doping (n=1.88) with
equal bandwidth and V=0.25. Full green (gray) lines: fillings
n+,n−, dashed dark lines: QP weights Z+, Z−. The vertical
dotted line marks the critical U at half filling.
Mott transitions in each layer, i.e., an OSMT scenario,
at which Z vanishes continuously. In the presence of a
non-zero V , this is replaced by a single discontinuous
transition for both orbitals. This is consistent with pre-
vious findings on the OSMT problem32.
We now consider the effect of finite doping away from
half filling. Fig. 9 displays the diagonal (Z11=Z22) ele-
ments as well as the now appearing Z12 element of the
QP weight matrix as a function of doping for U<Uc. Ad-
ditionally shown are the symmetry-adapted QP weights
Z+,− (occupations n+,−) which follow from diagonalizing
the Ẑ (∆̂(p)) matrix. Since Z12 is small in this case, the
Z+,− are rather similar to Z11=Z22 and merge with the
latter at half-filling. On the other hand, the polarization
of the (+,−) bands is still increasing.
As it is seen in Fig. 10 the off-diagonal component Z12
becomes increasingly important for larger U (>Uc) in the
doped case. It follows that in this regime the QP weights
Z+,− for the bonding/antibonding bands have rather dif-
ferent magnitude/behavior. Whereas Z− is monotoni-
cally decreasing, Z+ turns around and grows again (as
U = 0
U = 3.5
R Γ X M Γ
FIG. 11: (color online) QP bands of the doped bilayer model
(n=1.88) with equal bandwidths and V=0.25. The domi-
nately filled band is the bonding one, respectively.
also is the filling of the bonding band). Hence, this model
is a simple example in which a differentiation between QP
properties in different regions of the FS occur. Fig. 11
shows the QP (+,−)-bands in the noninteracting and in-
teracting case (U>Uc), exhibiting strong orbital polariza-
tion and different band narrowing close to the insulating
state. For very small doping and large U a transition to
a new metallic phase is found, which will be discussed in
detail in a forthcoming publication41.
Finally, we have also investigated a case in which
the interlayer (interorbital) hopping does not have a lo-
cal component (V=0), but does have a non-local one
V=t12 6=0, treated in the band term of the hamiltonian.
Hence the corresponding energy matrix reads here
ε(k) = −2
t11 t12
t12 t22
µ=xyz
cos(kµa) , (63)
with the choice t11=t22=0.5 and t12=0.25, as well as
a=1. In that case, a continuous Mott transition within
0 1 2
U/(2t)
FIG. 12: (color online) Half-filled bilayer, with equal band-
widths and V=0, but with a non-local interlayer hybridization
an OSMT scenario can be recovered, with, interestingly,
a sizeable value of the off-diagonal Z12 (Fig. 12). At the
transition Z11=Z22=Z12≡Zc holds, i.e., the Ẑ-matrix has
a zero eigenvalue, associated with the (antibonding) in-
sulating band. Note however that no net orbital polar-
ization appears with V being purely non-local.
C. Application to the momentum-dependence of
the quasiparticle weight within cluster extensions of
In this section, we finally consider the implications of
the rotationally invariant SB technique for the Mott tran-
sition and the momentum-dependence of the QP weight,
in the framework of cluster extensions of DMFT.
For simplicity, we consider a CDMFT approach to the
two-dimensional Hubbard model with nearest-neighbor
hopping t and a next-nearest neighbor hopping t′, based
on clusters consisting of two sites (dimers), arranged in
a columnar way on the square lattice (see Fig. 13). The
“local” hamiltonian on each dimer is formally identical
to the one introduced in the previous section for the bi-
layer model, i.e. Eq. (62), with the value V=−t of the
inter-‘orbital’ hybridization. The interdimer kinetic en-
ergy matrix reads (we set again a=1):
ε11(K) = ε22(K) = −2t cosKy (64)
ε12(K) = ε
21(K) = −t ei 2Kx − 2t′
1 + ei 2Kx
cosKy
in which K denotes a momentum in the reduced Bril-
louin zone (BZ) of the superlattice: Kx ∈ [−π/2,+π/2],
Ky ∈ [−π,+π]. Note again that in SB calculations, the
intradimer t has to be treated separately from the rest
of the kinetic energy within Hloc. It is easy to check
that when putting back −t into the offdiagonal elements
of the above kinetic-energy matrix, the eigenvalues just
correspond to the one of a single band:
ε(k) = −2t(cos kx + cos ky)− 4t′ cos kx cos ky (65)
PSfrag replacements 1 2
FIG. 13: Square lattice in the 2-site CDMFT picture.
0 0.5 1.0 1.5 2.0
U/(2t)
0 0.5 1.0 1.5 2.0
FIG. 14: (color online) Half-filled two-dimensional Hubbard
model within 2-site CDMFT. Left: QP weights and band
fillings, right: static dimer self-energy Σ.
in the full BZ of the original lattice kx,y ∈ [−π,+π].
In the following, we set t=0.25 and consider successively
t′=0 and t′= −0.3 t (a value appropriate to hole-doped
cuprates). Note that, in this article, we do not consider
a bigger cluster than the 2-site dimer, even in the pres-
ence of t′. Hence, the cluster self-energy will only con-
tain Σ11 and Σ12 components, i.e., has a spatial range
limited to the dimer. As a result, no renormalization of
the effective t′ is taken into account. This is of course
an oversimplification (particularly in view of the demon-
strated physical importance11,12 of Σ13 close to the Mott
transition). Larger clusters will be considered within the
present SBMFT in a further publication. The goal of the
present (simplified) study is to make a point of princi-
ple, namely that the SB formalism can indeed produce
a momentum-dependent Z(k). As the cluster symme-
try of the problem at hand is identical to the bilayer
model from the last section, the number of nonzero SBs
amounts again to 36. Figures 14, 15 and 16 summa-
rize the main findings, at half-filling and as a function of
doping, respectively. Let us first concentrate on the case
0.75 0.80 0.85 0.90 0.95 1.00
filling per site
FIG. 15: (color online) Doped two-dimensional Hubbard
model (t′=0) within 2-site CDMFT, for U=2.195 (Uc∼2.197).
0.0 0.5 1.0 1.5 2.0 2.5
U/(2t)
0.5 1.0 1.5 2.0 2.5
-0.02
-0.01
FIG. 16: (color online) Two-dimensional Hubbard model
within 2-site CDMFT at fixed doping n=0.94. Left: t-only
model, right: t-t′ model (t′=−0.3t). The vertical lines denote
the critical U at half filling, respectively.
t′=0. Obviously, the Mott transition at half-filling in this
case occurs in a manner which is very similar to the bi-
layer model with a finite interlayer hybridization studied
in the previous section (Fig. 7): a first-order transition is
found. The static part of the self-energy Σ11=Σ22 equals
U/2, while Σ12 (which has no frequency dependence at
half-filling within SBMFT) has a more complicated nega-
tive amplitude close to the transition. Note that we only
discuss the paramagnetic solution, though, of course, the
system is in principle unstable against antiferromagnetic
order for any U . Upon doping, a finite value of Z12 is
generated. The behavior of Z12 is rather similar to the
case of the bilayer model, except for the change of sign
(Figs. 15, 16). Hence its amplitude is again significantly
enhanced for U>Uc, i.e., the Z+,− values tend to mani-
festly deviate from each other. This is therefore signalling
an increasingly nonlocal component of Z(k) as the Mott
insulating state is approached at strong coupling. Again,
further studies in the latter regime at small doping will
be published soon41.
Including the effect of a nonzero nearest-neighbor hop-
ping t′ 6=0 turns out to lead to significant differences. Al-
though the first-order character of the transition remains
stable, the critical U is significantly lower (Fig. 14). The
static components of the self-energy behave rather simi-
larly to the t-only case, with some minor quantitative dif-
ferences. There is a small negative Z12 with a maximum
amplitude ∼0.02, remaining nonzero also at the Mott
transition (∼0.01). The main difference in comparison
to t′=0 is that here, in the doped case, Z12 changes sign
from negative to positive close to the insulating regime
for U>Uc+δ (with δ>0) (see Fig. 16). Thus the degree of
correlation of the effective (bonding-antibonding) bands
is inverted. These differences have to be interpreted with
caution however, since again the hopping range on the
lattice is larger than our cluster size, and definitive con-
clusions will have to be drawn from a study involving Σ13
as well.
Nonetheless, keeping with the simplified treatment
based on a 2-site cluster, we now describe the resulting
momentum dependence of the QP weight Z(k) for the t-t′
model. The matrix elements Σ11 and Σ12 of the cluster
(physical) self-energy matrix Σc are obtained from the
SB amplitudes at saddle point according to (55). The
self-energy is then periodized on the whole lattice, in the
form42,43:
Σlat(k, ω) = Σ11(ω) +
Σ12(ω)(cos kx + cos ky) . (66)
The interacting FS is defined as follows
µ− ε(k)− Σlat(k, ω=0) = 0 . (67)
For our case, using Eqs. (65, 66), this reads
µ− Σ11(0) +
Σ12(0)
(cos kx + cos ky)
+4t′ cos kx cos ky = 0 . (68)
Hence the FS deforms in a nontrivial way in the presence
of Σlat, when including t
′ in the present 2-site CDMFT
description. The QP weight Z(k) can be derived from
Σlat according to:
Z(k) =
Σlat(k, ω)
]−1∣∣∣∣∣
, (69)
which leads here to:
Z(k) =
−1]11 +
−1]12(cos kx + cos ky)
Z211 − Z212
Z11 −
Z12(cos kx + cos ky)
A contour plot of this function is displayed in Fig. 17.
Note that it varies only according to (cos kx+cos ky). Be-
cause the interacting FS involves both t and t′, and hence
 0.195
 0.196
 0.197
 0.198
 0.199
 0.201
 0.202
 0.203
 0.204
 0.205
PSfrag replacements
(0,0)
(π, π)
0.278
0.280
0.282
0.284
0.286
FIG. 17: (color online) Interacting Fermi surface (solid lines)
for the CDMFT treatment of the 2D t-t′ Hubbard model with
t′=−0.3t and U=2.5 at n=0.94 (per site). The color contours
show the variation of Z(k) (smallest at antinodes).
both lattice harmonics (cos kx + cos ky) and cos kx cos ky
(for t′ 6=0), it cuts through different contour lines of Z(k).
This results in a QP weight which varies on the FS. Fig-
ure 17 shows Z(k) for k close to the interacting FS.
Albeit the momentum variation is quantitatively quite
small, the key qualitative effect of Z being different on
different part of the FS is indeed found. It is seen that
the QPs along the nodal direction, i.e., along (0,0)-(π, π),
have slightly larger Z than the ones in the the antinodal
direction ((0,0)-(0, π)). Hence these results are indeed
in qualitative agreement with ARPES measurements on
cuprates. Note that to get nodal points to be more co-
herent that antinodal ones in this 2-site scheme, Z12>0
is actually crucial.
Our results provide, to our knowledge, the first ex-
ample of a SB calculation which can address the issue of
the momentum dependence of the QP weight. We believe
that the too small variation of Z along the FS found here
is due to the oversimplified 2-site description in which Σ13
is neglected. We intend to consider improvements on this
issue using the present SBMFT in a forthcoming work.
Finally, let us make contact with previous work on
the two-dimensional Hubbard model. Of course, this
model has been intensively studied with a variety of
methods such as: quantum Monte-Carlo44,45, exact diag-
onalization46,47, path-integral renormalization group48,
functional renormalization group49 and various quan-
tum cluster methods (dynamical cluster approxima-
tion50, cluster extensions of dynamical mean-field the-
ory7, and variational cluster perturbation theory8). We
shall not attempt here a detailed comparison between
the rotationally-invariant SB method (which anyhow is
a mean-field technique tailored to address low-energy is-
sues) with the results of these numerical methods over
the whole phase diagram (note in particular that we
have not yet investigated long-range ordered phases, such
as antiferromagnetism or superconductivity). Rather,
we would like to point out that some recent numerical
studies using the above methods10–12,49 have indeed re-
vealed the emergence of momentum-space differentiation
in the two-dimensional Hubbard model. We hope that
the rotationally-invariant SB method will help under-
stand qualitatively the low-energy physics emerging from
these results.
IV. CONCLUSION AND PERSPECTIVES
In this paper, we extended and generalized the
rotationally invariant formulation of the slave-boson
method13,14.
Our formulation achieves two goals: (i) extending the
slave-boson method in order to accommodate the most
general crystal fields, interactions and multiplet struc-
tures and (ii) the development of a technique which can
describe QP weights and Fermi liquid parameters which
vary along the Fermi surface.
The key aspect of the formalism is to introduce slave
boson fields which form a matrix with entries labeled by
a pair of a physical state and a QP state (within an ar-
bitrary choice of basis set). As a result, a density matrix
is constructed instead of just a probability amplitude for
each state.
While the first objective (i) could also be achieved
by generalizing appropriately the Gutzwiller approxima-
tion15,16, we find the slave-boson approach to be some-
what more flexible, in the sense that it is a mean field
theory which can in principle be improved by computing
fluctuations around the saddle point. Our application to
the two band model seems promising. While further work
is needed to benchmark the accuracy of the rotationally-
invariant slave-boson method against exact quantum im-
purity solvers, it is clear that already in the single site
multiorbital DMFT setting, our method has numerous
advantages. It obeys the Luttinger theorem even in the
presence of multiplets, and can accommodate full atomic
physics information. Furthermore, the off-diagonal el-
ements of the matrix of QP weights can be calculated
within this method, while the standard slave-boson or
Gutzwiller approximations (using probability amplitudes
instead of a density matrix) cannot achieve this goal.
Our technique achieves the second objective (ii) via a
detour, namely the use of cluster extensions of dynami-
cal mean-field theory in order to reduce the lattice to a
multisite (molecular) impurity problem, to which we ap-
ply our rotationally-invariant slave-boson method as an
impurity solver. Because the intersite matrix elements
of the QP weight can be calculated, it leads on the lat-
tice to a momentum-dependence of the QP residue Z(k).
We successfully demonstrated this point, in the frame-
work of a 2-site CDMFT study of the single-band 2D
Hubbard model. We did find that the QP weight at
the nodes is somewhat larger than at the antinodes, al-
though the magnitude of this effect is expected to in-
crease within a more realistic study involving a larger
cluster (e.g., a square plaquette), which is left for fu-
ture work. A major challenge is the direct extension of
our slave-boson approach to the lattice, without resorting
to the cluster-DMFT detour. In this context, we men-
tion that other slave-boson techniques, which introduce
magnetic correlations through the use of link variables to
decouple the superexchange J term1, can be interpreted
in terms of a k dependent self-energy. However, within
such schemes, the derivative of the self-energy with re-
spect to frequency is momentum independent (in contrast
to the static part), yielding a k independent QP residue.
Hence, our approach goes beyond these methods, at least
in conjunction with the cluster-DMFT approach. We
hope that having an economic impurity solver based on
SBs will allow us to study larger cluster sizes than feasible
with other methods, and most importantly help us un-
derstanding the low-energy physics emerging from these
cluster dynamical mean-field theories.
Finally, we limited our study to slave bosons which
do not mix the particle number. The extension to
full charge-rotational invariance and superconductivity
is possible (see Refs. [14,22] in the single-orbital case),
and will be presented in a separate paper. In this con-
text, the slave boson method will incorporate the SU(2)
charge symmetry and its extension away from half fill-
ing considered by Wen and Lee51 and the rotationally-
invariant slave-boson formalism can serve as a power-
ful tool for interpreting the low-energy physics emerging
from plaquette-CDMFT studies of this issue.
Acknowledgments
We are grateful to Pablo Cornaglia and Michel Ferrero
for very useful discussions and remarks. As this work
was being completed, we learned of a parallel effort by
Michele Fabrizio52, in the framework of the Gutzwiller
approximation. In particular, the form of the constraints
advocated in this work matches our constraints (28,29)
in SB language. A.G. also thanks him for discussions.
This work has been supported by the “Chaire Blaise Pas-
cal” (région Ile de France and Fondation de l’Ecole Nor-
male Supérieure), the European Union (under contract
“Psi-k f-electrons” HPRN-CT-2002-00295), the CNRS
and Ecole Polytechnique. G.K. is supported by the NSF
under Grant No. DMR 0528969.
APPENDIX A: SINGLE-ORBITAL CASE AND
CONNECTION WITH PREVIOUS WORK
Here, we briefly consider the single-orbital case
(M=2), which also allows to make contact with Refs. [13,
14]. These authors introduced in this case a rotationally-
invariant formalism, with the calculation of response
functions associated with the saddle point as their main
motivation. For M=2, the following local basis set can
be considered (whether or not Hloc is diagonal in this
basis):
N = 0 : |0〉 ,
N = 1 : |σ〉 = d †σ |0〉 (A1)
N = 2 : |D〉 = d †↑d
↓ |0〉 .
Hence, we introduce the following bosons (not mixing
sectors with different particle numbers, i.e., not consid-
ering superconducting states):
φ00 ≡ φE , φσσ′ , φ↑↓ ≡ φD . (A2)
Up to normalizations, the bosons p
σσ′ introduced in
Ref. 13 correspond to φ
σσ′ . In contrast, the stan-
dard Kotliar-Ruckenstein17 formalism introduces only
two bosons p†σ in the one-particle sector. The represen-
tatives (27) of the physical states read here:
|0〉 = φ†E |vac〉
|σ〉 =
σ′ |vac〉 (A3)
|D〉 = φ†Dd
↓ |vac〉 ,
and the constraints (28,29) read:
1 = φ
EφE +
σσ′φσσ′ + φ
DφD (A4)
f †αfα = φ
DφD +
φ†σαφσα (A5)
↑f↓ =
σ↓φσ↑ (A6)
↓f↑ =
σ↑φσ↓ . (A7)
Not including, for simplicity, the square-root normaliza-
tions in (37), needed however in order to insure a correct
U=0 limit at saddle point, the ‘simplest’ expression (32)
of the electron creation operators read:
↑βφE + (−1)
Dφ↓β ] f
β (A8)
↓βφE − (−1)
Dφ↑β ] f
β . (A9)
Apart from the motivations of Ref. [13] (associated with
fluctuations and response functions), the usefulness of the
rotationally-invariant scheme in the single-orbital case
can be demonstrated on a toy model consisting of a one-
band Hubbard model with a magnetic field, purposely
written in the Sx direction (i.e., in the form h d
↑d↓+h.c,
analogous to a hybridization). Although the direction of
the field should not matter, a direct application of the
standard Kotliar-Ruckenstein formalism is impossible in
that case. The rotationally-invariant formalism can be
shown to lead to the correct saddle point, independently
of the spin-quantization axis.
APPENDIX B: DERIVATION OF EQ. (29)
In this section, we show that the physical states of the
form (27) are exactly those selected by the constraints
(29) and (28). First, it is easy to check that states of the
form (27) do satisfy these constraints. Indeed, let us act
on the state |C〉 ≡ 1√
Cm|vac〉 ⊗ |m〉f and with
(29). The l.h.s leads to:
f †α fα′ |C〉 =
Cm|vac〉 ⊗ f
α fα′ |m〉f
〈m′|f †α fα′ |m〉φ
Cm|vac〉 ⊗ |m
′〉f .
When acting with the r.h.s, only the term A=C and n=m
gives a non-vanishing contribution, hence:
An′φAn 〈n|f
αfα′ |n′〉|C〉 =
〈n|f †αfα′ |n′〉φ
Cn′ |vac〉 ⊗ |n〉f .
We now prove that (29) are sufficient conditions, which
is a bit more difficult. Since (28) excludes states with
more than one boson, it is enough to consider a general
state of the form:
|C;W 〉 ≡
Wpq φ
Cp|vac〉 ⊗ |q〉f . (B2)
and to show that (29) implies Wpq ∝ δpq. Acting on this
state with each term in the constraint (29) yields for the
l.h.s:
f †α fα′ |C;W 〉 =
Wpq φ
Cp|vac〉 ⊗ f
αfα′ |q〉f
Cp|vac〉 ⊗ |r〉f
Wpq〈r|f †αfα′ |q〉 .(B3)
Let us now act with the r.h.s. Only the terms with A=C
and n=p contribute, leading to:
An′φAn 〈n|f
αfα′ |n′〉|C;W 〉
Wpq φ
Cn′ |vac〉 ⊗ |q〉f 〈p|f
αfα′ |n′〉
Cp|vac〉 ⊗ |r〉f
Wqr〈q|f †αfα′ |p〉 ,(B4)
where the last expression comes from a change of indices
n′ → p, q → r, p→ q.
We see that the constraint is satisfied provided that
the following identity holds, for all orbital indices αα′
and all states p, r:
Wpq〈r|f †αfα′ |q〉 =
Wqr〈q|f †αfα′ |p〉 . (B5)
Let us first look at the case α=α′, which reads:
rαWpr = pαWpr . (B6)
Hence Wpr=0 unless pα=rα for all α, so that:
Wpq = wp δpq . (B7)
Substituting this into (B5), we obtain:
wp 〈r|f †αfα′ |p〉 = wr 〈r|f †αfα′ |p〉 . (B8)
Thus, if r and p are related by a move of a QP from
one state to another (a transposition of two occupation
numbers), then wp=wr. Moreover, two Fock states in
the same sector HN of the Hilbert space are related by
a permutation of the occupied states, which can be de-
composed in a product of transpositions. Hence, wp is a
constant for p ∈ HN , and Wpq∝δpq as claimed.
APPENDIX C: PHYSICAL CREATION
OPERATOR
1. Proximate expression
First let us note that there is a systematic route to
find the expression for d, which consists in writing the
operator as :
d†α =
〈A|d †α|B〉 |A〉〈B|
n∈HA,m∈HB
〈A|d †α|B〉φ
AnφBmX
nm , (C1)
withXfnm=|n〉f 〈m|f in usual Hubbard notations. Xfnm is
obviously not just a one-particle operator f †, even when
restricted to the sectors of interest in the above formula,
since the states n and m can differ in many places. How-
ever, because any transposition of two QPs, when acting
on a physical state, can be replaced by a corresponding
operation on bosons using the constraint (29), and be-
cause any product of bosonic operators which cannot be
reduced to a quadratic form will produce a state which is
out of the physical subspace, the physical operator must
in the end take the form:
d†α =
CAnBm(α, β) φ
AnφBm f
β . (C2)
One can solve for the coefficients CAnBm(α, β), requesting
proper action on the physical states.
In this section however, we restrict ourselves to proving
that Eq. (32) does the job, i.e., that
d†α =
β,AB,nm
〈A|d †α|B〉〈n|f
β |m〉√
NA(M −NB)
AnφBm f
β (C3)
satisfies
d†α |B〉 =
〈A|d †α|B〉 |A〉 . (C4)
We start by proving the formula :
〈n|f †β|p〉 f
β |p〉 = (N + 1) |n〉 , (C5)
where the sum over p runs over the basis of the subspace
HN (states with N QPs) of the Fock space. First, we
have in general :
〈n|f †β |p〉 f
β|p〉 =
n′∈HN+1
an′ |n′〉 , (C6)
but, because f † is the creation operator (it connects one
basis state to only one another) :
an′ =
β,p∈HN
〈n|f †β|p〉〈n
′|f †β|p〉 ∝ δnn′
β,p∈HN
∣∣∣〈n|f †β |p〉
We now use the fact that the tensor is invariant
(it has the same expression in every basis) and use
the notations introduced for Eq. (40): U is a uni-
tary transformation of the one QP states and U is
the corresponding transformation in the Fock states.
〈n|f †β |m〉=Uββ′U∗nn′〈n′|f
β′ |m′〉 Umm′ and we have :
∣∣∣〈n|f †β |m〉
β,m∈HN
β′m′∈HN
n′∈HN+1
β′′m′′∈HN
n′′∈HN+1
Uββ′U
ββ′′U∗nn′Unn′′
×Umm′U∗mm′′〈n′|f
′〉〈n′′|f †β′′ |m
n′,n′′∈HN+1
U∗nn′Unn′′〈n′|f
β′ |p〉〈n
′′|f †β′ |p〉
β,p∈HN
∣∣∣〈(Un)|f †β |m〉
Moreover, any couple of elements of the basis of the Fock
state can be connected by a U transformation (with a U
that permutes the one QP basis state), therefore an≡a is
a independent of n. a can then be determined by sum-
ming (C5) over n :
n∈HN+1
n∈HN+1
∣∣∣〈n|f †β |p〉
n∈HN+1
βfβ|n〉
n∈HN+1
(N + 1) , (C9)
leading to a=N + 1. This completes the proof of (C5).
It is now simple to compute the action of (32) : acting
on |C〉≡ 1√
Cp|vac〉⊗ |p〉f with this operator, only
the term B=C,m=p contributes, and we get :
d†α|C〉 =
DC(NC + 1)(M −NC)
A,n∈HNC+1
〈A|d †α|C〉φ
An|vac〉 ⊗
p∈HNC
〈n|f †β |p〉 f
NC + 1
DC(M −NC)
A,n∈HNC+1
〈A|d †α|C〉φ
An|vac〉 ⊗ |n〉
〈A|d †α|C〉 |C〉 , (C10)
which is identical to (31).
2. Improved expression
In this section we present arguments for the improved
formula used in this paper. First, it is useful to define the
“natural orbitals” (NO) basis as the basis which diago-
nalizes the quasiparticle and quasihole density matrices
corresponding to the average constraint, which is given
αβ [φ] ≡
φ∗AnφAm〈m|f †αfβ |n〉 (C11a)
αβ [φ] ≡
φ∗AnφAm〈m|fβf †α|n〉 (C11b)
φ∗AnφAn − ∆̂
αβ [φ] .
Let us denote by ξλ, |λ〉 the eigenvalues and eigenvectors
of those matrices :
ξλ〈α|λ〉〈λ|β〉 , |λ〉 =
〈α|λ〉 |α〉 , (C12)
which is equivalent to use the NO quasiparticle operator
λ such that :
〈λ|α〉f †α , 〈ψ
λψµ〉 = δλµξλ . (C13)
To be fully explicit, we can consider the particular basis
transformation (in the notations of the section above)
f †α=Uαλψ
λ which rotates to the NOs, and the corre-
sponding rotation on the bosons: φAn=U(U)nn′ΩAn′ .
The rotation matrix is:
Uαλ = 〈α|λ〉 , (C14)
and in the NO basis:
Ω∗AnΩBm〈m|ψ
λψµ|n〉 = δλµ
Ω∗AnΩAnnλ
= δλµξλ({ΩAn}) , (C15)
αβ [φ] =
Uαλ ξλ [U
†]λβ . (C16)
The idea is to generalize the Kotliar-Ruckenstein nor-
malization factor in the NO basis, where the QP density
being diagonal, its probabilistic interpretation is more
transparent. Hence, the improved expression of d reads:
d†α =
λ,AB,nm
〈A|d †α|B〉 〈n|ψ
λ|m〉√
ξλ({ΩAn})(1− ξλ({ΩAn}))
AnΩBm ψ
(C17)
Note that the formal square-root normalisation, i.e.,
NA(M −NB), does not appear in this representa-
tion. We can now rotate back to the generic basis we
started from and use the gauge invariance, leading to :
d†α =
AB,nm,βγ
CAnBm(α, β)φ
AnφBm
〈β|λ〉〈λ|γ〉
ξλ(1 − ξλ)
AB,nm,βγ
CAnBm(α, β)φ
AnφBm 〈β|[∆̂
(p)∆̂(h)]−
2 |γ〉 f †γ ,(C18)
with CAnBm(α, β) =
A|d†α|B
n|f †β|m
.(C19)
Hence this yields the following form for the R-matrix:
R[φ]∗αβ =
AB,nm,δ
CAnBm(α, δ)φ
AnφBm 〈δ|[∆̂
(p)∆̂(h)]−
2 |β〉 .
(C20)
In the actual implementation of the saddle point calcu-
lations, the explicit use of both the quasiparticle and the
quasihole density matrices has been utilized (i.e., not us-
ing their relation). Although at convergence the different
representations yield the same values, writing the equa-
tions via both, particle and hole density matrix, appears
to be necessary within the minimization cycle. This is
due to the fact that the derivatives with respect to the
slave bosons have to be symmetric, however when using
only the particle density matrix (or its eigensystem de-
composition) for instance the derivative with respect to
the empty boson vanish, although this one exists in the
KR case. In the end an even more symmetrized form,
i.e., 1
(∆̂(p)∆̂(h)+∆̂(h)∆̂(p)) was used for the square root
in eq. (C18). Thus defining the following matrix:
Mγβ =
∣∣∣∣∣
(∆̂(p)∆̂(h) + ∆̂(h)∆̂(p))
∣∣∣∣∣
, (C21)
the electron operators are written as
d†α =
CAnBm(α, γ)φ
AnφBmMγβ f
R∗αβf
β (C22)
CAnBm(α, γ)φ
BmφAnMβγ fβ
Rαβfβ , (C23)
and correspondingly independently written the elements
of the R,R†-matrices read
R[φ]αβ =
AB,nm,γ
CAnBm(α, γ)φ
BmφAn M̂βγ (C24)
R†[φ]αβ ≡ R[φ]∗βα =
AB,nm,γ
CAnBm(β, γ)φ
AnφBm M̂γα .(C25)
APPENDIX D: DETAILS ON THE
SADDLE-POINT EQUATIONS AND THEIR
NUMERICAL SOLUTION
The saddle-point equations for T=0 are obtained by
performing the partial derivatives with respect to all the
variables, i.e., condensed slave-boson amplitudes and La-
grange multipliers:
AnϕAn (D1)
= 〈f †αfβ〉 −
A,nn′
n|d†αdβ |n′
ϕAn (D2)
Am + λ0ϕ
m|d†αdβ |n′
Cn′ (D3)
Bm + λ0ϕCm
n|d†αdβ |m
ϕCn (D4)
with 〈f †αfβ〉 =
f̃kj 〈α|νkj〉 〈νkj |β〉 . (D5)
The εkj are the eigenvalues (with band index j) of the QP
matrix (R†(ϕ)ε(k)R(ϕ) +Λ) with corresponding eigen-
vector |νkj〉, while f̃kj denotes the occupation number of
the state |νkj〉 for a given total number of particles, to
be evaluated by standard k-integration techniques (e.g.
tetrahedron method, Gaussian smearing, etc.).
1. Some slave-boson derivatives
a. Eigenvalues and R matrices. The derivatives of
the eigenvalues with respect to the slave bosons, i.e.,
may be performed pertubatively:
ε(k)R + Λ
)∣∣∣∣ νkj
ε(k)R +R†ε(k)
∣∣∣∣ νkj
∣∣∣∣∣∣
〈β|ε(k)R
ε(k)|α〉∂R̂αβ
∣∣∣∣∣∣
〈νkj |α〉
〈β|ε(k)R|νkj〉
+〈νkj |R†ε(k)|α〉
∂R̂αβ
〈β|νkj〉
.(D6)
The therefore needed explicit expressions for the deriva-
tives of the R,R†-matrices read as follows (using
eqs.(C24,C25)):
∂R̂αβ
AB,nn′,γ
CAnBn′(α, γ)ϕ
δCmAn M̂βγ + ϕAn
∂M̂βγ
∂R̂αβ
AB,nn′,γ
CAnBn′(α, γ)ϕAn
δCmBn′M̂βγ + ϕ
∂M̂βγ
analogous for [R̂†]αβ
b. The M matrix. As it is seen, the derivatives in-
volve the derivative of the M matrix (C21). This deriva-
tive is computed as follows. Lets first write M as
M̂γβ = 〈γ|K−1/2|β〉 . (D7)
What we are looking for is the derivative of K−1/2 with
respect to the SBs. In order to get access to this quantity
we use the identity
−1/2)K−1/2 +K−1/2(∂ϕK
−1/2) = ∂ϕK
−1 (D8)
⇔ XK−1/2 +K−1/2X = Y , (D9)
with X=∂ϕK
−1/2 and Y=∂ϕK
−1. We then apply P
which transforms K to its eigensystem. This yields
L + LX′ = Y′ , (D10)
where the prime denotes that the quantities de-
fined above are expressed in that eigensystem, and
L=P†K−1/2P. Since in the eigensystem K−1/2 is di-
agonal, i.e., L is, the last equation can be written in
components and X′ determined:
X ′ijLj + LiX
ij = A
ij ⇔ X ′ij =
Li + Lj
. (D11)
Backtransforming to X=PX′P†=∂ϕK
−1/2 yields the
desired derivative of K and subsequently of M.
To perform the described computation we need to
know ∂ϕK
−1 in eq. (D8), however this quan-
tity may be straightforwardly calculated when start-
ing from the identity KK−1/2K−1/2=1, resulting in
−1=−K−1(∂ϕK)K−1.
2. Mixing
In order to solve the saddle-point equations, a method
to deal with a system of nonlinear equations F as a func-
tion of the variables (slave bosons, Lagrange multipliers)
x has to be utilized:
F(x) = 0 (D12)
In the present work we tested several quasi-Newton tech-
niques (e.g. Broyden53, modified Broyden54, etc.) to
handle this numerically. Thereby from a starting guess
for x the variables are updated via
x(m+1) = x(m) + J
F(m) , (D13)
since we want F(m+1) to be zero to linear order. The
jacobian J is here defined as follows
Jij ≡ −
(D14)
and is not calculated exactly (this would involve sec-
ond derivatives and would lead to the Newton-Raphson
method) but is computed at each step m via formulae
which dictate several constraints on how J should evolve.
In our numerical implementation we found the modified
Broyden scheme to be well suited for the so far investi-
gated applications Note that there is usually no need for
explicitly fixing the gauge for the numerical solution of
the saddle-point equations. The initial amplitudes of the
variational parameters, i.e., slave bosons and lagrange
multipliers, at the start of the iteration, together with
the choice of the atomic basis |A〉, always ensured proper
convergence to one of the family of solutions within our
implementation.
∗ Electronic address: Frank.Lechermann@physnet.
uni-hamburg.de
1 G. Kotliar, in Strongly Interacting Fermions and High-Tc
superconductivity, edited by B. Doucot and J. Zinn-Justin
(Elsevier, 1995), Les Houches, Session LVI, p. 197.
2 M. Imada, A. Fujimori, and Y. Tokura, Rev. Mod. Phys.
70, 1039 (1998).
3 A. Damascelli, Z. Hussain, and Z.-X. Shen, Rev. Mod.
Phys. 75, 473 (2003), cond-mat/0208504.
4 M. R. Norman et al, Nature 392, 157 (1998).
5 M. Le Tacon, A. Sacuto, A. Georges, G. Kotliar, Y. Gallais,
D. Colson, and A. Forget, Nature Physics 2, 537 (2006),
arXiv:cond-mat/0603392.
6 T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev.
Mod. Phys. 77, 1027 (2005).
7 G. Kotliar, S. Y. Savrasov, K. Haule, V. S. Oudovenko,
O. Parcollet, and C. A. Marianetti, Rev. Mod. Phys. 78,
865 (2006).
8 A. M. S. Tremblay, B. Kyung, and D. Senechal, J. Low
Temp. Physics 32, 424 (2006).
9 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,
Reviews of Modern Physics 68, 13 (1996).
10 D. Sénéchal and A.-M. S. Tremblay, Phys. Rev. Lett. 92,
126401 (2004).
11 M. Civelli, M. Capone, S. S. Kancharla, O. Parcollet, and
G. Kotliar, Phys. Rev. Lett. 95, 106402 (2005).
12 O. Parcollet, G. Biroli, and G. Kotliar, Phys. Rev. Lett.
92, 226402 (2004).
13 T. Li, P. Wölfle, and P. J. Hirschfeld, Phys. Rev. B. 40,
6817 (1989).
14 R. Frésard and P. Wölfle, Int. J. Mod. Phys. B 6, 685
(1992).
15 C. Attaccalite and M. Fabrizio, Phys. Rev. B 68, 155117
(2003).
16 M. Ferrero, Ph.D. thesis, SISSA-Trieste (2006).
17 G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57,
1362 (1986).
18 R. Frésard and G. Kotliar, Phys. Rev. B 56, 12909 (1997).
19 J. Bünemann, W. Weber, and F. Gebhard, Phys. Rev. B
67, 6896 (1998).
20 B. R. Trees, A. J. Fedro, and M. R. Norman, Phys. Rev.
B 51, 6167 (1995).
21 X. Dai, G. Kotliar, and Z. Fang, cond-mat/0611075 (2006).
22 B. R. Bulka and S. Robaszkiewicz, Phys. Rev.B 54, 13138
(1996).
23 J. H. Shim, K. Haule, and G. Kotliar, Nature 446, 513
(2007).
24 Y. Ono, M. Potthoff, and R. Bulla, Phys. Rev. B 67,
035119 (2003).
25 M. Rozenberg, Phys. Rev. B 55, R4855 (1997).
26 T. Pruschke and R. Bulla, Eur. Phys. J. B 44, 217 (2005).
27 K. Inaba and A. Koga, Phys. Rev. B 73, 155106 (2006).
28 J. E. Han, M. Jarrell, and D. L. Cox, Phys. Rev. B 58,
R4199 (1999).
29 Y. Song and L.-J. Zou, Phys. Rev. B 72, 085114 (2005).
30 A. Koga, N. Kawakami, T. M. Rice, and M. Sigrist, Phys.
Rev. Lett. 92, 216402 (2004).
31 A. Rüegg, M. Indergaard, S. Pilgram, and M. Sigrist, Eur.
Phys. J. B 48, 55 (2005).
32 L. de’ Medici, A. Georges, and S. Biermann, Phys. Rev. B
72 (2005), cond-mat/0503764.
33 V. I. Anisimov, I. Nekrasov, D. Kondakov, T. M. Rice, and
M. Sigrist, Eur. Phys. J. B 25, 191 (2002).
34 A. Liebsch, Phys. Rev. B 70, 165103 (2004).
35 M. Ferrero, F. Becca, M. Fabrizio, and M. Capone, Phys.
Rev. B 72, 205126 (2005).
36 H. Monien, N. Elstner, and A. J. Millis, cond-mat/9707051
(1997).
37 G. Moeller, V. Dobrosavljević, and A. E. Ruckenstein,
Phys. Rev. B 59, 6846 (1999).
38 A. Fuhrmann, D. Heilmann, and H. Monien, Phys. Rev. B
73, 245118 (2006).
39 S. S. Kancharla and S. Okamoto, cond-mat/0703728
(2007).
40 M. Capone, private communication (2007).
41 M. Ferrero, P. S. Cornaglia, and et al., unpublished (2007).
42 G. Biroli, O. Parcollet, and G. Kotliar, Phys. Rev. B 69,
205108 (2004).
43 G. Kotliar, S. Y. Savrasov, G. Pálson, and G. Biroli, Phys.
Rev. Lett. 87, 186401 (2001).
44 J. E. Hirsch, Phys. Rev. B 31, 4403 (1985).
45 N. Furukawa and M. Imada, J. Phys. Soc. Jpn. 61, 3331
(1992).
46 G. Fano, F. Ortolani, and A. Parola, Phys. Rev. B 42,
6877 (1990).
47 E. Dagotto, A. Moreo, F. Ortolani, D. Poilblanc, and J. Ri-
era, Phys. Rev. B 45, 10741 (1992).
48 T. Kashima and M. Imada, J. Phys. Soc. Jpn. 70, 2287
(2001).
49 D. Rohe and W. Metzner, Phys. Rev. B 71, 115116 (2005).
50 T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev.
Mod. Phys. 77, 1027 (2005).
51 X.-G. Wen and P. A. Lee, Phys. Rev. Lett. 76, 503 (1996).
52 M. Fabrizio, unpublished (2007).
53 C. G. Broyden, Math. Comp. 19, 577 (1965).
54 D. Vanderbilt and S. Louie, Phys. Rev. B 30, 6118 (1984).
55 In this paper, the chemical potential is usually included in
the one-body part ε0αβ of the local hamiltonian.
56 Note however that both sides of (24) have identical matrix
elements between physical states
