...

Exact results for strongly correlated fermions: Hubbard and Falicov-Kimball models

by user

on
Category: Documents
15

views

Report

Comments

Transcript

Exact results for strongly correlated fermions: Hubbard and Falicov-Kimball models
Exact results for strongly correlated fermions: Hubbard and
Falicov-Kimball models
Pedro Goldbaum
A Dissertation
Presented to the Faculty
of Princeton University
in Candidacy for the Degree
of Doctor of Philosophy
Recommended for Acceptance
by the Department of
Physics
May, 2005
c Copyright 2005 by Pedro Goldbaum.
All rights reserved.
Abstract
In this work, we study two models of strongly correlated fermions: Hubbard and FalicovKimball. We present a proof of the existence of real and ordered solutions to the nested
Bethe Ansatz equations for the one dimensional Hubbard model on a finite lattice, with
periodic boundary conditions. The existence of a continuous set of solutions extending from
any U > 0 to U = ∞ is also shown, where U is the strength of the interaction between
fermions at the same site. We use this continuity property, combined with the proof that
the norm of the wavefunction obtained with the generalized Bethe Ansatz is not zero, to
prove that the solution gives us the ground state of the finite system, as assumed by Lieb
and Wu. For the absolute ground state at half-filling, we show that the solution converges
to a distribution in the thermodynamic limit. This limit distribution satisfies the integral
equations that led to the Lieb-Wu solution of the 1D Hubbard model.
We also obtain a lower bound for the ground state energy of the Falicov-Kimball model.
This bound is given by a bulk term, plus a term proportional to the boundary. A numerical
value for the coefficient of the boundary term is determined. The explicit derivation is
important in the proof of the conjecture of segregation of the two kinds of fermions in the
Falicov-Kimball model, for sufficiently large interactions.
iii
Acknowledgements
I would first like to thank my advisor Elliott Lieb, for the suggestion of a great PhD
problem and all the guidance throughout the process. Elliott’s research contains seminal
contributions to the three topics covered in this thesis: Hubbard model, Bethe Ansatz and
Falicov-Kimball model. It has been a true privilege to have him as my advisor.
I am also indebted to the mathematical physics community in Princeton, including
Michael Aizenman, Robert Seiringer, Thomas Chen, Robert Sims, Simone Warzel, Vassilios
Papathanakos and others no longer at Princeton, for the stimulating conversations, and
also for the subsidized Beijing duck after the seminars. My graduate student friends have
contributed equally to my academic and social life at Princeton, and I would like to recognize
some of them: Kumar Raman, Srinivas Raghu, Subroto Mukerjee, Alexander Baitine and
Mike Leung.
Finally, I would like to thank my family, for all their support, in particular my wife,
Mariela Steiervalt. And my son, Felipe Goldbaum, for reminding me to enjoy the simple
things in life.
iv
Contents
Abstract
iii
Acknowledgements
iv
Contents
v
List of Figures
viii
1 Introduction
1
2 The Hubbard model
4
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4
2.2
Symmetries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6
2.3
Ferromagnetism and the Hubbard model . . . . . . . . . . . . . . . . . . . .
7
2.3.1
Lieb’s theorems and ferrimagnetism . . . . . . . . . . . . . . . . . .
8
2.3.2
Saturated ferromagnetism . . . . . . . . . . . . . . . . . . . . . . . .
10
Attractive case: monotonicity of the energy with spin . . . . . . . . . . . .
12
2.4
3 The Bethe Ansatz
16
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2
Heisenberg model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
17
3.3
1D Hubbard model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21
3.3.1
23
Checking for consistency at the boundaries . . . . . . . . . . . . . .
v
3.3.2
Periodic boundary condition
. . . . . . . . . . . . . . . . . . . . . .
27
3.3.3
M=1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
30
3.3.4
M=2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
33
3.3.5
General solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
36
4 Existence of solution to the Bethe Ansatz equations
38
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
38
4.2
Existence of a solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.2.1
Defining the map . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
42
4.2.2
The fixed point theorem . . . . . . . . . . . . . . . . . . . . . . . . .
43
4.2.3
The fixed point is a solution to the Bethe Ansatz equations . . . . .
45
4.3
Existence of a continuous curve of solutions . . . . . . . . . . . . . . . . . .
46
4.4
Non-vanishing norm of the wavefunction . . . . . . . . . . . . . . . . . . . .
50
4.5
Extensions to excited states . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
5 The thermodynamic limit
60
5.1
Defining the densities
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
60
5.2
Convergence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
62
5.3
Integral equations
64
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
6 Beyond Hubbard: Falicov-Kimball model
69
6.1
From asymmetric Hubbard to Falicov-Kimball . . . . . . . . . . . . . . . . .
69
6.2
Crystalline and segregated phases . . . . . . . . . . . . . . . . . . . . . . . .
72
6.2.1
Crystalline phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
72
6.2.2
The segregated phase . . . . . . . . . . . . . . . . . . . . . . . . . .
74
The one-dimensional model . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
6.3.1
Neutral case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
6.3.2
Charged case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
78
Perturbative results in two dimensions . . . . . . . . . . . . . . . . . . . . .
78
6.3
6.4
vi
7 Lower bound for the segregation energy in the FK model
83
7.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
83
7.2
The boundary term for n = 1/2 . . . . . . . . . . . . . . . . . . . . . . . . .
85
7.2.1
Projection of the boundary vector . . . . . . . . . . . . . . . . . . .
85
7.2.2
A bound for (hΛ − 2d)bk 2 . . . . . . . . . . . . . . . . . . . . . . .
86
7.2.3
A bound for |∇j f (k)| . . . . . . . . . . . . . . . . . . . . . . . . . .
89
7.2.4
The lower bound . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
7.3
d=2: A better result . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
91
7.4
The result for n < 1/2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
92
7.5
Finite U . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
94
7.6
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
95
A Groebner bases
97
References
101
vii
List of Figures
2.1
Bipartite lattice with |B| = 2|A|. . . . . . . . . . . . . . . . . . . . . . . . .
10
2.2
3-dimensional example with |B| = 3|A|, or C = 1/2. . . . . . . . . . . . . .
10
2.3
Kagome lattice. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
11
2.4
Tasaki’s model for near flat-band ferromagnetism.
. . . . . . . . . . . . . .
12
3.1
Consistency relations for exchanges of two and three indices. . . . . . . . .
25
3.2
An example with exchanges of four indices. . . . . . . . . . . . . . . . . . .
26
3.3
Reidemeister moves. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
26
3.4
Symmetry types and Young’s diagrams for N=10, M=4. . . . . . . . . . . .
30
4.1
Two-dimensional vector fields with zeros of index +1 and −1. . . . . . . . .
47
4.2
Sketch of the curve of zeros connecting the basis of the cylinder. . . . . . .
48
6.1
Ground state for ne = nc = 1/2.
. . . . . . . . . . . . . . . . . . . . . . . .
73
6.2
Construction of whom for Nc = 7 and N = 11. . . . . . . . . . . . . . . . . .
77
6.3
Ground states for nc = 1/5 and nc = 1/4. . . . . . . . . . . . . . . . . . . .
80
6.4
Three by three blocks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81
7.1
Diagram of second nearest neighbors. . . . . . . . . . . . . . . . . . . . . . .
87
7.2
Fermi surface εk = 4 and extended regions of integration. . . . . . . . . . .
91
7.3
Diagram for (0, k). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
94
7.4
Diagram for (π/2, π/2 − y) and the sum over set S. . . . . . . . . . . . . . .
94
viii
Chapter 1
Introduction
This thesis is a collection of results in the field of exactly solved quantum many-body
problems. We study two hamiltonians, Hubbard and Falicov-Kimball, which have been
used as simplifying models to explain non trivial phenomena in many-body physics, such
as ferromagnetism, superconductivity, metal-insulator transition and formation of crystals.
For the Hubbard model, we restrict to the one-dimensional case, which is still the only one
with an exact solution available.
Most original results in this thesis are presented in chapters 4, 5 and 7, so the reader
with a solid background in exactly solved models can in principle focus on these chapters.
The remaining chapters are brief surveys of exact results for the models and techniques
considered, with the exception of a partial proof of magnetic ordering in chapter 2.
Chapter 2 presents an introduction to the Hubbard model, stating a few of the known
exact results on the model. We start by describing the hamiltonian and its origin, and
devote special attention to its SO(4) symmetry. The magnetic properties of the ground
state are also studied, followed by a partial proof of magnetic ordering for the attractive
Hubbard model. We hope this proof will be a small contribution in validating the conjecture
that the ground state energy E(S) on a given sector with fixed spin S is strictly increasing
with S. The result is know to be valid in one dimension, by a general theorem by Lieb and
Mattis [LM].
1
2
In Chapter 3, we discuss the Bethe Ansatz, starting with its application to the solution
of the one-dimensional Heisenberg model. We briefly mention applications of the technique
in the solution of other integrable models, and devote a section to the derivation of the
Bethe Ansatz equations for the Hubbard model. At this point, the reader is exposed to all
necessary background to understand the proof of the existence of a solution to the Bethe
Ansatz equations.
Chapter 4 presents the proof that the ground state of the Hubbard model on a finite
lattice can be obtained by the Bethe Ansatz. This assumption was the starting point of
the celebrated solution of the model by Lieb and Wu [LW1]. We first show existence of
solution to the Lieb-Wu equations. We prove continuity of the solution with respect to the
interaction strength and that the resulting wavefunction is a proper (non zero) eigenstate.
To conclude the proof, we connect the solution for finite U with the solution at U = ∞. A
few technical details concerning the norm of the wavefunction are deferred to the appendix
on Groebner bases. Existence of solutions for excited states are also discussed.
The thermodynamic limit of this solution is studied in chapter 5. By proving the
existence of a thermodynamic limit, we can connect the solution of the previous chapter
with the explicit Lieb-Wu solution at half-filling. The results from chapters 4 and 5 can
also be found in cond-mat/0403736 and have been recently accepted for publication in
Communications in Mathematical Physics [Go1].
In chapter 6, we study the Falicov-Kimball model, which is introduced as a limiting
case of the asymmetric Hubbard model. Again, we survey the history and exact results
on the model. We focus on the crystalline and segregated phases, as a preparation for
the presentation of a lower bound for the segregation energy in the Falicov-Kimball model.
Although most of our original work in chapter 7 is valid in any dimension, the one and
two-dimensional case are also discussed, since a more complete rigorous description of the
phase diagram is available.
Chapter 7 presents our lower bound for the boundary term on the ground state energy
for the Falicov-Kimball model at large U , also called the segregation energy. Our motivation
3
was the original proof of segregation given by Freericks, Lieb and Ueltschi, where the authors
show that the ground state energy is bounded above and below by a bulk term plus a
boundary term. A coefficient for the boundary term in the lower bound was only available
for very low densities. In this part of the thesis, we show how to obtain a coefficient for
all densities. The essence of this chapter is the paper published in Journal of Physics A:
Mathematical and General [Go2] or cond-mat/0211183.
Chapter 2
The Hubbard model
2.1
Introduction
The Hubbard model was introduced in 1963, independently by Gutzwiller [Gu] and
Hubbard [H], in an attempt to create the simplest possible model that could describe the
effect of correlations on d-band electrons, and hence explain ferromagnetism in transition
metals. It was originally conceived as a model for spin-half fermions, and we will restrict
our analysis to this case, although we should point out that the analog for bosons (known
as the Bose-Hubbard model) is also widely studied.
In this simplified model, the positions of the electrons are constrained to be in the
sites of a graph, and the electrons hop between any two sites connected by an edge. The
hamiltonian consists of the usual hopping term of the non-interacting tight binding model,
plus an interaction term for every pair of fermions occupying the same site. If we denote
our graph by Λ, and index the sites by integers between 1 and |Λ|, the hamiltonian can be
written
H=−
|Λ|
tij c†iσ cjσ
i,j=1 σ=↑,↓
+
|Λ|
Ui ni↑ ni↓ ,
(2.1)
i=1
where tij are the elements of the hermitian hopping matrix T , c†iσ and ciσ are the usual
creation an annihilation operators for electrons at site i and spin σ, and niσ = c†iσ ciσ is the
number operator at site i. The number of particles is conserved, since H commutes with
4
5
N=
i (ni↑
+ ni↓ ). If the strength of the interaction is constant Ui = U for every site, we
can add a constant term −U N/2 to (2.1) and obtain a hamiltonian that is symmetric with
respect to particles and holes,
H=−
|Λ|
i,j=1 σ=↑,↓
tij c†iσ cjσ
+U
|Λ| i=1
ni↑ −
1 1
ni↓ −
.
2
2
(2.2)
In what follows, we will be working with the symmetric form of the hamiltonian (2.2). From
chapter 3 on, we will fix the filling factor, so both are equivalent and we will use (2.1) for
simplicity. Both the repulsive (U > 0) and repulsive case (U < 0) will be studied. Note
that we assumed the hopping matrix to be independent of the spin. Although this is the
most intuitive case, we shall generalize the model later to include spin dependent hopping
matrices, obtaining the asymmetric Hubbard model. A limiting case of the asymmetric
Hubbard model is known as the Falicov-Kimball model, which we will discuss in chapters 6
and 7.
One of the main goals in studying the Hubbard model is to understand the relation
between the interactions in a system of correlated fermions and its magnetic properties,
particularly at the ground state. We have devoted a section to describe recent results on
the search for conditions under which the model exhibits ferromagnetic behavior in the
ground state.
It is widely believed that the Hubbard model is the right framework to study highTc superconductors. The nearly half-filled Hubbard model with large U is actually the
starting point of Anderson’s resonating valence bond (RVB) theory for high temperature
superconductivity in cuprates [A], which revived the importance of the model in condensed
matter physics.
In higher dimensions, the model is also used to investigate Mott insulator-conductor
transition, since the U = 0 limit reduces to a non-interacting conductor, contrasting with
the case U = ∞ at half filling, where every site is occupied and the system is a Mott
insulator. In one dimension, no metal-insulator transition can occur for T > 0, due to the
Mermin-Wagner theorem. The exact solution for the ground state in 1D also shows the
6
absence of this transition even at T = 0, for any U > 0 [LW1].
In this chapter, we will describe some of the exact results known for the Hubbard
model. The one-dimensional model remains as the only one with an exact solution, which
is the topic of chapter 3. The existence of an exact solution, very unusual in interacting
models, makes the Hubbard model a paradigm to study non-trivial behavior in many-body
problems. It has also led to many new applications of the model, including the study of the
ring-shaped aromatic molecules, as benzene [HL], [LN].
We should point out that a few good reviews on the Hubbard model exist in the literature, and they are much more complete than this brief introduction. In particular, [Mo]
contains a collection of reprints of the major contributions to the model, until the early
90’s. Another collection of reprints, but on the more general topic of exactly solved models
in one dimension can be found in [Ma]. For a review on the known rigorous results on the
Hubbard model, as well as suggestions of open problems of recent interest, we refer the
reader to [L1].
2.2
Symmetries
Let us restrict now to bipartite lattices, which are lattices of the kind Λ = A ∪ B, with
A and B disjoint, and where all the edges connect one element of A with one element of
B, but do not connect elements within A or B. The square lattice is a good example of a
bipartite lattice, as is the linear chain with an even number of sites.
We know that our matrix T has to be hermitian, in order for the hamiltonian to be
self-adjoint. In the general case, tij will be complex numbers, but we can assume T is real
and symmetric in the absence of a magnetic field. In this case, there is a convenient unitary
transformation that takes
c†i↓ → (−1)x ci↓
,
ci↓ → (−1)x c†i↓ ,
where x is even in the sublattice A and odd in B. The operators c†i↑ and ci↑ are left
unchanged. If we denote the hamiltonian in (2.2) by H(U ), the unitary transformation takes
7
H(U ) → H(−U ), but it also changes the number of spin down electrons from N↓ → Na −N↓ .
Therefore, if we solve Schrödinger’s equation for the repulsive case, we get the solution of
the attractive case as well, with a different number of down spins.
If we apply a similar transformation to the spins up, we will map H(−U ) back into
H(U ), and now we have Na − N↑ spins up and Na − N↓ spins down. When we prove
the existence of solution to the Bethe Ansatz equations, we will restrict the number of
electrons to N ≤ Na . There is no loss of generality in this restriction, since we can apply
the sequence of unitary transformations above to map any case where N > Na (above halffilling) to N = 2Na − N < Na . We will also impose the further restriction N↓ ≤ N↑ with
no loss of generality, given the symmetric roles played by spins up and down.
We have the usual SU (2) invariance, with the generators given by the spin operators
Sz =
1
(N↑ − N↓ ) ,
2
S+ =
c†x↑ cx↓
,
S− =
Λ
c†x↓ cx↑ .
Λ
But the particle hole symmetry described above induces another SU (2) symmetry. The
generators are given by transforming back the spin operators of the hamiltonian H(−U ) to
our original basis,
1
Ŝ z = (N↑ + N↓ − Na )
2
,
Ŝ + =
(−1)x c†x↑ c†x↓
,
Ŝ − =
(−1)x cx↑ cx↓ .
Λ
Λ
The operators above are known as the pseudospin operators. Application of the operators S ± and Ŝ ± are necessary to obtain a complete set of eigenstates of the Hubbard model
from the Bethe Ansatz states.
Combining this two SU (2) symmetries, Yang and Zhang [YZ] showed that the group
symmetry of the Hubbard model is SO(4) = SU (2) × SU (2)/Z2 .
2.3
Ferromagnetism and the Hubbard model
Now we will analyze the magnetic properties of eigenstates of the Hubbard model, particularly in the ground state. It is remarkable that we can observe non-trivial magnetic
8
properties in very simple lattices , and slightly more complicated lattices can lead to saturated ferromagnetism. This suggests that the Hubbard model could indeed shed some light
into the theory of ferromagnetism for a system of itinerant electrons.
2.3.1
Lieb’s theorems and ferrimagnetism
If we consider systems of itinerant electrons interacting with a finite potential, the first
example where the spin of the ground state is shown to be proportional to the number of
electrons was given by Lieb [L2], in a paper where two important theorems on the magnetic
properties of the ground state of the Hubbard model are proved.
The first theorem applies for the attractive case U ≤ 0, when the number of electrons N
is even. When U = −∞, we know that the electrons pair up in the ground state, resulting
in total spin S = 0. Also for U = 0 (non-interacting case), we can show that there is a
spin zero ground state. Using spin space reflection positivity, Lieb showed that this result is
actually true for any U ≤ 0, and among the ground states of the system, there is one with
spin S = 0. If U < 0, the S = 0 ground state is unique. The conditions for the theorem
to hold are slightly more general, as Ux is allowed to be a function of the site, provided
Ux ≤ 0.
As a consequence of this theorem, it can be shown that the ground state in a sector
with Sz fixed has minimal pseudospin, that is,
|Λ| − N Ŝ = |Sˆz | = 2
when the number of electrons is even and U < 0. Notice that in principle we could have
any Ŝ between |(|Λ| − N )/2| and |Λ|/2, so the ground state in this sector has the lowest
total pseudospin that is compatible with the value of Ŝz .
To prove that, we notice that for U = −∞, the ground state energy is strictly decreasing
as a function of the number of particles, E(2) > E(4) > · · · > E(2|Λ|). Let the ground state
for N particles be ψ(N ). Suppose there is a level crossings between E(N ) and E(N + 2),
with N +2 ≤ |Λ|, as we reduce the interaction strength to finite negative values. That would
imply a degeneracy in the ground state, since Ŝ + ψ(N ) and ψ(N + 2) would be orthogonal
9
to each other and have the same energy, which can’t be, due to Lieb’s theorem. Therefore,
there can be no crossings, and the energy is ordered by number of particles, for N ≤ |Λ|
and negative U .
Now if we apply S + to ψ(N ), N ≤ |Λ|, we obtain either 0 or a state of N − 2 electrons. But since that state would have an energy lower than ψ(N − 2) we necessarily have
S + ψ(N ) = 0, which proves our claim that the ground state has minimal pseudospin.
The same argument can be used to show absence of crossing between E(N ) and E(N +2)
for N ≥ |Λ|, since it would imply that Ŝ − ψ(N + 2) and ψ(N ) are orthogonal and have the
same energy. Again, it follows that Ŝ − ψ(N ) = 0, which proves the claim for N ≥ |Λ|.
The spin space reflection positivity introduced by Lieb to prove this theorem turned
out to be a powerful technique to derive further rigorous results on the Hubbard model.
In particular, Tian [Ti1], [Ti2] used it to prove that the charged gap of the repulsive case
U > 0, defined by
∆c = E(S = 0, Ŝ = 1) − E(S = 0, Ŝ = 0),
is always larger than the spin excitation gap
∆s = E(S = 1, Ŝ = 0) − E(S = 0, Ŝ = 0),
where E(S, Ŝ) is the ground state in the subspace of fixed spin S and pseudospin Ŝ. Tian’s
result applies also to the periodic Anderson model and the Kondo lattice model.
For the repulsive case U > 0 on a bipartite lattice Λ = A ∪ B, Lieb proved that the
ground state of the half-filled band (N = |Λ|) is (2S + 1)-fold degenerate and has spin
S=
|B| − |A|
,
2
where we assume |B| ≥ |A|. This remarkable theorem not only proves the long assumed fact
that the ground state has spin zero when |B| = |A|, but also shows that for certain simple
lattices, the total spin is extensive, S = CSmax , where Smax = N/2 is the maximum possible
spin for a system of N electrons. In this case, 0 < C < 1, hence the name ferrimagnetism.
Figure 2.1 shows an example with C = 1/3.
10
A
B
Figure 2.1: Bipartite lattice with |B| = 2|A|.
The next figure shows an example in three dimensions, where |B| = 3|A| and C = 1/2.
It is clear that the corresponding d dimensional lattice would have |B| = d|A| and C =
(d − 1)/(d + 1), which approaches 1 when d is large. We will discuss examples of saturated
magnetism (C = 1) next.
Figure 2.2: 3-dimensional example with |B| = 3|A|, or C = 1/2.
2.3.2
Saturated ferromagnetism
The first example of saturated ferromagnetism in the Hubbard model precedes Lieb’s
results, but requires infinite Coulomb repulsion U = ∞. In this case, for a very general
11
class of lattices, Nagaoka [N] and Thouless [T] proved that when N = |Λ| − 1, the ground
state has S = N/2.
A more generic case of ferromagnetism was given by Mielke [Mi1], [Mi2], in an example
of what we call flat-band ferromagnetism, since the models need to have a single-electron
band that is flat at a particular filling. For example, it was shown that the Hubbard model
in the kagome lattice (figure 2.3) exhibits saturated ferromagnetism when N = |Λ/3| [Mi3].
Shortly after Mielke, other models of flat-band ferromagnetism were proposed by Tasaki
[Ta1], Tasaki and Mielke [MT].
Figure 2.3: Kagome lattice.
Although these models constitute a major achievement in the quest for ferromagnetism
in the Hubbard model, they are probably not totally satisfactory in order to describe real
ferromagnets. We expect a good model of ferromagnetism to arise from a subtle interplay
between quantum dynamics and Coulomb repulsion. In the examples of flat-band ferromagnetism, as in Lieb’s unsaturated magnetism, the ground state of the free electron problem
at the specified filling is highly degenerate, having states of many different total spins. The
Coulomb interaction then chooses among those the one with the highest spin, for any positive U . We expect real ferromagnetism to be stable only for U > Uc for some critical value
Uc , whereas for U << Uc paramagnetism should prevail.
12
It is believed that a necessary condition for ferromagnetism is to have U ρ large, where
ρ is the density of states at the Fermi level. Lieb, Tasaki and Mielke achieve that by having
ρ = ∞, whereas Thouless and Nagaoka have U = ∞. It would be interesting (and more
realistic) to have an example where both U and ρ were large but finite. This was first
achieved by Tasaki [Ta2], by adding a few small hopping terms to his own model of flatband ferromagnetism, as it is shown in figure 2.4. The hopping term is −ν 2 s between the
upper sites, ν 2 t between the lower sites, and ν(t+s) in the diagonal bonds connecting them.
At quarter-filling (N = |Λ|/2), this model exhibits saturated ferromagnetism for sufficiently
large U/s and t/s, for ν > 0. Tasaki also proved near flat-band ferromagnetism in other
dimensions [Ta3].
Figure 2.4: Tasaki’s model for near flat-band ferromagnetism.
We conclude this section by pointing out that the search for realistic lattices that explain
magnetic properties of real materials is still very active, and the models we mentioned
constitute a great progress in that direction. They have also motivated the design of new
ferromagnetic materials, such as organic polymers in [ASKA]. All the examples mentioned
above suggest that the Hubbard model is indeed a good framework to study ferromagnetism.
2.4
Attractive case: monotonicity of the energy with spin
Lieb’s theorem for the U < 0 Hubbard model shows us that the ground state has spin
zero for any even number of electrons. In one dimension, a general theorem of Lieb and
Mattis [LM] holds, so we can make an even stronger claim. Let E(S) be the lowest energy
over all states with a given total spin S, we have
E(N/2) > E(N/2 − 1) > · · · > E(0)
(2.3)
13
We should remind the reader that the conditions for the theorem are much more general,
but include the lattice model with nearest neighbor hopping and arbitrary many body
potentials as a particular case.
A natural question that comes up is whether or not (2.3) is valid for any dimension.
The result is certainly true for |U | large, since interaction favors double occupancy of sites.
Therefore, if we consider the ground state with M spins down (M ≤ N/2), or Sz = (N −
2M )/2, only the remaining (N − 2M )/2 spins will be unpaired, resulting in spin S = Sz .
Although we are not able to give a proof that (2.3) is true in any dimension, we will now
prove that at least part of the inequalities have to be satisfied, in the translation invariant
case.
Theorem 2.4.1. Consider the hamiltonian (2.1), with U < 0, in a hypercubic lattice with
|Λ| sites. If the hopping matrix T = (txy ) is translation invariant, the following sequence of
inequalities holds at half filling (N = |Λ|):
E(N/2) > E(N/2 − 1) > · · · > E(N/2 − k),
where k is the largest integer satisfying k < |Λ|1/2 .
Proof. Let us start with a subspace of fixed Sz = N/2 − M , so we have M spins down and
|Λ| − M spins up. Let us consider M < |Λ|1/2 . We can define a basis {ψ α } for one species
of M spinless fermions by
ψα =
c†x |0 >,
x∈α
where α is a set of M distinct sites, and |0 > is the vacuum. Analogously, we define a basis
{φβ } for the N − M spinless fermions, indexed by sets β of N − M distinct sites.
We need to prove that all the lowest energy states on this subspace of fixed Sz have
S = Sz . Suppose this is not true, and there is a lowest energy state with S > Sz . We can
write this state as
Ψ=
Mαβ ψα ⊗ φβ
α,β
where
ψα
corresponds to the state of the spins down and φβ of the spins up. S > Sz implies
that there is no term in this sum where the spins down are all on top of the spins up.
14
Therefore U Ψ > (m − 1)U , where U Ψ is the expectation value of the potential term
|Λ|
i=1 Ui ni↑ ni↓ in this state.
Now consider a state Ψ (T ) obtained by translating all spins down along the lattice by
a translation T ,
Ψ (T ) =
Mαβ ψT α ⊗ φβ ,
α,β
where the translation of the set α is defined in the usual sense, and we are assuming periodic
boundary conditions. If we consider the set of all possible translations, and calculate the
average of U over all such states, we obtain
U T = U
M (|Λ| − M )
,
|Λ|
since the M spins down will see an uniform distribution (in average) of spins up with density
(|Λ| − M )/|Λ|.
Using M < |Λ|1/2 we get
U T < (M − 1)U,
which implies the existence of a state Ψ (T ) with lower potential energy than Ψ. Since f
and f (T ) differ only by a translation, they have the same kinetic energy, and Ψ (T ) will
have lower total energy, which is in contradiction with the fact that Ψ is a ground state in
this subspace.
This result can also be proved by taking the ground state of the kinetic term and
calculating its potential energy. This will be an upper bound for the potential energy of the
true ground state, which can be shown to be less than (M − 1)U if M < |Λ|1/2 .
In particular, for every hypercubic lattice we have
E(N/2) > E(N/2 − 1) > E(N/2 − 2).
This result can easily be extended away from half-filling. In general, the average potential energy over all translations is
U
M (N − M )
,
|Λ|
15
and the condition to be fulfilled is
M (N − M )
> M − 1.
|Λ|
Chapter 3
The Bethe Ansatz
3.1
Introduction
In this chapter, we will discuss the Bethe Ansatz and its applications. Over the past
74 years, the Ansatz has been used to solve many models in one dimension, including the
Heisenberg model, continuum models of bosons and fermions with delta interactions and the
Hubbard model, just to mention a few of the known integrable systems of physical interest.
It has also been used to calculate the eigenvalues of the transfer matrix of the six vertex
model, respecting the condition of neutrality or “ice rule”, which eventually led to an exact
calculation of the residual entropy of two-dimensional ice [L3],[L4], and to the solution of
the F model [L5] and KDP model [L6].
More recently, the six vertex solution has attracted renewed attention due to its connections with the stochastic evolution of the probabilistic cellular automaton (PCA) [GS].
That gives us a measure of the wide range of problems that can be solved used the Bethe
Ansatz, which makes it one of the most important techniques for the study of exactly
solvable models.
Although the original results in this thesis follow from the application of the Bethe
Ansatz to the Hubbard model, we will begin by introducing the solution of the Heisenberg
model. This will serve as a good historical introduction to the topic, since its solution
16
17
motivated the original Ansatz by Bethe [Be]. The goal is also to familiarize the reader with
the technique, without getting distracted by the additional complications that come up in
the solution of a model of fermions.
The Hubbard model will be studied in the sequence. The solution is considerably more
complicated than in the Heisenberg case, since it requires the usual Bethe Ansatz for the
spacial part of the wavefunction, followed by a second Ansatz for the spin part.
3.2
Heisenberg model
The one-dimensional Heisenberg model describes localized spins in interaction with its
neighbors. It can also be interpreted as a model of lattice bosons. We will restrict our
analysis to the spin-1/2 isotropic antiferromagnetic case, but the technique applies also to
ferromagnetic and anisotropic cases.
The hamiltonian for a chain of Na spins with periodic boundary conditions is
H=J
Na
i=1
σi σj = J
Na
(σx,i σx,i+1 + σy,i σy,i+1 + σz,i σz,i+1 ),
(3.1)
i=1
where σx,i , σy,i , σz,i denote the x, y, z components of the Pauli spin matrix σi , multiplied by
a constant to make σx2 = σy2 = σz2 = I. From now on, we will set J = 1. The boundary
condition can be expressed as Na + 1 ≡ 1.
Our goal now is to find the eigenstates of (3.1). Since the hamiltonian commutes with
the z projection of the spin, Sz = i σz,i , we can look for eigenstates with a definite number
of spins up and down. Let m be the number of down spins and Na − m the number of up
spins so that the eigenstate has Sz = Na − 2m. Due to the symmetric roles of spins up and
down, we can restrict to m ≤ Na /2 with no loss of generality.
We can express every eigenstate in a basis {|x1 , ..., xm >}, where 1 ≤ x1 < x2 · · · <
xm ≤ Na specify the position of the down spins. In this basis, the eigenstates are written
ψ = ψ(x1 , . . . , xm )|x1 , ...xm > .
(3.2)
18
Bethe [Be] suggested the following Ansatz for the amplitudes:
ψ(x1 , . . . , xm ) =
AP exp i[kP 1 x1 + · · · + kP m xm ],
(3.3)
P
where k1 , · · · , km are m unequal real numbers between −π and π, and the sum is taken over
all the permutations of m elements. AP are coefficients depending only on the permutation
P . Schrödinger’s equation and the boundary conditions will impose conditions on these
coefficients, which can be expressed as a system of linear equations in the m! variables AP .
The consistency of those will imply constraints on the values of ki , which are called the
Bethe Ansatz equations for the Heisenberg model. Let us now go through the derivation of
these equations.
First we need to determine the energy of each eigenstate. For simplicity, let’s consider the
case m = 2, so that ψ(x1 , x2 ) = A1 ei(k1 x1 +k2 x2 ) + A2 ei(k2 x1 +k1 x2 ) . Consider a configuration
in which the two down spins are not neighbors, something of the form
· · · ↑↓↑ · · · ↑↓↑ · · · .
This configuration is represented by a state |xa , xb >, with xb > xa +1. We want to compute
< xa , xb |Hψ >, for a given eigenstate ψ. By Schrödinger’s equations, this should be equal
a
to Eψ(xa , xb ). Writing the hamiltonian as Na + N
i=1 [σx,i σx,i+1 + σy,i σy,i+1 + σz,i σz,i+1 − 1],
we see that the operator in brackets applied to a configuration where the spins i and i + 1
are either both up or both down is equal to zero. If we denote this operator by Oi , that
means Oi |xa , xb >= 0, except for i = xa − 1, xa , xb − 1, xb . Keeping track of the remaining
terms, we have
< xa , xb |Hψ >= [Na − 8 + 8 cos k1 + 8 cos k2 ]ψ(xa , xb ).
So the energy is given by E(k1 , k2 ) = Na − 8 + 8 cos k1 + 8 cos k2 . It can be easily generalized
to m down spins, and the same argument applied to a configuration with no down spins as
nearest neighbors gives the energy
E(k1 , · · · , km ) = Na − 4
m
(1 − cos ki ).
i=1
(3.4)
19
Now let us consider a configuration with exactly one pair of nearest neighbors down
spins. Back to the m = 2 case, that corresponds to a state |xa , xa + 1 >, representing a
configuration like
· · · ↑↓↓↑ · · · .
For this configuration, we have Oi |xa , xa + 1 >= 0 except at i = xa − 1 or i = xa + 1.
If the energy is given by (3.4), Schrödinger’s equation implies
A1
2eik1 − 1 − ei(k1 +k2 )
= ik
.
A2
2e 2 − 1 − ei(k1 +k2 )
In the arbitrary m case, we can generalize to
2eip − 1 − ei(p+q)
AP
= iq
,
AQ
2e − 1 − ei(p+q)
where P and Q are two permutations satisfying P j = Q(j + 1), P (j + 1) = Qj, P i = Qi
otherwise, and kP j = p, kQj = q. (Q is obtained from P by interchanging j and j + 1).
Defining the function
Θ(p, q) = −2 tan−1
tan (p/2) − tan (q/2) 2
,
we have
AP
= −e−iΘ(p,q) .
AQ
(3.5)
Finally, if P0 is the identity permutation P i = i, and we denote its coefficients by A0 , we
can calculate all the other coefficients by
AP
= (−1)P exp −i
Θ(kj , kl ),
A0
(3.6)
where the sum is taken over all pairs of indices that need to be exchanged to go from the
permutation P to the identity, by permutation of nearest neighbors only. This follows from
repeated applications of (3.5).
Lastly, we need to consider the periodic boundary condition. If P is a permutation
obtained by P from a cyclic permutation, P = (P 2, P 3, · · · , P 1), this condition can be
expressed as
AP = AP eikP 1 Na .
20
Using (3.6), we see that this set of conditions is equivalent to
eikj Na = (−1)m−1 exp −i
m
Θ(kj , kl ),
j = 1, . . . , m.
l=1
Taking the logarithm, we see that one set of solutions is given by
Na kj = 2πIj −
m
Θ(kj , kl ),
j = 1, . . . , m,
l=1
where Ij is an increasing sequence given by
Ij = j −
m+1
.
2
The equations above are known as the Bethe Ansatz equations for the Heisenberg model.
Although we derived them for the isotropic antiferromagnetic case, they can be generalized
to the anisotropic and ferromagnetic hamiltonians, given by
H =−
Na
(σx,i σx,i+1 + σy,i σy,i+1 + ∆σz,i σz,i+1 ),
i=1
∆ ≤ 1. For ∆ = −1, this is equivalent to the hamiltonian (3.1) (just differs by a unitary
transformation), whereas ∆ = 1 corresponds to the isotropic ferromagnetic case. A rigorous
proof that the Bethe Ansatz does indeed give us the ground state for these hamiltonians
can be found in [YY].
We just concluded the presentation of the essence of the Bethe Ansatz. Although the
technique became more sophisticated over the years, as it started to be applied to the
solution of more complicated problems (especially for fermions), all Bethe Ansatz solutions
have the trial wavefunction (3.3) as common feature. For fermions, we mentioned that a
second Ansatz is needed for the spin part, which is not going to be a sum of plane waves
like (3.3). However, the trial solution for the spin wave function will still be of the “Bethe
Ansatz type”, in the sense that it is a product of “single particle states”. We will now
show how the technique can be applied to the solution of the Hubbard model, which is the
subject of most of the original work in this thesis.
21
3.3
1D Hubbard model
The success of the Bethe Ansatz was not restricted to the Heisenberg model; it turned
out to solve continuum models as well. The solution of the boson gas with delta interaction
was first obtained by Lieb and Liniger in [LLi] and Lieb in [L7], where the first paper treats
the general solution for the eigenfunctions and analyzes the ground state in detail, whereas
the latter describes the excitation spectrum. Orthogonality and completeness of the Bethe
Ansatz states were shown by Dorlas [D].
The conjugate symmetry between the spatial part of the wavefunction and the spin part
makes the problem more complicated for fermions. For a system of N spin half particles,
we can think of energy eigenstates with a fixed number of spins up (N − M ) and down (M ).
The case M = 0 is trivial since antisymmetry implies absence of interaction. M = 1 was
solved by McGuire in [McG1] and [McG2]. A major breakthrough in this problem occurred
in 1967, with the solution of the case M = 2 by Flicker and Lieb [FL], which shortly after
was extended to the solution of the general problem by Yang [Y] and Gaudin [G].
The Hubbard model was solved in 1968, by Lieb and Wu [LW1], using a technique
similar to the one used in the solution for fermions with delta interaction. Besides deriving
the Bethe Ansatz equations, they were able to solve them in the thermodynamic limit for
the absolute ground state. We will get back to this solution in the next two chapters. Now
we shall focus on the finite lattice case, and how to obtain the equations that led to the
celebrated solution of the model. The existence of solution to these equations is the most
important original result presented in this thesis, along with the proof that it does indeed
give us the ground state of the system, as Lieb and Wu suggested.
Let us start with the one dimensional Hubbard hamiltonian:
H=−
Na Na
(c†i+1,σ ci,σ + h.c.) + U
ni↑ ni↓ .
i=1
σ
i=1
We are using here periodic boundary conditions (Na + 1 ≡ 1), so we are actually working
on a ring with Na sites. Particle-hole symmetry permits us to restrict to the case N ≤ Na .
Na
a
Since the hamiltonian commutes with N
i=1 ni↑ and
i=1 ni↓ , we can look for energy
22
eigenstates with a fixed number of spins up and down. Let us denote these states by
|M, M , where M is the number of spin-down, and M the number of spin-up electrons.
Spin-up and spin-down symmetry allows us to restrict to the case M ≤ M .
We can write the space part of |M, M in a basis of states of localized electrons,
|M, M space =
f (x1 , . . . , xN )|x1 , . . . , xN ,
(3.7)
1≤xi ≤Na
where |x1 , . . . , xN denotes the state where the electrons occupy the sites x1 , . . . , xN . The
final eigenstate will be a product of this space part with the spin part, which we shall discuss
later.
Our states are defined in the region R = {1 ≤ xi ≤ Na , i = 1, . . . , N }. Given any
permutation Q : {1, . . . , N } → {Q1, . . . , QN }, we can define the region
RQ = {1 ≤ xQ1 ≤ xQ2 ≤ · · · ≤ xQN ≤ Na } ⊂ R.
It is clear that
∪Q RQ = R,
where the union is taken over the N ! possible permutations.
The solution of the 1D Hubbard model is based on a generalized version of the BetheAnsatz [Be]. We know that the energy eigenstates are given by (3.7). The Bethe Ansatz
consists of making the following assumption for the amplitudes f (x1 , . . . , xN ) in each region
RQ :
fQ (x1 , . . . , xN ) =
P
N
AQ {P } exp i
kP j xQj ,
j=1
where k1 < k2 < · · · < kN is a sequence of unequal and ordered numbers and the sum is
taken over all permutations P . The (N !)2 coefficients AQ {P } have to be determined in
order for the state to be well defined in the boundary of two regions RQ and RQ . Also, we
should make sure that the resulting state is antisymmetric with respect to exchange of two
identical particles, and that it satisfies the periodic boundary conditions.
Now suppose for a moment that the density of electrons is low enough to permit the
existence of configurations where no particles are nearest neighbors. If we apply the hamil-
23
tonian to our trial wavefunction, and keep track of the coefficient of such configuration, we
obtain the energy as a function of the ki
E(k1 , . . . , kN ) = −2
N
cos kj ,
j=1
which is the expected result, since within each sector the particles behave as free electrons.
We will require this to be true also for higher densities.
3.3.1
Checking for consistency at the boundaries
The first difficulty comes from checking that the wavefunction is well defined in the
boundary of two sectors. Let us pick an integer 1 ≤ i < N and let j = i + 1. If Q and Q
are two permutations such that Qi = Q j, Qj = Q i and Qk = Q k otherwise, the domains
RQ and RQ have a boundary at xQi = xQj = x. Hence the wavefunction can be given by
either fQ or fQ , and that implies
AQ {P }e[i(···+(kP i +kP j )x+··· )] =
P
Now if
P
AQ {P }e[i(···+(kP i +kP j )x+··· )] .
P
is obtained from P by interchanging P i and P j (just like for Q and Q ), we
can write the condition above as
(AQ {P } + AQ {P })e[i(···+(kP i +kP j )x+··· )] =
(AQ {P } + AQ {P })e[i(···+(kP i +kP j )x+··· )] ,
P
P
and we see that a sufficient condition is
AQ {P } + AQ {P } = AQ {P } + AQ {P }.
(3.8)
Still considering the case of two particles occupying the same site x = xQi = xQj , we
see that Schrödinger’s equation is satisfied if
−(AQ {P } + AQ {P })(e−ikP i + eikP j ) − (AQ {P } + AQ {P })(e−ikP i + eikP j )
+U (AQ {P } + AQ {P }) = −2(cos kP i + cos kP j )(AQ {P } + AQ {P }),
or
AQ {P }(U + eikP i + e−ikP j ) + AQ {P }(U + e−ikP i + eikP j ) =
AQ {P }(e−ikP i + eikP j ) + AQ {P }(eikP i + e−ikP j ).
(3.9)
24
Equations (3.8) and (3.9) can be rearranged to
i(sin kP i − sin kP j )AQ {P } − U AQ {P }/2
U/2 − i(sin kP i − sin kP j )
i(sin kP i − sin kP j )AQ {P } − U AQ {P }/2
AQ {P } = −
U/2 − i(sin kP i − sin kP j )
AQ {P } = −
We can write this equations in a more compact form if we define the N !-dimensional vector
A{P }, whose components are given by AQ {P },
For instance, if P and P are identical, except for P a = i = P b and P b = j = P a,
b = a + 1, the equations above become equivalent to the matrix equations
A{. . . ij . . . } = Yjiab A{· · · ji . . . }.
(3.10)
where
Yijab = −
U/2 + i(sin ki − sin kj )Pab
,
U/2 − i(sin ki − sin kj )
(3.11)
and Pab acts on the vector A{P } by interchanging Qa and Qb.
There are N !(N − 1) equations of this type (N ! permutations times (N − 1) choices of
pairs ab). They give us a prescription to obtain every A{P } from the coefficients for the
identity permutation. But since we can obtain the permutation P by different sequences
of operations on the identity, the first question that comes up is: are the equations (3.10)
consistent?
Figure 3.1 shows two examples of sequences of exchanges that result in the same permutation. This representation in terms of diagrams is presented in Sutherland’s lectures on
the Bethe Ansatz for the Hubbard model [S], from where we borrowed some of the following
diagrams. Consistency of (3.10) with the first diagram implies Yijab Yjiab = 1; for the second,
ab Y bc Y ab = Y bc Y ab Y ab .
we need Yjk
ij ik jk
ik ij
Indeed,
Yijab Yjiab =
(U/2)2 + (sin ki − sin kj )2 (Pab )2
= 1,
(U/2)2 + (sin ki − sin kj )2
2 = 1. Similarly, one can check that the second relation also holds.
since Pab
At first, the task of checking for all consistency relations might seem overwhelming.
However, all equivalent interchanges of four or more indices are consistent if the elementary
25
1
2
1
2
1
2
1
2
3
2
1
3
2
1
1
2
3
1
2
3
Figure 3.1: Consistency relations for exchanges of two and three indices.
two-indices and three-indices consistency relations are satisfied. For example, we can think
of two equivalent ways of obtaining (1432) from (1234), hence we have at least two ways
34 Y 23 Y 12 Y 34 Y 12 A{1234} or Y 23 Y 34 Y 23 A{1234}. Consistency
of calculating A{1432}: Y23
24 21 34 12
34 24 23
implies
34 23 12 34 12
23 34 23
Y24 Y21 Y34 Y12 = Y34
Y24 Y23
Y23
12 and Y 34 commute, we can reverse the order of the two operations,
But if we notice that Y21
34
which corresponds to shrinking part of the diagram, as in figure 3.2, to derive
34 23 12 34 12
34 23 34 12 12
34 23 34
23 34 23
Y24 Y21 Y34 Y12 = Y23
Y24 Y34 Y21 Y12 = Y23
Y24 Y34 = Y34
Y24 Y23 ,
Y23
and consistency is checked. Note that all we used was the consistency of the diagrams of
figure 3.1.
This heuristic argument can be made precise by using a theorem from knot theory due
to Reidemeister ([R], see also [Ka]). The theorem is a result for objects known as links,
which are defined as embeddings of a collection of circles into three dimensional space. Two
26
1
4
3
2
1
4
3
2
1
4
3
2
1
4
3
2
1
2
3
4
1
2
3
4
1
2
3
4
1
2
3
4
Figure 3.2: An example with exchanges of four indices.
links L0 and L1 are called ambient isotopic if there is a continuous family of embeddings
parametrized by L(t) such that L(0) = L0 and L(1) = L1 . A diagram is a projection of the
link in two dimensions, where no more than two points are mapped to the same point in
the plane. If we restrict to links that can be represented by a piecewise linear embedding
composed of straight line segments, also known as tame links, we have the following result:
Theorem 3.3.1 (Reidemeister). Two diagrams represent ambient isotopic links if and
only if one diagram can be transformed to the other by a finite sequence of the following
moves
(a)
(b)
(c)
Figure 3.3: Reidemeister moves.
The moves (a), (b) and (c) above are known as the Reidemeister moves.
Notice that the Reidemeister moves are local, and in each step the rest of the diagram
27
is left unchanged. All segments in figure 3.3 are continuous, and the discontinuity at the
crossings is just a way of showing which segment is on top.
We can create links by representing our diagrams of permutations in three dimensional
space. For example, for the diagram in figure 3.2, we draw the path corresponding to number
1 in the plane z = 0. The second path is drawn in the plane z = 1, and so on. Strictly
speaking, to obtain links we would have to connect the endpoints, but this is equivalent
to keeping them fixed, so we will abuse notation and call this objects with fixed endpoints
links. The diagrams drawn in 3.2 are the two dimensional projection of these links, which
are tame by construction.
It is clear that two diagrams connecting the same initial and final permutations will be
ambient isotopic, since each line can move independently of each other in different planes,
and only the endpoints are fixed. By Reidemeister’s theorem,we can obtain one from the
other by a finite sequence of moves like in 3.3. In our case, it is clear we don’t need to
consider moves of the type (a). Hence, equivalence of the diagrams transformed into each
ab Y bc Y ab = Y bc Y ab Y ab )
other by moves (b) and (c) (which follows from Yijab Yjiab = 1 and Yjk
ij ik jk
ik ij
implies equivalence of any pair diagrams corresponding to ambient isotopic links.
With that, we conclude the proof of consistency of A{P } = Yjiab A{P }. These equations
give us a prescription to calculate every A{P } from the vector corresponding to the identity
permutation, A{I}. Now it is time to consider the boundary conditions.
3.3.2
Periodic boundary condition
If Q = {Q1, Q2, . . . , QM } and Q = {QN, Q1, . . . , Q(n − 1)} are permutations that
differ only by a cyclic permutation, and x1 ≤ · · · ≤ xN ≤ x1 + Na , the periodic boundary
conditions can be written as
fQ ({x2 , x3 , . . . , xN , x1 + Na }) = fQ ({x1 , x2 , . . . , xN −1 , xN }),
where fQ ({x1 , x2 , . . . , xN }) is just a notation for fQ (xn1 , xn2 , . . . , xNN ), with Qni = i. A
sufficient condition is given by
AQ {P }ei[kP 1 x2 +···+kP N (x1 +Na )] = AQ {P }ei[kP 1 x1 +···+kP N xN ] ,
28
or
AQ {P } = AQ {P }eiNa kP N ,
(3.12)
where P = {P N, P 1, . . . , P (N − 1)}. Defining new vectors A {P } = P −1 A{P }, we can
write (3.12) as
A {P } = A{P }eiNa kP N .
(3.13)
Let us first solve (3.12) for P = I = {1, . . . , N } and P = {N, 1, . . . , N − 1}.
We can define a sequence of permutations interpolating between I and P :
P (0) = I = {1, . . . , N },
P (1) = PN −1,N P (0) = {1, . . . , N − 2, N, N − 1},
..
.
P (N ) = P1,2 P (N −1) = {N, 1, . . . , N − 1} = P .
Since P (1) differs from the identity just by an elementary permutation of the last two
indices, we have
−1,N
A{I} = XN −1,N A{I},
A {P (1) } = PN −1,N A{P (1) } = PN −1,N YNN−1,N
where Xij ≡ Pij Yijij .
We also have
−2,N −1
A{P (1) } =
A {P (2) } = PN −1,N PN −2,N −1 A{P (2) } = PN −1,N PN −2,N −1 YNN−2,N
−2,N −1
PN −1,N A {P (1) }.
= PN −1,N PN −2,N −1 YNN−2,N
Using (3.11) for Y , and noting that PN −1,N PN −2,N −1 PN −1,N = PN −2,N , we obtain
−2,N A {P (1) } = XN −2,N A {P (1) }.
A {P (2) } = PN −2,N YNN−2,N
Notice that this is the first time we introduce the elementary permutation between nonnearest neighbors PN −2,N .
Proceeding by induction, the following recursion relation is obtained:
A {P (n) } = XN −n,N A {P (n−1) } ,
n = 1, 2, . . . , N
,
29
which implies
−1
X
A{I} = eiNa kN A{I}.
A {P } = ΠN
lN
l=1
Repeating the same argument when adding Na to each of the other variables x2 , . . . , xN ,
we see that the periodic boundary condition is satisfied if the following identities hold:
ΠN
l=j+1 Xlj
Πj−1
l=1 Xlj A{I} = λj A{I}
,
j = 1, . . . , N
,
(3.14)
where
Xlj = Plj Yljlj = −
U Plj /2 + i(sin kl − sin kj )
,
U/2 − i(sin kl − sin kj )
and λj = eiNa kj . We need to solve the N eigenvalue equations (3.14), with the same
eigenvector A{I}.
To solve these equations, we need to determine the symmetry type of the eigenfunction.
In general, the operators Pij define a (N !)2 dimensional representation of the group SN . For
every irreducible representation R that determines the symmetry type of the wavefunction
we will get different values of λj , so the eigenvalues are functions of the representation R,
λj = λj (R).
For example, the boson case can be obtained by taking R as the identity representation
[N ], where by this notation we mean the wavefunction is symmetric with respect to all N
variables. This corresponds to a Young’s diagram of one row of length N . The equations
have to be solved with Pij = 1. On the other extreme, if we take the total antisymmetric
representation [1N ], Pij = −1 (one column Young’s diagram), we see that Xij = 1 and
exp ikj Na = 1. This gives us the solution of the one component problem, when all the spins
are up.
If we fix the symmetry part of the spatial part of the wavefunction, the spin part will need
to have the conjugate symmetry to result in the antisymmetry of the total wavefunction.
The equivalent of (3.14) for the spin part will then be
ΠN
l=j+1 Xlj
Πj−1
l=1 Xlj Φ = µj Φ
,
j = 1, . . . , N
,
(3.15)
30
where Xlj is obtained from Xlj by inverting all signs in Plj , to account for the conjugate
symmetry:
Xlj =
U Plj /2 − i(sin kl − sin kj )
.
U/2 − i(sin kl − sin kj )
The eigenvalues λj and µj are related by µj (R) = λj (R̂), where R and R̂ are conjugate
representations (if R = [1N ], R̂ = [N ] and so on).
Suppose we are looking for an eigenstate with a given value for Sz , the z projection
of the total spin. That means we fix the number of electrons with spin down as M , with
the remaining N − M spins up, resulting in Sz = (N − 2M )/2. Now we also fix the total
spin as S = Sz , since we can obtain all other states by applying spin lowering operators.
The Young’s diagrams corresponding to these representations are shown in figure 3.4. Fixing S and Sz is equivalent to solve (3.15) for the symmetry [N − M, M ]. The conjugate
representation for the spatial part is [2M 1N −2M ].
N=10 , M=4
M N−2M
[2 1
]
[N−M,M]
Figure 3.4: Symmetry types and Young’s diagrams for N=10, M=4.
Obtaining the solutions to (3.15) is not that simple, so let’s try to use the M = 1 and
M = 2 case to guide us to the solution of the general problem.
3.3.3
M=1
In the case M = 1, the states are labeled by the position of the single down spin.
Given the representation in which we are working, it makes sense to write the eigenstate
as Φ = i φi |i >, where i is the position of the down spin in the string of spins. The first
31
:
operator to act on Φ is X(j−1)j
Φ = · · · + (tj−1,j φj−1 + rj−1,j φj )|j − 1 > +(tj−1,j φj + rj−1,j φj−1 )|j > + · · · ,
X(j−1)j
where
tij =
−i(sin ki − sin kj )
U/2 − i(sin ki − sin kj )
,
rij =
U/2
.
U/2 − i(sin ki − sin kj )
Only the amplitudes of |j − 1 > and |j > were affected so far, all the others remain
constant due to the symmetry. Also note that the amplitude of |j − 1 > will not be affected
by the other operators, so the eigenvalue equation implies
µj φj−1 = tj−1,j φj−1 + rj−1,j φj .
(1)
We will define φj
as the new amplitude of the state |j >,
(1)
φj = tj−1,j φj + rj−1,j φj−1 .
to X(j−1)j
Φ results in
Now applying the operator X(j−2)j
(1)
(1)
· · · + (tj−2,j φj−2 + rj−2,j φj )|j − 2 > +µj φj−1 |j − 1 > +(tj−2,j φj + rj−2,j φj−2 )|j > + · · · .
As we mentioned before, the amplitude for |j − 1 > is already the correct one. But since
the amplitude for |j − 2 > will be no longer changed, we should also have
(1)
µj = tj−2,j + rj−2,j
(2)
We also define φj
φj
φj−2
.
as the new amplitude of |j >.
If we continue this process, after the k-th iteration, which corresponds to the application
(with the convention that 0 ≡ N ), we will have the recursion
of the operator X(j−k)j
relations
(k−1)
µj = tj−k,j + rj−2,j
(k)
φj
(k−1)
= tj−k,j φj
φj
,
(3.16)
+ rj−k,j φj−k
(3.17)
φj−k
32
(k)
Using (3.16) (one iteration further) to calculate φj , and substituting into (3.17), we obtain
the ratio between consecutive amplitudes
−U + U µj /(µj − 1) − 2i(sin kj−k − sin kj )
φj−k−1
.
=
φj−k
U µj /(µj − 1) − 2i(sin kj−k−1 − sin kj )
(3.18)
Since the eigenvectors are the same for all the N eigenvalues µj , this ratio should be independent of j. This is only possible if U µj /(µj − 1) + 2i sin kj is a constant. We define a
constant Λ by
U µj
+ i sin kj = U/4 + iΛ,
2(µj − 1)
and solving for the eigenvalue we have
µj =
i sin kj − iΛ − U/4
≡ σj (Λ).
i sin kj − iΛ + U/4
(3.19)
For those familiar with the Bethe Ansatz equations, this is exactly the point where the
rapidities Λ come up.
The ratio (3.18) becomes
i sin kj−k − iΛ + U/4
φj−k−1
=
φj−k
i sin kj−k − iΛ − U/4
from where we can deduce the expression for all φj in terms of φ1 . Normalizing the state
to have φ1 = 1, we obtain
Fj (Λ) ≡
j−1
l=1
i sin kl − iΛ − U/4
.
i sin kl+1 − iΛ + U/4
(3.20)
The following relation will be useful:
Fj (Λ) =
i sin k1 − iΛ − U/4
σ1 · · · σj−1 .
i sin kj − iΛ + U/4
(3.21)
If we start with the eigenstate given by · · · + Fj (Λ)|j >, the iterated amplitude of the
state |j > can be calculated using (3.16),
(k)
Fj (Λ) =
i sin kj−k−1 − iΛ − U/4
Fj (Λ)
Fj−k−1 =
,
i sin kj − iΛ + U/4
σj−1 · · · σj−k
where for the second equality, we make use of (3.21).
33
Let us consider now the case where j = N . By the end of the N − 1 iterations, all
amplitudes are forced by construction to satisfy the eigenvalue equation, except for the
amplitude of |j >, that will be given by
(N −1)
FN
(Λ) =
FN (Λ)
.
σ1 · · · σN −1
By the eigenvalue equation, the above expression should be equal to σN FN (Λ), which results
in an additional constraint:
N
j=1
N
i sin kj − iΛ − U/4
= 1.
σj =
i sin kj − iΛ + U/4
j=1
This equation, combined with the eigenvalue equations
eikj Na =
i sin kj − iΛ − U/4
i sin kj − iΛ + U/4
,
j = 1, . . . , N,
for a set of N + 1 coupled equations on the N + 1 variables kj and Λ. Those are the Bethe
Ansatz equations for the Hubbard model, when M = 1.
Before moving on to the case M = 2, let us first rewrite the recursion relations (3.16)
and (3.17) in terms of Fj (Λ):
rj−k,j
Fj
= σj Fj−k ,
σj−1 · · · σj−k+1
Fj
+ rj−k Fj−k =
.
σj−1 · · · σj−k
tj−k Fj−k + rj−k,j
(3.22)
Fj
σj−1 · · · σj−k+1
(3.23)
As we mentioned before, equation (3.22) shows how the |j − k > component has the right
amplitude after the k-th iterations. Equation (3.23) shows how the amplitude of |j > is
prepared for the (k + 1)-th iteration, in order to “correct” the amplitude of |j − k − 1 >.
All this operations have to be done in perfect synchrony, and we just showed this can only
happen if the Bethe Ansatz equations are satisfied.
3.3.4
M=2
Let’s consider two spins flipped, or M = 2. Now the states are labeled by two indices,
one for each spin down. We will try a solution of the form
Φ=
i,j
φij |ij >,
34
with
φij = A1 Fi (Λ1 )Fj (Λ2 ) + A2 Fi (Λ2 )Fj (Λ1 ) = A1 Fi (Λ1 )Fj (Λ2 ) + (1 ↔ 2),
where A1 , A2 , Λ1 , Λ2 are constants to be determined, and Fj is defined as in the M = 1
case. The notation (1 ↔ 2) just means we exchange the indices 1 and 2 (in Λ and A only)
in the preceding expression, and that will just save us some space.
It is clear that this trial solution is similar to the Bethe Ansatz, where we try a wavefunction that is the product of one-particle eigenstates, summed over all permutation. The
only difference here is that the function F (Λ) plays the role of the free particle state eikx .
Let’s consider again the case where j = M , but keeping the index j in the equations,
acts only on the states
just to simplify the generalization later on. The operator Xj−1,j
where one of the spins is at j − 1 or j, but not both. In particular, the amplitude of
|j − 1, j > is unchanged. If there is only a down spin at j − 1, equation shows us that
Fj (Λ2 ) is multiplied by σj (Λ2 ); if there is a down spin at j, equation shows us that Fj (Λ2 )
is divided by σj−1 (Λ2 ).
It works exactly like for M = 1, but the analogy doesn’t hold for |j, j − 1 >, since its
amplitude is unchanged and equal to
A1 Fj−1 (Λ1 )Fj (Λ2 ) + (1 ↔ 2).
This is a problem, because the amplitude of the state |j − 2, j − 1 > has the new value
A1 Fj−2 (Λ1 )Fj−1 (Λ2 )σj (Λ2 ) + (1 ↔ 2),
and the application of Xj−2,j will result in the wrong amplitude for this state, unless the
amplitude of |j − 1, j > were given by
A1
Fj (Λ1 )
Fj−1 (Λ2 )σj (Λ2 ) + (1 ↔ 2).
σj−1 (Λ1 )
We have to hope that A1 and A2 are related in such a way that the following identity
holds:
A1 Fj−2 (Λ1 )Fj−1 (Λ2 )σj (Λ2 ) + (1 ↔ 2) =
A1 Fj (Λ1 )Fj−1 (Λ2 )σj (Λ2 )
+ (1 ↔ 2).
σj−1 (Λ1 )
(3.24)
35
If this is the case, it follows from (3.23) that the final amplitude for |j − 2, j − 1 > will be
given by
σj (Λ1 )σj (Λ2 )φ(j − 2, j − 1),
and we identify the eigenvalue
µj = σj (Λ1 )σj (Λ2 ).
So far, so good. However, that only solves the problem of the amplitude |j − 1, j >.
What about the amplitude of |j − k, j > for k > 1? After the k-th iteration, it is given by
A1 Fj−k (Λ1 )Fj (Λ2 )
+ (1 ↔ 2),
σj−1 (Λ2 ) · · · σj−k+1 (Λ2 )
whereas the correct value to go into the eigenvalue equation for the amplitude of |j − k −
1, j − k > in (k + 1)-th iteration would be
A1 Fj (Λ1 )Fj−k (Λ2 )σj (Λ2 )
+ (1 ↔ 2).
σj−1 (Λ1 ) · · · σj−k (Λ1 )
For our scheme to work, we need the following identity valid for k = 1, . . . , j − 1:
A1 Fj (Λ1 )Fj−k (Λ2 )σj (Λ2 )
A1 Fj−k (Λ1 )Fj (Λ2 )
+ (1 ↔ 2) =
+ (1 ↔ 2).
σj−1 (Λ2 ) · · · σj−k+1 (Λ2 )
σj−1 (Λ1 ) · · · σj−k (Λ1 )
(3.25)
Notice that (3.25) is all we need, since (3.24) is just a particular case when k = 1. The
ratio between A1 and A2 has to be
Fj−k (Λ2 )Fj (Λ1 )/σj−1 (Λ1 ) · · · σj−k+1 (Λ1 ) − Fj (Λ2 )Fj−k (Λ1 )σj (Λ1 )/σj−1 (Λ2 ) · · · σj−k (Λ2 )
A1
=−
A2
(1 ↔ 2)
After using (3.20), (3.19), and canceling common factors between the denominator and
numerator, we end up with
i(Λ1 − Λ2 ) + U/2
A1
.
=
A2
i(Λ1 − Λ2 ) − U/2
(3.26)
This expression is independent of k, so the equation above implies that (3.25) holds for all
k. That ensures that all amplitudes satisfy the eigenvalue equation, except for those of the
form |j, j − k >, which we need to check. The condition for these amplitudes is
A1 FM (Λ1 )FM −k (Λ2 )σM (Λ2 )
+ (1 ↔ 2) = σM (Λ1 )σM (Λ2 )[A1 FM −k (Λ1 )FM (Λ2 ) + (1 ↔ 2)].
σM −1 (Λ1 ) · · · σ1 (Λ1 )
36
One sufficient condition is
A1 FM (Λ1 )FM −k (Λ2 )σM (Λ2 )
= σM (Λ1 )σM (Λ2 )A2 FM −k (Λ2 )FM (Λ1 ),
σM −1 (Λ1 ) · · · σ1 (Λ1 )
or
M
A1
=
σj
A2
(3.27)
j=1
Consistency between this ratio and (3.26) implies
N
2
−iΛβ + iΛα + U/2
i sin kj − iΛα − U/4
=−
i sin kj − iΛα + U/4
−iΛβ + iΛα − U/2
j=1
,
β = 1, 2
,
β=1
which combined with the eigenvalue equations
ikj Na
µj = e
M
i sin kj − iΛβ − U/4
=
i sin kj − iΛβ + U/4
,
j = 1, . . . , N
,
β=1
are called the Bethe Ansatz equations, for M=2.
3.3.5
General solution
At this point, we can make the generalization for an arbitrary number M of flipped
spins. The solution of equations (3.15) is of the form
Φ = φ(y1 , . . . , yM )|y1 , . . . , yM ,
where y1 < y2 < · · · < yM are the coordinates of the M down spins along the chain of spins,
and φ(y1 , . . . , yM ) is given by the second Ansatz
φ(y1 , . . . , yM ) =
AP ΠM
i=1 F (ΛP i , yi ),
P
where the sum is taken over all permutations of M numbers. The function F is given by
F (Λ, y) =
y−1
j=1
i sin kj − iΛ − U/4
.
i sin kj+1 − iΛ + U/4
If Q = {Q1, . . . , Qi, Qj, . . . , QM } and Q = {Q1, . . . , Qj, Qi, . . . , QM } are two permutations that differ only by the exchange Qi ↔ Qj, their coefficients are related by the
analogous of (3.26),
i(ΛQi − ΛQj ) + U/2
AQ
.
=
AQ
i(ΛQi − ΛQj ) − U/2
37
Now if Q = {Q2, . . . , QM, Q1} is obtained from Q by a cyclic permutation, consecutive
applications of the identity above imply
i(ΛQ1 − ΛQβ ) + U/2
AQ
,
=−
AQ
i(ΛQ1 − ΛQβ ) − U/2
M
β=1
which should also be equal to
N
σj (ΛQ1 ),
j=1
by the same argument used to derive (3.27). The eigenvalues are again given by
M
β=1 σj (Λβ ),
and we can finally write the set of equations
eiNa kj
=
M
σj (Λβ )
,
β=1
N
j=1
σj (Λα ) = −
M
−iΛβ + iΛα + U/2
−iΛβ + iΛα − U/2
.
β=1
The N + M coupled equations above are known as the Bethe Ansatz equations for the
Hubbard model. The goal of this section was to give a somewhat didactic derivation of these
equations, in order to prepare the reader for the presentation of our original results in the
following chapters. In Chapter 4 we prove the existence of real and ordered solutions for the
Bethe Ansatz equations, and show that the Ansatz can be used to determine the absolute
ground state of the system. In Chapter 5, we show the existence of a thermodynamic limit
to the distribution of k-s and Λ-s, connecting the finite lattice solution with the celebrated
solution of the model in the thermodynamic limit.
Chapter 4
Existence of solution to the Bethe
Ansatz equations
4.1
Introduction
In 1968, Lieb and Wu [LW1] solved the one-dimensional Hubbard model [H], and showed
the absence of Mott transition in the ground state. Many technical details of the exact
solution for the ground state were left out of the original paper and were recently published
[LW2], motivated by the increasing importance of the 1D Hubbard model in Condensed
Matter Physics.
Despite being one of the most studied models to describe interacting electron systems,
many questions concerning the Hubbard model remain as open problems, and the 1D Hubbard model is the only one for which an exact solution was found. Still, many results on
the Bethe Ansatz to the Hubbard model rely on unproved hypotheses. The purpose of
this paper is to answer some of the questions raised in the recent paper by Lieb and Wu,
concerning the existence of a solution to the Bethe Ansatz equations for the 1D Hubbard
model on a finite lattice, and whether or not this solution gives us the ground state. We
will also discuss the solution in the thermodynamic limit.
Let us briefly summarize the technique presented in the last chapter to obtain the
38
39
eigenstates of the Hubbard hamiltonian. We look for eigenstates of the form
|M, M space =
f (x1 , . . . , xN )|x1 , . . . , xN ,
(4.1)
1≤xi ≤Na
where the amplitudes f (x1 , . . . , xN ) are given in each region RQ by the Bethe Ansatz
fQ (x1 , . . . , xN ) =
N
AQ {P } exp i
kPj xQj ,
(4.2)
j=1
P
We need to make sure that the resulting wavefunction satisfies Schröedinger’s equation,
the periodic boundary condition, and that it is antisymmetric. Those restrictions result
in a system of linear and homogeneous equations for the coefficients AQ {P }. The validity
of the Bethe Ansatz is determined by the consistency of the equations. Fortunately, these
equation were studied in detail by Gaudin [G] and Yang [Y], in the context of a system
of fermions with delta interaction. The only difference here is that we are working on a
lattice, and not in the continuum, so the kj that appear on the system of equations have to
be substituted by sin kj , which does not interfere much with the algebraic analysis.
The analysis performed in [G] and [Y] to verify the existence of a non-trivial solution
for the coefficients AQ {P } was reproduced in the last chapter for fermions in a lattice. It
results in a set of conditions that can be written as the generalized Bethe Ansatz equations
ikj Na
e
M
i sin kj − iΛβ − U/4
=
i sin kj − iΛβ + U/4
,
j = 1, . . . , N
,
(4.3)
β=1
N
M
−iΛβ + iΛα + U/2
i sin kj − iΛα − U/4
=−
i sin kj − iΛα + U/4
−iΛβ + iΛα − U/2
j=1
,
α = 1, . . . , M
, (4.4)
β=1
M
in terms of our original set {kj }N
j=1 , and a set of auxiliary parameters {Λα }α=1 .
Defining
θ(x) = −2 tan−1 (2x/U )
and taking the logarithm of the equations above, we obtain the Lieb-Wu equations
Na kj = 2πIj +
j
θ(2 sin kj − 2Λβ ),
β
θ(2 sin kj − 2Λα ) = 2πJα −
β
θ(Λα − Λβ ),
j = 1, . . . , N ,
(4.5)
α = 1, . . . , M .
(4.6)
40
The coefficients Ij are integers if M is even, and half-integers otherwise. Similarly, Jα are
integers if N − M = M is odd and half-integers otherwise. We will restrict to the case
where M is odd, and N is even, since in this case the ground state is unique for every U .
For the ground state, we can choose the Ij and Jα that will give us the correct solution in
the limit U → ∞. In this limit, the equations decouple, and the choice
Ij = j −
N +1
2
,
Jα = α −
M +1
2
(4.7)
implies a solution for fQ of the form of a Slater determinant of plane waves with wavenumbers kj = 2πIj /Na , which is the unique ground state (provided N is even). Also note that
in order to satisfy (4.6), each Λα has to be proportional to U in this limit.
In the thermodynamic limit, taking Na , N, M, M → ∞, keeping their ratios fixed, (4.5)
and (4.6) become integral equations. The analysis of the resulting equations led to the well
known solution of the 1D Hubbard model [LW1], which includes a rigorous expression for
the ground state at half filling. Some excited states can be obtained by different choices of
integers/half-integers for Ij and Jα . It is also known that not all states are of the Bethe
Ansatz form. However, Essler, Korepin and Schoutens [EKS] showed that the remaining
states can be obtained from the Bethe Ansatz states, by using the SO(4) symmetry of the
model.
Some questions concerning the finite lattice case remained as open problems and the
purpose of this paper is to solve them rigorously. In particular, we are going to show that the
equations (4.5), (4.6) do indeed have a solution, and that it is ordered in j and α (kj+1 > kj ,
Λα+1 > Λα ). The ordering is important in the derivation of the integral equations in the
thermodynamic limit. We also show that there is a continuous curve of solutions extending
from U = ∞ to any U > 0. We are ready to state our main result:
Theorem 4.1.1. The equations (4.5) and (4.6) have at least one real solution satisfying
−π ≤ k1 < · · · < kN ≤ π,
Λ1 < · · · < ΛM .
41
Furthermore, there is a continuous curve of solutions defined by γ(t) : [0, 1] → RN +M +1 ,
where
γ(t) = (k1 (t), . . . , kn (t), Λ1 (t), . . . , ΛM (t), U (t)),
extending from any U > 0 to U = ∞.
We know that the state obtained by the solution to (4.5), (4.6) corresponds to the true
ground state in the limit U = ∞, since we chose the coefficients Ij and Jα in order to match
these states. To prove that the ground state is given by the Bethe Ansatz solution for any
positive U , we need a few additional facts. Since M and M are odd, the ground state
is non-degenerate for any U [LW2], and our result shows that starting with the solution
at U = ∞, and going along the curve γ(t), the state is indeed the ground state, provided
f (x1 , . . . , xN ) defined by (4.2) does not vanish for all (x1 , . . . , xN ).
In section 4, we prove that the generalized Bethe Ansatz equations have a solution that
is algebraic in U . Since the Bethe Ansatz wavefunction is a rational function in the variables
zj = eikj , Λα and U (see [W]), its norm will have at most finitely many zeros and poles as a
function of U . We show that we can redefine our wavefunction at these problematic points,
in order to obtain the physical solution for all U . We use the continuity of the ground state
energy and the fact that there is always a finite gap to the next excited state. With this
we conclude the study of the ground state of the finite lattice.
Once the existence of a solution that corresponds to the ground state of our finite system
is shown, our next goal is to determine what happens with the set of k-s and Λ-s in the
thermodynamic limit. As we increase the lattice size, a sequence of distributions measuring
the density of k-s and Λ-s can be defined. In the thermodynamic limit, we show that a
subsequence of these distributions converge to the solution of the integral equations derived
by Lieb and Wu.
We are particularly interested in the absolute ground state at half-filling, in which case
an exact solution can be found, including an explicit expression for the ground state energy.
In this case, we show that the whole sequence converges to the explicit solution obtained
in [LW1].
42
We did not prove here the uniqueness of the solution at the thermodynamic limit away
from half-filling, and leave this as an open problem. This would imply that the whole
sequence converges, as is the case at half-filling.
4.2
Existence of a solution
The goal in this section is to prove the existence of solutions to the Lieb-Wu equations
that are ordered in j and α (kj+1 > kj , Λα+1 > Λα ). The proof will follow from Brouwer’s
fixed point theorem, but first we need to define a convenient map to apply the theorem.
4.2.1
Defining the map
Because of the symmetry of the coefficients Ij and Jα with respect to the origin, we can
look for symmetric solutions of the form kj = −kN −j+1 , and Λα = −ΛM −α+1 . For instance,
if kj = k solves (4.5) for a given Ij = I, kj = −k would also solve it for Ij = −I. Therefore,
our set of k-s and Λ-s are totally defined by
k = {kj } ,
N
+ 1, . . . , N ,
2
M +3
,...,M .
α=
2
j=
Λ = {Λα } ,
The other components of k and Λ are obtained by the symmetry condition. In particular,
Λ(M +1)/2 = 0.
We can define the map φ : R 2 +
N
M −1
2
→R2+
N
M −1
2
by
φ(k, Λ) = (k , Λ ).
where Λ is defined by
θ(2 sin kj − 2Λα ) = 2πJα −
j
for (M + 3)/2 ≤ α ≤ M , and k by
β
θ(Λα − Λβ ),
(4.8)
43
kN +1 = max{0,
2
whereas for i >
N
2
1
(2πI N +1 +
θ(2 sin k N +1 − 2Λβ ))},
2
2
Na
β
+1
,
ki = max{ki−1
1
(2πIi +
θ(2 sin ki − 2Λβ ))}.
Na
β
The function
Λ (k, Λ)
is well defined since the left side of (4.8) is continuous and strictly
increasing in Λα . It is also clear that k (k, Λ) is continuous, so the map φ(k, Λ) = (k , Λ )
is continuous.
It is also important to observe that
Λα < CN,U ,
where CN,U is a constant for each value of N and U . Indeed, from (4.8) we have
θ(2 sin kj − 2Λα ) ≤
j
which implies
π
π 4
4Λα
< tan
−
+ ,
U
2
2N
U
or
Λα <
4.2.2
(M − 1)
2π + M π ≤ (N − 1)π,
2
π
π U
tan
−
+ 1 ≡ CN,U .
4
2 2N
The fixed point theorem
Now we are ready to state the result that allows us to use Brouwer’s fixed point theorem.
Lemma 4.2.1. If we define the domain Ω ⊂ R 2 +
N
M −1
2
by
0 ≤ k N +1 ≤ . . . ≤ kN ≤ π,
2
0 ≤ Λ M +3 ≤ . . . ≤ ΛM ≤ CN,U ,
2
we have
φ(Ω) ⊂ Ω.
44
Proof. By definition,
.
0 ≤ kN +1 ≤ · · · ≤ kN
2
Also,
≤
kN
Regarding Λα , we have
since Jα > 0 and −
β
N
1 N − 1 2π
<π
≤ π.
Na
2
Na
θ(2 sin kj − 2Λα ) > 0,
j
θ(Λ − Λβ ) is strictly increasing with Λ (and equal to 0 at Λ = 0).
Since the left side is also increasing with Λα and vanishes for Λα = 0, we have
Λα > 0.
Furthermore, the right side of (4.8) is increasing in Λα and Jα+1 > Jα , whereas the left
side is increasing in Λα . Restricting the map φ to Ω, this implies
Λα+1 ≥ Λα .
Combining these results with ΛM < CN,U we conclude the proof of the lemma.
Since Ω is compact and convex (in particular, it has no “holes”), there is a continuous
and invertible map
ψ :Ω→D2+
N
where D 2 +
N
M −1
2
M −1
2
,
is a disk in (N + M − 1)/2 dimensions. The map
ψ ◦ φ ◦ ψ −1 : D 2 +
N
M −1
2
→D2+
N
M −1
2
is continuous and takes the disk into the disk, and by Brouwer’s fixed point theorem, it
must have a fixed point. Since ψ is one-to-one, the fixed point of ψ corresponds to a fixed
point of φ.
45
4.2.3
The fixed point is a solution to the Bethe Ansatz equations
Since we artificially introduced the maximum in the definition of k , we need to show
that the fixed point (k∗ , Λ∗ ), such that
φ(k∗ , Λ∗ ) = (k∗ , Λ∗ ),
is indeed a solution to the original equations.
In other words, we need to show that
∗
.
0 < k∗N +1 < · · · < kN
2
Taking kj = 0, we have
kj =
2πIj
>0
Na
⇒
k∗N +1 > 0.
2
That also implies
k∗N +1 =
2
1 θ(2 sin k∗N +1 − 2Λβ ) .
2πI N +1 +
2
Na
2
β
Let us consider now
1 θ(2 sin k∗N +3 − 2Λβ ) .
2πI N +3 +
k∗N +3 = max k∗N +1 ,
2
Na
2
2
2
β
Assuming that k∗N +3 = k∗N +1 , we have
2
2
2π > k∗N +1 ,
k∗N +3 = max k∗N +1 , k∗N +1 +
Na
2
2
2
2
which is a contradiction. Therefore, k∗N +3 > k∗N +1 , and proceeding by induction we get
2
2
∗
.
k∗N +1 < k∗N +3 < · · · < kN
2
2
We have therefore concluded the proof that the fixed point (k∗ , Λ∗ ) satisfies (4.5) and
(4.6).
46
4.3
Existence of a continuous curve of solutions
We want to analyze now what happens to our solutions as U changes. Our goal is to
show that there is a connected set of solutions extending from any U > 0 to U = ∞.
From what we have seen in the last section, we know that on the part of the boundary
of Ω defined by k N +1 = 0 or ki+1 = ki , the vector field defined by
2
vU (k, Λ) = φU (k, Λ) − (k, Λ)
has a normal component that always points inward. Also, for Λα = 0 we have Λα > 0, and
Λα = Λα+1 implies Λα < Λα+1 . Therefore, the vector field actually points inward in all the
boundary of Ω.
All we need to prove then is the following:
Proposition 4.3.1. Given any family of continuous vector fields
vt : Dn → Rn
,
t ∈ [0, 1],
such that the normal component of vt points inward on the boundary of the disk Dn , there
is a connected subset of Dn × [0, 1] in which vt (k, Λ) = 0, extending all the way from t = 0
to t = 1.
Proof. The proof of this statement is trivial if for each t ∈ [0, 1] the vector field v is nondegenerate. (Here, we consider a vector field to be non-degenerate if, for every point x such
that v(x) = 0, the determinant of the jacobian |∂j vi | does not vanish.) In this case, since
the zeros of the vector field are always isolated, each zero at t = 0 will follow a smooth
curve as we increase t from 0 to 1, by the implicit function theorem.
In the general case, however, zeros can collide with each other and disappear, as we
change t, so the proof is not that simple. Let us assume now that we have a general family
of vector fields, possibly degenerate.
We can always approximate vt by a polynomial vector field, in the sense that for every
> 0 there is a vector field p(x) = (p1 (x), . . . , pn (x)), with pi polynomials of x, such that
vt − p < . Here the norm is defined by v2 = Dn |v(x)|2 dx, with |v(x)|2 = i (vi (x))2 .
47
If vt is a polynomial, it will have a finite number of sets of connected zeros. We can
define an index for every isolated set of zeros. For an isolated zero x0 , we define its index
[GP] by the degree of the map φ : S → S n−1 defined by
φ(x) =
vt (x)
,
|vt (x)|
where S is a sphere center at x0 with radius , and is sufficiently small so that x0 is the
only zero of vt in its interior. So the index is just a measure of how many times we map
the sphere S into the sphere S n−1 . In the two dimensional case, the index is simply the
number of times the vector field rotates (clockwise) completely as we move around the circle
clockwise (counterclockwise rotations count as −1). If, for instance, the image of the map
is not the whole sphere S n−1 , we say that the index is zero. Figure 4.1 shows examples of
zeros with indices +1 and −1.
Index = +1
Index = −1
Figure 4.1: Two-dimensional vector fields with zeros of index +1 and −1.
For any connected set of zeros, we can define the index analogously, by isolating the set
from other zeros by a surface that can be smoothly deformed into a sphere. Alternatively,
we could define the index by adding a small perturbation to remove the degeneracy, and
adding the index of all resulting isolated zeros. A general definition of the index for vector
fields can be found in [GS].
Since vt always points inwards in the boundary of the disk, the sum of the indices of all
connected sets of zeros is equal to one, for all t. As t varies from 0 to 1, the zeros describe a
continuous trajectory. A given zero of index +1 cannot simply disappear, unless it collides
48
with a zero of index −1. Every connected set of zeros can at some point break up into
other sets, provided the sum of indices is preserved. Also, new sets of zeros can be created,
provided their total index is 0.
Let Z ⊂ D n × [0, 1] be the set of all zeros of v, where we now drop the index t, and
consider the vector field defined on the cylinder Dn × [0, 1]. If a zero of v0 with index +1 is
not connected to a zero of v1 by Z, then it has to be connected to another zero of v0 with
index −1. Since the total index is always equal to +1, there are not enough zeros of index
−1 to annihilate all the +1 zeros, and there will be at least one connected set of zeros from
t = 0 to t = 1. That concludes the proof for a polynomial vector field.
+1
+1
E
−1
+1
C
B
+1
A
−1
D
F
t
Figure 4.2: Sketch of the curve of zeros connecting the basis of the cylinder.
Figure 4.2 sketches one possible curve of zeros connecting the two basis of the cylinder,
D0 and D1 . In D0 we have three zeros, and the total index is +1 + 1 − 1 = +1. As we move
forward in t, these zeros describe continuous trajectories, until we reach A, where two zeros
of opposite index annihilate each other. In B we have a creation of a pair +1, −1 of zeros.
The −1 zero eventually collides with a +1 zero at C, but the remaining +1 zero goes all
the way until D1 . Note that we can go from D0 to B to C to D1 along a curve of zeros of
the vector field. In D we have a creation of another pair, which is eventually annihilated in
F , whereas the pair created in E reaches D1 .
49
We can find a sequence of vector fields {wn }, such that v + wn is a polynomial vector
field and limn→∞ wn = 0. For each n, we have a continuous curve of zeros of v + wn ,
from t = 0 to t = 1. What we need to do is to prove that the curve cannot be disrupted as
we take the limit wn → 0. Since the vector field is continuous and defined on a compact
domain, we can also assume that (v + wn )(x) converges uniformly to v(x).
Let us also define the set
Z = D0 ∪ Z ∪ D1
where D0 is the disk at t = 0, and D1 the disk at t = 1. Consider the topology τ , where
the open sets are defined by the intersection of the open sets of the usual metric topology
with D × [0, 1].
Let us assume for a moment that there is no subset of Z connecting D0 to D1 . In this
case, we can find two disjoint open sets O1 , O2 ∈ τ such that
D0 ⊂ O0
,
D1 ⊂ O1 .
But the set (O1 ∩ O2 )c is nonempty and closed, so
|v(x)| = 0
,
∀x ∈ (O1 ∩ O2 )c ,
|v(x)| > δ
,
∀x ∈ (O1 ∩ O2 )c .
implies
Hence, for n sufficiently large,
|(v + wn )(x)| > δ/2,
and there is no curve of zeros of v + wn connecting D0 to D1 . This is in contradiction with
the existence of a sequence of polynomial vector fields converging to v.
Therefore, we conclude the proof of the existence of a connected set of solutions of the
Bethe-Ansatz equations, extending from any U > 0 to arbitrarily large values. Along this
curve, the energy given by
E(U (t)) = −2
N
j=1
is continuous in t.
cos kj (t)
50
4.4
Non-vanishing norm of the wavefunction
As we discussed before, the existence of a continuous curve of solutions implies that
the right side of (4.1) is indeed the true ground state, provided that its norm is not zero.
However, since we are working on a lattice, and not in the continuum, the wavefunction
(4.1) could in principle vanish (f (x1 , . . . , xN ) ≡ 0) , even though the coefficients AQ {P }
are not all zero. Therefore, to complete the proof that (4.1) is the true ground state, we
need to prove that f (x1 , . . . , xN ) does not vanish identically. We should also point out that
the norm of the Bethe Ansatz wavefuncion was also studied in [GK], where a determinant
formula for the norm is conjectured. This conjecture is motivated by the existence of a
similar determinant formula for a large class of models that can be solved with the Bethe
Ansatz.
We will follow here the same strategy used by C. N. Yang and C. P. Yang [YY] to prove
that the Bethe Ansatz does indeed give the true ground state for the anisotropic Heisenberg
model.
Let us first define the variables zj by
zj =
⎧
⎨ eikj
⎩ Λ
j−N
,
j = 1, . . . , N,
,
j = N + 1, . . . , N + M.
We can write the equations (4.3), (4.4) as polynomial equations in the variables zj and
U . So we have a system of N + M polynomial equations, in N + M + 1 variables,
pi (zj , U ) = 0
,
i = 1, . . . , N + M,
(4.9)
and our goal is to find solutions given by algebraic functions zj (U ).
Let us remind ourselves of some basic facts concerning polynomial equations in several
variables. A more detailed background on the subject can be found in [CLO]. Consider the
space of polynomials in n complex variables, with complex coefficients, which we denote by
C[x1 , . . . , xn ]. Given a set of polynomials {fi } ⊂ C[x1 , . . . , xn ], where i = 1, . . . , m, we can
51
define an ideal I(f1 , . . . , fm ) ⊂ C[x1 , . . . , xn ] by
m
hi fi , h1 , . . . , hm ∈ C[x1 , . . . , xn ]}.
I(f1 , . . . , fm ) = {
i=1
We can also define an affine variety V (I) ⊂ Cn by
V (I) = {(a1 , . . . , an ) ∈ Cn : f (a1 , . . . , an ) = 0 , ∀f ∈ I}.
V (I) is the set of solutions to the original polynomial equations
f1 (a1 , . . . , an ) = f2 (a1 , . . . , an ) = · · · = fm (a1 , . . . , an ) = 0.
(4.10)
We say that the set of polynomials {fj } form a basis to the ideal I. But an ideal can be
defined by many bases. A particularly convenient basis for an ideal is the Groebner basis
(see appendix A or [CLO] for details). The Groebner basis is a set of polynomials with the
same roots as the original set, with the nice property that it contains polynomials where
some of the variables are eliminated. We will be using a Groebner basis for the ideal defined
by the Bethe Ansatz equations shortly, but we still need to define the projection of an affine
variety.
Given V (I), we can define π (k) : V (I) ⊂ Cn → C by
π (k) (a1 , . . . , an ) = ak
The projection π (k) (V (I)) is the set of values of xk for which the system of polynomial
equations has at least one solution.
We are ready to state the technical result that will assist us in our proof.
Lemma 4.4.1. Let V ⊂ Cn be an affine variety. Then one of the following has to be true:
• π (k) (V ) consists of finitely many points (possibly zero),
• C − π (k) (V ) consists of finitely many points (possibly zero).
Proof. If we define the ideal I (k) = I ∩ C[xk ], given by the polynomials in I that depend
only on the variable xk , we have
V (I (k) ) − W ⊂ π (k) (V ) ⊂ V (I (k) ),
52
where W is a proper subset of V (I (k) ). The right side of this relation is trivial, since
(a1 , . . . , ak , . . . , an ) ∈ V (I) implies that ak is a root of any polynomial in I (k) . The left side
follows from [CLO] (Theorem 3, p.123).
But V (I (k) ) and W are affine varieties in C, which can either be C or a set with finitely
many points (roots of a polynomial equation p(xk ) = 0). If V (I (k) ) = C, W is a finite set
and C − π (k) (V ) consists of finitely many points. Otherwise, π (k) (V ) is finite.
As a first application of this lemma, we see that our system of polynomial equations will
actually have a solution for any complex U , except for at most finitely many values, since
we proved that it has at least one solution for U real and positive. If we choose the ordering
z1 > z2 > · · · > zN +M > U , our Groebner basis will not have any polynomial depending
only on U . Therefore, our last polynomial in the basis should be of the form
p(zM +N , U ) = 0.
There is nothing special about the ordering we chose, so can also construct a Groebner
basis with the ordering z1 > · · · > zk > U , for any k between 1 and N + M . Therefore, the
solutions of (4.9) will have to satisfy
p(1) (z1 , U ) = 0,
p(2) (z2 , U ) = 0,
..
.
p(M +N ) (zM +N , U ) = 0.
Let us assume for a moment that none of these polynomials vanishes identically. Then,
we can factorize each of them into irreducible polynomials
(k)
p(k) (zk , U ) = Πlj=1 pl (zk , U ).
(k)
Every solution will be given by roots of some combination of pl . Since there are finitely
many combinations, and infinite solutions, there will be one particular combination that
will give us solutions for all but finitely many values of U (we are again using Lemma 4.1).
53
If one of these polynomials vanishes identically, it means that our system of equations
is degenerate, and we can add equations to it, reducing the set of solutions, in such a way
that the previous argument will work for the extended set of polynomial equations.
Therefore, each of our variables zk will be an algebraic function of U defined by
(k)
pl (zk , U ) = 0.
(4.11)
Since we know that our wavefunction is a rational function in zk and U , if we go from
U = ∞ to any finite U > 0 along the set of solutions defined by (4.11), our wavefunction
will have finitely many zeros and poles. Since those are isolated, we can redefine the
wavefunction at these points by the limit of the normalized wavefunctions as we approach
them, since the energy of those states varies continuously and is equal to the ground state
energy for all U . Again, this idea goes back to [YY], in the context of a spin system. We
can summarize the previous results in the following theorem:
Theorem 4.4.2. There is a solution to the Bethe Ansatz equations in which all zi are
algebraic functions of U .
If N/Na ≤ 2/3, we could actually prove a slightly stronger version of this result, in
which we show that the wavefunction that results from the solution of (4.5) and (4.6) is not
zero for all U > 0.
Theorem 4.4.3. If the density of electrons per site satisfies
2
N
≤ ,
Na
3
the norm of the state given by the Bethe-Ansatz method is strictly positive:
f (x1 , . . . , xN )|x1 , . . . , xN > > 0.
1≤xi ≤Na
Therefore, the Bethe-Ansatz gives us the true and unique ground state of the system.
Proof. Let us start by considering the case M = M = N/2. Since we know that for a given
solution of (4.5) and (4.6), the coefficients AQ {P } are not all equal to zero, we just need to
54
show that f (x1 , . . . , xN ) ≡ 0 in the region R implies
AQ {P } = 0,
for any P and Q.
Let us consider the region R0 defined by
R0 = {1 ≤ x1 ≤ x2 ≤ · · · ≤ xN ≤ Na } ⊂ R.
Gaudin [G] showed that f (x1 , . . . , xN ) in R is totally determined by its value in the
region R0 , and A0 {P } = 0, for any P, implies AQ {P } = 0, for all P, Q. Therefore we just
need to show that f (x1 , . . . , xN ) ≡ 0 implies A0 {P } = 0, for any P. We will denote the
restriction of f (x1 , . . . , xN ) to R0 by
f0 (x1 , . . . , xN ) =
N
A0 {P } exp i
kPj xj .
P
j=1
Also in [G], it is shown that for any two permutations P and Q, such that Q only acts
on variables corresponding to one kind of spin, we have
A0 {P Q} = I(Q)A0 {P },
where I(Q) is the sign of the permutation Q. This can be proved by combining the boundary
conditions at the different regions RQ with the antisymmetry of the wavefunction with
respect to exchange of identical particles. In particular, if P and P differ only by the
exchange of coordinates of two identical particles,
A0 {P } = −A0 {P }
Therefore, we have
f0 (x2 , x1 , . . . , xN ) =
−
A0 {P } exp i(kP 1 x2 + kP 2 x1 + · · · + kP N xN ) =
P
A0 {P } exp i(kP 1 x1 + kP 2 x2 + · · · + kP N xN ) = −f0 (x1 , x2 , . . . , xN ).
P
Notice that this is not a trivial consequence of the antisymmetry alone, since f0 is not
the true wavefunction when x2 < x1 .
55
As a result, we have that f (x1 , . . . , xN ) ≡ 0 implies f0 (x1 , . . . , xN ) = 0 in the region
{1 ≤ x1 , . . . , xM ≤ xM +1 , . . . , xN ≤ Na } ⊂ R
If Na ≥ 3N/2, that includes
S = {1 ≤ x1 , . . . , xM ≤ N , N + 1 ≤ xM +1 , . . . , xN ≤ 3N/2} ⊂ R
If f0 is zero in S, we can fix x2 , . . . , xN and consider the system of equations
f0 (x1 , x2 , . . . , xN ) = 0 ,
1 ≤ x1 ≤ N,
or
eik1 x1 C (1) (x2 , . . . , xN ) + · · · + eikN x1 C (N ) (x2 , . . . , xN ) = 0
,
1 ≤ x1 ≤ N,
(4.12)
where each C (n) is given by a sum over permutations such that P (1) = n:
C (n) =
A0 {P } exp i(kP 2 x2 + · · · + kP N xN ) .
P :P 1=n
Therefore, (4.12) is a homogeneous system of N linear equations in the variables C (N ) .
The determinant
⎛
ik1
⎜e
⎜
⎜ e2ik1
⎜
⎜ .
⎜ ..
⎜
⎝
eN ik1
eik2
···
e2ik2
..
.
···
..
.
eN ik2
···
⎞
eikN ⎟
⎟
e2ikN ⎟
⎟ iki
(e − eikj ) = 0
=
.. ⎟
⎟
. ⎟ j<i
⎠
eN ikN
does not vanish, since −π < k1 < · · · < kN < π, and we have
C (n) = 0 ,
1 ≤ n ≤ N.
We can fix x3 , . . . , xN , and write
C (n) =
eikm x2 C (n,m) (x3 , . . . , xN ) = 0 ,
1 ≤ x2 ≤ N − 1,
m=n
where
C (n,m) =
P :P 1=n,P 2=m
A0 {P } exp i(kP 3 x3 + · · · + kP N xN ) .
56
Again, the determinant does not vanish and
C (n,m) = 0.
After repeating this argument for all spin-down variables, we will have only M = N/2
pseudo-momenta ki in each term. That is why it suffices to have N + 1 ≤ xM +1 , . . . , xN ≤
3N/2. By the time we get to the variable xN , we will have
A0 {P } = 0
for all P .
If the set {ki } is given by a solution of the equations (4.5) and (4.6), they correspond
to a non trivial set of coefficients A0 {P }, and by the result above, f (x1 , . . . , xN ) = 0 for
some (x1 , . . . , xN ) ∈ R, and the state given by the right side of (4.1) has a strictly positive
norm. That concludes the proof of the theorem.
In general, M = M , and the argument above still applies provided Na ≥ N + M (which
includes Na ≥ 3N/2 as a particular case).
4.5
Extensions to excited states
Our method to show existence of solution to the Bethe Ansatz equations can be extended
to other excited states with k-s and Λ-s real. In our proof so far, we have considered a very
particular set of coefficients Ij and Jα given by (4.7), which guarantees that our solution
connects to the ground state at U = ∞. To obtain the excited states, we just need to change
these coefficients, and check if the proof still works. We will also consider the general case
where M and N can be odd or even.
It is clear that it is not going to work for any choice of Ij and Jα , as we have a few
constraints to satisfy. First, we would like a strictly ordered solution, so we will restrict to
a set of coefficients that is also strictly ordered. We would also like to have all k-s within
the interval [−π, π]. If we go back to the Lieb-Wu equations (4.5), we have
2πIj − M π < Na kj < 2πIj + M π.
57
To make sure that −π ≤ kj ≤ π we need
−
Na − M
Na − M
≤ Ij ≤
,
2
2
for all j = 1, . . . , N .
The fixed point theorem used requires a bounded domain, so the Λ-s have to be bounded.
Again, using the Lieb-Wu equations (4.6) we obtain
2πJα − M π <
θ(2 sin kj − 2Λα ) < 2πJα + M π
j
If all Jα are within the interval
−
N −M
N −M
+ < Jα <
− ,
2
2
for some > 0, the Λ-s will be bounded by
|Λα | ≤
π
π U
∗
tan
−
+ 1 ≡ CN,U
4
2 2N
for α = 1, . . . , M .
Let us know state the generalized result of existence of solution to the Lieb- Wu equations.
Theorem 4.5.1. Let M , N and Na be integers with Na ≥ N ≥ 2M . If {Ij }N
j=1 and
{Jα }M
α=1 are sequences of real numbers satisfying
Na − M
Na − M
≤ I1 < · · · < IM ≤
,
2
2
N −M
N −M
< J1 < · · · < JM <
,
−
2
2
−
then the Lieb-Wu equations (4.5), (4.6) have a solution that is real and ordered in j and α.
Proof. Since J1 > −(N − M )/2 and JM < (N − M )/2, we can take
= min{J1 −
and define
∗
=
CN,U
N −M N −M
,
− JM } > 0
2
2
π
π U
tan
−
+ 1.
4
2 2N
58
Our set of k-s and Λ-s are now represented by the vectors
k = {kj } ,
Λ = {Λα }
j = 1, . . . , N
,
α = 1, . . . , M
,
.
Analogously to what we did in the beginning of the chapter, we can define a map
φ : RN +M → RN +M by
φ(k, Λ) = (k , Λ ).
where Λ is given by
θ(2 sin kj − 2Λα ) = 2πJα −
j
θ(Λα − Λβ ),
(4.13)
β
for 1 ≤ α ≤ M , and k by
k1 =
1
(2πI1 +
θ(2 sin k N +1 − 2Λβ )),
2
Na
β
whereas for 2 ≤ i ≤ N
,
ki = max{ki−1
1
(2πIi +
θ(2 sin ki − 2Λβ ))}.
Na
β
We are now going to show that if Ω ⊂
RN +M
is given by
−π ≤ k N +1 ≤ . . . ≤ kN ≤ π,
2
∗
∗
≤ Λ1 ≤ . . . ≤ ΛM ≤ CN,U
,
−CN,U
we have
φ(Ω) ⊂ Ω.
(4.14)
By definition,
∗
|Λα | < CN,U
,
|kj | ≤ π,
so we just need to show that the map does not interfere with the ordering. For the k-s, this
is trivial, since the k -s are defined to be ordered. In the case of the Λ -s, we notice that
59
the right side of (4.13) is monotone increasing as a function of Λα , and also as a function
of Jα . Since the left side is also increasing in Λα , the ordering is preserved. That concludes
the proof of (4.14).
By Brouwer’s fixed point theorem, the map φ has a fixed point in Ω. To show that this
fixed point is a solution of the Lieb-Wu equations, we just need to make sure that it does
not occur at a boundary ki = ki+1 , which follows from the exact same argument used in for
the ground state. Therefore we conclude the proof of the theorem.
We showed the existence of ordered solutions to the Lieb-Wu equations for a very general
set of coefficients Ij and Jα . Not any choice of ordered coefficients will give us a solution
of the original Bethe Ansatz equations though. If we look carefully at (4.3), (4.4), we see
that if M is even (odd), Ij have to be integers (half-integers), and if N − M is even (odd),
Jα have to be half-integers (integers). Those are the only cases of physical interest, and we
stated the theorem in the most general form just to avoid making the distinction between
these cases.
For any choice of Ij and Jα that differs from (4.7), the energy of the resulting eigenstate
is
E = −2
N
cos kj
i=1
which is larger than the ground state energy. When 2M < N < Na we obtain many excited
states by this procedure.
In the next section, we will discuss the thermodynamic limit of this solution. We will
focus primarily on the ground state, as there is no explicit solution for excited states in the
thermodynamic limit.
Chapter 5
The thermodynamic limit
5.1
Defining the densities
Let us now analyze the behavior of the solutions {kj , Λα } as we take the thermodynamic
limit, by taking Na → ∞, keeping the ratios n = N/Na and m = M/Na fixed.
For a chain of length Na = i, we can define
ρi (k) =
1
i|kj+1 − kj |
,
kj ≤ k < kj+1 ,
and ρi (k) = 0 for k < k1 or k ≥ kN . The above expression defines a sequence of functions
ρi ∈ L1 ([−π, π]), i ∈ N. We should keep in mind that every element in the sequence
corresponds to a different number of sites, with the conditions Na = i, n = N/i and
m = M/i, with n and m fixed for the whole sequence.
Analogously, we can define σi (Λ) ∈ L1 (R) by
σi (Λ) =
1
i|Λβ+1 − Λβ |
,
Λβ ≤ Λ < Λβ+1 ,
and σi (Λ) = 0 for Λ < Λ1 or Λ ≥ ΛM . We can now state our main result in this section,
which connects our finite lattice results with the solution of the model in the thermodynamic
limit:
Theorem 5.1.1. As we take the thermodynamic limit i → ∞, there exist subsequences of
60
61
{ρi } and {σi } such that
ρni → ρ ∈ L2 ([−π, π])
,
σni → σ ∈ L2 (R),
where the convergence is in L2 -norm, and ρ, σ are solutions of the Lieb-Wu integral equations (defined by (5.6), (5.7)). For the absolute ground state at half-filling (M = N/2,
N = Na ), all the sequence converges:
ρi → ρ ∈ L2 ([−π, π])
,
σi → σ ∈ L2 (R).
Proof. From the definition, the L1 norm of ρi and σi satisfy
ρi L1 =
1
N −1
=n−
i
i
,
σi L1 =
1
M −1
=m− .
i
i
In particular, these functions are uniformly bounded as we take the thermodynamic limit.
Let us also define
fi (k) = ik −
θ(2 sin k − 2Λβ ).
(5.1)
β
From the Lieb-Wu equations, we have fi (kj ) = 2πIj . Also,
fi (k) ≤ i +
which implies
8m 8
M ≤i 1+
,
U
U
8m ,
fi (kj+1 ) − fi (kj ) = 2π ≤ |kj+1 − kj |i 1 +
U
or
ρi (k) ≤
(1 + 8m/U )
.
2π
Since ρi (k) is bounded, ρi ∈ Lp ([−π, π]), for 1 ≤ p ≤ ∞. The sequence {ρi } is actually
uniformly bounded in Lp .
With a similar goal, we define
gi (Λ) =
j
θ(2 sin kj − 2Λ) +
θ(Λ − Λβ ),
β
such that gi (Λα ) = 2πJα . Its derivative is bounded by
gi (Λ) ≤
8 8in
1
.
4 sin kj −4Λ 2 ≤
U
U
j 1+
U
(5.2)
62
Once more, that implies
gi (Λα+1 ) − gi (Λα ) = 2π ≤ |Λα+1 − Λα |i
8n
,
U
or
σi (Λ) ≤
4n
.
πU
The first inequality on (5.2) shows that all gi are bounded by a function in L2 (R). Therefore,
{σi } is uniformly bounded in L2 .
5.2
Convergence
Since both {ρi } and {σi } are uniformly bounded, we can use Banach-Alaoglu theorem
to show that there is a subsequence converging weakly in L2 (see [LL]):
ρni ρw
,
σn i σw .
We will be working now with this subsequence, but we will keep the index as i just to
simplify the notation. Our goal is to show that this subsequence also converges pointwise.
If we take equation (4.5) with index j, subtract it from the same equation with index
j + 1, and divide the result by 2πi|kj+1 − kj |, we obtain
ρi (k) =
1 θ(2 sin kj+1 − 2Λβ ) − θ(2 sin kj − 2Λβ )
1
−
σi (Λ)|Λβ+1 − Λβ |,
2π 2π
kj+1 − kj
(5.3)
β
where kj ≤ k < kj+1 and Λβ ≤ Λ < Λβ+1 . We can fix k now and consider (5.3) as an
equation for ρi (k). Note that once we do that, the index j becomes a function or k and N
given by j(k, N ) = max{i = 1, . . . , N : ki ≤ k}.
To show the existence of a pointwise limit, we just need to show that the right side
of (5.3) converges for all k. We know that the fraction on the right side converges to the
derivative of θ, if |kj+1 − kj | → 0. Since we did not prove that this happens for every k, let
us divide the domain of ρ into Q ∈ [−π, π], defined by Q = {k ∈ [−π, π] : |kj+1 − kj | → 0},
and its complement Qc .
63
In Qc , the distance between the surrounding kj does not converge to zero. Therefore,
we can find a subsequence for which ρ(k) → 0 pointwise, for all k ∈ Qc , as i → ∞. Also,
we should note that the weak limit of ρi restricted to Qc is also zero, since the integral of
the characteristic function of any subset of Qc times ρi converges to zero.
Therefore, the subsequence {ρi } has a pointwise limit in Qc that agrees with the weak
limit. We will show that the same holds true in Q. A similar analysis can be done for
the sequence {σi }. We define the set B ∈ R by B = {Λ ∈ R : |Λα+1 − Λα | → 0}. In its
complement B c , the subsequence has a zero pointwise limit that agrees with the weak limit.
Now if k ∈ Q, the last term of (5.3) is
cos k (1)
[θ (2 sin k − 2Λβ ) + Ei (k, Λβ )]σi (Λ)|Λβ+1 − Λβ |,
π
(5.4)
β
(1)
where Ei (k, Λβ ) is the error in approximating the ratio in (5.3) by the derivative of θ.
Since θ is bounded, we have
(1)
Ei (k, Λβ )σi (Λ)|Λβ+1 − Λβ | ≤
β
C|kj+1 − kj |
i
β
= Cm|kj+1 − kj | −→ 0,
as i → ∞.
We need to show now that (5.4) converges to an integral. We first notice that
θ (2 sin k − 2Λβ )σi (Λ)|Λβ+1 − Λβ | =
β
+
θ (2 sin k − 2Λβ )σi (Λ)dΛ +
(2 sin k − 2Λβ ) − θ (2 sin k − 2Λβ )]σi (Λ)dΛ,
[θtr
(5.5)
is the truncated version of θ , defined by a piecewise constant function, whose value
where θtr
on the interval [Λβ , Λβ+1 ) is equal to θ (2 sin k − 2Λβ ). Again, since the second derivative
of θ is bounded uniformly by a function in L1 , we can write
(2 sin k − 2Λβ ) − θ (2 sin k − 2Λβ )| ≤ g(Λ)|Λβ+1 − Λβ |,
|θtr
where g(Λ) ∈ L1 (R). It follows that
lim
i→∞
(2 sin k
|θtr
− 2Λβ ) − θ (2 sin k − 2Λβ )|σi (Λ)dΛ ≤ lim
i→∞
g(Λ)
dΛ = 0.
i
64
Since {σi } has a weak limit, the first integral on the right side of (5.5) converges, and
we have
cos k
1
−
lim ρi (k) =
i→∞
2π
π
θ (2 sin k − 2Λ)σw (Λ)dΛ,
B
for all k ∈ Q. A similar analysis for σi (Λ) yields
lim σi (Λ) = −
i→∞
1
π
θ (2 sin k − 2Λ)ρw (k)dk +
B
1
2π
θ (Λ − Λ )σw (Λ)dΛ,
B
for all Λ ∈ B.
5.3
Integral equations
We just showed that the subsequences have a pointwise limit. But since {ρi }, {σi } are
uniformly bounded, dominated convergence implies the existence of a strong limit, which
equals the pointwise limit and the weak limit almost everywhere. Therefore, as we go to
the thermodynamic limit on this subsequence, our k-s and Λ-s converge to distributions
satisfying the integral equations
1
+ cos k
K(sin k − Λ)σ(Λ)dΛ,
ρ(k) =
2π
B
K(sin k − Λ)ρ(k)dk −
K 2 (Λ − Λ )σ(Λ)dΛ,
σ(Λ) =
Q
(5.6)
(5.7)
B
where we define
1
K(Λ − Λ ) = − θ (2Λ − 2Λ ),
π
∞
1 2
K(Λ − x)K(x − Λ )dx.
K (Λ − Λ ) = − θ (Λ − Λ ) =
2π
−∞
It was shown in [LW2] that it is convenient to restrict the integral in k to a domain
where the function sin k is monotone. Hence, we divide the range of integration Q into
Q− = Q ∩ [−π/2, π/2] and Q+ = Q ∩ ([−π, −π/2) ∪ (π/2, π]). The first thing to show is
that for every k = (π/2 + x) ∈ Q+ , there is a k = (π/2 − x) ∈ Q− .
Since the derivative of fi (k), as defined in (5.1) satisfies fi (k) ≥ i for k ∈ [−π/2, π/2],
the set Q− consists of a single interval [−k∗ , k∗ ] (possibly [−π/2, π/2]). Since f (k) ≥
65
f (k∗ ) + i(k − k∗ ) for any k ∈ [k∗ , π − k∗ ], and the kj -s are ordered, we conclude that the
existence of a gap in the interval [k∗ , π/2] implies that Q+ ∩ [π/2, π − k∗ ] = ∅, and for every
positive k in Q+ there is a corresponding point in Q− given by the reflection with respect
to k = π/2. The same holds for every negative k in Q+ , where now we need to take the
reflection with respect to −π/2.
Therefore, if we define the set Qr+ by the reflection of Q+ with respect to ±π/2,
Qr+ = {±(π − k) , k ∈ Q+ ∩ [π/2, π]},
and we will have Qr+ ⊂ Q− . Finally, we can avoid the integration in Q+ by noting that
K(sin k − Λ) is an even function of k − π/2, which implies by (5.6) that ρ(k) − 1/2π is odd
as a function of k − π/2. Therefore, we can rewrite part of the first integral in (5.7) as
+
Q+
Qr+
K(sin k − Λ)ρ(k)dk =
+
Q+
Qr+
2
K(sin k − Λ)
=
2π
2π
Qr+
K(sin k − Λ).
Now sin k is monotone in the domain of the integral in k, and we can make the change
in variables
t(x) =
(1 − x2 )−1/2
,
2π
f (x) = (1 − x2 )−1/2 ρ(sin−1 x),
for −1 ≤ x ≤ 1. If A, D are the images of Q− , Qr+ under sin k, the integral equations (5.6),
(5.7) can now be written
∞
K(x − x )B(x )σ(x )dx , x ∈ A,
f (x) = t(x) +
−∞
∞
K(x − x )(A(x )f (x ) + D(x )t(x ))dx
σ(x) =
−∞
∞
−
K 2 (x − x )B(x )σ(x )dx , x ∈ B,
(5.8)
(5.9)
−∞
where A(x) = χA\D , B(x) = χB , and D(x) = 2χD (χX denotes the characteristic function
of the set X).
These equations have to be solved for x in the sets defined above. However, we can let
the right side of (5.8) and (5.9) define f (x) and σ(x) for any x ∈ R, by setting t(x) ≡ 0
66
for |x| > 1. Since the integrals only depend on f (x) and σ(x) at the original intervals,
there is a one to one correspondence between the solutions to the extended equations and
the solutions of the original equations. In particular, the new equations will have a unique
solution if and only if the solution to (5.8) and (5.9) is unique.
It is convenient to write the integral equations in the operator form
f
= t + K̂ B̂σ,
σ = K̂ Âf + K̂ D̂t − K̂ 2 B̂σ,
(5.10)
(5.11)
where K̂ is the convolution with K and Â, B̂, D̂ are the multiplication by χA , χB and χD .
The proof that (5.10) and (5.11) have a unique solution for given sets Q and B was
done by Lieb and Wu in [LW2]. We did not prove here that the sets Q and B are uniquely
defined for every n = N/Na and m = M/Na , although we believe they are (and we also
believe that they are given by intervals Q = [−kmax , kmax ], B = [−Λmax , Λmax ]). However,
this is not a problem for the absolute ground state of the half filled band, in which case
an explicit solution for (5.10) and (5.11) can be found. In this case, we have n = 1 and
m = 1/2, and Lieb and Wu proved that Q = [−π, π] and B = R. Their proof assumes that
Q and B are intervals, but the generalization to other subsets is immediate. Let us restrict
now to the absolute ground state.
Since the solution is unique, and every subsequence of {ρi }, {σi } contains a further
subsequence that converges to the same ρ and σ, we can actually prove that the whole
sequence converges. Otherwise, we would be able to find a subsequence such that ρni −ρ >
or σni − σ > , which is a contradiction, and we conclude the proof of our theorem.
For the absolute ground state at half-filling, the limit distributions in the thermodynamic
limit for k-s and Λ-s are then given by the solution of (5.6), (5.7) with Q = [−π, π] and
B = R. This solution was obtained in [LW1],
ρ(k) =
σ(Λ) =
cos k ∞ cos (ω sin k)J0 (ω)
1
+
dω,
2π
π
1 + eωU/2
0
∞
J0 (ω) cos (ωΛ)
1
dω,
2π 0
cosh (ωU/4)
(5.12)
(5.13)
67
which yield the ground state energy
E(Na /2, Na /2) = −2Na
π
−π
ρ(k) cos kdk = −4Na
∞
0
J0 (ω)J1 (ω)
dω,
ω(1 + eωU/2 )
where Jn is the Bessel function of order n.
One of the main questions addressed in the original paper [LW1] is whether or not there
is a Mott transition at some critical U = Uc . We know that the system is a conductor for
U = 0 and an insulator for U = ∞ at half-filling, so the question is whether the system
is still a conductor for U small in the thermodynamic limit. By general arguments, we
know that there can be no such transition at positive temperatures for the one-dimensional
model, but that does not rule out the possibility of a transition occurring at T = 0 (ground
state).
If there is such a transition, it is characterized by a discontinuity in the chemical potential, which is defined by
µ(U ) =
dE(N/2, N/2)
.
dN
The chemical potential would not be well defined at the transition point, but convexity
of the ground state energy imply existence of the right and left derivatives µ+ and µ− . An
insulator is thus characterized by µ− < µ+ .
Since the ground state energy is smooth in N for all N < Na , µ+ = µ− and there can
be no transition below (or above) half filling. At half-filling, the situation is a bit more
complex, since the Lieb-Wu solution was obtained for N ≤ Na , and the energy for N > Na
has to be obtained using particle hole symmetry on both particles. This symmetry implies
E
Na − δ Na − δ Na + δ Na + δ ,
=E
,
+ δU,
2
2
2
2
and taking the derivative with respect to δ results in
µ+ + µ− = U.
Lieb and Wu calculated the left derivative
µ− = 2 − 4
∞
0
J1 (ω)
dω,
ω[1 + exp (ωU/2)]
68
and showed that µ− < U/2 for every positive U , which implies
µ− = µ+ .
Therefore the system is an insulator for any positive U , and there is no Mott transition.
Chapter 6
Beyond Hubbard: Falicov-Kimball
model
6.1
From asymmetric Hubbard to Falicov-Kimball
In the very beginning of our discussion on the Hubbard model, we assumed that the
hopping matrix was the same for spins up and down, which is a reasonable assumption,
given the symmetry between the two spin states in the absence of a magnetic field. However,
instead of considering a Hubbard model with spin half fermions, we can consider the same
model for two different kinds of “spinless” fermions, interacting with each other by the
Coulomb term. The hamiltonian for this model is given by (2.1), but where the indices ↑
and ↓ now denote one of the possible kinds of particles.
With this interpretation in mind, we have no reason to assume that the hopping terms
will be the same, especially since the fermions might have different masses. If we restrict
to the case were the hopping tij is constant for each kind of fermion, and Ui = U , we can
write the hamiltonian as
H=−
<xy>
c†x1 cy1 − t
<xy>
c†x2 cy2 + U
|Ω|
nx1 nx2 ,
(6.1)
x=1
where we set the hopping constant to be 1 for the fermions of the first kind and t for the
69
70
fermions of the second kind. The first sum is over all pairs of sites < xy > connected by a
bond in the lattice Ω. The Hubbard model we’ve been considering so far is just the limit
of (6.1) when t = 1.
In this chapter we will be discussing another limit of (6.1), namely when t = 0. That
corresponds to the case where one of the fermions is much more massive than the other.
In this limit, the number operator nx2 becomes a classical variable taking values 0 or 1,
depending on whether or not the site i is occupied by a heavy fermion. We will redefine
this variable as w(x), and denote by Λ the subset of Ω that is not occupied by classical
particles. In other words, w(x) = 0, if x ∈ Λ and w(x) = 1 otherwise. With this notation,
the Falicov-Kimball hamiltonian is
H=−
c†x cy + U
<xy>
|Ω|
nx w(x),
(6.2)
x=1
where we dropped the index for the light fermion. We will occasionally refer to the light
particles as electrons, and to the heavy particles as ions, in order to stress one of the
main applications of the model, which is the description of a solid composed of fixed nuclei
interacting with itinerant electrons. One might think that considering spinless electrons is
unrealistic, but we can simply add terms to (6.2) to get
H=−
σ <xy>
c†x cy
+U
|Ω|
(nx↑ + nx↓ )w(x),
x=1
This new hamiltonian is more realistic than (6.2), but it presents no additional mathematical
complexity, so we might as well restrict our analysis to (6.2).
From now on, we will denote by Ne the number of electrons, and Nc =
|Ω|
x=1 ω(x)
the
number of ions. The respective densities will be denoted ne = Ne /|Ω| and nc = Nc /|Ω|.
The Falicov-Kimball (FK) model, described by the hamiltonian above, was introduced
to investigate metal-insulator transitions in mixed valence compounds of rare earth and
transition metal oxides [FK]. These materials seem to be well described by a model of two
spinless fermions with different effective masses, one much larger than the other, since the
phase transition is supposedly due to interaction between conduction band electrons and
71
electrons that are tightly bound to the metallic ions in the periodic structure. Later, the
model was again considered to describe order in mixed valence systems and binary alloys.
A review of exact results can be found in [GM].
Fixed Λ, we determine a configuration of classical particles, and the problem reduces
to a trivial one particle problem, where we obtain the ground state by occupying the Ne
lowest eigenstates of (6.2). The problem becomes highly non-trivial when we allow the
configuration Λ to change, and try to determine which configuration has the lowest energy
ground state.
The approach we used to derive the Falicov-Kimball hamiltonian as a limiting case of a
Hubbard model was not in the original works of Falicov and Kimball, but is presented in
[KL], where the crystalline behavior of the absolute ground state of the model at half-filling
is proved. This result attracted renewed attention to the model as a paradigm to address
fundamental questions in condensed matter physics, like the formation of crystals. It turns
out this behavior is in contrast with a segregated phase away from half-filling, as will be
discussed in the next session.
In this chapter we will describe what is known about the crystalline and segregated
phases of the model, preparing the reader to our original result on the lower bound for the
segregation energy away from half- filling, which will be presented in chapter 7.
Before starting, let us just point out that the model has particle-hole symmetry with
respect to the two kinds of particles, analogously to the Hubbard model. In a bipartite
lattice, that permits us to map the repulsive case (U, ne , nc ),with U > 0 to the attractive
cases (−U, 1 − ne , nc ) or (−U, ne , 1 − nc ). In particular, the neutral attractive case (ne = nc )
corresponds to the repulsive case at half-filling (ne + nc = 1). By using both symmetries at
a time, we can also map (U, ne , nc ) into (U, 1 − ne , 1 − nc ). Therefore, we can always assume
that we have ne + nc ≤ 1.
72
6.2
Crystalline and segregated phases
Starting with the pioneering work of Kennedy and Lieb [KL], a significant progress
has been made in determining the conditions under which the ground state of the FalicovKimball model exhibits crystalline long range order. For certain densities of fermions in
the attractive case it was shown that when the number of electrons and ions are equal to
each other (neutral case), the arrangement of ions that minimizes the ground state energy
occurs at a periodic subset of the original lattice. It can be thought as the most homogeneous
configuration possible, for the given density. Using particle-hole symmetry, we can see that
the result also applies for the repulsive case, when the density of both fermions add to 1.
This crystalline behavior is in contrast to the segregated phase that occurs in the repulsive case when ne + nc < 1, for U sufficiently large. In this case, in order to minimize the
energy, the classical particles tend to group together, and the lattice splits in two compact
regions, one occupied by the classical particles, and its complement Λ where the electron
wavefunctions have its highest absolute values.
This section is devoted to discussing rigorous results on these phases for the ground
state in all dimensions, and also extensions to finite temperatures. Results specific to lower
dimensions will be discussed later.
6.2.1
Crystalline phase
For a bipartite lattice Ω = A ∪ B, T. Kennedy and E. H. Lieb [KL] proved the following
result:
Theorem 6.2.1. If U > 0, the restriction of H to the subspace Ne + Nc ≥ |Ω| has two
ground states:
1. Ne = |B|, Nc = |A| and Λ = B
2. Ne = |A|, Nc = |B| and Λ = A
If U < 0 and Ne + Nc ≤ 2|A| there is a ground state given by Ne = |A|, Nc = |A| and
Λ = B. This ground state is unique unless |A| = |B|, in which case it is doubly degenerate.
73
Restricting to the hypercubic lattice, we see that the ions occupy a sublattice that is
√
itself a hypercubic lattice of spacing 2. The theorem is more general, and applies also
when the sublattices have different number of sites, provided the lattice is bipartite.
Figure 6.1: Ground state for ne = nc = 1/2.
For d ≥ 2, the result was also extended for finite temperatures, and they could prove
that for sufficiently low temperatures the system still displays long range order, whereas for
sufficiently high temperatures no long range order persists. Here, low and high temperatures
are defined by β > β1 (U ) and β < β2 (U ), where β1 (U ) and β2 (U ) are curves that limit the
crystalline phase and the high temperature phase in the (β, U ) plane. Since in general these
curves are distinct, there is an intermediate region where there is no rigorous result in the
literature. An exception should be made to the limit d → ∞, where it was shown that the
two curves coincide [BM1], [BM2], [VV]. At higher temperatures, the ordering disappears,
which implies the existence of a phase transition at positive temperatures, for d ≥ 2.
A natural question is what happens in the case where nc = ne = 1/2 in the hypercubic
lattice. The same question can be rephrased in the context of the repulsive case, when ne +
nc = 1. Is there still crystallization, with the ions arranging themselves in a homogeneous
configuration?
Working in the grand canonical ensemble for the hypercubic lattice, the average densities
are equal to half when the chemical potentials satisfy µe = µc = U/2, and the Kennedy-Lieb
theorem can be applied to prove crystallization. An important generalization was obtained
by Lebowitz and Macris [LM], when they proved that the long range order actually occurs
74
in a strip in the (µe , µc ) plane surrounding the symmetry point (U/2, U/2). Therefore long
range order is preserved when ne and nc are both close to half, in the thermodynamic limit.
The result holds for U sufficiently large, but not too large (since the strip width goes as
U −1 ), and for sufficiently low temperatures.
6.2.2
The segregated phase
As we distance ourselves from half-filling, the behavior of the system turns out to be
quite distinct from the crystalline long range order phase. When nc + ne = 1 and U is
sufficiently large, numerical evidence led Freericks and Falicov [FF] to conjecture that the
absolute ground state would occur when the ions occupy a single compact domain on the
lattice, characterizing what is called the segregated case. This conjecture is intuitive, since
for large U , the ground state energy is equal to the kinetic energy to first order in 1/U .
To minimize the kinetic energy of the electrons, we should put them in a large box with a
small boundary (|∂Λ| ≈ |Λ|(1−1/d) ). This is in contrast with the crystalline phase showed
in 6.2.1, where the boundary of Λ has size |∂Λ| = |Λ|.
The proof of this conjecture for arbitrary dimensions was finally given in [FLU]. Freericks, Lieb and Ueltschi showed that the ground state energy is bounded above and below
by a bulk term plus a boundary term, and hence minimizing the energy is to some extent
equivalent to minimizing the boundary. Their result led to the conjecture of existence of a
phase transition as the chemical potential is varied. The calculation of a lower bound for the
boundary coefficient is presented in chapter 7, and can also be found in [Go2]. We will also
defer until chapter 7 to present the details of the solution in [FLU], since the calculations
in [Go2] are linked to their solution.
We strongly believe in the importance of the Falicov-Kimball model as a source of new
ideas for the analysis of the more complicated Hubbard model. Keeping this in mind, we
should point out that the connection between the two models is given by (6.1), increasing t
from 0 to 1. Segregation in the asymmetric Hubbard models is particularly important, since
in the context of the original Hubbard model, it corresponds to formation of regions that
75
are fully magnetized. Ueltschi [U] was able to extend segregation results for the hamiltonian
(6.1), provided t is sufficiently small.
The results presented so far are valid in any dimension. Let us now describe the progress
that has been made in one and two dimensions. We are not aware of any such results that
are only valid in three dimensions.
6.3
The one-dimensional model
The hamiltonian in one dimension can be written as a function of the configuration of
classical particles w by
H(w) = −
N
c†i+1 ci + U
i=1
N
w(x)c†i ci ,
i=1
where N is the number of sites. We will consider from now on periodic boundary conditions
(N + 1 ≡ 1).
Working with the attractive case, Lemberger [Lem] proved that for |U | sufficiently large,
the ground state is indeed the most homogeneous one. In other words, fixed N , nc and ne ,
there exists Uchom , such that for |U | > Uchom the minimum of the ground state energy
E(Ne , w) occurs at w = whom , where whom is the most homogeneous configuration of
classical particles with density nc . The notion of the most homogeneous configuration will
be made more precise in the following section, when we discuss this result in more detail.
Lemberger also obtained a sufficient condition for segregation in the non-neutral case:
given Λ, nc and ne = nc there exists Ucseg such that the ground state is the segregated
configuration wseg . This configuration consist of a single cluster of classical particles, which
is clearly unique up to an overall translation.
Let us describe in more detail each of these results individually.
6.3.1
Neutral case
In order to state Lemberger’s result more precisely, we need first to define what is meant
by “the most homogeneous configuration”. It is intuitive to think of this configuration as
76
the one in which the ions are nearly equally spaced. So if nc = 1/d, with d integer, it is
easy to construct the whom , just by putting each ion at a distant d from the previous ion
on the sequence. The configurations in which the ions distant from each other by a fixed
number d are to be defined as regular configurations.
However, in general we cannot construct a regular configuration for a given nc = p/q.
The best we can do then is to force the spacing between occupied sites to be either d or
d + 1, where
1
1
≤ nc ≤ .
d+1
d
Therefore, given a configuration w of Nc occupied sites, we say that the configuration is
homogeneous if the distance between two consecutive occupied sites xj and xj + 1, given
by dj = |xj+1 − xj | satisfies
dj ∈ {d, d + 1}
,
j = 1, . . . , Nc
Now the problem is that given Nc there may be more than one homogeneous configuration.
We need to somehow decide which of those is the most homogeneous. The intuitive picture
is to determine in which configuration the sequence of dj alternates between d and d + 1 in
the most regular way. To formalize this concept, it is convenient to introduce the definition
of derivative of a configuration. From now on, we will also denote a configuration w by
N
the sequence {w(j)}N
j=1 ∈ {0, 1} , where w(j) = 1 if j is an occupied site and w(j) = 0
otherwise.
Given a homogeneous, non-regular configuration w, we define its derivative w(1) ∈
{0, 1}Nc to be a new configuration with Nc sites, defined by w (j) = 1 if dj = d and
w(1) (j) = 0 if dj = d + 1, for j = 1, . . . , Nc .
Given Nc and N , we define whom as the homogeneous configuration for which the n-th
derivative w(n) is well defined and regular for some n. Notice that being well defined means
that all derivatives up to order n − 1 are homogeneous. Let us prove the existence of such
configuration by a constructive argument.
For each Nc and N , we can always construct a homogeneous configuration. If this
configuration is regular, we can take it to be the whom . If not, we take the derivative, and
77
we will get a sequence with Nc (w ) occupied sites, with Nc (w(1) ) < Nc (w). For the new
Nc (w(1) ) and N (w(1) ) we can also find a homogeneous configuration. If we continue with
this procedure, we will get a decreasing sequence Nc (w), Nc (w(1) ), Nc (w(2) ), . . . , which
means that for some integer n the only possible homogeneous configuration w(n) will be
regular. We now define whom to be the configuration obtained by integrating w(n) n times,
and reconstructing the sequences Nc (w(n) ) and N (w(n) ) up to Nc and N . Clearly there
exists a unique whom for given Nc and N , up to a translation.
Let us consider one example, just to clarify this idea, by taking Nc (w) = 7 and N (w) =
11. In this case, d = 1, Nc (w(1) ) = 3 and N (w(1) ) = 7. Proceeding with the argument,
we have d(1) = 2, Nc (w(2) ) = 2 and N (w(2) ) = 3. Finally, d(2) = 1, Nc (w(3) ) = 1 and
N (w(3) ) = 2. Clearly, w(3) is regular, and we can obtain whom by integrating it three times.
The following figure illustrates the procedure.
w(3)
w(2)
w(1)
whom
Figure 6.2: Construction of whom for Nc = 7 and N = 11.
After proving the existence and uniqueness of the most homogeneous configuration,
Lemberger showed that for every configuration w, there is a sequence {wj }kj=1 interpolating
between w1 = w to wk = whom , such that E(Ne , wj ) > E(Ne , wj+1 ), for U sufficiently
large.
Theorem 6.3.1 (Lemberger I). For Λ finite, and given Nc , the most homogeneous configuration whom (Nc ) is unique, up to a overall translation. If Ne = Nc , and nc = p/q, with
p and q relatively primes, there exists C > 0, independent of Λ, p, q and Ne such that for
|U | > C4q , the ground state energy E(Ne , w) takes its minimum at w = whom .
78
6.3.2
Charged case
When U = −∞ and ne < nc , the ground state will have all electrons on top of the ions.
In order to minimize their kinetic energy, the best possible configuration is wseg , in which
all ions form a single cluster of size Nc . Clearly wseg is unique up to an overall translation.
It minimizes the boundary cost to the kinetic energy of the electrons, and hence the total
energy. It is natural to expect that this remains true for large but finite |U |.
Using degenerate perturbation theory, Lemberger proved the following result:
Theorem 6.3.2 (Lemberger II). If Ne < Nc , there exists a function Ucseg (ne /nc ) such
that for U > Ucseg (ne /nc ), E(Ne , w) takes its unique minimum at w = wseg , where wseg is
defined as the configuration in which all ions form one unique cluster of size Nc .
The argument above cannot work when the densities are equal, since we showed that
in the neutral case the ground state has crystalline structure. Indeed, in the ground state
when ne = nc and U = −∞, all ions are coupled to an electron, hence the kinetic energy is
zero for any configuration Λ. The two theorems proved by Lemberger are consistent in the
limit when ne = nc , since both wseg and whom have the same ground state energy.
6.4
Perturbative results in two dimensions
For a given configuration, the problem of finding the spectrum of the FK hamiltonian
reduces to the one particle problem with hΛ = T + 2U W , where T and W are the Λ × Λ
matrices defined by Txy = −txy and Wxy = w(x)δxy (notice that we scale U by a factor
of 2). It is convenient to define the spin variables s(x) = 2w(x) − 1, in such a way that
s(x) = 1 if the site is occupied by a classical particle and s(x) = −1 otherwise. If we define
a matrix S = 2W − I, our hamiltonian becomes
hΛ = T + U S,
where we disregard the terms that do not depend on the configuration.
79
The diagonal matrix elements of hΛ are either U or −U , whereas for the off diagonal
terms we have x=y |txy | ≤ zt, where t = max{|txy |}. Therefore, for U < −zt, the number
of negative eigenvalues is given by Nc , and we can write the energy for N = Nc electrons
(neutral case) as
Eneutral
z U dz
Tr
=−
sx +
,
2
z − hΛ
γ 2πi
x∈Λ
where γ is a closed positively oriented curve in the complex plane enclosing only the negative
eigenvalues of hΛ . The recursive application of the identity
1
1
1
1
−
T
=
z − hΛ
z − U S z − U S z − hΛ
results in an expansion of the ground state energy in powers of U −1 . This expansion is the
basis of most of the perturbative results for large U in the FK model. It can be proved
that the series converge only when γ encloses all negative eigenvalues, so the method only
applies to the neutral case (or half filled repulsive case).
Gruber, Jedrzehewski and Lemberger [GJL] calculated explicitly the first few terms on
this expansion for the square lattice. Up to order U −3 , they obtained
E0 (s) =
1
8
U −1 −
9 −3 3
U
sx sy +
32
32
|x−y|=1
+
√
|x−y|= 2
sx sy
1 −3 5
U
sx sy + U −3
s(P ),
16
16
|x−y|=2
P
where the last term is the sum of s(P ) = sx sy sw sz over all the unit squares P in the lattice,
x, y, z, w ∈ P .
Based on calculations derived from this expansion, they suggested that the ground state
configuration for the neutral case would be the most homogeneous one for densities n = 1/3,
1/4 and 1/5. Figure 6.4 illustrates the ground state proposed for n = 1/4 and n = 1/5.
Shortly after those ground states were suggested, Kennedy [K1] gave a rigorous proof
that they were indeed correct, for periodic boundary conditions and U sufficiently large.
Using the same truncation to order U −3 , he wrote the energy as the sum of three by three
80
Figure 6.3: Ground states for nc = 1/5 and nc = 1/4.
blocks (squares with 9 lattice sites), E0 (s) =
B
EB , where the term corresponding to each
block is given by
388EB = [8U −1 − 18U −3 ]
sx sy + 9U −3
x,y∈B:|x−y|=1
+8U −3
√
x,y∈B:|x−y|= 2
sx sy + 30U −3
sx sy
s(P ).
(6.3)
P ⊂B
x,y∈B:|x−y|=2
If a configuration can be found such that it minimizes all the EB , than it will clearly be
the ground state for E0 . In general, there is no such configuration, but the following lemma
turns out to be sufficient.
i , i = 1, 2, 3, such that
Lemma 6.4.1. There exist KB
B
i depends only on the total
KB
i is minimized by the
density of fermions (which is assumed to be constant), and EB + KB
following configurations of three by three blocks (figure 6.4), up to rotations and reflections:
• I and II for i = 1.
• I, III and VI for i = 2.
• V and VI for i = 3.
One can easily check that the homogeneous configuration suggested for density 1/5 is
given by a combination of blocks I and II. Analogously for ρi = 1/4, we have a combination
81
I
II
III
IV
V
VI
Figure 6.4: Three by three blocks.
of I, III and IV , and for ρi = 1/3 we have blocks V and VI. Therefore the lemma proves
that the crystalline configurations of figure 6.4 minimize the truncated expansion 6.3. To
prove that they minimize the energy, one needs to work on some bounds for the remaining
terms. The interested reader is referred to the original proof [K1], but the basic idea is to
show that if we make slight changes to the ground state proposed, the possible negative
contribution from the terms U −n , n > 3 can be made smaller than the positive contribution
given by the U −3 terms, by taking U sufficiently large.
In the same work, Kennedy answered negatively to the question of whether or not it is
true for every density that the ground state configuration is the most homogeneous one (in
the neutral case). His answer is based on the proof that for 1/3 ≤ nC ≤ 1/2, and under
the assumption that the ground state is periodic, then either for m = 1 or m = −1, the
line with slope m is totally occupied by nuclei or totally empty of nuclei. Therefore, if
we consider n = 1/2 − , with small, the ground state will definitely not have the most
homogeneous configuration of nuclei.
Although there is no hope of finding a result analogous to Lemberger’s in two dimensions, the half-filled repulsive case exhibits periodicity for several other densities ne =
2/5, 2/9, 2/11, 1/6 [K2], [Ha1]. Haller [Ha2] also showed periodicity for any density of the
form ne = 1/[n2 + (n + 1)2 ], n integer. Haller and Kennedy [HK] proved it for any rational number between 1/3 and 2/5, in which case the ground state configuration is given
82
by constructing Lemberger’s whom for the given density on a line, and extending it to the
second dimension by making every 45 degree line either totally occupied or empty. For
2
) ∪ ( 15 , 29 ) ∪ ( 29 , 14 ), coexistence of two periodic phases can be shown [K2], [Ha1].
ne ∈ ( 16 , 11
All these results are valid for U sufficiently large.
Chapter 7
Lower bound for the segregation
energy in the FK model
7.1
Introduction
As we described in the last chapter, the Falicov-Kimball model assumes two kinds of
fermions in the lattice Ω: classical (infinitely massive) ‘ions’ with density nc = Nc /|Ω| and
electrons with density ne = Ne /|Ω|. For simplicity, the particles are assumed to be spinless
(without loss of generality, the spin variable can be introduced later). The Falicov-Kimball
hamiltonian can be written in the second quantized form
H=−
txy c†x cy + U
x,y∈Ω
nx w(x),
x∈Ω
where c†x and cx are the fermion creation and annihilation operators for the electrons in x,
and nx = c†x cx . The variable w(x) can be either 1 or 0, according to whether the site x is
occupied by a classical particle or not. We will assume here Ω ⊂ Zd
For a bipartite lattice Ω = A ∪ B, Kennedy and Lieb [KL] proved that the ground
state displays crystalline long range order at half filling (nc + ne = 1). This result illustrates the relevance of the model in fundamental problems in condensed matter physics,
like understanding the formation of crystals and molecules. Also, it is expected that the
83
84
better understanding of the Falicov-Kimball model will provide new insights to the Hubbard model. And in the context of the Hubbard model, other fundamental questions can be
addressed, like the existence of ferromagnetism in a system in which the spins are itinerant
(not localized).
A long standing conjecture for the Falicov-Kimball model [FF] was that, for sufficiently
strong interactions, the two kinds of particles should segregate away from half-filling. This
conjecture was proved in [FLU] where it is shown that the total ground state energy is
bounded above and below by a bulk term, plus a second term which is proportional to the
boundary of the region Λ devoid of classical particles. If EΛ,N is the ground state energy
for N electrons,
e(n)|Λ| + α1 (n)|∂Λ| ≤ EΛ,N ≤ e(n)|Λ| + α2 (n)|∂Λ|,
where e(n) is the energy per site for a density n = N/|Λ| of free electrons in the infinite
lattice Zd . Therefore, given that the bulk term is fixed for all configurations, lowering the
energy requires minimizing the boundary, which is accomplished by segregating the two
species of fermions from each other.
Also in [FLU], explicit coefficients α1 (n) are obtained for low densities of electrons
n ≤ |Sd |/(4π)d , where Λ is the domain devoid of classical particles and |Sd | is the volume of
the d-dimensional sphere, whereas α2 (n) is determined for all densities. The lower bound
is obtained by considering first the U = ∞ case. Taking txy ≡ 1, the hamiltonian acting on
a function ϕ(x) ∈ L2 (Λ) can be written
[hΛ ϕ](x) = 2dϕ(x) −
ϕ(x + e),
e:x+e∈Λ
where the sum is over the edges of the lattice. The eigenvalue equations are hΛ ϕj = ej ϕj ,
for j = 1, . . . , |Λ|. Their lower bound is derived from the inequality
1
1
−
ε
)
|(bk , ϕj )|2 dk,
(ε
EΛ,N − |Λ|e(n) ≥
F
k
d
(2π)
(4d)2
j:ej >eN
where εk = 2d − 2 i cos ki . Also, the concept of the boundary vector
bk (x) = χ∂Λ (x)e−ikx
e:x+e∈Λ
/
e−ike
(7.1)
85
is introduced. Therefore, the problem reduces to showing that the boundary vector has a
projection in the subspace spanned by the largest eigenvalues. The mathematical results are
bounds for the sum of the lowest eigenvalues of the Laplace operator. For the continuous
Laplace operator, one should refer to [LY].
This will be the starting point of our study here. Our goal is to obtain an explicit
coefficient for the boundary term for intermediate densities |Sd |/(4π)d < n < 1− |Sd|/(4π)d .
We are going to obtain results for U = ∞, from which the results for finite interaction can
be derived (see section 7.5 and [FLU]). Our main result in this limit is:
Theorem 7.1.1. For d = 2 and density of electrons n = 1/2, the ground state energy of
the Falicov-Kimball model is bounded below by
EΛ,N − |Λ|e(1/2) ≥ α1 (1/2)|∂Λ|,
where α1 (1/2) = 10−13 .
7.2
The boundary term for n = 1/2
7.2.1
Projection of the boundary vector
The goal is to prove that the boundary vector bk has a projection in the subspace
spanned by the eigenfunctions {ϕj }, ej > 2d (hΛ ϕj = ej ϕj ). If we can prove that this
projection is proportional to the boundary for a subset (of non-zero measure) of the region
in k-space limited by the fermi surface of n = 1/2 (εF = 2d), the boundary term can be
calculated.
Expanding bk in terms of the eigenfunctions of hΛ we have
−
|(ϕj , bk )|2 (ej − 2d) + 2
j
j:ej >2d
=
|(ϕj , bk )|2 (ej − 2d) =
|(ϕj , bk )|2 |ej − 2d|
j
j
|(ϕj , bk )|2
(hΛ − 2d)bk 2
(ej − 2d)2
≥
.
|ej − 2d|
2d
Therefore
j:ej >2d
|(ϕj , bk )|2 (ej − 2d) ≥
(hΛ − 2d)bk 2 (bk , (hΛ − 2d)bk )
+
,
4d
2
86
and
|(ϕj , bk )|2 ≥
j:ej >2d
(hΛ − 2d)bk 2 (bk , (hΛ − 2d)bk )
≡ f (k).
+
8d2
4d
(7.2)
Suppose we can find k such that εk = 2d and (hΛ − 2d)bk 2 ≥ α|∂Λ|, for some constant
α. For such k, we have to consider the two possible cases:
• (bk , (hΛ − 2d)bk ) ≥ 0
• (bk , (hΛ − 2d)bk ) < 0
We should only be concerned with the second case, where the negative contribution from the
second term could cancel out the boundary term. We claim that for k = k + (π, π, . . . , π),
εk = 2d, (hΛ − 2d)bk = (hΛ − 2d)bk and (bk , (hΛ − 2d)bk ) = −(bk , (hΛ − 2d)bk ) ≥ 0.
|Λ|
|Λ|
Indeed, if we consider the expansions bk (x) = j=1 cj ϕj (x) and bk (x) = j=1 dj ϕj (x),
and observing that bk (x) = (−1)|x|+1 bk (x) and ϕ|Λ|−j (x) = (−1)|x| ϕj (x), we have
dj = (ϕj , bk ) =
ϕ∗j (x)(−1)|x|+1 bk (x) = −
x
ϕ∗|Λ|−j (x)bk (x) = −c|Λ|−j .
x
Therefore
(hΛ − 2d)bk 2 =
|dj |2 (ej − 2d)2 =
j
|c|Λ|−j |2 (2d − e|Λ|−j )2 = (hΛ − 2d)bk 2 ,
j
and
(bk , (hΛ − 2d)bk ) =
|dj |2 (ej − 2d) =
j
|c|Λ|−j |2 (2d − e|Λ|−j ) = −(bk , (hΛ − 2d)bk ).
j
where the identity ej + e|Λ|−j = 4d was used.
7.2.2
A bound for (hΛ − 2d)bk 2
Now, it remains to show that the first term in the right side of (7.2) can not vanish for
all k in the fermi surface. First, let us consider the case d=2. Since
[(hΛ − 2d)bk ](x) = −
e
bk (x + e) = −e−ik·x
e:x+e∈∂Λ x+e+e ∈Λ
/
e−ik·(e+e ) ,
87
the absolute value will be given by a sum of exponentials over some of the second nearest
neighbors x + e + e . The following diagrams illustrate the real part of the terms associated
with each site, for particular values of k. Note that the configuration defines which terms
will be in the sum. If we can prove, for suitable values of k, that all the terms have positive
(or negative) real part, we conclude that they are not canceled out by each other, and a
lower bound is obtained.
−1
+1 +1
−1
x
−1 +1
−1 −1 −1
−1
x
−1 −1
+1 −1
+1
k = ( π2 , π2 )
k = (0, ±π)
+1
Figure 7.1: Diagram of second nearest neighbors.
/ Λ}. If
For x ∈ Λ, let [Qx ]ij = qx,ij = #{(e, e ) : e i, e j, x + e ∈ ∂Λ, x + e + e ∈
trQx = 0
2d
2d
1 2
1 2
|[(h
−
2d)b
](x)|
≥
[(h
−
2d)b
](x)
≥ 1,
Λ
Λ
ki
ki
2d
2d
i=1
where the sum is take over ki ∈
i=1
π
π
π
{(± 2 , ± 2 ), (± 2 , ∓ π2 )}.
(7.3)
Therefore we can conclude
(hΛ − 2d)bki 2 ≥ #{x ∈ Λ, trQx = 0},
(7.4)
for ki = (± π2 , ± π2 ) or (± π2 , ∓ π2 ). The same kind of argument makes (7.4) valid for d = 3
when the sum is taken over all ki in the subset
π π π π π π π π π π π π .
± ,± ,± , ± ,± ,∓ , ± ,∓ ,± , ∓ ,± ,±
2
2
2
2
2
2
2
2
2
2
2
2
On the other hand, if trQx = 0 and Qx = 0, |(hΛ − 2d)bk (x)|2 ≥ 1 and
(hΛ − 2d)bk 2 ≥ #{x ∈ Λ : trQx = 0 and Qx = 0},
for k = (0, ±π). For d = 3, the analogous result would be
(hΛ − 2d)bk 2 ≥
1
#{x ∈ Λ : trQx = 0 and Qx = 0},
3
(7.5)
88
for ki = (0, π2 , π) or some vector obtained by permutation of the coordinates.
If there are no isolated sites in Λ,
#{x ∈ Λ : Qx = 0} = α|∂Λ| , α ≥
1
.
2d
The reason that we need not consider the case where some of the sites in Λ are isolated lies in the fact that there is always a configuration obtained by joining this site to a
larger cluster, preserving the boundary. We only need to show that the energy of the new
configuration (Λ ) is lower than the original one.
If we have a cluster and a disjoint site, the hamiltonian can be written as a direct sum
⎞
⎛
T
h1 0
⎠,
hΛ = ⎝
0 2d
where h1 is the hamiltonian for the cluster. Consider now the perturbed hamiltonian
⎞
⎛
h1 v(λ)T
⎠,
h(λ) = ⎝
v(λ)
2d
where v(λ) = (0, . . . , 0, −λ, 0, . . . , 0), such that hΛ = h(0) and hΛ = h(1). We know that
the sum of first N eigenvalues (EN ) is a concave function of the perturbation λ. Also, there
is a unitary transformation that takes λ → −λ, which implies that each eigenvalue is an
even function of λ. Combining these two results, we see that the sum of the eigenvalues is
a decreasing function of λ, and EΛ,N ≥ EΛ ,N . So, from this point on we can consider Λ as
a single cluster.
For simplicity, the remainder of the proof will be presented for d = 2, but the method is
clearly general for arbitrary d. The only difference lies in the choice of the vectors ki . We
are going to consider the two following possible cases:
• Case I: #{x ∈ Λ : trQx = 0} ≥ α2 |∂Λ|
(hΛ − 2d)bki 2 ≥ α2 |∂Λ| for ki = (± π2 , ± π2 ) or (± π2 , ∓ π2 ).
• Case II: #{x ∈ Λ : Qx = 0 and trQx = 0} ≥ α2 |∂Λ|
(hΛ − 2d)bki 2 ≥ α2 |∂Λ| for ki = (0, ±π).
89
For each of the two cases we have
f (ki ) ≥
α
|∂Λ|.
24 d2
Now we need to know how rapidly can f (k) vary.
7.2.3
A bound for |∇j f (k)|
Lemma 7.2.1. For f (k) defined by (7.2), the j-component of the gradient is bounded by
f (k) ≤ 10αd3 .
∇j
|∂Λ|
Proof. First, we should write
∂
|(hΛ − 2d)bk (x)|2 = −i
(e1 + e2 − e3 − e4 )j e−ik(e1 +e2 −e3 −e4 ) ,
∂kj
e ,e ,e ,e
1
2
3
4
/
where the sum is taken over the edges such that x+e1 , x+e2 ∈ ∂Λ and x+e1 +e3 , x+e2 +e4 ∈
Λ. We can bound the expression in parenthesis by 4, and the number of terms by (2d)4 .
Also, the number of sites where |(hΛ − 2d)bk (x)|2 does not vanish is limited by α|∂Λ|.
Therefore
∂
(hΛ − 2d)bk 2 ≤ α26 d4 |∂Λ|.
∂kj
(7.6)
The same kind of argument leads to
∂
(bk , (hΛ − 2d)bk ) ≤ 6αd3 |∂Λ|,
∂kj
(7.7)
and combining the two results we conclude the proof of the lemma.
So, for
|k − ki | ≤
1
α
1
·
=
,
10αd3 25 d2
10 · 24 d5
(7.8)
we have
α
1
f (k)
≥ 5 2 ≥ 6 3.
|∂Λ|
2 d
2 d
(7.9)
The lemma used to bound the gradient of f (k) is useful in determining a result for any
dimension. However, if we focus on determining a better coefficient for d = 2, for instance,
90
we should improve inequality (7.6). Instead of using the bound (e1 + e2 − e3 − e4 )j ≤ 4, we
can sum over all possible vectors, using the real value of the expression in parenthesis. The
same can be done for (7.7). The lower bound obtained for the j-component of the gradient
of f (k) is 7α|∂Λ|. Therefore, it turns out that (7.9) is valid for the extended region
|k − ki | ≤
α
1
1
· 5 2 =
.
7α 2 d
7 · 27
(7.10)
We should make a remark concerning the fact that we don’t know in principle which value
of ki is the right one for Case I. But since εk is invariant under inversion of coordinates,
the result will be the same, regardless of the choice between the neighborhood around
ki = (± π2 , ± π2 ) or ki = (± π2 , ∓ π2 )
7.2.4
The lower bound
We are ready now to calculate the boundary term on the lower bound for the ground
state energy. First, we should recall that
1
1
)
|(bk , ϕj )|2 dk.
(2d
−
ε
EΛ,N − |Λ|e(1/2) ≥
k
(2π)d
(4d)2
j:ej >2d
But we proved that
j:ej >2d
|(bk , ϕj )|2 ≥ f (k)|∂Λ| ≥
1
|∂Λ|
26 d3
in the neighborhood of ki . We are ready now to state the preliminary result for d = 2.
Proposition 7.2.2. Let region I be the neighborhood of ki = (π/2, π/2) defined by εk <
εF (1/2) and (7.10), and region II be the neighborhood of ki = (0, ±π) defined in a similar
way. A lower bound for the ground state energy at n = 1/2 is given by
|∂Λ|
1
·
min
(2d − εk )dk = α1 (1/2)|∂Λ|,
EΛ,N − |Λ|e(1/2) ≥
(2π)d 28 d4 I,II I,II
where α1 (1/2) > 10−17 .
A similar result holds for d = 3. The regions I and II will be defined as the vicinities of
the vectors ki presented in the last section. Also, we should include a factor of 1/3, to take
into account (7.5).
91
7.3
d=2: A better result
Considering a diagram like in figure 1 for the vector ki = (k, π − k), we see that if
Qx = 0, trQx = 0 and cos 2k > 0,
|(hΛ − 2d)bki (x)| ≥ cos 2k,
whereas if trQx = 0 and cos 2k < 0
|(hΛ − 2d)bki (x)| ≥ −4 cos 2k,
ki
where the sum is taken over ki ∈ {(k, π−k), (−k, π−k), (k, −π+k), (−k, −π+k)}. Therefore
we can extend the region of integration, as shown in the figure below. The shape of the
internal boundary curve is defined by the cos2 (2k) dependence. The new lower bound will
be given by
EΛ,N
|∂Λ|
1
− |Λ|e(1/2) ≥
·
min
(2π)d 29 d4 I,II
I,II
(2d − εk ) cos2 (2kx )dk = α1 (1/2)|∂Λ|.
Calculating the integral we get α1 (1/2) > 10−13 , which proves our main result.
I
II
I
1
0
ky
0
1
0
1
0
1
II
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0
1
0II
1
0
1
0
1
0
1
I
II
kx
I
Figure 7.2: Fermi surface εk = 4 and extended regions of integration.
92
7.4
The result for n < 1/2
For simplicity, we presented the detailed proof for n = 1/2. We should stress, however,
that the method is quite general, and can be used to obtain the lower bound for the boundary
term for any density n. Taking n = N/|Λ|, we have an inequality which is equivalent to
(7.2):
|(ϕj , bk )|2 ≥
j:ej >eN
(hΛ − eN )bk 2 (bk , (hΛ − eN )bk )
+
.
8d2
4d
Again, if (bk , (hΛ − eN )bk ) < 0, we take the application k = k + (π, π, . . . , π), and we
have εk = 4d − εk ,
(hΛ − eN )bk 2 = (hΛ − e|Λ|−N )bk 2 ,
and
(bk , (hΛ − e|Λ|−N )bk ) = −(bk , (hΛ − eN )bk ).
Suppose we can prove (hΛ − eN )bk 2 ≥ α|∂Λ| for some k in the fermi surface. Then,
either
|(ϕj , bk )|2 ≥ α |∂Λ|,
j:ej >eN
or
|(ϕj , bk )|2 ≥ α |∂Λ|,
j:ej >e|Λ|−N
with α = α/8d2 , which means that the boundary contribution can be calculated for density
n or 1 − n. But due to particle-hole symmetry, the boundary term should be the same for
the two densities. Therefore, the problem reduces to proving
(hΛ − eN )bk 2 ≥ α|∂Λ|.
This observation, when combined with the lower bound for the minimum of (hΛ −
eN )bk 2 over the fermi surface obtained in [FLU] is enough to determine the lower bound for
the boundary term. To obtain a better and explicit coefficient, however, one should proceed
like in the last sections, find vectors ki such that the function f (k) cannot vanish for all
of them in arbitrary configurations, and take the integrals over the neighborhoods of these
93
points in k-space. The difference know is that we also need to perform the integration over
neighborhoods of the vectors ki , situated in the image fermi surfaces of density 1 − n. Then
we take the minimum over all of these integrals, to determine the lower bound. Therefore,
the choice of the vectors ki might depend on the density, but the method is quite general
otherwise.
It is not hard to find a vector ki with the properties mentioned above. In two dimensions,
all we need to do is find a pair (k, k ) in the fermi surface such that n1 cos 2k + n2 cos 2k +
n3 cos (k + k ) + n4 cos (k − k ) = 0 for any choice of ni = 0, 1, 2. Since cos k is irrational for
almost every k, only a set of measure 0 in the fermi surface will violate this condition.
However, it is still useful to have a systematic way of choosing the vectors, independent
of the density. Now we will show how to get these vectors for certain ranges of densities in
two dimensions.
Fixing the density is equivalent to fixing the sum of cosines in the fermi surface. If the
√
density is low enough such that 2 ≤ cos k1 + cos k2 < 2, we can pick the vector (k, k),
where k is determined by the density. We certainly have k < π/4, so if we draw a diagram
like in figure 7.1, all the terms will have positive real part, and the sum we get a lower
bound.
On the other hand, for higher densities, when 0 ≤ cos k1 + cos k2 < 1, the natural choice
follows the exact same procedure as the one described for n = 1/2. First, we need to find
a vector that gives us a bound when trQx = 0. Inspired by the choice (0, π) when n = 1/2,
we will choose a vector of the kind (0, k). Its diagram is shown in the figure below, and it
is clear that it will give us a bound provided k > π/2.
When trQx = 0, we try ki = (π/2, π/2 − y), where now y is determined by the filling.
Figure 7.4 (a) shows the diagram for this choice of vector. To eliminate the terms cos y
with alternating sign we take the sum
[(hΛ − 2d)bki ](x),
i
94
cos 2k
cos k 1 cos k x
cos k
1
cos k
cos 2k
(0, k)
Figure 7.3: Diagram for (0, k).
in analogy to the second term in (7.3), but now summing over the vectors in the set
S=
π π
π
π π
π π
π , −y , − , −y ,
− y, ,
− y, −
.
2 2
2 2
2
2
2
2
The resulting diagram is shown if 7.4 (b), where ∆ = −2 − 2 cos y, and it is clear that we
will get a bound provided 0 ≤ x < π/2, since ∆ < 0.
− cos 2y
+ cos y −1 − cos y x
∆
− cos y
0 −1
∆
+ cos y
0 − cos 2y
x
0
∆
0
∆
(a)
(b)
Figure 7.4: Diagram for (π/2, π/2 − y) and the sum over set S.
For intermediate densities (1 ≤ cos k1 + cos k2 ≤
√
2), similar vectors can be found. It
would be good to have a convenient choice that worked for all densities. We believe that a
good choice of vectors can result in a boundary coefficient a few orders of magnitude greater
than the one obtained here.
7.5
Finite U
We have so far considered only the U = ∞ case. The results can be generalized to large
but finite U by applying a result obtained in [FLU], where it was shown that the difference
95
between the ground state energy with infinite repulsion (EΛ,N ) and the ground state energy
U ) is of the order of the boundary. More precisely, it was shown that
with finite U (EΛ,N
U
≥ EΛ,N − γ(U )|∂Λ|,
EΛ,N
where
γ(U ) =
(U − 2d)2d
8d2
+ 2d+2 d
−1 .
d
U − 2d
(U (U − 4d))
Notice that γ(U ) → 0 as U → ∞.
Combining this result with our lower bound, we get
U
≥ e(n)|Λ| + (α1 (n) − γ(U ))|∂Λ|,
EΛ,N
and segregation will still hold provided α1 (n) − γ(U ) > 0. The value of U obtained by
setting α1 (n) = γ(Uc∗ ) should be interpreted as an upper bound for the critical value of
U = Uc above which segregation occurs. We probably have Uc∗ >> Uc , since our lower
bound is far from being optimal.
Extensions to positive temperatures depend on generalizing the ground state inequalities to the free energy at low temperatures. This is also discussed in [FLU], where it is
conjectured that there is a first order phase transition when varying the chemical potential,
for d ≥ 2. We refer the interested reader to their paper.
7.6
Conclusions
We showed here how to derive the lower bound for the boundary term of the ground
state energy of the Falicov-Kimball model. The existence of the boundary term is important
since the system will try to minimize the boundary (to some extent) in order to minimize
energy. Therefore, a segregated phase, where electrons and classical particles try to occupy
distinct regions of the lattice, is obtained. When contrasted to the half-filling case, where
crystalline long range order is observed, it suggests that the model has a first order phase
transition when varying the chemical potentials.
96
Our coefficient for intermediate densities is small when compared to the upper bound
obtained in [FLU]. This means that our energy is not very sensitive with respect to the
boundary size. However, the strength of the method is that it provides an explicit coefficient.
Also, α1 (n) can indeed be much smaller than α2 (n), since the upper bound is saturated by
configurations with isolated sites, whereas the lower bound is not.
The author is indebted to E. H. Lieb and D. Ueltschi for helpful comments and suggestions.
Appendix A
Groebner bases
A monomial in n variables is a product
xα1 1 xα2 2 · · · xαnn ,
where the exponents α1 , . . . , αn are non-negative integers. If α = (α1 , . . . , αn ), we can also
use the short notation
xα = xα1 1 · · · xαnn .
A polynomial p in n variables with coefficients in a field k is a finite linear combination of
monomials,
p=
cα xα ,
α
where cα ∈ k.
It is useful to introduce the concept of the leading term of a polynomial. For that we
need a prescription to order monomials in n variables. Any total ordering relation between
n-tuples of integers can be used to define an ordering between monomials. By total ordering
relation we mean any relation > between pairs of n-tuples such that:
• For any pair α, β, with α = β, either α > β or β > α .
• α > β and β > γ ⇒ α > γ.
97
98
The corresponding ordering relation between monomials is such that xα > xβ if α > β.
Since we want the ordering to preserve the algebraic structure of polynomial rings, we will
restrict to relations that satisfy xα xγ > xβ xγ if xα > xβ .
For example, given two distinct monomials m1 = xα and m2 = xβ , we can choose an
ordering such that m1 > m2 when αk > βk , where k is the smallest integer i such that
αi = βi . We denote this ordering by x1 > x2 > · · · > xn , since the exponent of x1 is the
most relevant, followed the exponent of x2 and so on. If we want to order the polynomials
x1 x2 x23 , x1 x2 x3 , x21 x2 , x1 x22 according to this rule, we get
x21 x2 > x1 x22 > x1 x2 x23 > x1 x2 x3 .
The leading term of a polynomial f is its largest monomial (times its coefficient), according
to the specified ordering, and we denote it by LT (f ). As an example, if f = x21 x2 + 3x1 x22
we have LT (f ) = x21 x2 .
We say that the set G = {gk }lk=1 is a Groebner basis for an ideal I(f1 , . . . , fm ) if
I(LT (g1 ), . . . , LT (gn )) = I(LT (I(f1 , . . . , fm ))),
where LT (I(f1 , . . . , fm )) is the set of leading terms of all polynomials in I(f1 , . . . , fm ). In
other words, G is a Groebner basis if the leading term of any element in I(f1 , . . . , fm ) is
divisible by the leading term of at least one gk .
Every ideal I(f1 , . . . , fm ) generated by the polynomials {fk }m
k=1 = {0} has a Groebner
basis. The general idea in constructing a Groebner basis is to eliminate one by one the
variables of the polynomials fk , obtaining a set of polynomials G = {gk }lk=1 that generate
the same ideal as {fk }m
k=1 . One of its main applications is to solve systems of polynomial
equations, since reducing the number of variables simplifies the task of finding solutions.
This process can be thought of as a generalization of the process of solving a system of
linear equations by elimination of variables. We will present here the essential facts about
Groebner bases in order to understand the results of section 4. The algorithm used to
obtain a Groebner basis given a set of polynomials can be found in [CLO].
99
Given an ideal I ⊂ C[x1 , . . . , xn ], we can define the k-elimination ideal by
Ik = I ∩ C[xk+1 , . . . , xn ].
Ik is the set of all polynomials in I that do not depend on the first k variables xk . What
makes the Groebner basis special is that
Gk = G ∩ C[xk+1 , . . . , xn ]
is also a basis for Ik .
As a consequence of this property, if the variety defined by (4.10) consists of finitely
many points, Gn−1 = {0}, and we can solve the corresponding equations for xn , to find
V (In−1 ). Then, we can move on to Gn−2 , and solve it for xn−2 , finding V (In−2 ), and so
on. In the end, we will find V (I) just by solving polynomial equations in one variable. We
are ordering the variables here by x1 > x2 > · · · > xn . This ordering determines which
variables will be eliminated first. Of course we can choose any other ordering, and we will
obtain a different Groebner basis, but the set of solutions should be the same
As one example, consider the set of polynomial equations
f1 = x2 + y 2 + z − 1 = 0,
f2 = x + y + z − 1 = 0,
f3 = x + y − z 2 − 1 = 0.
Fixing the ordering x > y > z, the ideal I(f1 , f2 , f3 ) has
G = {x + y + z − 1, y 2 − y + yz − z, z 2 + z}
as a possible Groebner basis. Therefore, the solution to the original equations will satisfy
x + y + z − 1 = 0,
y 2 − y + yz − z = 0,
z 2 + z = 0,
which can be easily solved to get V (I) = {(0, 1, 0), (1, 0, 0), (1, 1, −1)}.
100
Clearly our solutions V (I) could be more complicated than isolated points. It could
contain subsets of higher dimension. In this case, we would have a slightly different basis.
Consider for instance the set of solutions of f1 = f2 = 0. It clearly has infinitely many
solutions, and
G = {x + y + z − 1, 2y 2 − 2y + 2yz + z 2 − z}
is a Groebner basis. We see that G2 = {0} and V (I2 ) = C. In this case, we could not
eliminate x and y, since we have solutions for infinitely many values of z. In summary, the
number of variables we can actually eliminate by calculating the Groebner basis depends
on the dimension of V (I).
References
[A] P. W. Anderson, The resonating valence bond state in La2 CuO4 and superconductivity,
Science 235, 1196 (1987)
[ASKA] R. Arita, Y. Suwa, K. Kuroki, and H. Aoki, Gate-induced band ferromagnetism in
an organic polymer, Phys. Rev. Lett. 88, 127202 (2002).
[Be] H. A. Bethe, Zur Theorie der Metalle: I. Eigenwerte und Eigenfunktionen der linearen
Atom Kette, Zeits. f. Physik 71, 205 (1931)
[BM1] U. Brandt and C. Mielsh, Thermodynamics and correlation functions fo the FalicovKimball model in large dimensions, Z. Phys. B 75, 365 (1989)
[BM2] U. Brandt and C. Mielsh, Thermodynamics of the Falicov-Kimball model in large
dimensions II, Z. Phys. B 79, 295 (1990)
[CLO] D. Cox, J. Little and D. O’Shea, Ideals, Varieties, and Algorithms, Springer-Verlag,
New York, USA, 1992
[D] T. C. Dorlas, Orthogonality and completeness of the Bethe Ansatz eigenstates of the
nonlinear Schroedinger model, Commun. Math. Phys. 154, 347 (1993)
[EKS] F. H. L. Essler, V. E. Korepin and K. Schoutens, Completeness of the SO(4) extended
Bethe Ansatz for the one-dimensional Hubbard model, Nucl. Phys. B 384, 431 (1992)
[FF] J. K. Freericks and L. M. Falicov, Two-state one-dimensional spinless Fermi gas,
Phys. Rev. B 41, 2163 (1990)
101
102
[FK] L. M. Falicov and J. C. Kimball, Simple model for semiconductor-metal transitions:
SmB6 and transition-metal oxides, Phys. Rev. Lett. 22, 997 (1969)
[FL] E. H. Lieb and M. Flicker, Delta function Fermi gas with two spin deviates, Phys.
Rev. 161, 179-188 (1967)
[FLU] J. K. Freericks, E. H. Lieb and D. Ueltschi, Segregation in the Falicov-Kimball Model,
Commun. Math. Phys. 227, 243–279 (2002)
[G] M. Gaudin, Travaux de Michel Gaudin: Modèles exactement résolus, Les Éditions de
Physique, Paris & Cambridge, USA, 1995
[GJL] C. Gruber, J. Jedrzejewski, P. Lemberger, Ground states of the spinless FalicovKimball model II, J. Stat. Phys. 76, 913 (1992)
[GK] F. Göhmann and V. E. Korepin, The Hubbard chain: Lieb-Wu equations and norm
of the eigenfunctions, Phys. Lett. A 263, 293 (1999)
[GM] Ch. Gruber and N. Macris, The Falicov-Kimball model: a review of exact results and
extensions, Helv. Phys. Acta 76, 913 (1996)
[Go1] P. S. Goldbaum, Existence of solutions to the Bethe Ansatz equations for the 1D
Hubbard model: finite lattice and thermodynamic limit, cond-mat/0403736, to appear
in Commun. Math. Phys..
[Go2] P. S. Goldbaum, Lower bound for the segregation energy in the Falicov-Kimball Model,
J. Phys. A 36, 2227 (2003)
[GS] L. Gua and H. Spohn, Six-vertex model, roughened surfaces, and an asymmetric spin
hamiltonian, PRL 68, 725-728 (1992)
[GP] V. Guillemin and A. Pollack, Differential Topology, Prentice-Hall Inc., Englewood
Cliffs, NJ, USA, 1974
[GS] D. H. Gottlieb and G. Samaranayake, The index of discontinuous vector fields, New
York J. Math. 1, 130 (1995)
103
[Gu] M. C. Gutzwiller, Effect of correlation on the ferromagnetism of transition metals,
Phys. Rev. Lett. 10, 159 (1963)
[H] J. Hubbard, Electron correlation in narrow energy bands, Proc. Roy. Soc. (London) A
276, 238-257 (1963)
[Ha1] K. Haller, Ground state properties of the neutral Falicov-Kimball model, PhD Thesis,
University of Arizona, 1998
[Ha2] K. Haller, Ground state properties of the neutral Falicov-Kimball model, Commun.
Math. Phys. 210, 703 (2000)
[HK] K. Haller and T. Kennedy, Periodic ground states in the neutral Falicov-Kimball model
in two dimensions, J. Stat. Phys. 102, 15 (2001)
[HL] O. J. Heilmann and E. H. Lieb, Violation of the non-crossing rule: The Hubbard
Hamiltonian for benzene, Trans. N.Y. Acad. Sci. 33, 116 (1970)
[K1] T. Kennedy, Some rigorous results on the ground states of the Falicov-Kimball Model,
Rev. Math. Phys. 6, 901 (1994)
[K2] T. Kennedy, Phase separation in the neutral Falicov-Kimball model, J. Stat. Phys.
91, 829 (1998)
[Ka] L. H. Kauffman, Knots and statistical mechanics, in The interface of knots and physics,
Proceedings of Symposia in Applied Mathematics 51, 1 (1995)
[KL] T. Kennedy and E. H. Lieb, An itinerant electron model with crystalline or magnetic
long range order, Physica A 138, 320 (1986)
[L1] E. H. Lieb, The Hubbard model: Some rigorous results and open problems, in Proceedings of the XIth International Congress of Mathematical Physics, Paris, 1994,
International Press (1995)
[L2] E. H. Lieb, Two theorems on the Hubbard model, Phys. Rev. Lett. 62, 1201 (1989)
104
[L3] E. H. Lieb, Exact solution of the problem of the entropy of two-dimensional Ice, Phys.
Rev. Lett. 18, 692 (1967)
[L4] E. H. Lieb, Residual entropy of square ice, Phys. Rev. 162, 162 (1967)
[L5] E. H. Lieb, Exact solution of the F-model of an antiferroelectric, Phys. Rev. Lett. 18,
1046 (1967)
[L6] E. H. Lieb, Exact solution of the two-dimensional KDP model of a ferroelectric, Phys.
Rev. Lett. 19, 108 (1967)
[L7] E. H. Lieb, Exact analysis of an interacting bose gas. II. The excitation spectrum,
Phys. Rev. 130, 1616 (1963)
[Lem] P. Lemberger Segregation in the Falicov-Kimball model, J. Phys. A 25, 715 (1992)
[LeM] J. Lebowitz and N. Macris, Long range order in the Falicov-Kimball model: extension
of Kennedy-Lieb Theorem, Rev. Math. Phys. 6, 927–46 (1994)
[LL] E. H. Lieb and M. Loss, Analysis, 2nd edition, Amer. Math. Soc., 2001.
[LLi] E. H. Lieb and W. Liniger, Exact analysis of an interacting bose gas. I. The general
solution and the ground state, Phys. Rev. 130, 1605 (1963).
[LM] E. H. Lieb and D. C. Mattis, Theory of ferromagnetism and the ordering of electronic
energy levels, Phys. Rev. 125, 164 (1962)
[LN] E. H. Lieb and B. Nachtergaele, The stability of the Peierls instability for ring shaped
molecules, Phys. Rev. B 51, 4777 (1995)
[LW1] E. H. Lieb and F. Y. Wu, Absence of Mott transition in an exact solution of the
short-range one-band model in one dimension, Phys. Rev. Lett. 20, 1445-1448 (1968)
[LW2] E. H. Lieb and F. Y. Wu, The one-dimensional Hubbard model: A reminiscence,
Physica A 321, 1-27 (2003)
105
[LY] P. Li and S.-T. Yau, On the Schrödinger equation and the eigenvalue problem, Commun. Math. Phys. 88, 309–18 (1983)
[Ma] D. C. Mattis, The many-body problem, World Scientific, 1993
[Mi1] A. Mielke, Ferromagnetic ground states for the Hubbard model on line graphs, J.
Phys. A 24, L73 (91)
[Mi2] A. Mielke, Ferromagnetism in the Hubbard model on line graphs and further considerations, J. Phys. A 24, 3311 (91)
[Mi3] A. Mielke, Exact ground states for the Hubbard model on the kagome lattice, J. Phys.
A 25, 4335 (1992)
[Mo] A. Montorsi, The Hubbard model: a reprint volume, World Scientific, 1992
[McG1] J. B. McGuire, Interacting fermions in one dimension, I. repulsive potential, J.
Math. Phys. 6, 432 (1965)
[McG2] J. B. McGuire, Interacting fermions in one dimension, II. attractive potential, J.
Math. Phys. 7, 123 (1966)
[MT] A. Mielke and H. Tasaki, Ferromagnetism in the Hubbard model. Examples from
models with degenerate single-electron ground states, Commun. Math. Phys. 158, 341
(1993)
[N] Y. Nagaoka, Ferromagnetism in a narrow, almost half-filled s band, Phys. Rev. 147,
392 (1966)
[R] K. Reidemeister, Knotentheorie, Chealsea Publishing Co., New York (1948)
[S] B. Sutherland, Introduction to Bethe’s ansatz, , in Exactly Solvable Models in Condensed Matter and Relativistic Field Theory, Springer-Verlag Lecture Notes in Physics,
242.
106
[T] D. J. Thouless, Exchange in solid 3 He and the Heisenberg Hamiltonian, Proc. Phys.
Soc. London 86, 893 (1965)
[Ta1] H. Tasaki, Ferromagnetism in the Hubbard models with degenerate single-electron
ground states, Phys. Rev. Lett. 69, 1608 (1992)
[Ta2] H. Tasaki, Ferromagnetism in Hubbard models, Phys. Rev. Lett. 75, 4678 (1995)
[Ta3] H. Tasaki, Ferromagnetism in the Hubbard model: a constructive approach, Commun.
Math. Phys. 242, 445 (2003)
[Ti1] G. Tian, Charged and spin-excitation gaps in half-filled strongly correlated electron
systems: A rigorous result, Phys. Rev. B 58, 7612 (1998)
[Ti2] G. Tian, Lieb’s spin-reflection-positivity method and its applications to strongly correlated electron systems, J. Stat. Phys. 116, 629 (2004)
[U] D. Ueltschi, Segregation in the asymmetric Hubbard model, J. Stat. Phys. 116, 681
(2004)
[VV] P. G. Van Dongen and D. Vollhardt, Exact mean-field hamiltonian for fermionic lattice
models in high dimensions, Phys. Rev. Lett. 65, 1663 (1990)
[W] F. Woynarovich, Excitations with complex wavenumbers in a Hubbard chain: I. States
with one pair of complex wavenumbers, J. Phys. C 15, 85-96 (1982)
[Y] C. N. Yang Some exact results for the many-body problem in one dimension with repulsive delta-function interaction, Phys. Rev. Lett. 19, 1312-1314 (1967)
[YY] C. N. Yang and C. P. Yang, One-dimensional chain of anisotropic spin-spin interactions. I. Proof of Bethe’s hypothesis for ground state in a finite system, Phys. Rev.
150, 321-327 (1966)
[YZ] C. N. Yang and S. C. Zhang, SO4 symmetry in a Hubbard model, Mod. Phys. Lett.
B4, 759 (1990)
Fly UP