Gutzwiller description of non-magnetic Mott insulators: a dimer lattice model.
Michele Fabrizio,1,2
1 International School for Advanced Studies (SISSA) and CNR-INFM-Democritos
National Simulation Centre, Via Beirut 2-4, I-34014 Trieste, Italy
2 The Abdus Salam International Centre for Theoretical Physics (ICTP), P.O. Box 586, I-34014 Trieste, Italy
(Dated: September 12, 2021)
We introduce a novel extension of the Gutzwiller variational wavefunction able to deal with
insulators that escape any mean-field like description, as for instance non-magnetic insulators. As
an application, we study the Mott transition from a paramagnetic metal into a non-magnetic Peierls,
or valence-bond, Mott insulator. We analyze this model by means of our Gutzwiller wavefunction
analytically in the limit of large coordination lattices, where we find that: (1) the Mott transition
is first order; (2) the Peierls gap is large in the Mott insulator, although it is mainly contributed by
the electron repulsion; (3) singlet-superconductivity arises around the transition.
PACS numbers: 71.10.-w, 71.10.Fd, 71.30.+h
I. INTRODUCTION
Among the theoretical tools devised to deal with
strongly correlated metals close to a Mott metal-to-
insulator transition (MIT), the simplest one likely is
the variational approach introduced in the 60ths’ by
Gutzwiller1,2 to describe itinerant ferromagnetism and
narrow band conductors. In its original version, the
Gutzwiller variational wavefunction has the form
|ΨG〉 = P |φ〉 =
PR |φ〉, (1)
where |φ〉 is an uncorrelated wavefunction for which
Wick’s theorem holds, PR an operator at site R, and
both |φ〉 and PR have to be determined by minimizing
the variational energy. The role of the operator PR is to
modify, according to the on-site interaction, the weights
of the local electronic configurations with respect to their
values in the uncorrelated wavefunction.
In spite of its simplicity, the Gutzwiller wavefunction
is quite effective in capturing physical properties that
supposedly identify strongly correlated metals, as for in-
stance the large increase of the effective mass.3 However,
since the dependence upon the distance |R−R′| of inter-
site correlations are still determined by the uncorrelated
wavefunction, while the local operators PR just affect the
amplitudes, the Gutzwiller wavefunction can describe a
Mott insulator either if PR suppresses completely charge
fluctuations, that provides a very poor description of an
insulator, or if |φ〉 itself is insulating. The latter case
can be stabilized within the original Gutzwiller approach
only when |φ〉 is an admissible Hartree-Fock solution of
the Hamiltonian. As an example let us consider a sin-
gle band model at half-filling, for instance the Hubbard
model
H = −
RR′,σ
tRR′ c
nR,↑ nR,↓,
where c
and c
creates and annihilates, respectively,
an electron with spin σ =↑, ↓ at site R and nR,σ =
. This Hamiltonian admits at the mean-field
level two possible phases, one paramagnetic, 〈nR,↑〉 =
〈nR,↓〉, and the other magnetic, 〈nR,↑〉 6= 〈nR,↓〉. The
latter is the only one that can eventually describe an
insulator. In the Hubbard model the action of the oper-
ator P is to increase the weight of singly occupied sites
at expenses of doubly occupied and empty sites, in or-
der to minimize the Coulomb repulsion U . Evidently,
even when the repulsion is very strong, hence the model
is a Mott insulator, a realistic wavefunction should still
allow for charge fluctuations responsible for the super-
exchange, that survives even deep inside the Mott phase.
However, any paramagnetic uncorrelated wavefunction,
for instance the Fermi sea, is unable to generate any
super-exchange and necessarily leads to a non-realistic
Mott insulator where configurations with empty or dou-
bly occupied sites are fully suppressed.3 The only way
to generate super-exchange is to assume a magnetically
ordered |φ〉, which is also the only insulating wavefunc-
tion accessible within Hartree-Fock. However a magnetic
state might not always be the right choice, especially if
magnetism is sufficiently frustrated.
Recently, an improved version of the Gutzwiller wave-
function has been proposed,4 in which additional inter-
site correlations are provided by density-density Jastrow
factors, namely
|ΨG〉 → exp
vR,R′ nR nR′
|ΨG〉, (2)
where nR is the site R occupation number and vR,R′
variational parameters. This novel class of wavefunctions
has the capability to disentangle charge from other de-
grees of freedom, hence is more suitable to capture Mott
localization, as it has indeed been shown.4,5,6 However,
unlike the conventional Gutzwiller wavefunction (1), the
Gutzwiller-Jastrow wavefunction (2) can only be dealt
with numerically by variational Monte Carlo, which is in-
herently limited to finite-size systems, albeit quite large.7
An alternative approach, that is closely related to
recently proposed extensions of Dynamical Mean Field
Theory (DMFT) from the original single-site formula-
http://arxiv.org/abs/0704.1389v1
tion8 to a cluster one9,10,11,12,13, is to consider a vari-
ational wavefunction of the same form as (1) but de-
fined on a lattice with non-primitive unit cells. In this
case, the operator PR acts on all the available electronic
configurations of the lattice sites belonging to the non-
primitive cell. The advantage is that in this way one may
include additional short-range correlations without losing
the property of the wavefunction to be analytically man-
ageable, at least in infinite-coordination lattices. The ob-
vious disadvantage is that this wavefunction could bias
the variational solution towards translational-symmetry
breaking.
Within this scheme, the variational problem becomes
generically equivalent to optimize a Gutzwiller wave-
function for a multi-band Hamiltonian. There have
been recently an amount of attempts to extend the
Gutzwiller wavefunction to multi-orbital models that in-
clude further complications like for instance Coulomb ex-
change14,15,16,17. In this paper we introduce a further
extension that is capable to generate inter-site correla-
tions as the super-exchange for paramagnetic wavefunc-
tions, otherwise missed by the conventional Gutzwiller
approach. This novel class of wavefunctions also allows
to explore new kinds of variational solutions. Specifi-
cally, there are interesting examples of correlated models
where the Mott insulating phase escapes any Hartree-
Fock mean-field treatment, in other words can not be
represented by a single Slater determinant. A very sim-
ple case, that we will explicitly consider throughout this
work, is a Peierls insulator, namely a short-range valence-
bond crystal, in which pairs of nearest neighbor sites are
strongly bound into a singlet configuration, leading to a
state that is simply a collection of spin-singlets. Such a
Mott insulating state is not accessible by Hartree-Fock
theory, just because each singlet is itself not express-
ible as a Slater determinant, nor by the conventional
Gutzwiller approach, which, as mentioned, gives a poor
description of paramagnetic insulators.
The paper is organized as follows. In Section II we
present the variational wavefunction and discuss under
which conditions it can be deal with analytically. In Sec-
tion III we discuss how to build up the wavefunction in
the case in which the basic unit of the lattice model is
a dimer. Next, in Section IV, we solve the variational
problem for a specific lattice model of dimers. Conclu-
sions are given in Section V.
II. THE VARIATIONAL WAVEFUNCTION
In this Section, we introduce an extension of the
Gutzwiller wavefunction (1) which is particularly con-
venient to perform analytical calculation in the limit of
infinite-coordination lattices.14 Let us consider a generic
multi-band Hamiltonian. Each lattice site R contains
several orbitals that give rise to a bunch of electronic con-
figurations which we denote individually as |Γ;R〉. The
most general operator PR can be chosen of the form:
λ(R)ΓΓ′ |Γ;R〉〈Γ′;R|, (3)
where λ(R)ΓΓ′ are variational parameters. In general PR
needs not to be hermitean, namely for Γ 6= Γ′ it is not
required that λ(R)∗ΓΓ′ = λ(R)Γ′Γ. Indeed, as we shall
see, the non-hermitean character plays a very important
role. We further assume that the Wick’s theorem holds
for the uncorrelated wavefunction, hence that |φ〉 is either
a Slater determinant or a BCS wavefunction.
It was realized by Bünemann, Weber and Gebhard14
that average values of operators on the Gutzwiller wave-
function (1) can be analytically computed in infinite co-
ordination lattices provided the following two constraints
are imposed on PR:
〈φ| P†
|φ〉 = 〈φ|φ〉 = 1, (4)
〈φ| P†
CR |φ〉 = 〈φ| CR |φ〉, (5)
where CR is the local single-particle density-matrix op-
erator, with elements c
R,α cR,β and c
R,α c
, α labeling
single-particle states, while c
and c
create and an-
nihilate, respectively, an electron at site R in state α.
The first constraint, Eq. (4), does not actually limit
the variational freedom, since PR is defined up to a nor-
malization factor. On the contrary, the latter constraint,
Eq. (5), may reduce the variational freedom, although it
seems not in a relevant manner, at least in all cases that
we have so far investigated. We notice that Eq. (5) is not
the same as imposing
〈φ| P†
CRPR |φ〉 = 〈φ| CR |φ〉, (6)
unless P
commutes with CR, which is a further con-
straint to be imposed on PR. This actually is the only
case that has been hitherto considered, see e.g. Refs. 14
and 15. However, as we shall see, there are interesting
models which force to abandon the supplementary con-
dition (6), which is anyway unnecessary.17
By means of Wick’s theorem, the left-hand side of (5)
includes a disconnected term
〈φ| P†
|φ〉 〈φ| CR |φ〉 = 〈φ| CR |φ〉,
where the right-hand side follows from (4), plus con-
nected terms that are obtained by selecting in all possible
ways a pair of single-fermion operators from P†
, av-
eraging on |φ〉 what remains, and finally averaging the
two single-fermion operators with those of CR. There-
fore, imposing (5) means that the sum of all connected
terms vanishes, whatever is the element of the single-
particle density-matrix. In other words, the operator
that is left after taking out from P†
any pair of single-
fermion operators has null average on |φ〉. In turns, this
also implies that, when averaging on |φ〉 P†
multi-particle operators at different sites, the only con-
nected terms that survive are those that involve four or
R’+ + . . .
=R’R R
FIG. 1: Graphical representation of the average on |φ〉 of
PR, drawn as a box, times a generic multi-particle oper-
ator at site R′, drawn as a circle. Lines that join the two
operators represent the average of two single-fermion opera-
tors, one at R and the other at R′. The dots include all terms
where the two sites are joined by more than four lines. The
important thing to notice is the absence of terms in which the
two sites are connected by two lines.
more single-fermion operators of P†
, that are repre-
sented graphically in Fig. 1 as lines coming out of P†
This property of PR turns out to be extremely useful
in infinite-coordination lattices. In this limit, the con-
tribution to the average value on |φ〉 of terms in which
more than two fermionic lines come out of P†
can be
shown to vanish14, which simplifies considerably all cal-
culations. For instance, the average value on (1) of any
local operator OR becomes
〈φ| P†OR P |φ〉 = 〈φ| P†R OR PR|φ〉, (7)
which also implies, taking OR = 1, that the variational
wavefunction (1) is normalized. In addition, the average
value of the inter-site single-particle density matrix turns
out to be
〈φ| P† c†
R,α cR′,β P |φ〉
= 〈φ| P†
Z(R)αγ Z(R′)
〈φ| c†
R,γ cR′,δ |φ〉
Z(R)αγ ∆(R′)
〈φ| c†
∆(R)αγ Z(R′)
〈φ| c
R,γ cR′,δ |φ〉
∆(R)αγ ∆(R′)
〈φ| c
|φ〉, (8)
〈φ| P† c†
P |φ〉
= 〈φ| P†
R,α PRP
Z(R)αγ Z(R′)βδ 〈φ| c†R,γ c
Z(R)αγ ∆(R′)βδ 〈φ| c†R,γ cR′,δ |φ〉
∆(R)αγ Z(R′)βδ 〈φ| cR,γ c
∆(R)αγ ∆(R′)βδ 〈φ| cR,γ cR′,δ |φ〉, (9)
where the matrices Z and ∆ are determined by inverting
the following set of equations
〈φ| P†
R,α PR cR,β |φ〉 (10)
Z(R)αγ 〈φ| c†R,γ cR,β |φ〉
∆(R)αγ 〈φ| cR,γ cR,β |φ〉, (11)
〈φ| P†
|φ〉 (12)
Z(R)αγ 〈φ| c†R,γ c
∆(R)αγ 〈φ| cR,γ c
|φ〉, . (13)
Näıvely speaking, it is as if, when calculating the inter-
site density matrix, a fermionic operator transforms ef-
fectively into
Z(R)αβ c
∆(R)αβ cR,β, (14)
namely that a particle turns into a particle or a hole
with probabilities Z and ∆, respectively. Although all
the above expressions are strictly valid only in infinite-
coordination lattices, it is quite common to use the same
formulas also to evaluate average values on the Gutzwiller
wavefunction in finite-coordination lattices. This ap-
proximation is refereed to as the Gutzwiller approxima-
tion1,2,18,19,20, and is known to be equivalent to the sad-
dle point solution within the slave-boson technique.21
We conclude by noting that the constraint (5) turns
out to be useful also when the variational wavefunction
(1) is applied to Anderson impurity models. In this case
the operator P acts only on the electronic configurations
|Γ〉 of the impurity, namely
λΓΓ′ |Γ〉〈Γ′|.
If we impose
〈φ| P† P |φ〉 = 1, (15)
〈φ| P† P Cimp |φ〉 = 〈φ| Cimp |φ〉, (16)
where Cimp is the single-particle density matrix of the
impurity, then, for any operator of the conduction bath,
Obath, and because of (4) and (5), the following result
holds
〈φ| P† Obath P |φ〉 = 〈φ| Obath |φ〉. (17)
A. Some formal definitions
In order to perform actual calculations, it is convenient
to introduce some notations. We define a matrix Fα with
elements
(Fα)Γ1Γ2 = 〈Γ1;R| cR,α |Γ2;R〉,
as well as its hermitean conjugate, F †α, where we assumed
that the definition of the local configurations is the same
for all sites. It follows that
Fα = δαβ I,
Fα Fβ + Fβ Fα = 0,
where I is the identity. Next, we introduce the uncor-
related occupation-probability matrix, P0(R), with ele-
ments
(P0(R))Γ1Γ2 = 〈φ|Γ1;R〉 〈Γ2;R|φ〉, (18)
that satisfies
1 = Tr (P0(R)) ,
〈φ| c†
R,α cR,β |φ〉 = Tr
P0(R)F
〈φ| c†
|φ〉 = Tr
P0(R)F
Analogously, the variational parameters that define P
λ(R)Γ1Γ2 , are interpreted as elements of a matrix λ(R).
With these definitions, Eqs. (4) and (5) become
〈φ| P†
|φ〉 = Tr
P0(R)λ(R)
† λ(R)
〈φ| P†
|φ〉 = Tr
P0(R)λ(R)
† λ(R)F †α Fβ
= 〈φ| c†
R,α cR,β |φ〉,
〈φ| P†
R,α c
|φ〉 = Tr
P0(R)λ(R)
† λ(R)F †α F
= 〈φ| c†
that suggests to introduce a variational occupation-
probability matrix P (R) = P0(R)λ(R)
† λ(R) with ma-
trix elements
(P (R))Γ1Γ2 =
(P0(R))Γ1Γ3 λ(R)
λ(R)Γ4Γ2 ,
that must satisfy
Tr (P (R)) = 1, (20)
P (R)F †α Fβ
= 〈φ| c†
R,α cR,β |φ〉, (21)
P (R)F †α F
= 〈φ| c†
|φ〉, , (22)
Eqs. (20), (21) and (22) replace the constraints (4) and
(5). With these definitions, the matrices Z and ∆, see
Eqs. (11) and (13), are obtained by solving
P0(R)λ(R)
† F †α λ(R)Fβ
Z(R)αγ Tr
P0(R)F
∆(R)αγ Tr
P0(R)Fγ Fβ
, (23)
P0(R)λ(R)
† F †α λ(R)F
Z(R)αγ Tr
P0(R)F
∆(R)αγ Tr
P0(R)Fγ F
. (24)
The above equations simplify if one uses the natural
basis, namely the single-particle basis that diagonalizes
the density-matrix,
〈φ| c†
R,α cR,β|φ〉 = n(R)α δαβ ,
〈φ| c†
R,α c
|φ〉 = 0.
In this case
Z(R)αβ =
P0(R)λ(R)
† F †α λ(R)Fβ
n(R)β
, (25)
∆(R)αβ =
P0(R)λ(R)
† F †α λ(R)F
1− n(R)β
. (26)
Moreover, if one constructs the states |Γ;R〉 so that
P0(R) is diagonal
(P0(R))ΓΓ′ = δΓΓ′ P0(R; Γ),
(P (R))Γ1Γ2 = P0(R; Γ1)
λ(R)Γ3Γ2 , (27)
B. The variational problem
We are now in position to settle up the variational
problem. We consider a generic tight-binding Hamilto-
H = −
R,αcR′,β
E(R)ΓΓ′ |Γ;R〉〈Γ′;R|, (28)
where α and β stem for spin, orbital and lattice site in
the chosen unit cell, and the hermitean matrix E(R) with
elements E(R)ΓΓ′ may be also unit-cell dependent. The
average value of this Hamiltonian on the Gutzwiller wave-
function (1) in the limit of infinite coordination lattices
or, in finite coordination ones, within the Gutzwiller ap-
proximation, is
Evar = −
Z(R)αγ Z(R′)
〈φ| c†
Z(R)αγ ∆(R′)
〈φ| c†
∆(R)αγ Z(R′)
〈φ| c
R,γ cR′,δ |φ〉
∆(R)αγ ∆(R′)
〈φ| c
P0(R)λ(R)
† E(R)λ(R)
≡ Ehop + Eint. (29)
The last term depends only on the local properties of the
uncorrelated wavefunction |φ〉, specifically on the occupa-
tion probabilities P0(R). Therefore, for any given choice
of P0(R), the optimal |φ〉 that minimizes the variational
energy is the ground state of the Hamiltonian
Hvar = −
Z(R)αγ Z(R′)
Z(R)αγ ∆(R′)
∆(R)αγ Z(R′)
R,γ cR′,δ
∆(R)αγ ∆(R′)
µ(R)αβ c
R,αcR,β
ν(R)αβ c
+H.c.
, (30)
where the parameters µ(R)αβ and ν(R)αβ are Lagrange
multipliers to be determined by imposing that the ground
state has indeed the chosen P0(R). The last task is to
find the values of the variational parameters λ(R)ΓΓ′ as
well as of P0(R) for which the variational energy (29) is
minimum. We note that the variational Hamiltonian (30)
that has to be solved may include also inter-site pairing
terms, which are absent in the original Hamiltonian (28).
Analogously to other more conventional variational ap-
proaches, like Hartree-Fock theory, it is common to in-
terpret the single-particle spectrum of the variational
Hamiltonian (30) as an approximation of the true co-
herent spectrum of quasi-particles.22
III. THE MODEL
Let us now apply the variational wavefunction to spe-
cific models that are inspired by the valence-bond crys-
tal example we mentioned in the introduction, and where
the off-diagonal elements of the operator PR as well as
its non-hermitean character do play an important role.
Since the operator PR is built out of purely local prop-
erties, namely the available on-site electronic configura-
tions plus a variational guess for the uncorrelated on-site
single-particle density-matrix, a lot of preliminary results
can be extracted without even specifying how lattice-sites
are coupled together. Therefore we start our analysis
from defining some local properties and later we will con-
sider a specific lattice model.
A. The isolated dimer
The basic unit of the model we are going to investigate
consists of a dimer with Hamiltonian
Hdimer = −t⊥
1σc2σ +H.c.
(ni − 1)2
≡ H⊥ +HU , (31)
where 1 and 2 refer to the two sites of the dimer and ni,
i = 1, 2, is the on-site occupation number.
It is more convenient to work in the basis of the even
(bonding) and odd (anti-bonding) combinations, defined
through
ceσ =
(c1σ + c2σ) , coσ =
(c1σ − c2σ) .
and use this basis to built the available electronic con-
figurations, which we will denote as |n,Γ〉, with n that
refers to the number of electrons. The empty and the
fourfold occupied dimer states are denoted as |0〉 and
|4〉, respectively, while the singly-occupied states as
|1, e(o), σ〉 = c†
e(o)σ
and the states with 3 electrons as
|3, e(o), σ〉 = c†
e(o)σ
o(e)↑ c
o(e)↓ |0〉.
There are six doubly-occupied configurations. Two are
spin-singlets with two electrons in the even or in the odd
orbital, |2, e〉 and |2, o〉, respectively. When each orbital
is singly occupied, the two electrons form either a spin
triplet, |2, 1, Sz〉 with Sz = −1, 0, 1, or a spin singlet,
|2, 0〉. Since we are not going to consider variational so-
lutions that break spin-SU(2) symmetry, it is convenient
to define the projector operators
|1, e(o)〉〈1, e(o)| =
|1, e(o), σ〉〈1, e(o), σ|,
|3, e(o)〉〈3, e(o)| =
|3, e(o), σ〉〈3, e(o), σ|,
|2, 1〉〈2, 1| =
Sz=−1
|2, 1, Sz〉〈2, 1, Sz|.
The isolated-dimer ground state in the subspace with
two electrons is
|Ψ〉 = cos θ√
1↑c2↓ + c
2↑c1↓
sin θ√
1↑c1↓ + c
2↑c2↓
(cos θ + sin θ) |2, e〉
(cos θ − sin θ) |2, o〉,
where tan 2θ = 4t⊥/U and has energy
+ 4t2⊥. (32)
|Ψ〉 can be always rewritten in the form of a Gutzwiller
wavefunction. First of all, we needs to choose an uncor-
related wavefunction |φ〉. A natural choice might be the
ground state at U = 0, namely |2, e〉. Indeed |Ψ〉 can be
written as
|Ψ〉 = P |2, e〉,
where
P = |Ψ〉〈2, e| = 1√
(cos θ + sin θ) |2, e〉〈2, e|
(cos θ − sin θ) |2, o〉〈2, e|. (33)
and obviously satisfies both (4) and (5).
Another possibility, that we are also going to con-
sider in what follows, is to use an uncorrelated wavefunc-
tion that corresponds to a dimer in which the two sites
are only coupled by an intersite singlet-Cooper pairing,
namely with 〈c†1↑c
2↓〉 = 〈c
1↓〉 6= 0. In this case
|φ〉 = 1
|0〉+ |2, e〉 − |2, o〉 − |4〉
and, once again, the true ground state can be written as
|Ψ〉 = |Ψ〉〈φ| |φ〉 ≡ P |φ〉.
Already at this stage one can appreciate how important
is the role of the off-diagonal elements in P , especially
for large U/t⊥.
B. The non-isolated dimer: variational density
matrix
When the dimer is coupled to the rest of the system, in
order to built the operator P we need to specify an un-
correlated local single-particle density matrix based on a
variational guess of the uncorrelated wavefunction |φ〉. A
simple guess would be a magnetic wavefunction in which
the two sites of each dimer have opposite magnetization.
This choice is also the only one admitted by an Hartree-
Fock decomposition of the interaction term HU . How-
ever, a magnetic wavefunction is not the most suitable
choice to reproduce the limit of isolated dimers, which is
a collection of singlets.
Alternatively, one can consider a paramagnetic |φ〉 that
has built in the tendency of each dimer to lock into a
spin-singlet. This can be accomplished in two ways that
do not exclude each other. The first is to assume an
uncorrelated wavefunction with a huge splitting between
even and odd orbitals, namely with
〈φ| c†eσceσ |φ〉 ≫ no =
〈φ| c†oσcoσ |φ〉.
This implies that, among the doubly-occupied configura-
tions of each dimer, mainly the spin-singlet |2, e〉 survives
in the uncorrelated wavefunction. The latter can then be
turned into the isolated dimer configuration by an ap-
propriate Gutzwiller operator P , as shown before. The
other possibility is to include Cooper pairing correlations
in the singlet channel
∆SC = 〈φ| c†1↑c
2↓ + c
1↓ |φ〉.
In this case, the isolated dimer can be recovered by as-
suming a very strong pairing ∆SC ≃ 1 and suppressing,
through P , configurations with none or two singlet-pairs.
Note that both ne − no and ∆SC do not appear by a
mean-field decoupling of HU , so that a variational wave-
function with such correlations built in can not be sta-
bilized within Hartree-Fock theory. Here the role of P
becomes crucial.
Therefore, let us assume for |φ〉 a BCS-wavefunction
defined such that
〈φ| c†1σc1σ |φ〉 = 〈φ| c
2σc2σ |φ〉 =
, (34)
〈φ| c†1σc2σ |φ〉 = 〈φ| c
2σc1σ |φ〉 =
, (35)
〈φ| c†1↑c
2↓ |φ〉 = 〈φ| c
1↓ |φ〉 =
, (36)
with real ∆SC . In the even/odd basis this translates into
〈φ| c†eσceσ |φ〉 =
, (37)
〈φ| c†oσcoσ |φ〉 =
, (38)
〈φ| c†
e↓ |φ〉 = −〈φ| c
o↓ |φ〉 =
, (39)
where ne + no = n. As previously mentioned, the cal-
culations simplify considerably in the natural basis, that
is derived in the Appendix for this particular choice of
density matrices.
As a particular application, we assume hereafter that
the model is half-filled, namely ne +n0 = 2. The density
matrix of the operators in the natural basis, d
e(o)σ
e(o)σ
is, by Eqs. (A.3) and (A.4),
〈φ| d†
e(o)σ
e(o)σ
|φ〉 = 1
where
δ2 +∆2
. (40)
The two angles θe and θo, that are defined by Eq. (A.2)
and identify the unitary transformation from the original
to the natural basis, are given by θe = θ and θo = θ−π/2,
where
tan 2θ =
. (41)
We note that, for q → 1/2, the uncorrelated wavefunc-
tion describes an insulator where charge fluctuations are
completely suppressed since each natural orbital is fully
occupied. It is obvious that, if our choice of the varia-
tional wavefunction is correct, then the optimal uncor-
related wavefunction must asymptotically acquire
q = 1/2 for U → ∞.
The expression in the natural basis of the hopping,
Eq. (A.6), and interaction, Eq. (A.7), operators can be
derived through (A.5) and have a relatively simple ex-
pression at half-filling:
H⊥ = −t⊥
2 cos 2θ
|4̃〉〈4̃| − |0̃〉〈0̃|
+ cos 2θ
|3̃〉〈3̃| − |1̃〉〈1̃|
− sin 2θ
|1̃〉〈3̃|+H.c.
2 sin 2θ
|0̃〉〈2̃,+|+ |4̃〉〈2̃,+|+H.c.
|0̃〉〈0̃|+ |4̃〉〈4̃| − |0̃〉〈4̃| − |4̃〉〈0̃|
|2̃,+〉〈2̃,+|+ |2̃,−〉〈2̃,−|+ |2̃, 0〉〈2̃, 0|
|1̃〉〈1̃|+ |3̃〉〈3̃|
, (43)
where we have defined
|1̃(3̃)〉〈1̃(3̃)| = |1̃(3̃), e〉〈1̃(3̃), e|+ |1̃(3̃), o〉〈1̃(3̃), o|,
|1̃〉〈3̃| = |1̃, e〉〈3̃, e|+ |1̃, o〉〈3̃, o|,
|2̃,±〉 = 1√
|2̃, e〉 ± |2̃, o〉
and denoted the local configurations in the natural basis
as |ñ,Γ〉 to distinguish them from the analogous ones in
the original representation.
C. The Gutzwiller operator P
The most general Gutzwiller operator P should include
at least all the projectors |ñ,Γ〉〈ñ,Γ| as well as all the
off-diagonal operators |ñ,Γ〉〈ñ′,Γ′| that appear in the lo-
cal Hamiltonian, Eqs. (42) and (43). As we mentioned
before, our expectation is that the uncorrelated wave-
function which better connects to the large-U Mott in-
sulator should be identified by q → 1/2, in which locally
only the configurations |3̃〉 and |4̃〉 are occupied with non-
negligible probability. This suggests that P must include
at least those off-diagonal operators that would turn |4̃〉
into the isolated dimer ground state, namely |0̃〉〈4̃| and
|2̃,+〉〈4̃|. The latter forces to include also |1̃〉〈3̃|, as we
shall see.
Therefore we assume for P the following variational
ansatz:
λnΓ |ñ,Γ〉〈ñ,Γ|+ λ13 |1̃〉〈3̃|
+λ04 |0̃〉〈4̃|+ λ2+ 4 |2̃,+〉〈4̃|, (44)
with real λ’s. We define
P (n,Γ) = λ2nΓ P0(ñ,Γ), (45)
for all n 6= 3, 4, while, for n = 3, 4,
P (3) =
λ23 + λ
P0(3̃), (46)
P (4) =
λ24 + λ
04 + λ
P0(4̃). (47)
Then the conditions Eqs. (4) and (5) read
P (n,Γ) = 1, (48)
nP (n,Γ) = 2 + 4q, (49)
λ13 λ1
P0(3̃)P0(1̃) =
2λ2+ 4 λ2+
P0(4̃)P0(2̃,+). (50)
Here P0(ñ,Γ) are the occupation probabilities in the nat-
ural basis of the uncorrelated wavefunction. Specifically
P0(ñ,Γ) = gen,Γ
)en (
)4−en
where gen,Γ is the degeneracy of the configuration.
Eq. (50) guarantees that the anomalous averages
〈φ| P† P d†
e(o)↑d
e(o)↓ |φ〉
vanish in the natural basis, and explains why we have
included |1̃〉〈3̃| in (44). It is convenient to rewrite
P (4)
P0(4̃)
ui for i = 4, 04, 2 + 4,
P (3)
P0(3̃)
ui for i = 3, 13,
where u23 + u
13 = 1, which can be satisfied by choosing
u3 = cosψ and u13 = sinψ, and u
2+4 = 1. The
latter parameters can be expressed by means of another
unit vector v = (v1, v2, v3), through
(u4 + u04) ,
cos 2θ√
(u4 − u04)− sin 2θ u2+4,
sin 2θ√
(u4 − u04) + cos 2θ u2+4,
In terms of all the variational parameters, the P (n,Γ)’s,
θ, q, ψ and v, the average values per dimer of the in-
teraction, HU , and intra-dimer hopping, H⊥, are readily
found to be
EU = 〈φ| P†HU P |φ〉 =
P (3) + P (1)
P (0) + U
P (2,+) + P (2,−) + P (2, 0)
v22 + v
P (4), (52)
E⊥ = 〈φ| P†H⊥ P |φ〉 = −2 t⊥ δ∗, (53)
where the actual correlated values of the hybridization
and of the anomalous average are
2 δ∗ = 〈φ| P (ne − no) P |φ〉
= 4 v1 v2 P (4) + cos (2θ + 2ψ) P (3)
−2 cos 2θ P (0)− cos 2θ P (1), (54)
2∆∗ = 〈φ| P
e↓ − c
o↓ +H.c.
P |φ〉
= 4 v1 v3 P (4) + sin (2θ + 2ψ) P (3)
−2 sin2θ P (0)− sin 2θ P (1), (55)
We note that δ∗ and ∆∗ are mutually exclusive, namely
the choice of parameters that maximizes one of the two,
makes the other vanishing.
Upon the action of P , the single fermion operators
in the Nambu spinor representation transform effectively
e(o)↑
e(o)↓
Z +∆ e−i β τ2
e(o)↑
e(o)↓
, (57)
where τi, i = 1, 2, 3, are the Pauli matrices that act on
the Nambu spinor components,
tanβ =
. (58)
and, finally,
〈φ| P† d†
e(o)σ
e(o)σ
1− 4q2
P (0)P (1) +
P (1)P (2,+)
P (1)P (2,−) + 1
P (1)P (2, 0)
P (1)P (2, 1) + cosψ
P (3)P (2, 1)
+ cosψ
P (3)P (2, 0) + cosψ
P (3)P (2,−)
+ cosψ
P (3)P (2,+)
v1 cosψ + v2 cos (2θ + ψ)
+ v3 sin (2θ + ψ)
P (4)P (3)
, (59)
〈φ| P† d†
e(o)↑ P d
e(o)↓ |φ〉
1− 4q2
P (3)P (2,−)
P (3)P (2,+)
P (3)P (2, 1)
P (3)P (2, 0)
v1 sinψ − v2 sin (2θ + ψ)
+v3 cos (2θ + ψ)
P (4)P (3)
, (60)
with real
Z and
∆. Therefore, if the dimers are cou-
pled one to another by the single particle hopping term
R 6=R′
i,j=e,o
τ3 ΨR′,j , (61)
where
R,i↑, cR,i↓),
and Ψ its hermitean conjugate, the uncorrelated wave
function |φ〉 minimizes the effective hopping
Tvar = (Z +∆)
R 6=R′
i,j=e,o
−2 i β τ2 Ψ
under the condition that the local density matrix satisfies
Eqs.(37)-(39). One can readily show that this amounts to
find the ground state |φ〉 of the variational Hamiltonian
Hvar = (Z +∆)
R 6=R′
i,j=e,o
τ3 ΨR′,j
R,e τ3 ΨR,e −Ψ
R,o τ3 ΨR,o
R,e τ1 ΨR,e −Ψ
R,o τ1 ΨR,o
, (63)
with µ3 and µ1 such that
〈φ| Hvar |φ〉 + 4q µ3 cos(2θ + 2β) + 4q µ1 sin(2θ + 2β),
is maximum.
Before we consider specific lattice models, it is worth
re-deriving within this variational scheme the isolated-
dimer ground-state energy (32) at half-filling. For that
purpose, we take all P (n,Γ) zero but P (4) = 1. The
variational energy is simply
Evar = EU + E⊥ = U
v22 + v
− 4 t⊥ v1 v2.
The minimum under the constraint v · v = 1 is obtained
for v3 = 0 and exactly reproduces (32). We note that
the minimum energy is independent on θ, namely there
exists a continuous family of variational solutions with
equal energy parametrized by θ. However, in spite of the
fact that the uncorrelated wavefunction may describe a
superconductor, the actual value of the anomalous aver-
age ∆∗ = 0.
IV. A LATTICE MODEL OF DIMERS
As a particular application, let us consider the follow-
ing lattice model
H = −
tRR′ c
R,iσcR′,iσ +H.c.
(nR,i − 1)2
<RR′>
R,1σcR,2σ +H.c.
(ǫk − t⊥) c†k,eσck,eσ + (ǫk + t⊥) c
(nR,i − 1)2 , (64)
where nR,i =
and ǫk is the band dis-
persion induced by tRR′ , with half-bandwidth D. The
Hamiltonian (64) represents two Hubbard models cou-
pled by a single-particle hopping t⊥, each model being
defined on a lattice with coordination number z. As we
mentioned, the variational results that we have so far
derived are rigorous strictly speaking only if z → ∞,
although, in the spirit of the Gutzwiller approximation,
they can be used for generic z as well.
If U ≫ D, t⊥, (64) describes at half-filling a Mott insu-
lator which may be magnetic at t⊥ ≪ D, but is certainly
non-magnetic at t⊥ ≫ D, where the ground state reduces
essentially to a collection of singlets. For instance, in the
case of a Bethe lattice with nearest neighbor hopping,
the transition is at t⊥ = D/
8, value that is going to
decrease if frustration is included. If U is small and the
Fermi surface is not nested, then the model is metallic
for t⊥ ≤ D and is a band insulator otherwise. In fact,
in the absence of nesting there is generically a finite win-
dow of t⊥ values in which, upon increasing U , the model
undergoes a transition from a paramagnetic metal into a
non-magnetic Mott insulator, and this is just the case we
are going to consider in what follows. The same model
have been recently studied by Fuhrmann, Heilmann and
Monien using DMFT23 and by Kancharla and Okamoto24
using DMFT and cluster DMFT, respectively, that gives
us the opportunity to directly check the accuracy of our
wavefunction.
The variational Hamiltonian (63) of the model (64) has
a very simple expression,
Hvar =
ǫk∗ − µ3
τ3 − µ1 τ1
ǫk∗ + µ3
τ3 + µ1 τ1
, (65)
where
ǫk∗ = (Z +∆) ǫk.
The variational single-particle spectrum has the conven-
tional BCS form with eigenvalues
Eek =
(ǫk∗ − µ3)2 + µ21, Eok =
(ǫk∗ + µ3)
+ µ21,
hence, for any µ1 6= 0, has a gap equal to 2µ1. On the
contrary, when µ1 = 0, the spectrum is gapless for |µ3| ≤
D, otherwise is gaped. The Lagrange multipliers µ1 and
µ2 are obtaining by maximizing
Ehop = −
Eek + Eok
+ 4q µ3 cos(2θ + 2β) + 4q µ1 sin(2θ + 2β). (66)
In terms of (66), (52) and (53) the variational energy is
Evar = Ehop + EU + E⊥, (67)
and depends on eight independent variational parame-
ters.
We have solved numerically the variational problem at
fixed t⊥/D = 0.5 as function of U/D. To simplify cal-
culations, we have assumed for the band dispersion ǫk
either a flat or a semi-circular density of states, although
both would give rise to nesting that could stabilize mag-
netic phases, which we do not take into account. How-
ever, from the point of view of the paramagnetic-metal
to paramagnetic-insulator transition, this choice is not
influential.
We find that the variational solution displays a first
order phase transition at Uc ≃ 2.05 D for a flat density
of states, as shown by the behavior of the variational en-
ergy in Fig. 2. This result agrees almost quantitatively
with the DMFT calculation23 obtained with a semicir-
cular density of states, that also predicts a first order
transition with a coexistence region between U ≃ 1.5 D
and 1.8 D at the same value of t⊥ = 0.5 D. We note that
the energy is everywhere finite and vanishes like 1/U for
large U , see the inset of Fig. 2. The asymptotic behavior
UEvar/D
2 ∼ −7/6 is compatible with second order per-
turbation theory in t and t⊥ using as zeroth-order state
a collection of dimers, as explained below. In Fig. 3 we
show the behavior across the transition of the three con-
tribution to the energy, namely EU , E⊥ and Ehop. We
FIG. 2: Variational energy Evar in units of D and for a flat
non-interacting density of states, as function of U/D for t⊥ =
D/2. At Uc ≃ 2.05 D a first order transition occurs. The
inset shows the asymptotic value of UEvar/D
FIG. 3: The different contributions to the variational energy,
EU , E⊥ and Ehop.
find that the transition is accompanied by an energy loss
in Ehop, but a gain in both E⊥ and EU .
In order to characterize physically the two phases, in
Fig. 4 we plot the values of µ1 and µ3 across the tran-
sition. Since µ1 = 0, within our numerical accuracy,
and |µ3| < D, the phase at U < Uc is gapless hence
metallic, see the behavior of the density of states (DOS)
drawn in Fig. 5. On the contrary, on the U > Uc side
FIG. 4: The behavior of the parameters µ1 and µ3 in units of
D as function of U/D, see Eq. (65).
of the transition, µ1 6= 0, that implies a finite gap in the
single-particle variational spectrum, see Fig. 5. In the
gaped phase at U > Uc the spectrum looks like the one
of a Peierls insulator with a very large hybridization gap,
not consistent with the bare value of t⊥. In reality this
gap is, more properly, the Mott-Hubbard gap. Indeed
the DOS has weight both below and above the chemical
potential, suggestive of asymmetric Mott-Hubbard side-
bands. Moreover, as we are going to discuss below, the
actual difference between the occupations of the bond-
ing and anti-bonding bands, which we denoted as 2δ∗ in
Eq. (54), decreases with U , unlike the single-particle gap,
see Fig. 7. This behavior is reminiscent of what has been
found by Biermann et al.25 as an attempt to understand
the physics of VO2.
The other quantities that identify the variational spec-
trum are Z and ∆, shown in Fig. 6. We see that Z is
decreasing with U but reaches a finite value Z = 1/4 for
U → ∞. On the contrary, ∆ = 0 for U < Uc, while
∆ 6= 0 for U > Uc and increases monotonically to reach
asymptotically the same value 1/4 for large U .
Further insights can be gained by the average values of
the intra-dimer hopping and pairing, Eqs. (54) and (55),
drawn in Fig. 7. As we mentioned, the intra-dimer hy-
bridization 2δ∗ is monotonically decreasing with U , apart
from the jump at the first order transition. More inter-
estingly, around the transition the variational solution
is characterized by a sizeable BCS order parameter ∆∗,
which does not follow the behavior of the BCS coupling
µ1 present in the variational Hamiltonian. Indeed, while
µ1 is zero within our numerical accuracy for U < Uc,
yet a non negligible ∆∗ develops just before the transi-
tion, see Fig. 7. Moreover, although µ1 starts already
U/D = 2.0
U/D = 0
U/D = 2.05
U/D = 3.0
−2  0.5
FIG. 5: The variational single-particle spectrum for the even,
i.e. bonding, band, solid lines, and odd, i.e. anti-bonding,
one, dashed lines, across the transition for a non-interacting
semi-circular density of states.
large for U > Uc and increases monotonically with U ,
see Fig. 4, the actual order parameter ∆∗ is apprecia-
ble only near the transition and fastly decreases with
U to very tiny values. Also interesting is that, besides
the intra-dimer superconducting order parameter, also an
inter-dimer one arises, that can be for instance defined
through
∆̃∗ =
tRR′ 〈φ| P† Ψ†R,i τ1 ΨR′,i P |φ〉
FIG. 6: The behavior of the parameters Z and ∆.
FIG. 7: The average values of the intra-dimer hopping, δ∗,
and pairing, ∆∗.
4ǫk=0
where V is the number of sites. We find that ∆̃∗ has
actually the opposite sign of ∆∗, and both closely follow
each other, rapidly decreasing with U , see Fig. 8. Even
though ∆∗ and ∆̃∗ are everywhere finite for U > Uc, sug-
gestive of a superconducting phase that survives up to
very large U , we believe that superconductivity, it it oc-
curs at all, may appear only very close to the transition,
where the value of the order parameter is larger. Indeed,
FIG. 8: The intra-dimer, circles, and, with reversed sign, the
inter-dimer, triangles, superconducting order parameters.
the optimal solution with finite ∆∗ and ∆̃∗ and another
solution in which both are forced to be zero are practi-
cally degenerate within our numerical accuracy for large
U . Moreover, since the variational values of the order
parameter are extremely small but close to Uc, the in-
clusion of quantum fluctuations, for instance in the form
of a Jastrow factor as in Eq. (2), would likely suppress
superconductivity leading to a bona fide insulating wave-
function.
Unfortunately, we can not compare this result with the
DMFT analyses of Refs. 23,24, where superconductivity
has not been looked for.26
A. Large U limit
In order to appreciate qualities and also single out
defects of the variational wavefunction, it is worth dis-
cussing the large-U solution. To leading order in 1/U ,
one can assume all P (n,Γ) = 0 but P (3) and P (4). The
two constraints Eqs. (48) and (49) can be solved by defin-
P (3) = 2− 4q ≡ 4d, P (4) = 1− 4d,
with d ≪ 1, namely q → 1/2, while Eq. (50) is already
satisfied since P (1) = P (2,+) = 0. Moreover, one readily
recognizes that the variational solution asymptotically
tends to acquire ψ ≃ β → π/4, θ → 0 and v1 ≫ v2, v3,
which is indeed what we find by numerical minimization.
It then follows that
∆ → 1
1− 4q2
4d(1− 4d) → 1
This implies that µ3 → 0, hence that µ1 is determined
by maximizing
Ehop = −
k∗ + µ
1 + 4q µ1,
which leads to µ1 =
ǫ2/4d and
Ehop = −4
ǫ2 d, (68)
where
At leading order, the variational energy per dimer is
therefore
Evar = −4
ǫ2 d− 4 t⊥ v1 v2 + U v22 + 2U d,
that is minimized by d = ǫ2/U
2, v1 ≃ 1 and v2 ≃ 2t⊥/U ,
and takes the value
Evar = −
= − 1
tRR′ tR′R
. (69)
In the case of a flat density of states
δ (ǫ− ǫk) =
θ (D − |ǫ|) ,
and with t⊥ = 0.5 D we recover the numerical result
Evar ≃ −7/6 D2/U , see Fig. 2.
We note that, in spite of the hybridization δ∗ ∼ t⊥/U
being small, the single-particle gap of the variational
spectrum 2µ1 ≃ U is large, as one should expect in a
Mott insulator.
Coming back to the large-U value of the variational
energy (69), one can readily realize that it coincides with
the second order correction in tRR′ to the energy of the
state
|Ψ〉 =
R,1↑c
R,2↓ + c
R,2↑c
which is just a collection of singlets. In other words, in
spite of being non-magnetic, our variational wavefunc-
tion is able to reproduce the correct super-exchange be-
tween dimers. This is a remarkable property that ac-
tually derives from the square-root dependence upon d
of Ehop, see Eq. (68). If we considered a more conven-
tional Gutzwiller operator P commuting with the single-
particle density matrix, that amounts to further impose
δ2∗ +∆
∗, we would find Ehop ∝ d, implying a
transition into an unrealistic insulator with d = 0 above
a critical U . The obvious defect of the wavefunction is
that, since it emphasizes strongly the role of individual
dimers, the hopping among dimers, although finite for
any U , is under-estimated with respect to the intra-dimer
one. Therefore we do not expect the wavefunction to be
particularly accurate for small t⊥/D.
anti−bonding
bonding
FIG. 9: The non-interacting density of states of the lattice
of dimers. The bonding and anti-bonding state of each dimer
give rise to two bands that overlap, leading to a metallic phase
in the absence of interaction.
V. CONCLUSIONS
In this work we have proposed an extension of the
Gutzwiller variational approach to account for correlated
models which display metal-insulator transitions into
Mott insulators that escape any simple single-particle
descriptions, like the Hartree-Fock approximation. The
wavefunction has still the same form as the conventional
Gutzwiller wavefunction,
|ΨG〉 = P |φ〉 =
PR |φ〉,
with R identifying unit cells that may also be non-
primitive ones, with the novel feature that the opera-
tor PR is non-hermitean and does not commute with the
local single-particle density matrix. In essence, this prop-
erty realizes a variational implementation of a Schrieffer-
Wolff transformation27, although only restricted to the
lattice sites within each unit cell. We have shown that, by
slightly reducing the variational freedom, this wavefunc-
tion, like the conventional Gutzwiller wavefunction,14 al-
lows for an extension of the Gutzwiller approximation
to evaluate average values, approximation that becomes
exact in the limit of infinite coordination lattices.
As an application, we have considered the Mott tran-
sition into a Peierls, or valence-bond, insulator, namely
an insulator that is adiabatically connected to a collec-
tion of independent dimers. Such an insulator can not
be described by Hartree-Fock simply because the singlet
configuration of each dimer is not a Slater determinant.
Specifically, we have considered the hypothetical situa-
tion shown in Fig. 9, where the splitting between the
bonding and antibonding orbitals of each dimer is as-
sumed not to be sufficient to lead, in the absence of in-
teraction, to a band insulator. When interaction is taken
into account, in the form of an on-site repulsion U , one
expects, above a critical U , a transition from the metal
into a Mott insulator. If magnetism is prevented, for
instance by a sufficiently large splitting between bond-
ing and anti-bonding orbitals and/or by frustration, the
Mott insulator is non-magnetic. We have shown that
our wavefunction overcomes the difficulties of Hartree-
Fock theory and allows to study, albeit variationally, this
transition. In particular we find that:
(i) at the variational level the Mott transition is first
order;
(ii) the variational spectrum inside the Mott insulator
looks similar to that of a Peierls insulator with a
large hybridization gap, namely a large splitting be-
tween bonding and anti-bonding bands. In reality
the gap is the Mott-Hubbard gap and the actual
difference between the occupations of the bonding
and anti-bonding bands, is small;
(iii) inter-site singlet-superconductivity appears around
the transition.
While (i) and (ii) are presumably true, as they have
been also found by more rigorous calculations23,25,28, the
emergence of superconductivity might be an artifact of
the variational wavefunction.26 Nevertheless, the possi-
ble occurrence of superconductivity is quite suggestive.
It is known for instance that two-leg Hubbard ladders
with nearest neighbor hopping display dominant super-
conducting fluctuations with the same symmetry that we
find variationally29, although at half-filling they always
describe non-magnetic spin-gaped insulators30,31 because
of nesting. Moreover, the uncorrelated wavefunction |φ〉
is quite similar to the wavefunctions used in Refs. 32,33
to simulate t-J ladders. It would be surprising and inter-
esting if this tendency towards superconductivity turned
into a true symmetry breaking instability in higher di-
mensionality, as suggested by our analysis, which we
think it is worth deserving further investigations.
Note added: During the completion of this work, we be-
came aware of a recent extension of slave-boson technique
whose saddle-point solution closely resembles our varia-
tional approach.34 Indeed the two conditions we impose
on the Gutzwiller operator, Eqs. (4) and (5), are in one-
to-one correspondence with the constraints identified in
Ref. 34 within the slave-boson formalism.
Acknowledgments
We are grateful to C. Castellani and E. Tosatti for their
helpful comments and suggestions. We also thanks A.
Georges for useful discussions in connection with Ref. 34.
APPENDIX: THE NATURAL BASIS
Let us assume that, in the Nambu-spinor representa-
the uncorrelated wavefunction has the following density
matrices
Ĉe =
ne/2 ∆SC/2
∆SC/2 1− ne/2
, Ĉo =
no/2 −∆SC/2
−∆SC/2 1− no/2
(A.1)
The natural orbitals are obtained by the unitary trans-
formation
e(o)↑ = cos θe(o) ce(o)↑ + sin θe(o) c
e(o)↓
e(o)↓ = cos θe(o) ce(o)↓ − sin θe(o) c
e(o)↑,
where
tan 2θe =
ne − 1
, tan 2θo =
no − 1
, (A.2)
and posses a diagonal density matrix with the non-
vanishing elements given by
〈φ| d†
e(o)σ
e(o)σ
|φ〉 = 1
+ qe(o), (A.3)
where
qe(o) =
(ne(o) − 1)2 +∆2SC . (A.4)
In the natural basis we introduce states that have the
same formal expression as in the original basis but are
built with d-operators, and denote them as |ñ,Γ〉. The
transformation rules from these states to the original ones
|0〉 = cos θe cos θo|0̃〉+ sin θe sin θo |4̃〉
+cos θe sin θo |2̃, o〉+ sin θe cos θo |2̃, e〉,
|1, e(o), σ〉 = cos θo(e) |1̃, e(o), σ〉 + sin θo(e) |3̃, e(o), σ〉,
|2, e(o)〉 = cos θe cos θo |2̃, e(o)〉+ cos θe(o) sin θo(e) |4̃〉
− sin θe(o) cos θo(e) |0̃〉
− sin θe(o) sin θo(e) |2̃, o(e)〉,
|2, 1, Sz〉 = |2̃, 1, Sz〉, (A.5)
|2, 0〉 = |2̃, 0〉,
|3, e(o), σ〉 = cos θo(e) |3̃, e(o)〉 − sin θo(e) |1̃, e(o)σ〉,
|4〉 = cos θe cos θo|4̃〉+ sin θe sin θo |0̃〉
− cos θe sin θo |2̃, e〉 − sin θe cos θo |2̃, o〉.
The inverse transformation is obtained by letting θe(o) →
−θe(o).
The hopping operator in the original representation
can be written as
H⊥ = −2t⊥
1σc2σ +H.c.
= −2t⊥
c†eσceσ − c†eσceσ
= −2t⊥
|3, o, σ〉〈3, o, σ| − |3, e, σ〉〈3, e, σ|
+|1, e, σ〉〈1, e, σ| − |1, o, σ〉〈1, o, σ|
+2 |2, e〉〈2, e| − 2 |2, o〉〈2, o|
, (A.6)
while the interaction operator as
(ni − 1)2 =
|0〉〈0|+ |4〉〈4|
n=1,3
|n, e, σ〉〈n, e, σ|+ |n, o, σ〉〈n, o, σ|
|2, e〉+ |2, o〉
〈2, e|+ 〈2, o|
+2 |2, 0〉〈2, 0|
. (A.7)
Their expression in the natural basis can be obtained by
the transformation rules (A.5).
1 M. C. Gutzwiller, Phys. Rev. 134, A923 (1964).
2 M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
3 W. F. Brinkman and T. M. Rice, Phys. Rev. B 2, 4302
(1970).
4 M. Capello, F. Becca, M. Fabrizio, S. Sorella, and
E. Tosatti, Phys. Rev. Lett. 94, 026406 (2005).
5 M. Capello, F. Becca, S. Yunoki, M. Fabrizio, and
S. Sorella, Phys. Rev. B 72, 085121 (2005).
6 M. Capello, F. Becca, M. Fabrizio, and S. Sorella, Super-
fluid to Mott-insulator transition in Bose-Hubbard models
(2006),
7 S. Sorella, Phys. Rev. B 71, 241103 (2005).
8 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,
Rev. Mod. Phys. 68, 13 (1996).
9 A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 62,
R9283 (2000).
10 D. Senechal, D. Perez, and M. Pioro-Ladriere, Phys. Rev.
Lett. 84, 522 (2000).
11 M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev.
Lett. 91, 206402 (2003).
12 T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev.
Mod. Phys. 77, 1027 (2005).
13 G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys.
Rev. Lett. 87, 186401 (2001).
14 J. Bünemann, W. Weber, and F. Gebhard, Phys. Rev. B
57, 6896 (1998).
15 C. Attaccalite and M. Fabrizio, Phys. Rev. B 68, 155117
(2003).
16 Q.-H. Wang, Z. D. Wang, Y. Chen, and F. C. Zhang, Phys-
ical Review B (Condensed Matter and Materials Physics)
73, 092507 (pages 4) (2006).
17 M. Ferrero, Ph.D. thesis, SISSA-Trieste (2006).
18 W. Metzner and D. Vollhardt, Phys. Rev. Lett. 59, 121
(1987).
19 W. Metzner and D. Vollhardt, Phys. Rev. B 37, 7382
(1988).
20 F. Gebhard, Phys. Rev. B 41, 9452 (1990).
21 G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57,
1362 (1986).
22 J. Bünemann, F. Gebhard, and R. Thul, Phys. Rev. B 67,
075103 (2003).
23 A. Fuhrmann, D. Heilmann, and H. Monien, Physical Re-
view B 73, 245118 (2006).
24 S. S. Kancharla and S. Okamoto, Band
insulator to Mott insulator transition in
a bilayer Hubbard model (2007), URL
http://www.citebase.org/abstract?id=oai:arXiv.org:cond-mat/0703728.
25 S. Biermann, A. Poteryaev, A. I. Lichtenstein, and
A. Georges, Physical Review Letters 94, 026404 (pages 4)
(2005).
26 There are however recent DMFT calculations that seem to
support the existence of a superconducting region around
the transition. [A. Privitera, M. Schiro’, M. Capone and
C. Castellani, unpublished (2007).].
27 J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491
(1966).
28 G. Moeller, V. Dobrosavljević, and A. E. Ruckenstein,
Phys. Rev. B 59, 6846 (1999).
29 M. Fabrizio, Phys. Rev. B 48, 15838 (1993).
30 S. P. Strong and A. J. Millis, Phys. Rev. B 50, 9911 (1994).
31 D. G. Shelton, A. A. Nersesyan, and A. M. Tsvelik, Phys.
Rev. B 53, 8521 (1996).
32 G. Sierra, M. A. Mart́ın-Delgado, J. Dukelsky, S. R. White,
and D. J. Scalapino, Phys. Rev. B 57, 11666 (1998).
33 S. Sorella, G. B. Martins, F. Becca, C. Gazza, L. Capriotti,
A. Parola, and E. Dagotto, Phys. Rev. Lett. 88, 117002
(2002).
34 F. Lechermann, A. Georges, G. Kotliar, and O. Parcollet,
unpublished (2007).
http://www.citebase.org/abstract?id=oai:arXiv.org:cond-mat/0703728

  We introduce a novel extension of the Gutzwiller variational wavefunction
able to deal with insulators that escape any mean-field like description, as
for instance non-magnetic insulators. As an application, we study the Mott
transition from a paramagnetic metal into a non-magnetic Peierls, or
valence-bond, Mott insulator. We analyze this model by means of our Gutzwiller
wavefunction analytically in the limit of large coordination lattices, where we
find that: (1) the Mott transition is first order; (2) the Peierls gap is large
in the Mott insulator, although it is mainly contributed by the electron
repulsion; (3) singlet-superconductivity arises around the transition.

Gutzwiller description of non-magnetic Mott insulators: a dimer lattice model.
Michele Fabrizio,1,2
1 International School for Advanced Studies (SISSA) and CNR-INFM-Democritos
National Simulation Centre, Via Beirut 2-4, I-34014 Trieste, Italy
2 The Abdus Salam International Centre for Theoretical Physics (ICTP), P.O. Box 586, I-34014 Trieste, Italy
(Dated: September 12, 2021)
We introduce a novel extension of the Gutzwiller variational wavefunction able to deal with
insulators that escape any mean-field like description, as for instance non-magnetic insulators. As
an application, we study the Mott transition from a paramagnetic metal into a non-magnetic Peierls,
or valence-bond, Mott insulator. We analyze this model by means of our Gutzwiller wavefunction
analytically in the limit of large coordination lattices, where we find that: (1) the Mott transition
is first order; (2) the Peierls gap is large in the Mott insulator, although it is mainly contributed by
the electron repulsion; (3) singlet-superconductivity arises around the transition.
PACS numbers: 71.10.-w, 71.10.Fd, 71.30.+h
I. INTRODUCTION
Among the theoretical tools devised to deal with
strongly correlated metals close to a Mott metal-to-
insulator transition (MIT), the simplest one likely is
the variational approach introduced in the 60ths’ by
Gutzwiller1,2 to describe itinerant ferromagnetism and
narrow band conductors. In its original version, the
Gutzwiller variational wavefunction has the form
|ΨG〉 = P |φ〉 =
PR |φ〉, (1)
where |φ〉 is an uncorrelated wavefunction for which
Wick’s theorem holds, PR an operator at site R, and
both |φ〉 and PR have to be determined by minimizing
the variational energy. The role of the operator PR is to
modify, according to the on-site interaction, the weights
of the local electronic configurations with respect to their
values in the uncorrelated wavefunction.
In spite of its simplicity, the Gutzwiller wavefunction
is quite effective in capturing physical properties that
supposedly identify strongly correlated metals, as for in-
stance the large increase of the effective mass.3 However,
since the dependence upon the distance |R−R′| of inter-
site correlations are still determined by the uncorrelated
wavefunction, while the local operators PR just affect the
amplitudes, the Gutzwiller wavefunction can describe a
Mott insulator either if PR suppresses completely charge
fluctuations, that provides a very poor description of an
insulator, or if |φ〉 itself is insulating. The latter case
can be stabilized within the original Gutzwiller approach
only when |φ〉 is an admissible Hartree-Fock solution of
the Hamiltonian. As an example let us consider a sin-
gle band model at half-filling, for instance the Hubbard
model
H = −
RR′,σ
tRR′ c
nR,↑ nR,↓,
where c
and c
creates and annihilates, respectively,
an electron with spin σ =↑, ↓ at site R and nR,σ =
. This Hamiltonian admits at the mean-field
level two possible phases, one paramagnetic, 〈nR,↑〉 =
〈nR,↓〉, and the other magnetic, 〈nR,↑〉 6= 〈nR,↓〉. The
latter is the only one that can eventually describe an
insulator. In the Hubbard model the action of the oper-
ator P is to increase the weight of singly occupied sites
at expenses of doubly occupied and empty sites, in or-
der to minimize the Coulomb repulsion U . Evidently,
even when the repulsion is very strong, hence the model
is a Mott insulator, a realistic wavefunction should still
allow for charge fluctuations responsible for the super-
exchange, that survives even deep inside the Mott phase.
However, any paramagnetic uncorrelated wavefunction,
for instance the Fermi sea, is unable to generate any
super-exchange and necessarily leads to a non-realistic
Mott insulator where configurations with empty or dou-
bly occupied sites are fully suppressed.3 The only way
to generate super-exchange is to assume a magnetically
ordered |φ〉, which is also the only insulating wavefunc-
tion accessible within Hartree-Fock. However a magnetic
state might not always be the right choice, especially if
magnetism is sufficiently frustrated.
Recently, an improved version of the Gutzwiller wave-
function has been proposed,4 in which additional inter-
site correlations are provided by density-density Jastrow
factors, namely
|ΨG〉 → exp
vR,R′ nR nR′
|ΨG〉, (2)
where nR is the site R occupation number and vR,R′
variational parameters. This novel class of wavefunctions
has the capability to disentangle charge from other de-
grees of freedom, hence is more suitable to capture Mott
localization, as it has indeed been shown.4,5,6 However,
unlike the conventional Gutzwiller wavefunction (1), the
Gutzwiller-Jastrow wavefunction (2) can only be dealt
with numerically by variational Monte Carlo, which is in-
herently limited to finite-size systems, albeit quite large.7
An alternative approach, that is closely related to
recently proposed extensions of Dynamical Mean Field
Theory (DMFT) from the original single-site formula-
http://arxiv.org/abs/0704.1389v1
tion8 to a cluster one9,10,11,12,13, is to consider a vari-
ational wavefunction of the same form as (1) but de-
fined on a lattice with non-primitive unit cells. In this
case, the operator PR acts on all the available electronic
configurations of the lattice sites belonging to the non-
primitive cell. The advantage is that in this way one may
include additional short-range correlations without losing
the property of the wavefunction to be analytically man-
ageable, at least in infinite-coordination lattices. The ob-
vious disadvantage is that this wavefunction could bias
the variational solution towards translational-symmetry
breaking.
Within this scheme, the variational problem becomes
generically equivalent to optimize a Gutzwiller wave-
function for a multi-band Hamiltonian. There have
been recently an amount of attempts to extend the
Gutzwiller wavefunction to multi-orbital models that in-
clude further complications like for instance Coulomb ex-
change14,15,16,17. In this paper we introduce a further
extension that is capable to generate inter-site correla-
tions as the super-exchange for paramagnetic wavefunc-
tions, otherwise missed by the conventional Gutzwiller
approach. This novel class of wavefunctions also allows
to explore new kinds of variational solutions. Specifi-
cally, there are interesting examples of correlated models
where the Mott insulating phase escapes any Hartree-
Fock mean-field treatment, in other words can not be
represented by a single Slater determinant. A very sim-
ple case, that we will explicitly consider throughout this
work, is a Peierls insulator, namely a short-range valence-
bond crystal, in which pairs of nearest neighbor sites are
strongly bound into a singlet configuration, leading to a
state that is simply a collection of spin-singlets. Such a
Mott insulating state is not accessible by Hartree-Fock
theory, just because each singlet is itself not express-
ible as a Slater determinant, nor by the conventional
Gutzwiller approach, which, as mentioned, gives a poor
description of paramagnetic insulators.
The paper is organized as follows. In Section II we
present the variational wavefunction and discuss under
which conditions it can be deal with analytically. In Sec-
tion III we discuss how to build up the wavefunction in
the case in which the basic unit of the lattice model is
a dimer. Next, in Section IV, we solve the variational
problem for a specific lattice model of dimers. Conclu-
sions are given in Section V.
II. THE VARIATIONAL WAVEFUNCTION
In this Section, we introduce an extension of the
Gutzwiller wavefunction (1) which is particularly con-
venient to perform analytical calculation in the limit of
infinite-coordination lattices.14 Let us consider a generic
multi-band Hamiltonian. Each lattice site R contains
several orbitals that give rise to a bunch of electronic con-
figurations which we denote individually as |Γ;R〉. The
most general operator PR can be chosen of the form:
λ(R)ΓΓ′ |Γ;R〉〈Γ′;R|, (3)
where λ(R)ΓΓ′ are variational parameters. In general PR
needs not to be hermitean, namely for Γ 6= Γ′ it is not
required that λ(R)∗ΓΓ′ = λ(R)Γ′Γ. Indeed, as we shall
see, the non-hermitean character plays a very important
role. We further assume that the Wick’s theorem holds
for the uncorrelated wavefunction, hence that |φ〉 is either
a Slater determinant or a BCS wavefunction.
It was realized by Bünemann, Weber and Gebhard14
that average values of operators on the Gutzwiller wave-
function (1) can be analytically computed in infinite co-
ordination lattices provided the following two constraints
are imposed on PR:
〈φ| P†
|φ〉 = 〈φ|φ〉 = 1, (4)
〈φ| P†
CR |φ〉 = 〈φ| CR |φ〉, (5)
where CR is the local single-particle density-matrix op-
erator, with elements c
R,α cR,β and c
R,α c
, α labeling
single-particle states, while c
and c
create and an-
nihilate, respectively, an electron at site R in state α.
The first constraint, Eq. (4), does not actually limit
the variational freedom, since PR is defined up to a nor-
malization factor. On the contrary, the latter constraint,
Eq. (5), may reduce the variational freedom, although it
seems not in a relevant manner, at least in all cases that
we have so far investigated. We notice that Eq. (5) is not
the same as imposing
〈φ| P†
CRPR |φ〉 = 〈φ| CR |φ〉, (6)
unless P
commutes with CR, which is a further con-
straint to be imposed on PR. This actually is the only
case that has been hitherto considered, see e.g. Refs. 14
and 15. However, as we shall see, there are interesting
models which force to abandon the supplementary con-
dition (6), which is anyway unnecessary.17
By means of Wick’s theorem, the left-hand side of (5)
includes a disconnected term
〈φ| P†
|φ〉 〈φ| CR |φ〉 = 〈φ| CR |φ〉,
where the right-hand side follows from (4), plus con-
nected terms that are obtained by selecting in all possible
ways a pair of single-fermion operators from P†
, av-
eraging on |φ〉 what remains, and finally averaging the
two single-fermion operators with those of CR. There-
fore, imposing (5) means that the sum of all connected
terms vanishes, whatever is the element of the single-
particle density-matrix. In other words, the operator
that is left after taking out from P†
any pair of single-
fermion operators has null average on |φ〉. In turns, this
also implies that, when averaging on |φ〉 P†
multi-particle operators at different sites, the only con-
nected terms that survive are those that involve four or
R’+ + . . .
=R’R R
FIG. 1: Graphical representation of the average on |φ〉 of
PR, drawn as a box, times a generic multi-particle oper-
ator at site R′, drawn as a circle. Lines that join the two
operators represent the average of two single-fermion opera-
tors, one at R and the other at R′. The dots include all terms
where the two sites are joined by more than four lines. The
important thing to notice is the absence of terms in which the
two sites are connected by two lines.
more single-fermion operators of P†
, that are repre-
sented graphically in Fig. 1 as lines coming out of P†
This property of PR turns out to be extremely useful
in infinite-coordination lattices. In this limit, the con-
tribution to the average value on |φ〉 of terms in which
more than two fermionic lines come out of P†
can be
shown to vanish14, which simplifies considerably all cal-
culations. For instance, the average value on (1) of any
local operator OR becomes
〈φ| P†OR P |φ〉 = 〈φ| P†R OR PR|φ〉, (7)
which also implies, taking OR = 1, that the variational
wavefunction (1) is normalized. In addition, the average
value of the inter-site single-particle density matrix turns
out to be
〈φ| P† c†
R,α cR′,β P |φ〉
= 〈φ| P†
Z(R)αγ Z(R′)
〈φ| c†
R,γ cR′,δ |φ〉
Z(R)αγ ∆(R′)
〈φ| c†
∆(R)αγ Z(R′)
〈φ| c
R,γ cR′,δ |φ〉
∆(R)αγ ∆(R′)
〈φ| c
|φ〉, (8)
〈φ| P† c†
P |φ〉
= 〈φ| P†
R,α PRP
Z(R)αγ Z(R′)βδ 〈φ| c†R,γ c
Z(R)αγ ∆(R′)βδ 〈φ| c†R,γ cR′,δ |φ〉
∆(R)αγ Z(R′)βδ 〈φ| cR,γ c
∆(R)αγ ∆(R′)βδ 〈φ| cR,γ cR′,δ |φ〉, (9)
where the matrices Z and ∆ are determined by inverting
the following set of equations
〈φ| P†
R,α PR cR,β |φ〉 (10)
Z(R)αγ 〈φ| c†R,γ cR,β |φ〉
∆(R)αγ 〈φ| cR,γ cR,β |φ〉, (11)
〈φ| P†
|φ〉 (12)
Z(R)αγ 〈φ| c†R,γ c
∆(R)αγ 〈φ| cR,γ c
|φ〉, . (13)
Näıvely speaking, it is as if, when calculating the inter-
site density matrix, a fermionic operator transforms ef-
fectively into
Z(R)αβ c
∆(R)αβ cR,β, (14)
namely that a particle turns into a particle or a hole
with probabilities Z and ∆, respectively. Although all
the above expressions are strictly valid only in infinite-
coordination lattices, it is quite common to use the same
formulas also to evaluate average values on the Gutzwiller
wavefunction in finite-coordination lattices. This ap-
proximation is refereed to as the Gutzwiller approxima-
tion1,2,18,19,20, and is known to be equivalent to the sad-
dle point solution within the slave-boson technique.21
We conclude by noting that the constraint (5) turns
out to be useful also when the variational wavefunction
(1) is applied to Anderson impurity models. In this case
the operator P acts only on the electronic configurations
|Γ〉 of the impurity, namely
λΓΓ′ |Γ〉〈Γ′|.
If we impose
〈φ| P† P |φ〉 = 1, (15)
〈φ| P† P Cimp |φ〉 = 〈φ| Cimp |φ〉, (16)
where Cimp is the single-particle density matrix of the
impurity, then, for any operator of the conduction bath,
Obath, and because of (4) and (5), the following result
holds
〈φ| P† Obath P |φ〉 = 〈φ| Obath |φ〉. (17)
A. Some formal definitions
In order to perform actual calculations, it is convenient
to introduce some notations. We define a matrix Fα with
elements
(Fα)Γ1Γ2 = 〈Γ1;R| cR,α |Γ2;R〉,
as well as its hermitean conjugate, F †α, where we assumed
that the definition of the local configurations is the same
for all sites. It follows that
Fα = δαβ I,
Fα Fβ + Fβ Fα = 0,
where I is the identity. Next, we introduce the uncor-
related occupation-probability matrix, P0(R), with ele-
ments
(P0(R))Γ1Γ2 = 〈φ|Γ1;R〉 〈Γ2;R|φ〉, (18)
that satisfies
1 = Tr (P0(R)) ,
〈φ| c†
R,α cR,β |φ〉 = Tr
P0(R)F
〈φ| c†
|φ〉 = Tr
P0(R)F
Analogously, the variational parameters that define P
λ(R)Γ1Γ2 , are interpreted as elements of a matrix λ(R).
With these definitions, Eqs. (4) and (5) become
〈φ| P†
|φ〉 = Tr
P0(R)λ(R)
† λ(R)
〈φ| P†
|φ〉 = Tr
P0(R)λ(R)
† λ(R)F †α Fβ
= 〈φ| c†
R,α cR,β |φ〉,
〈φ| P†
R,α c
|φ〉 = Tr
P0(R)λ(R)
† λ(R)F †α F
= 〈φ| c†
that suggests to introduce a variational occupation-
probability matrix P (R) = P0(R)λ(R)
† λ(R) with ma-
trix elements
(P (R))Γ1Γ2 =
(P0(R))Γ1Γ3 λ(R)
λ(R)Γ4Γ2 ,
that must satisfy
Tr (P (R)) = 1, (20)
P (R)F †α Fβ
= 〈φ| c†
R,α cR,β |φ〉, (21)
P (R)F †α F
= 〈φ| c†
|φ〉, , (22)
Eqs. (20), (21) and (22) replace the constraints (4) and
(5). With these definitions, the matrices Z and ∆, see
Eqs. (11) and (13), are obtained by solving
P0(R)λ(R)
† F †α λ(R)Fβ
Z(R)αγ Tr
P0(R)F
∆(R)αγ Tr
P0(R)Fγ Fβ
, (23)
P0(R)λ(R)
† F †α λ(R)F
Z(R)αγ Tr
P0(R)F
∆(R)αγ Tr
P0(R)Fγ F
. (24)
The above equations simplify if one uses the natural
basis, namely the single-particle basis that diagonalizes
the density-matrix,
〈φ| c†
R,α cR,β|φ〉 = n(R)α δαβ ,
〈φ| c†
R,α c
|φ〉 = 0.
In this case
Z(R)αβ =
P0(R)λ(R)
† F †α λ(R)Fβ
n(R)β
, (25)
∆(R)αβ =
P0(R)λ(R)
† F †α λ(R)F
1− n(R)β
. (26)
Moreover, if one constructs the states |Γ;R〉 so that
P0(R) is diagonal
(P0(R))ΓΓ′ = δΓΓ′ P0(R; Γ),
(P (R))Γ1Γ2 = P0(R; Γ1)
λ(R)Γ3Γ2 , (27)
B. The variational problem
We are now in position to settle up the variational
problem. We consider a generic tight-binding Hamilto-
H = −
R,αcR′,β
E(R)ΓΓ′ |Γ;R〉〈Γ′;R|, (28)
where α and β stem for spin, orbital and lattice site in
the chosen unit cell, and the hermitean matrix E(R) with
elements E(R)ΓΓ′ may be also unit-cell dependent. The
average value of this Hamiltonian on the Gutzwiller wave-
function (1) in the limit of infinite coordination lattices
or, in finite coordination ones, within the Gutzwiller ap-
proximation, is
Evar = −
Z(R)αγ Z(R′)
〈φ| c†
Z(R)αγ ∆(R′)
〈φ| c†
∆(R)αγ Z(R′)
〈φ| c
R,γ cR′,δ |φ〉
∆(R)αγ ∆(R′)
〈φ| c
P0(R)λ(R)
† E(R)λ(R)
≡ Ehop + Eint. (29)
The last term depends only on the local properties of the
uncorrelated wavefunction |φ〉, specifically on the occupa-
tion probabilities P0(R). Therefore, for any given choice
of P0(R), the optimal |φ〉 that minimizes the variational
energy is the ground state of the Hamiltonian
Hvar = −
Z(R)αγ Z(R′)
Z(R)αγ ∆(R′)
∆(R)αγ Z(R′)
R,γ cR′,δ
∆(R)αγ ∆(R′)
µ(R)αβ c
R,αcR,β
ν(R)αβ c
+H.c.
, (30)
where the parameters µ(R)αβ and ν(R)αβ are Lagrange
multipliers to be determined by imposing that the ground
state has indeed the chosen P0(R). The last task is to
find the values of the variational parameters λ(R)ΓΓ′ as
well as of P0(R) for which the variational energy (29) is
minimum. We note that the variational Hamiltonian (30)
that has to be solved may include also inter-site pairing
terms, which are absent in the original Hamiltonian (28).
Analogously to other more conventional variational ap-
proaches, like Hartree-Fock theory, it is common to in-
terpret the single-particle spectrum of the variational
Hamiltonian (30) as an approximation of the true co-
herent spectrum of quasi-particles.22
III. THE MODEL
Let us now apply the variational wavefunction to spe-
cific models that are inspired by the valence-bond crys-
tal example we mentioned in the introduction, and where
the off-diagonal elements of the operator PR as well as
its non-hermitean character do play an important role.
Since the operator PR is built out of purely local prop-
erties, namely the available on-site electronic configura-
tions plus a variational guess for the uncorrelated on-site
single-particle density-matrix, a lot of preliminary results
can be extracted without even specifying how lattice-sites
are coupled together. Therefore we start our analysis
from defining some local properties and later we will con-
sider a specific lattice model.
A. The isolated dimer
The basic unit of the model we are going to investigate
consists of a dimer with Hamiltonian
Hdimer = −t⊥
1σc2σ +H.c.
(ni − 1)2
≡ H⊥ +HU , (31)
where 1 and 2 refer to the two sites of the dimer and ni,
i = 1, 2, is the on-site occupation number.
It is more convenient to work in the basis of the even
(bonding) and odd (anti-bonding) combinations, defined
through
ceσ =
(c1σ + c2σ) , coσ =
(c1σ − c2σ) .
and use this basis to built the available electronic con-
figurations, which we will denote as |n,Γ〉, with n that
refers to the number of electrons. The empty and the
fourfold occupied dimer states are denoted as |0〉 and
|4〉, respectively, while the singly-occupied states as
|1, e(o), σ〉 = c†
e(o)σ
and the states with 3 electrons as
|3, e(o), σ〉 = c†
e(o)σ
o(e)↑ c
o(e)↓ |0〉.
There are six doubly-occupied configurations. Two are
spin-singlets with two electrons in the even or in the odd
orbital, |2, e〉 and |2, o〉, respectively. When each orbital
is singly occupied, the two electrons form either a spin
triplet, |2, 1, Sz〉 with Sz = −1, 0, 1, or a spin singlet,
|2, 0〉. Since we are not going to consider variational so-
lutions that break spin-SU(2) symmetry, it is convenient
to define the projector operators
|1, e(o)〉〈1, e(o)| =
|1, e(o), σ〉〈1, e(o), σ|,
|3, e(o)〉〈3, e(o)| =
|3, e(o), σ〉〈3, e(o), σ|,
|2, 1〉〈2, 1| =
Sz=−1
|2, 1, Sz〉〈2, 1, Sz|.
The isolated-dimer ground state in the subspace with
two electrons is
|Ψ〉 = cos θ√
1↑c2↓ + c
2↑c1↓
sin θ√
1↑c1↓ + c
2↑c2↓
(cos θ + sin θ) |2, e〉
(cos θ − sin θ) |2, o〉,
where tan 2θ = 4t⊥/U and has energy
+ 4t2⊥. (32)
|Ψ〉 can be always rewritten in the form of a Gutzwiller
wavefunction. First of all, we needs to choose an uncor-
related wavefunction |φ〉. A natural choice might be the
ground state at U = 0, namely |2, e〉. Indeed |Ψ〉 can be
written as
|Ψ〉 = P |2, e〉,
where
P = |Ψ〉〈2, e| = 1√
(cos θ + sin θ) |2, e〉〈2, e|
(cos θ − sin θ) |2, o〉〈2, e|. (33)
and obviously satisfies both (4) and (5).
Another possibility, that we are also going to con-
sider in what follows, is to use an uncorrelated wavefunc-
tion that corresponds to a dimer in which the two sites
are only coupled by an intersite singlet-Cooper pairing,
namely with 〈c†1↑c
2↓〉 = 〈c
1↓〉 6= 0. In this case
|φ〉 = 1
|0〉+ |2, e〉 − |2, o〉 − |4〉
and, once again, the true ground state can be written as
|Ψ〉 = |Ψ〉〈φ| |φ〉 ≡ P |φ〉.
Already at this stage one can appreciate how important
is the role of the off-diagonal elements in P , especially
for large U/t⊥.
B. The non-isolated dimer: variational density
matrix
When the dimer is coupled to the rest of the system, in
order to built the operator P we need to specify an un-
correlated local single-particle density matrix based on a
variational guess of the uncorrelated wavefunction |φ〉. A
simple guess would be a magnetic wavefunction in which
the two sites of each dimer have opposite magnetization.
This choice is also the only one admitted by an Hartree-
Fock decomposition of the interaction term HU . How-
ever, a magnetic wavefunction is not the most suitable
choice to reproduce the limit of isolated dimers, which is
a collection of singlets.
Alternatively, one can consider a paramagnetic |φ〉 that
has built in the tendency of each dimer to lock into a
spin-singlet. This can be accomplished in two ways that
do not exclude each other. The first is to assume an
uncorrelated wavefunction with a huge splitting between
even and odd orbitals, namely with
〈φ| c†eσceσ |φ〉 ≫ no =
〈φ| c†oσcoσ |φ〉.
This implies that, among the doubly-occupied configura-
tions of each dimer, mainly the spin-singlet |2, e〉 survives
in the uncorrelated wavefunction. The latter can then be
turned into the isolated dimer configuration by an ap-
propriate Gutzwiller operator P , as shown before. The
other possibility is to include Cooper pairing correlations
in the singlet channel
∆SC = 〈φ| c†1↑c
2↓ + c
1↓ |φ〉.
In this case, the isolated dimer can be recovered by as-
suming a very strong pairing ∆SC ≃ 1 and suppressing,
through P , configurations with none or two singlet-pairs.
Note that both ne − no and ∆SC do not appear by a
mean-field decoupling of HU , so that a variational wave-
function with such correlations built in can not be sta-
bilized within Hartree-Fock theory. Here the role of P
becomes crucial.
Therefore, let us assume for |φ〉 a BCS-wavefunction
defined such that
〈φ| c†1σc1σ |φ〉 = 〈φ| c
2σc2σ |φ〉 =
, (34)
〈φ| c†1σc2σ |φ〉 = 〈φ| c
2σc1σ |φ〉 =
, (35)
〈φ| c†1↑c
2↓ |φ〉 = 〈φ| c
1↓ |φ〉 =
, (36)
with real ∆SC . In the even/odd basis this translates into
〈φ| c†eσceσ |φ〉 =
, (37)
〈φ| c†oσcoσ |φ〉 =
, (38)
〈φ| c†
e↓ |φ〉 = −〈φ| c
o↓ |φ〉 =
, (39)
where ne + no = n. As previously mentioned, the cal-
culations simplify considerably in the natural basis, that
is derived in the Appendix for this particular choice of
density matrices.
As a particular application, we assume hereafter that
the model is half-filled, namely ne +n0 = 2. The density
matrix of the operators in the natural basis, d
e(o)σ
e(o)σ
is, by Eqs. (A.3) and (A.4),
〈φ| d†
e(o)σ
e(o)σ
|φ〉 = 1
where
δ2 +∆2
. (40)
The two angles θe and θo, that are defined by Eq. (A.2)
and identify the unitary transformation from the original
to the natural basis, are given by θe = θ and θo = θ−π/2,
where
tan 2θ =
. (41)
We note that, for q → 1/2, the uncorrelated wavefunc-
tion describes an insulator where charge fluctuations are
completely suppressed since each natural orbital is fully
occupied. It is obvious that, if our choice of the varia-
tional wavefunction is correct, then the optimal uncor-
related wavefunction must asymptotically acquire
q = 1/2 for U → ∞.
The expression in the natural basis of the hopping,
Eq. (A.6), and interaction, Eq. (A.7), operators can be
derived through (A.5) and have a relatively simple ex-
pression at half-filling:
H⊥ = −t⊥
2 cos 2θ
|4̃〉〈4̃| − |0̃〉〈0̃|
+ cos 2θ
|3̃〉〈3̃| − |1̃〉〈1̃|
− sin 2θ
|1̃〉〈3̃|+H.c.
2 sin 2θ
|0̃〉〈2̃,+|+ |4̃〉〈2̃,+|+H.c.
|0̃〉〈0̃|+ |4̃〉〈4̃| − |0̃〉〈4̃| − |4̃〉〈0̃|
|2̃,+〉〈2̃,+|+ |2̃,−〉〈2̃,−|+ |2̃, 0〉〈2̃, 0|
|1̃〉〈1̃|+ |3̃〉〈3̃|
, (43)
where we have defined
|1̃(3̃)〉〈1̃(3̃)| = |1̃(3̃), e〉〈1̃(3̃), e|+ |1̃(3̃), o〉〈1̃(3̃), o|,
|1̃〉〈3̃| = |1̃, e〉〈3̃, e|+ |1̃, o〉〈3̃, o|,
|2̃,±〉 = 1√
|2̃, e〉 ± |2̃, o〉
and denoted the local configurations in the natural basis
as |ñ,Γ〉 to distinguish them from the analogous ones in
the original representation.
C. The Gutzwiller operator P
The most general Gutzwiller operator P should include
at least all the projectors |ñ,Γ〉〈ñ,Γ| as well as all the
off-diagonal operators |ñ,Γ〉〈ñ′,Γ′| that appear in the lo-
cal Hamiltonian, Eqs. (42) and (43). As we mentioned
before, our expectation is that the uncorrelated wave-
function which better connects to the large-U Mott in-
sulator should be identified by q → 1/2, in which locally
only the configurations |3̃〉 and |4̃〉 are occupied with non-
negligible probability. This suggests that P must include
at least those off-diagonal operators that would turn |4̃〉
into the isolated dimer ground state, namely |0̃〉〈4̃| and
|2̃,+〉〈4̃|. The latter forces to include also |1̃〉〈3̃|, as we
shall see.
Therefore we assume for P the following variational
ansatz:
λnΓ |ñ,Γ〉〈ñ,Γ|+ λ13 |1̃〉〈3̃|
+λ04 |0̃〉〈4̃|+ λ2+ 4 |2̃,+〉〈4̃|, (44)
with real λ’s. We define
P (n,Γ) = λ2nΓ P0(ñ,Γ), (45)
for all n 6= 3, 4, while, for n = 3, 4,
P (3) =
λ23 + λ
P0(3̃), (46)
P (4) =
λ24 + λ
04 + λ
P0(4̃). (47)
Then the conditions Eqs. (4) and (5) read
P (n,Γ) = 1, (48)
nP (n,Γ) = 2 + 4q, (49)
λ13 λ1
P0(3̃)P0(1̃) =
2λ2+ 4 λ2+
P0(4̃)P0(2̃,+). (50)
Here P0(ñ,Γ) are the occupation probabilities in the nat-
ural basis of the uncorrelated wavefunction. Specifically
P0(ñ,Γ) = gen,Γ
)en (
)4−en
where gen,Γ is the degeneracy of the configuration.
Eq. (50) guarantees that the anomalous averages
〈φ| P† P d†
e(o)↑d
e(o)↓ |φ〉
vanish in the natural basis, and explains why we have
included |1̃〉〈3̃| in (44). It is convenient to rewrite
P (4)
P0(4̃)
ui for i = 4, 04, 2 + 4,
P (3)
P0(3̃)
ui for i = 3, 13,
where u23 + u
13 = 1, which can be satisfied by choosing
u3 = cosψ and u13 = sinψ, and u
2+4 = 1. The
latter parameters can be expressed by means of another
unit vector v = (v1, v2, v3), through
(u4 + u04) ,
cos 2θ√
(u4 − u04)− sin 2θ u2+4,
sin 2θ√
(u4 − u04) + cos 2θ u2+4,
In terms of all the variational parameters, the P (n,Γ)’s,
θ, q, ψ and v, the average values per dimer of the in-
teraction, HU , and intra-dimer hopping, H⊥, are readily
found to be
EU = 〈φ| P†HU P |φ〉 =
P (3) + P (1)
P (0) + U
P (2,+) + P (2,−) + P (2, 0)
v22 + v
P (4), (52)
E⊥ = 〈φ| P†H⊥ P |φ〉 = −2 t⊥ δ∗, (53)
where the actual correlated values of the hybridization
and of the anomalous average are
2 δ∗ = 〈φ| P (ne − no) P |φ〉
= 4 v1 v2 P (4) + cos (2θ + 2ψ) P (3)
−2 cos 2θ P (0)− cos 2θ P (1), (54)
2∆∗ = 〈φ| P
e↓ − c
o↓ +H.c.
P |φ〉
= 4 v1 v3 P (4) + sin (2θ + 2ψ) P (3)
−2 sin2θ P (0)− sin 2θ P (1), (55)
We note that δ∗ and ∆∗ are mutually exclusive, namely
the choice of parameters that maximizes one of the two,
makes the other vanishing.
Upon the action of P , the single fermion operators
in the Nambu spinor representation transform effectively
e(o)↑
e(o)↓
Z +∆ e−i β τ2
e(o)↑
e(o)↓
, (57)
where τi, i = 1, 2, 3, are the Pauli matrices that act on
the Nambu spinor components,
tanβ =
. (58)
and, finally,
〈φ| P† d†
e(o)σ
e(o)σ
1− 4q2
P (0)P (1) +
P (1)P (2,+)
P (1)P (2,−) + 1
P (1)P (2, 0)
P (1)P (2, 1) + cosψ
P (3)P (2, 1)
+ cosψ
P (3)P (2, 0) + cosψ
P (3)P (2,−)
+ cosψ
P (3)P (2,+)
v1 cosψ + v2 cos (2θ + ψ)
+ v3 sin (2θ + ψ)
P (4)P (3)
, (59)
〈φ| P† d†
e(o)↑ P d
e(o)↓ |φ〉
1− 4q2
P (3)P (2,−)
P (3)P (2,+)
P (3)P (2, 1)
P (3)P (2, 0)
v1 sinψ − v2 sin (2θ + ψ)
+v3 cos (2θ + ψ)
P (4)P (3)
, (60)
with real
Z and
∆. Therefore, if the dimers are cou-
pled one to another by the single particle hopping term
R 6=R′
i,j=e,o
τ3 ΨR′,j , (61)
where
R,i↑, cR,i↓),
and Ψ its hermitean conjugate, the uncorrelated wave
function |φ〉 minimizes the effective hopping
Tvar = (Z +∆)
R 6=R′
i,j=e,o
−2 i β τ2 Ψ
under the condition that the local density matrix satisfies
Eqs.(37)-(39). One can readily show that this amounts to
find the ground state |φ〉 of the variational Hamiltonian
Hvar = (Z +∆)
R 6=R′
i,j=e,o
τ3 ΨR′,j
R,e τ3 ΨR,e −Ψ
R,o τ3 ΨR,o
R,e τ1 ΨR,e −Ψ
R,o τ1 ΨR,o
, (63)
with µ3 and µ1 such that
〈φ| Hvar |φ〉 + 4q µ3 cos(2θ + 2β) + 4q µ1 sin(2θ + 2β),
is maximum.
Before we consider specific lattice models, it is worth
re-deriving within this variational scheme the isolated-
dimer ground-state energy (32) at half-filling. For that
purpose, we take all P (n,Γ) zero but P (4) = 1. The
variational energy is simply
Evar = EU + E⊥ = U
v22 + v
− 4 t⊥ v1 v2.
The minimum under the constraint v · v = 1 is obtained
for v3 = 0 and exactly reproduces (32). We note that
the minimum energy is independent on θ, namely there
exists a continuous family of variational solutions with
equal energy parametrized by θ. However, in spite of the
fact that the uncorrelated wavefunction may describe a
superconductor, the actual value of the anomalous aver-
age ∆∗ = 0.
IV. A LATTICE MODEL OF DIMERS
As a particular application, let us consider the follow-
ing lattice model
H = −
tRR′ c
R,iσcR′,iσ +H.c.
(nR,i − 1)2
<RR′>
R,1σcR,2σ +H.c.
(ǫk − t⊥) c†k,eσck,eσ + (ǫk + t⊥) c
(nR,i − 1)2 , (64)
where nR,i =
and ǫk is the band dis-
persion induced by tRR′ , with half-bandwidth D. The
Hamiltonian (64) represents two Hubbard models cou-
pled by a single-particle hopping t⊥, each model being
defined on a lattice with coordination number z. As we
mentioned, the variational results that we have so far
derived are rigorous strictly speaking only if z → ∞,
although, in the spirit of the Gutzwiller approximation,
they can be used for generic z as well.
If U ≫ D, t⊥, (64) describes at half-filling a Mott insu-
lator which may be magnetic at t⊥ ≪ D, but is certainly
non-magnetic at t⊥ ≫ D, where the ground state reduces
essentially to a collection of singlets. For instance, in the
case of a Bethe lattice with nearest neighbor hopping,
the transition is at t⊥ = D/
8, value that is going to
decrease if frustration is included. If U is small and the
Fermi surface is not nested, then the model is metallic
for t⊥ ≤ D and is a band insulator otherwise. In fact,
in the absence of nesting there is generically a finite win-
dow of t⊥ values in which, upon increasing U , the model
undergoes a transition from a paramagnetic metal into a
non-magnetic Mott insulator, and this is just the case we
are going to consider in what follows. The same model
have been recently studied by Fuhrmann, Heilmann and
Monien using DMFT23 and by Kancharla and Okamoto24
using DMFT and cluster DMFT, respectively, that gives
us the opportunity to directly check the accuracy of our
wavefunction.
The variational Hamiltonian (63) of the model (64) has
a very simple expression,
Hvar =
ǫk∗ − µ3
τ3 − µ1 τ1
ǫk∗ + µ3
τ3 + µ1 τ1
, (65)
where
ǫk∗ = (Z +∆) ǫk.
The variational single-particle spectrum has the conven-
tional BCS form with eigenvalues
Eek =
(ǫk∗ − µ3)2 + µ21, Eok =
(ǫk∗ + µ3)
+ µ21,
hence, for any µ1 6= 0, has a gap equal to 2µ1. On the
contrary, when µ1 = 0, the spectrum is gapless for |µ3| ≤
D, otherwise is gaped. The Lagrange multipliers µ1 and
µ2 are obtaining by maximizing
Ehop = −
Eek + Eok
+ 4q µ3 cos(2θ + 2β) + 4q µ1 sin(2θ + 2β). (66)
In terms of (66), (52) and (53) the variational energy is
Evar = Ehop + EU + E⊥, (67)
and depends on eight independent variational parame-
ters.
We have solved numerically the variational problem at
fixed t⊥/D = 0.5 as function of U/D. To simplify cal-
culations, we have assumed for the band dispersion ǫk
either a flat or a semi-circular density of states, although
both would give rise to nesting that could stabilize mag-
netic phases, which we do not take into account. How-
ever, from the point of view of the paramagnetic-metal
to paramagnetic-insulator transition, this choice is not
influential.
We find that the variational solution displays a first
order phase transition at Uc ≃ 2.05 D for a flat density
of states, as shown by the behavior of the variational en-
ergy in Fig. 2. This result agrees almost quantitatively
with the DMFT calculation23 obtained with a semicir-
cular density of states, that also predicts a first order
transition with a coexistence region between U ≃ 1.5 D
and 1.8 D at the same value of t⊥ = 0.5 D. We note that
the energy is everywhere finite and vanishes like 1/U for
large U , see the inset of Fig. 2. The asymptotic behavior
UEvar/D
2 ∼ −7/6 is compatible with second order per-
turbation theory in t and t⊥ using as zeroth-order state
a collection of dimers, as explained below. In Fig. 3 we
show the behavior across the transition of the three con-
tribution to the energy, namely EU , E⊥ and Ehop. We
FIG. 2: Variational energy Evar in units of D and for a flat
non-interacting density of states, as function of U/D for t⊥ =
D/2. At Uc ≃ 2.05 D a first order transition occurs. The
inset shows the asymptotic value of UEvar/D
FIG. 3: The different contributions to the variational energy,
EU , E⊥ and Ehop.
find that the transition is accompanied by an energy loss
in Ehop, but a gain in both E⊥ and EU .
In order to characterize physically the two phases, in
Fig. 4 we plot the values of µ1 and µ3 across the tran-
sition. Since µ1 = 0, within our numerical accuracy,
and |µ3| < D, the phase at U < Uc is gapless hence
metallic, see the behavior of the density of states (DOS)
drawn in Fig. 5. On the contrary, on the U > Uc side
FIG. 4: The behavior of the parameters µ1 and µ3 in units of
D as function of U/D, see Eq. (65).
of the transition, µ1 6= 0, that implies a finite gap in the
single-particle variational spectrum, see Fig. 5. In the
gaped phase at U > Uc the spectrum looks like the one
of a Peierls insulator with a very large hybridization gap,
not consistent with the bare value of t⊥. In reality this
gap is, more properly, the Mott-Hubbard gap. Indeed
the DOS has weight both below and above the chemical
potential, suggestive of asymmetric Mott-Hubbard side-
bands. Moreover, as we are going to discuss below, the
actual difference between the occupations of the bond-
ing and anti-bonding bands, which we denoted as 2δ∗ in
Eq. (54), decreases with U , unlike the single-particle gap,
see Fig. 7. This behavior is reminiscent of what has been
found by Biermann et al.25 as an attempt to understand
the physics of VO2.
The other quantities that identify the variational spec-
trum are Z and ∆, shown in Fig. 6. We see that Z is
decreasing with U but reaches a finite value Z = 1/4 for
U → ∞. On the contrary, ∆ = 0 for U < Uc, while
∆ 6= 0 for U > Uc and increases monotonically to reach
asymptotically the same value 1/4 for large U .
Further insights can be gained by the average values of
the intra-dimer hopping and pairing, Eqs. (54) and (55),
drawn in Fig. 7. As we mentioned, the intra-dimer hy-
bridization 2δ∗ is monotonically decreasing with U , apart
from the jump at the first order transition. More inter-
estingly, around the transition the variational solution
is characterized by a sizeable BCS order parameter ∆∗,
which does not follow the behavior of the BCS coupling
µ1 present in the variational Hamiltonian. Indeed, while
µ1 is zero within our numerical accuracy for U < Uc,
yet a non negligible ∆∗ develops just before the transi-
tion, see Fig. 7. Moreover, although µ1 starts already
U/D = 2.0
U/D = 0
U/D = 2.05
U/D = 3.0
−2  0.5
FIG. 5: The variational single-particle spectrum for the even,
i.e. bonding, band, solid lines, and odd, i.e. anti-bonding,
one, dashed lines, across the transition for a non-interacting
semi-circular density of states.
large for U > Uc and increases monotonically with U ,
see Fig. 4, the actual order parameter ∆∗ is apprecia-
ble only near the transition and fastly decreases with
U to very tiny values. Also interesting is that, besides
the intra-dimer superconducting order parameter, also an
inter-dimer one arises, that can be for instance defined
through
∆̃∗ =
tRR′ 〈φ| P† Ψ†R,i τ1 ΨR′,i P |φ〉
FIG. 6: The behavior of the parameters Z and ∆.
FIG. 7: The average values of the intra-dimer hopping, δ∗,
and pairing, ∆∗.
4ǫk=0
where V is the number of sites. We find that ∆̃∗ has
actually the opposite sign of ∆∗, and both closely follow
each other, rapidly decreasing with U , see Fig. 8. Even
though ∆∗ and ∆̃∗ are everywhere finite for U > Uc, sug-
gestive of a superconducting phase that survives up to
very large U , we believe that superconductivity, it it oc-
curs at all, may appear only very close to the transition,
where the value of the order parameter is larger. Indeed,
FIG. 8: The intra-dimer, circles, and, with reversed sign, the
inter-dimer, triangles, superconducting order parameters.
the optimal solution with finite ∆∗ and ∆̃∗ and another
solution in which both are forced to be zero are practi-
cally degenerate within our numerical accuracy for large
U . Moreover, since the variational values of the order
parameter are extremely small but close to Uc, the in-
clusion of quantum fluctuations, for instance in the form
of a Jastrow factor as in Eq. (2), would likely suppress
superconductivity leading to a bona fide insulating wave-
function.
Unfortunately, we can not compare this result with the
DMFT analyses of Refs. 23,24, where superconductivity
has not been looked for.26
A. Large U limit
In order to appreciate qualities and also single out
defects of the variational wavefunction, it is worth dis-
cussing the large-U solution. To leading order in 1/U ,
one can assume all P (n,Γ) = 0 but P (3) and P (4). The
two constraints Eqs. (48) and (49) can be solved by defin-
P (3) = 2− 4q ≡ 4d, P (4) = 1− 4d,
with d ≪ 1, namely q → 1/2, while Eq. (50) is already
satisfied since P (1) = P (2,+) = 0. Moreover, one readily
recognizes that the variational solution asymptotically
tends to acquire ψ ≃ β → π/4, θ → 0 and v1 ≫ v2, v3,
which is indeed what we find by numerical minimization.
It then follows that
∆ → 1
1− 4q2
4d(1− 4d) → 1
This implies that µ3 → 0, hence that µ1 is determined
by maximizing
Ehop = −
k∗ + µ
1 + 4q µ1,
which leads to µ1 =
ǫ2/4d and
Ehop = −4
ǫ2 d, (68)
where
At leading order, the variational energy per dimer is
therefore
Evar = −4
ǫ2 d− 4 t⊥ v1 v2 + U v22 + 2U d,
that is minimized by d = ǫ2/U
2, v1 ≃ 1 and v2 ≃ 2t⊥/U ,
and takes the value
Evar = −
= − 1
tRR′ tR′R
. (69)
In the case of a flat density of states
δ (ǫ− ǫk) =
θ (D − |ǫ|) ,
and with t⊥ = 0.5 D we recover the numerical result
Evar ≃ −7/6 D2/U , see Fig. 2.
We note that, in spite of the hybridization δ∗ ∼ t⊥/U
being small, the single-particle gap of the variational
spectrum 2µ1 ≃ U is large, as one should expect in a
Mott insulator.
Coming back to the large-U value of the variational
energy (69), one can readily realize that it coincides with
the second order correction in tRR′ to the energy of the
state
|Ψ〉 =
R,1↑c
R,2↓ + c
R,2↑c
which is just a collection of singlets. In other words, in
spite of being non-magnetic, our variational wavefunc-
tion is able to reproduce the correct super-exchange be-
tween dimers. This is a remarkable property that ac-
tually derives from the square-root dependence upon d
of Ehop, see Eq. (68). If we considered a more conven-
tional Gutzwiller operator P commuting with the single-
particle density matrix, that amounts to further impose
δ2∗ +∆
∗, we would find Ehop ∝ d, implying a
transition into an unrealistic insulator with d = 0 above
a critical U . The obvious defect of the wavefunction is
that, since it emphasizes strongly the role of individual
dimers, the hopping among dimers, although finite for
any U , is under-estimated with respect to the intra-dimer
one. Therefore we do not expect the wavefunction to be
particularly accurate for small t⊥/D.
anti−bonding
bonding
FIG. 9: The non-interacting density of states of the lattice
of dimers. The bonding and anti-bonding state of each dimer
give rise to two bands that overlap, leading to a metallic phase
in the absence of interaction.
V. CONCLUSIONS
In this work we have proposed an extension of the
Gutzwiller variational approach to account for correlated
models which display metal-insulator transitions into
Mott insulators that escape any simple single-particle
descriptions, like the Hartree-Fock approximation. The
wavefunction has still the same form as the conventional
Gutzwiller wavefunction,
|ΨG〉 = P |φ〉 =
PR |φ〉,
with R identifying unit cells that may also be non-
primitive ones, with the novel feature that the opera-
tor PR is non-hermitean and does not commute with the
local single-particle density matrix. In essence, this prop-
erty realizes a variational implementation of a Schrieffer-
Wolff transformation27, although only restricted to the
lattice sites within each unit cell. We have shown that, by
slightly reducing the variational freedom, this wavefunc-
tion, like the conventional Gutzwiller wavefunction,14 al-
lows for an extension of the Gutzwiller approximation
to evaluate average values, approximation that becomes
exact in the limit of infinite coordination lattices.
As an application, we have considered the Mott tran-
sition into a Peierls, or valence-bond, insulator, namely
an insulator that is adiabatically connected to a collec-
tion of independent dimers. Such an insulator can not
be described by Hartree-Fock simply because the singlet
configuration of each dimer is not a Slater determinant.
Specifically, we have considered the hypothetical situa-
tion shown in Fig. 9, where the splitting between the
bonding and antibonding orbitals of each dimer is as-
sumed not to be sufficient to lead, in the absence of in-
teraction, to a band insulator. When interaction is taken
into account, in the form of an on-site repulsion U , one
expects, above a critical U , a transition from the metal
into a Mott insulator. If magnetism is prevented, for
instance by a sufficiently large splitting between bond-
ing and anti-bonding orbitals and/or by frustration, the
Mott insulator is non-magnetic. We have shown that
our wavefunction overcomes the difficulties of Hartree-
Fock theory and allows to study, albeit variationally, this
transition. In particular we find that:
(i) at the variational level the Mott transition is first
order;
(ii) the variational spectrum inside the Mott insulator
looks similar to that of a Peierls insulator with a
large hybridization gap, namely a large splitting be-
tween bonding and anti-bonding bands. In reality
the gap is the Mott-Hubbard gap and the actual
difference between the occupations of the bonding
and anti-bonding bands, is small;
(iii) inter-site singlet-superconductivity appears around
the transition.
While (i) and (ii) are presumably true, as they have
been also found by more rigorous calculations23,25,28, the
emergence of superconductivity might be an artifact of
the variational wavefunction.26 Nevertheless, the possi-
ble occurrence of superconductivity is quite suggestive.
It is known for instance that two-leg Hubbard ladders
with nearest neighbor hopping display dominant super-
conducting fluctuations with the same symmetry that we
find variationally29, although at half-filling they always
describe non-magnetic spin-gaped insulators30,31 because
of nesting. Moreover, the uncorrelated wavefunction |φ〉
is quite similar to the wavefunctions used in Refs. 32,33
to simulate t-J ladders. It would be surprising and inter-
esting if this tendency towards superconductivity turned
into a true symmetry breaking instability in higher di-
mensionality, as suggested by our analysis, which we
think it is worth deserving further investigations.
Note added: During the completion of this work, we be-
came aware of a recent extension of slave-boson technique
whose saddle-point solution closely resembles our varia-
tional approach.34 Indeed the two conditions we impose
on the Gutzwiller operator, Eqs. (4) and (5), are in one-
to-one correspondence with the constraints identified in
Ref. 34 within the slave-boson formalism.
Acknowledgments
We are grateful to C. Castellani and E. Tosatti for their
helpful comments and suggestions. We also thanks A.
Georges for useful discussions in connection with Ref. 34.
APPENDIX: THE NATURAL BASIS
Let us assume that, in the Nambu-spinor representa-
the uncorrelated wavefunction has the following density
matrices
Ĉe =
ne/2 ∆SC/2
∆SC/2 1− ne/2
, Ĉo =
no/2 −∆SC/2
−∆SC/2 1− no/2
(A.1)
The natural orbitals are obtained by the unitary trans-
formation
e(o)↑ = cos θe(o) ce(o)↑ + sin θe(o) c
e(o)↓
e(o)↓ = cos θe(o) ce(o)↓ − sin θe(o) c
e(o)↑,
where
tan 2θe =
ne − 1
, tan 2θo =
no − 1
, (A.2)
and posses a diagonal density matrix with the non-
vanishing elements given by
〈φ| d†
e(o)σ
e(o)σ
|φ〉 = 1
+ qe(o), (A.3)
where
qe(o) =
(ne(o) − 1)2 +∆2SC . (A.4)
In the natural basis we introduce states that have the
same formal expression as in the original basis but are
built with d-operators, and denote them as |ñ,Γ〉. The
transformation rules from these states to the original ones
|0〉 = cos θe cos θo|0̃〉+ sin θe sin θo |4̃〉
+cos θe sin θo |2̃, o〉+ sin θe cos θo |2̃, e〉,
|1, e(o), σ〉 = cos θo(e) |1̃, e(o), σ〉 + sin θo(e) |3̃, e(o), σ〉,
|2, e(o)〉 = cos θe cos θo |2̃, e(o)〉+ cos θe(o) sin θo(e) |4̃〉
− sin θe(o) cos θo(e) |0̃〉
− sin θe(o) sin θo(e) |2̃, o(e)〉,
|2, 1, Sz〉 = |2̃, 1, Sz〉, (A.5)
|2, 0〉 = |2̃, 0〉,
|3, e(o), σ〉 = cos θo(e) |3̃, e(o)〉 − sin θo(e) |1̃, e(o)σ〉,
|4〉 = cos θe cos θo|4̃〉+ sin θe sin θo |0̃〉
− cos θe sin θo |2̃, e〉 − sin θe cos θo |2̃, o〉.
The inverse transformation is obtained by letting θe(o) →
−θe(o).
The hopping operator in the original representation
can be written as
H⊥ = −2t⊥
1σc2σ +H.c.
= −2t⊥
c†eσceσ − c†eσceσ
= −2t⊥
|3, o, σ〉〈3, o, σ| − |3, e, σ〉〈3, e, σ|
+|1, e, σ〉〈1, e, σ| − |1, o, σ〉〈1, o, σ|
+2 |2, e〉〈2, e| − 2 |2, o〉〈2, o|
, (A.6)
while the interaction operator as
(ni − 1)2 =
|0〉〈0|+ |4〉〈4|
n=1,3
|n, e, σ〉〈n, e, σ|+ |n, o, σ〉〈n, o, σ|
|2, e〉+ |2, o〉
〈2, e|+ 〈2, o|
+2 |2, 0〉〈2, 0|
. (A.7)
Their expression in the natural basis can be obtained by
the transformation rules (A.5).
1 M. C. Gutzwiller, Phys. Rev. 134, A923 (1964).
2 M. C. Gutzwiller, Phys. Rev. 137, A1726 (1965).
3 W. F. Brinkman and T. M. Rice, Phys. Rev. B 2, 4302
(1970).
4 M. Capello, F. Becca, M. Fabrizio, S. Sorella, and
E. Tosatti, Phys. Rev. Lett. 94, 026406 (2005).
5 M. Capello, F. Becca, S. Yunoki, M. Fabrizio, and
S. Sorella, Phys. Rev. B 72, 085121 (2005).
6 M. Capello, F. Becca, M. Fabrizio, and S. Sorella, Super-
fluid to Mott-insulator transition in Bose-Hubbard models
(2006),
7 S. Sorella, Phys. Rev. B 71, 241103 (2005).
8 A. Georges, G. Kotliar, W. Krauth, and M. J. Rozenberg,
Rev. Mod. Phys. 68, 13 (1996).
9 A. I. Lichtenstein and M. I. Katsnelson, Phys. Rev. B 62,
R9283 (2000).
10 D. Senechal, D. Perez, and M. Pioro-Ladriere, Phys. Rev.
Lett. 84, 522 (2000).
11 M. Potthoff, M. Aichhorn, and C. Dahnken, Phys. Rev.
Lett. 91, 206402 (2003).
12 T. Maier, M. Jarrell, T. Pruschke, and M. H. Hettler, Rev.
Mod. Phys. 77, 1027 (2005).
13 G. Kotliar, S. Y. Savrasov, G. Pálsson, and G. Biroli, Phys.
Rev. Lett. 87, 186401 (2001).
14 J. Bünemann, W. Weber, and F. Gebhard, Phys. Rev. B
57, 6896 (1998).
15 C. Attaccalite and M. Fabrizio, Phys. Rev. B 68, 155117
(2003).
16 Q.-H. Wang, Z. D. Wang, Y. Chen, and F. C. Zhang, Phys-
ical Review B (Condensed Matter and Materials Physics)
73, 092507 (pages 4) (2006).
17 M. Ferrero, Ph.D. thesis, SISSA-Trieste (2006).
18 W. Metzner and D. Vollhardt, Phys. Rev. Lett. 59, 121
(1987).
19 W. Metzner and D. Vollhardt, Phys. Rev. B 37, 7382
(1988).
20 F. Gebhard, Phys. Rev. B 41, 9452 (1990).
21 G. Kotliar and A. E. Ruckenstein, Phys. Rev. Lett. 57,
1362 (1986).
22 J. Bünemann, F. Gebhard, and R. Thul, Phys. Rev. B 67,
075103 (2003).
23 A. Fuhrmann, D. Heilmann, and H. Monien, Physical Re-
view B 73, 245118 (2006).
24 S. S. Kancharla and S. Okamoto, Band
insulator to Mott insulator transition in
a bilayer Hubbard model (2007), URL
http://www.citebase.org/abstract?id=oai:arXiv.org:cond-mat/0703728.
25 S. Biermann, A. Poteryaev, A. I. Lichtenstein, and
A. Georges, Physical Review Letters 94, 026404 (pages 4)
(2005).
26 There are however recent DMFT calculations that seem to
support the existence of a superconducting region around
the transition. [A. Privitera, M. Schiro’, M. Capone and
C. Castellani, unpublished (2007).].
27 J. R. Schrieffer and P. A. Wolff, Phys. Rev. 149, 491
(1966).
28 G. Moeller, V. Dobrosavljević, and A. E. Ruckenstein,
Phys. Rev. B 59, 6846 (1999).
29 M. Fabrizio, Phys. Rev. B 48, 15838 (1993).
30 S. P. Strong and A. J. Millis, Phys. Rev. B 50, 9911 (1994).
31 D. G. Shelton, A. A. Nersesyan, and A. M. Tsvelik, Phys.
Rev. B 53, 8521 (1996).
32 G. Sierra, M. A. Mart́ın-Delgado, J. Dukelsky, S. R. White,
and D. J. Scalapino, Phys. Rev. B 57, 11666 (1998).
33 S. Sorella, G. B. Martins, F. Becca, C. Gazza, L. Capriotti,
A. Parola, and E. Dagotto, Phys. Rev. Lett. 88, 117002
(2002).
34 F. Lechermann, A. Georges, G. Kotliar, and O. Parcollet,
unpublished (2007).
http://www.citebase.org/abstract?id=oai:arXiv.org:cond-mat/0703728
