Bachelor Thesis:

Comparison of
Mott-Hubbard-transition on the
Bethe-lattice and the 2D-square-lattice

Supervisors:
Associate Prof. Dr. Alessandro Toschi
Dipl.-Ing. Patrick Chalupa
Institute of Solid State Physics
Quantum many body physics group
Technische Universität Wien
by

Mathias Pelz
1429137
Abstract
Many intriguing phenomena are encountered in studies of strongly correlated
electron systems. One of the most prominent among them is arguably the
Mott-Hubbard-metal-insulator-transition (MIT), a phase transition between
a metal and an insulator, which is originated by the electronic interaction.
The presence of strong correlations poses however specific problems to the
forefront theoretical treatment. For instance a challenging aspect, that has
been recently discovered is the existence of multiple divergences of the irre-
ducible vertex functions, observed in several many-electron models. Among
them also the Hubbard model, solved by means of the dynamical mean field
theory (DMFT). For this case, depending on the lattice model studied, a
peculiar discrepancy in the low-temperature behavior of the vertex diver-
gences lines was found. To study, whether this discrepancy is rooted in the
different behavior of the MIT, the coexistence regions for the Bethe- and the
2D-square-lattice are analysed in the thesis. This comparison shows also the
influence of the shape of the density of states (DOS) on the MIT. For the nu-
merical DMFT calculations the w2dynamics package has been used, which
is an efficient implementation of a continuous-time quantum Monte Carlo
impurity solver. For both lattices the self-energy and the double occupancy
at different temperatures are calculated. From these quantities the critical
interaction values Uc1 and Uc2 as a function of temperature are obtained.

2
Kurzfassung
Viele faszinierende Phänomene treten bei Untersuchungen von stark kor-
relierten Elektronensystemen auf. Einer der bekanntesten unter ihnen ist
wohl der Mott-Hubbard-Metall-Isolator-Übergang (MIT), ein Phasenüber-
gang zwischen einem Metall und einem Isolator, der durch die elektronische
Wechselwirkung entsteht. Das Vorhandensein starker Korrelationen wirft
jedoch spezifische Probleme für die vorderste theoretische Behandlung auf.
Ein herausfordernder Aspekt, der kürzlich entdeckt wurde, ist die Existenz
multipler Divergenzen der irreduziblen Vertexfunktionen, die in mehreren
Vielelektronenmodellen beobachtet wurden. Darunter auch das Hubbard-
Modell, gelöst mit der Dynamic Mean Field Theory (DMFT). Für diesen
Fall wurde, abhängig von dem untersuchten Gittermodell, eine besondere
Diskrepanz im Tieftemperaturverhalten der Vertexdivergenzlinien gefunden.
Um zu untersuchen, ob diese Diskrepanz aus dem unterschiedlichen Verhal-
ten des MIT herrührt, werden in der Arbeit die Koexistenzregionen für das
Bethe- und das 2D-Quadratgitter analysiert. Dieser Vergleich zeigt auch
den Einfluss der Form der Zustandsdichte (DOS) auf die MIT. Für die nu-
merischen DMFT Rechnungen wurde das w2dynamics-paket verwendet, das
eine effiziente Implementierung eines zeit kontinuierlichen Quanten-Monte-
Carlo Verunreinigungslösers darstellt. Für beide Gitter wird die Selbsten-
ergie und die Doppelbesetzung bei unterschiedlichen Temperaturen berech-
net. Aus diesen Größen werden die kritischen Wechselwirkungswerte Uc1
und Uc2 als Funktion der Temperatur erhalten.

3
Contents

1 Introduction 5
1.1 Strongly correlated electron systems . . . . . . . . . . . . . . . 5
1.1.1 Fermi-liquid theory . . . . . . . . . . . . . . . . . . . . 5
1.1.2 Hubbard model . . . . . . . . . . . . . . . . . . . . . . 6
1.1.3 Dynamical mean-field theory . . . . . . . . . . . . . . . 8
1.1.4 Mott-Hubbard metal-insulator-transition . . . . . . . . 9
1.2 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.2.1 Density of states . . . . . . . . . . . . . . . . . . . . . 13

2 Methods and Calculation 18

2.1 CT-QMC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Self-energy and double occupancy . . . . . . . . . . . . . . . . 18
2.2.1 Self-enegry . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2.2 Double occupancy . . . . . . . . . . . . . . . . . . . . 23
2.3 Phase-diagram of the Bethe-lattice . . . . . . . . . . . . . . . 27
2.4 Phase-diagram of the 2D-square-lattice . . . . . . . . . . . . . 28

3 Comparison and Discussion 30

4 Conclusion 31

4
1 Introduction
1.1 Strongly correlated electron systems
Systems with strong interaction between the electrons are called strongly
correlated electron systems. The description of such systems is not possible
with models which neglect the interaction. Many materials show very inter-
esting properties, due to the correlation among the electrons. The reason for
the high interaction is the shape of the orbitals, namely a very narrow shape,
which increases the probability for two electrons to be close to each other.
So the confinement leads to a high Coulomb interaction, which cannot be
neglected. This effect is very significant for the 3d- and 4f-orbitals.

Figure 1: Shapes of the d- and f-orbitals with the elements in 3d- and 4f-
configurations, taken from [1]

1.1.1 Fermi-liquid theory

A higher Coulomb interaction affects the mobility of electrons. This can
be described by the Fermi-liquid theory, which uses quasi-particles instead
of electrons. For the change from electrons to quasi-particles the concept of
adiabatic continuity is used. The starting point is the ground state of the sys-
tem with one excited electron and no interaction. In an adiabatic manner the
interaction is turned on, thus the non interacting eigenstates are changing to

5
eigenstates of the interacting system. The quantum numbers can be assumed
to stay the same, where else the energy and the wave-function are changing.
However the consequence is a one-to-one correspondence between the free
electrons and the quasi-particles. The lifetime equals the Eq. (1), i.e. the
lifetime diverges at T = 0 at the Fermi surface and decreases quadratically for
finite temperatures per energy. In other words, this concept works only close
to the Fermi-surface and at low temperatures w.r.t. the Fermi-temperature.
1
∝ (ε~k − εF )2 (1)
τ~k

The change of the energy though the interaction is reflected in the mass of
the quasi-particles, which is called the effective mass m∗ . m∗ and also the
movement are connected to the dispersion relation ε~k given by Eq. (2).

  1 h̄~k
~v ~k = ∇~k ε~k = ∗ , |~k| ∼ kF (2)
h̄ m
Physically an electron „carries“ other electrons through the interaction. So
the electrons have a different speed and that can be described with a free
electron movement with a different mass. For a more detailed description
see [2]. In our case the self-energy Σ(iωn ) is the connection to the quasi-
particle properties (see section 2.2.1). As a result of the kinetic energy and
the interaction, strongly correlated electron systems are difficult to describe.
One model describing these systems is the Hubbard-model.

1.1.2 Hubbard model

The Hubbard model is used for the description of quantum many-body sys-
tems [3]. In the simplest realization, this model considers the electron hop-
ping between neighboring lattice sites and the Coulomb-interaction among
electrons on the same lattice site. There are four configurations on one lattice
site possible (|0i, | ↑i, | ↓i and | ↑↓i). The Hamiltonian of the Hubbard-model
is given in Eq. (3).

UX
c†iσ cjσ +
X
H = −t ni↑ ni↓ + h.c. (3)
hi,ji,σ
2 i

The first term in Eq. (3) describes the electron-hopping between lattice sites
i and j where t is the hopping integral, cjσ is the annihilation and c†iσ is the
creation operator. The kinetic energy is given by this term, which can be
diagonalized in the momentum space.

6
The second term describes the Coulomb-interaction between a spin-up and a
spin-down electron on lattice site i, which is diagonal in position space and is
defining the potential energy. The U in this term corresponds to the strength
of the screened Coulomb-interaction and the occupation number operators
niσ = c†iσ ciσ measures the number of electrons with spin σ on lattice site i.
1
The long range behavior of the Coulomb-potential r and also other long
range effects are neglected. Hence the so called Hubbard U is only local and
constant in space. The value of U can be estimated by a Coulomb-integral
in a Wannier-basis set and the strength strongly depends on the confinement
of the orbitals.

Figure 2: Schematic representation of the Hubbard-model to illustrate the

change of states over time. It shows electrons in the solid hopping between
lattice sites (hopping: t) and the double occupancy of up-electrons (↑) and
down-electrons (↓) on the same sites (interaction: U ). Picture taken from [4]

Besides U and t, the temperature, the dimension of the lattice and the filling
are parameters of the Hubbard model. Based on the filling, the density of
the electrons n = n↑ + n↓ can also be used as a parameter. From the Pauli
principle it follows that nσ ≤ 1 for σ =↑, ↓ and so n ≤ 2. A special case is the
half filled case, where n = 1.
Other special cases are U = 0 or t = 0. In the case of U = 0, the electrons
can move freely as Bloch-waves. The case t = 0 describes a system without
hopping which is called atomic limit. The atomic limit is reached for U >> t.
Exact solutions for the Hubbard-model are not available, except for the case
of one dimension and infinite dimensions. One of the most commonly used
methods to approach the model is the dynamical mean-field theory (DMFT).

7
1.1.3 Dynamical mean-field theory
The Hubbard model can be solved using the dynamical mean-field theory
(DMFT). In DMFT the lattice problem is mapped onto a local impurity
problem which is solved self-consistently. The idea behind DMFT comes
from mean-field treatments of classical problems. An example for that is
the Ising-Hamiltonian, which describes localized spins interacting with the
surrounding spins. To solve the Ising-Hamiltonian in three dimensions, the
surrounding spins are approximated with their average and described by a
global Weiss-mean-field. The problem is solved self-consistently until the
effective magnetization of the Weiss-field matches with the behavior of the
neighboring spins. This description becomes rigorously exact for infinite di-
mensions, respectively infinitely many neighbors.

Figure 3: Reduction of the Hubbard-model to a dynamical mean-field in the

limit of Z or d → ∞. Electrons can hop onto the single site and back. The
electrons can interact on the impurity site. G(ω) is the local propagator and
Σ(ω) is the self-energy of the dynamical mean field, which both are used in
the calculations. Picture taken from [4]

The principle behind DMFT is similar. All non-local correlations are aver-
aged out in the limit of infinite dimensions. The dynamics in DMFT comes
in because the mean-field in DMFT is a time dependent quantity, the effec-
tive Weiss-field Gσ = Gσ (t) or hef f = hef f (t). Thus, the possible states of the
single impurity site (|0i, | ↑i, | ↓i and | ↑↓i) are changing over time and so
the full dynamics of the problem can be captured. DMFT can describe tem-
poral correlations, but the spatial correlations are lost. The lattice problem
in dimension d → ∞ can be mapped onto a single interacting impurity site in

8
a bath of conduction electrons, hence an Anderson impurity model. To solve
this impurity problem several methods are available, in this work CT-QMC
(see section 2.1) was used.
The local propagator Gσ (iωn ) is given by a functional integral over Grass-
mann variables (4)
1Z Y †
Gσ (iωn ) = − Dcσ Dcσ [cσ (iωn )c†σ (iωn )]exp[−Sloc ] (4)
Z σ

where Z is the partition function:

Z Y
Z= Dc†σ Dcσ exp[−Sloc ] (5)
σ

and the local action Sloc is given in Eq. (6)

Z βZ β Z β
0
c†σ (τ 0 )Gσ−1 (τ 0 − τ )cσ (τ ) + U
X
Sloc = − dτ dτ n↑ (τ )n↓ (τ ) (6)
0 0 σ 0

Here Gσ is the effective local propagator respectively the bath Green function
or Weiss-mean-field.
After solving the AIM problem, the solution is used for the self-consistency
equation (7).

Σσ (iωn ) = Gσ−1 (iωn ) − G−1

σ (iωn ) (7)

The self-consistency equation yields the self-energy Σσ (iωn ), which is used

to calculate the lattice Greens function Gloc as shown in equation (8).
X 1
Gloc (iωn ) = (8)
k∈BZ
iωn − εk + µ − Σσ (iωn )

That lattice Green function is, then used to calculate a new Weiss-mean-field
(9)

Gσ−1 (iωn ) = G−1

loc (iωn ) + Σσ (iωn ) (9)

In DMFT, starting with an initial guess e.g. for Σσ (iωn ), this DMFT-cycle is
running until the requirement Gloc (iωn ) = Gσ (iωn ) is fulfilled. The DMFT-
cycle is mapped in Fig. 4 and for more details see [6] or [4].

1.1.4 Mott-Hubbard metal-insulator-transition

One of the most famous and intriguing phenomena in strongly correlated elec-
tron systems is the Mott-Hubbard-transition or metal-insulator-transition

9
 Figure 4: The DMFT cycle, taken from [7]

(MIT). It describes a first order phase transition between a paramagnetic

metal (PM) and a paramagnetic insulator (PI). This type of transition is
completely different from the transition between metals and band insulators,
which arises from the band filling of the electrons. The MIT is connected to
the correlation effects of strongly correlated electrons. It can be explained
through the decreasing probability of the double occupancy by an increas-
ing Coulomb-interaction U in the Hubbard-model. For higher interactions
the energy costs for two electrons to be on the same lattice site increases.
However the total energy of the system can be reduced by localizing the elec-
trons. Thus, after a critical interaction value is reached, the metal becomes
a paramagnetic insulator.
Fig. 5 shows the splitting of the density of states (DOS), i.e. the formation
of the Hubbard bands, and the emergence of the band gap, depending on
U scaled with the band-width W . The DOS of the non-interacting system
(a) converts to a three-peak-structure (c) with increasing U . This depends
on a splitting of the local propagator Gσ (iωn ) in an incoherent part, which
is responsible for the two peaks beside the Fermi-energy (EF ± U2 ), and a
coherent part at EF . The width of the coherent peak is proportional to the
quasi-particle weight Z and moreover Z ∝ m1∗ . For larger U the peak disap-
pears (d) thereby Z = 0 and m∗ → ∞. Thus the mobility of the electrons
vanishes and the system shows an insulating behavior.
A typical phase-diagram is shown in Fig. 6. The region inside the borders of
the MIT phases is called the coexistence region, where both phases (PM and
PI) can exist. The coexistence region is located at low temperatures between

10
Figure 5: Change of the density of states depending on the Hubbard U , taken
from [8]

the critical values on the metallic side Uc1 and on the insulating side Uc2 . The
point where Uc1 and Uc2 are crossing is the critical end point. At this point
the system has a second order transition. The real first order transition line
(dashed line) in the coexistence region is obtained by analyzing the minimum
of the free energy. It proceeds from the critical end point downwards and
bends slightly towards the Uc2 line, because the Uc2 line has a second order
transition at T = 0.

1.2 Motivation
Recent studies have shown that irreducible vertex divergences occur in sev-
eral models for strongly correlated electrons, such as the half-filled two-

11
 Figure 6: Typical structure of a MIT phase diagram, taken from [8]

dimensional Hubbard model, solved by means of DMFT [9], [10]. For this
case the corresponding phase-diagram and the irreducible divergence lines
are shown in the left panel of Fig. 7.
One finds multiple divergence lines surrounding the MIT, where the differ-
ent colors of divergence lines corresponds to the two distinct kinds of di-
vergences appearing in specific physical scattering channels (red: charge Γc ,
orange: charge Γc and particle-particle up-down Γpp,↑↓ simultaneously). In
the region of high T and large U the divergences in the Hubbard model are
approximated well by those of the Hubbard atom. On the other side the lines
are bending rightwards when the metallic behavior of the system increases.
Because of this the lines get for T → 0 finite values of U , which are dependent
on the density of states. The reason of this behavior at low temperatures is
not known yet. In particular the right panel of Fig. 7 shows the first diver-
gence line at low temperature for the Bethe-lattice and the 2D-square-lattice.
The motivation for this work is to compare the phase-diagrams for these two
lattices in particular the behavior of the corresponding MIT’s. Due to the
huge difference in the density of states of the two lattices a distinct behavior
of the low temperature coexistence regions of the corresponding MIT’s can be
expected. In this case the difference in the low T behavior of the divergence
lines might be attributed to that.

12
 (a) DMFT phase-diagram with the lines (b) First divergence line for the 2D square
surrounding the MIT in the half-filled lattice (solid line) and the Bethe lattice
square-lattice Hubbard model (dashed line) at low temperature.

Figure 7: Irreducible vertex divergences in the Hubbard model solved by

DMFT, pictures taken from [9]

1.2.1 Density of states

In this thesis, the Hubbard model is solved in DMFT using CT-QMC meth-
ods implemented in the w2dynamics package (see section 2.1). As a starting
point the non-interacting density of states (DOS) is needed. The following
sections show the non-interacting DOS for the two lattices.

Bethe-lattice DOS
The Bethe-lattice DOS has a simple form, allowing sometimes analytic cal-
culations, which is why it is widely used in theoretical solid-state physics.
This DOS has a semi-elliptic shape (see Fig. 8) and is given in Eq. (10),
[11]. A band-width of W = 2 implies a hopping integral value of t = 12 .
1 q 2 2
4 q 2
Dbethe (ε) = 4t − ε = W − 4ε2 (10)
2πt2 πW 2

Square-lattice DOS
The DOS of the square-lattice is constructed by the eigenvalues of the non-
interacting Hubbard-Hamilton (Eq. (3)). To get the eigenvalues we consider
the case of an one-dimensional tight-binding lattice. The non-interacting
part of the Hubbard-Hamilton, considering only the hopping between the
nearest neighbors is given in Eq. (11).
X †
cj−1,σ cj,σ + c†j+1,σ cj,σ

H = −t (11)
j,σ

13
Figure 8: The semi-elliptic Bethe-lattice density of states given by Eq. (10)

Here the lattice site index is j = 0, 1, 2, , N and the location of the lattice site
is given by xj = ja, with the lattice constant a. The total length of the lattice
is L = N a and the basis transformation to the momentum space is given by
Eq. (13).
1 X ikxj †
c†j,σ = √ e ck,σ (12)
L k,σ
1 X −ik0 xj
cj,σ = √ e ck0 ,σ (13)
L k0 ,σ

The periodic boundary condition x0 = xN implies k = 2π 0

L j. Assuming k =
k + G and j 0 = j + 1, k + G = L (j + 1) takes effect, G is here a reciprocal
2π

lattice vector. So the annihilation and creation operators for each lattice site
can be reduced to the first Brillouin zone.
1 X −ik0 xj 1 X −i(k+G)xj
cj,σ = √ e ck0 ,σ = √ e ck+G,σ =
L k0 ,σ L k+G,σ
1 X −ikxj iGxj
=√ e e ck+G,σ =
L k+G,σ (14)
 
= Gxj = jGa = j2π , ck+G,σ ≡ ck,σ  =
 

1 X −ikxj
=√ e ck,σ
L k,σ

14
Thus we get the Hamiltonian (Eq. (11)) in the momentum space and thereby
the eigenvalues εk , where c†k,σ ck,σ is the occupation number operator nk,σ .

L
!
X 1 X ika(j−1) −ikaj 1 X ika(j+1) −ikaj †
Hk,σ = −t e e + e e ck,σ ck,σ =
xj ,σ L k,σ L k,σ
X    (15)
−ika ika
X
= −t e +e nk,σ = −2t cos ka nk,σ
k,σ k,σ
 
⇒ ε~k = −2tcos ka (16)
The eigenvalues of the square-lattice are identical, up to the additional di-
mension, which yields
h    i
ε~k = −2t cos kx a + cos ky a (17)

To calculate the DOS of the square-lattice Eq. (18) is used.

1 X
D(ε) = δ(ε − ε~k ) (18)
Ld ~
k

This sum over the delta distributions with the eigenvalues of the square-
lattice is only numerically solvable. So the delta distribution must be re-
placed by a numerically usable approximation, in our case the Lorentzian
form δ∆ (x) = π1 x2 +∆
∆
2 , with lim δ∆ (x) = δ(x). In this case, the expression
∆→0
for the DOS reads as follows:
)−1
N X
N i 2
(" #
1 X ∆ h   
Dsquare (ε) ' 2
ε + 2t cos kx a + cos ky a + ∆2 (19)
(N a) kx ky π

The number of wave vectors N and the width of the Lorentzian form ∆ are
parameters of the numerical solution. The precision increases for increasing
N and decreasing ∆. As a first guess for ∆ the value of the imaginary part
of a self-energy for the Bethe-lattice at iωn = 0 and temperature β = 200 is
used. Polynomials in several intervals and of different orders are used for
the fitting of Im[Σ(iωn )]. These polynomials are given in the following and
plotted in Fig. 9.

• p1 (x) = a1 + b1 x with x ∈ [0, 0.05]

• p2 (x) = a2 + b2 x + c2 x2 with x ∈ [0, 0.1]

• p3 (x) = a3 + b3 x + c3 x2 with x ∈ [0, 0.2]

15
 • p4 (x) = a4 + b4 x + c4 x2 + d4 x3 with x ∈ [0, 0.2]

• p5 (x) = a5 + b5 x + c5 x2 + d5 x3 with x ∈ [0, 0.3]

• p6 (x) = a6 + b6 x + c6 x2 + d6 x3 with x ∈ [0, 0.4]

• p7 (x) = a7 + b7 x + c7 x2 + d7 x3 + e7 x4 with x ∈ [0, 0.2]

• p8 (x) = a8 + b8 x + c8 x2 + d8 x3 + e8 x4 with x ∈ [0, 0.3]

• p9 (x) = a9 + b9 x + c9 x2 + d9 x3 + e9 x4 with x ∈ [0, 0.4]

Figure 9: Imaginary part of the self-energy Σ(iωn ) of the Bethe-lattice at

β = 200 approximated by polynomials for an estimation of the scattering
rate γ, which is needed for the square-lattice DOS calculation.

The insert of Fig. 9 shows the value of Im[Σ(iωn = 0)] ∼ 0.025, which is a
guess for the maximal ∆. The value of N is set to a value, such that the
numerical solution for ∆ < 0.025 results in a smooth function. The square-
lattice DOS calculated with Eq. (19) is shown in Fig. 10, for two sets of
parameters ∆ and N , showing the big difference between the two DOS. In
contrast to the semi-elliptic shape of the Bethe-lattice DOS, the square-lattice
DOS has a Van-Hove-singularity, which is located where the dispersion rela-
tion ε~k has an extremum, i.e. where |∇ε ~ ~ | = 0.
k
In our case we set N = 1000, to save computing time whereas the value of

16
Figure 10: Square-lattice density of states numerically calculated by Eq.
(19), showing a logarithmic divergence at ε = 0. (a) uses ∆ = 0.001 and
N = 3000, (b) ∆ = 0.01 and N = 1000

Figure 11: Square-DOS calculated by w2dynamics, and exploited for the

DMFT calculations. Note the presence of the central peak, which would be
associated for infinite resolution (N → ∞) to a Van Hove singularity.

∆ = 10−10 . The resulting DOS of the square lattice is shown in Fig. 11.

17
2 Methods and Calculation
This chapter gives a short overview of the methods exploited to obtain the
Uc1 (T ) and Uc2 (T ) lines for both lattices discussed in this thesis. In partic-
ular, we will present results for the self-energy and the double occupancy.
Eventually, the corresponding phase-diagrams of the two lattices will be
shown.

2.1 CT-QMC
CT-QMC (Continuous-time quantum Monte Carlo) is a numerical method to
solve the impurity problem. That is the bottleneck of any DMFT cycle. The
Hamiltonian of the problem is separated into two terms, where one is exactly
solvable and the other is expanded into the interaction part (CT-INT) or
the hybridisation part (CT-HYB). The different expansions need different
algorithmic descriptions (refer to recent literature [4] [6] and [11]).
In this Bachelor thesis the CT-HYB implementation of the w2dynamics [12],
[13] has been used. At the center of this work is a single-band Hubbard model
with a finite bandwidth set to W = 2. In these two cases this also ensures,
that the variance of the non-interacting DOS is the same. To enforce a half-
filled band, which will be considered in this thesis, the chemical potential is
set to µ = U2 .

2.2 Self-energy and double occupancy

The main work of this bachelor thesis is the calculation of the DMFT phase-
diagrams for the two different lattice models considered. To achieve this goal
several calculations were performed, in particular of the self-energy Σ(iωn )
and the double occupancy d = hn↑ n↓ i, because their properties identify the
corresponding phase of the system.

2.2.1 Self-enegry
The self-energy Σ(iωn ) is defined by the Dyson Eq. (7), for more details see
section 1.1.3. The imaginary part of Σ(iωn ) indicates whether the system
is in the metallic or the insulating phase. This information is encoded in
the behavior of positive, low Matsubara frequencies. Fig. 12 shows the two
qualitatively different behaviors of the self-energy.
Physically, the self-energy is related to the mobility of the quasi-particles of
the strongly correlated electron system. Through an expansion of the self-
energy one obtains quantities, which are related to physical properties, e.g.

18
 (a) Metallic behavior (b) Insulating behavior

Figure 12: Behavior of the self-energy for the different phases (left: metallic;
right: Mott-insulating) of the Hubbard model, pictures taken from [1]

for the metallic case Im[Σ(iωn )] ∼ −γ − αωn . The value γ represents the
scattering rate and is related to the inverse lifetime τ of the quasi-particles.
The slope α in the second term defines the effective mass m∗ (see section
1.1.1). For an insulating behavior Im[Σ(iωn )] ∼ − ω1n holds, which yields
vanishing lifetime τ = 0 and thus, no free quasi-particles. For more details
see [11].
The procedure to obtain a converged result for the self-energy is described
briefly in the following. First, the convergence of the DMFT cycle is reached
by performing the calculations with a low relatively statistical precision (Fig.
13).
Once the DMFT convergence is found, only a few iterations are made with
a higher precision to reduce the statistical error. This procedure can be re-
peated for different interaction values U , taking into account a previous con-
verged calculation with low statistical error at the same T and a very close
interaction value. In this way a large range of interaction values is scanned
for a given temperature. Starting from low interaction values, a metallic
behavior is observed, where parameters γ and α are increasing with increas-
ing U. This reflects that the system becomes, gradually, a poorer conductor.
When U = Uc2 (T ) is reached, the behavior of Σ(iωn ) changes dramatically:
This is shown in Fig. 14 for β = 200.
When an interaction value U > Uc2 is reached the self-energy shows an in-
sulating behavior (orange points in Fig. 14). The self-energies for U values
around the critical interaction Uc2 for β = 200 are shown in Fig. 15. The
phase-diagram can be obtained by analyzing several temperatures T thus
way.
Fig. 16 shows analogously to Fig. 14 and Fig. 15 the self-energy for different

19
Figure 13: Imaginary part of the self-energy Σ(iωn ) for several DMFT iter-
ations, showing the convergence of the DMFT cycle on the metallic side at
U = 2.0 and β = 50 for the Bethe-lattice. Inset: Converged iterations with
small offsets between the calculations

interaction values, however the interaction is now decreasing. This evokes

that the system changes from an insulator to a metal, which happens at the
critical interaction value U = Uc1 (T ). Note that because of the first order
nature of the phase transition it applies Uc1 (T ) 6= Uc2 (T ) (with the exception
of the critical endpoints). From an algorithmic perspective, for obtaining the
insulating solution in the coexistence region, it is crucial to take into account
previous converged calculation, because the stability of the metallic phase is
typically greater. 1

1 Notethat, the value of the mixing parameter for the DMFT self-consistency is bigger.
(mixing parameter for metal: 0.6 and insulator: 0.95)

20
Figure 14: DMFT self-energy Σ(iωn ) for the Bethe-lattice for β = 200 and dif-
ferent interactions values U , corresponding to metallic (U = 2.20, 2.50, 2.66)
and insulating (U = 2.67) solutions.

21
Figure 15: Same as Fig. 14 with a finer resolution in U around Uc2 (∼ 2.6675)
for the Bethe-lattice at β = 200.

(a) U corresponding to metallic (U = 2.355) (b) U around Uc1

and insulating (U = 2.360, 2.440, 2.600) so-
lutions

Figure 16: DMFT self-energy Σ(iωn ) for the Bethe-lattice for β = 200 and
different interactions values U

22
2.2.2 Double occupancy
The other parameter analyzed in this work is the double occupancy d =
hn↑ n↓ i, which is relevant for the identification of the different phases around
the MIT and in the coexistence region. It is the expectation value of find-
ing a spin-up and spin-down electron on the same lattice site. This mea-
sures the average number of double occupied sites associated to an energy
cost equal to U : d is proportional to the potential energy of the system
Epot = hHpot i = U hn↑ n↓ i = U d. For a high interaction U , double occupan-
cies are energetically unfavorable and therefore suppressed. In our case, we
consider lattices at half-filling, i.e. n = n↑ + n↓ = 1. This implies that the
maximum of the double occupancy is hn↑ n↓ i = hn↑ ihn↓ i = 14 at U = 0, where
the system is uncorrelated. On the other hand, for the limit U → ∞ the
double occupancy approaches 0. Physically, it is easy to understand that
in the metallic phase the double occupancy is higher than for the insulating
phase, because of the larger mobility of the electrons, which corresponds in
average to a higher possibility of finding two electrons on the same lattice
site.

Figure 17: Expectation value of the double occupancy d = hn↑ n↓ i as a func-

tion of the interaction U at β = 75 for the Bethe-lattice

The DMFT results for the double occupancy as a function of the interaction
U is shown in Fig. 17 for the Bethe-lattice at temperature β = 75, for two
data sets. The purple set shows the calculations, where U was increased

23
starting from a metallic solution. Opposed to that the green set represents
the results for decreasing interaction values, starting from an insulating so-
lution. The area where the results of the two data sets are different is the
coexistence region. There, the metallic and the insulating phase are found
simultaneously, reflecting the typical hysteresis of the first order MIT. The
corresponding critical interaction values Uc1 and Uc2 are shown as well, taken
from the analysis of the behavior of the self-energy (see section 2.2.1). In
general, we observe a significant difference between the two data sets. More
in detail: for the purple data sets the double occupancy drops continuously
with the decreasing U and at Uc2 the double occupancy is quite close to
double occupancy on the other data set. For the green data set the dou-
ble occupancy increases only very weakly for decreasing U until an abrupt
change is observed at Uc1 .

Figure 18: Expectation value of the DMFT double occupancy d = hn↑ n↓ i as

function of the interaction U at β = 75 for the square-lattice

The calculations for the square-lattice, shown in Fig. 18 are qualitatively

similar to the ones of the Bethe-lattice. A quantitative difference is found,
however, for the location in parameter space of the coexistence region. For
the square-lattice the critical interactions Uc1 and Uc2 are slightly smaller
although the overall width of the coexistence region remains nearly equal.
This happens for all temperatures considered in this work. Fig. 19 shows
the comparison for β = 40, close to the critical endpoint and Fig. 20 for the

24
 (a) bethe-lattice (b) square-lattice

Figure 19: Expectation value of the double occupancy d = hn↑ n↓ i for the two
lattice models and different values of the interaction U at β = 40

lowest temperature calculated (β = 200).

As a general trend, we observe that the coexistence region gets continuously
broader with decreasing temperature. At the temperature of the critical end-
point (β ∼ 33), the data sets of the double occupancy merge, displaying a
kink at the critical U , while for the temperatures above it d is a continuous
function for all U .

Let us note that the CT-QMC calculations do not always preserve particle-
hole symmetry. This can be ascribed to a non-ergodicily of the DMFT cal-
culation for the insulating solution. The violation is stronger for lower tem-
peratures. As an example in Fig. 21, G(τ ) is composed with itself mirrored
around τ = β2 , which should not yield any difference. This has, however, no
influence on the results for the self-energy or the double occupancy (at least
for the precision level required by our analysis) meaning that it does not
affect the conclusions presented in this thesis.

25
 (a) bethe-lattice (b) square-lattice

Figure 20: Expectation value of the double occupancy d = hn↑ n↓ i for the two
lattice models and different values of the interaction U for calculations at the
lowest temperature β = 200

Figure 21: Local propagator G(τ ) for U = 2.55 at β = 100 on the insulating
side. Inset: Same plot on a logarithmic scale.

26
2.3 Phase-diagram of the Bethe-lattice
We can summarizing the results of our DMFT calculations described in sec-
tion 2.2, by plotting the corresponding phase-diagram of the Bethe-lattice,
shown in Fig.22.

Figure 22: DMFT phase-diagram of the half-filled Hubbard model on a

Bethe-lattice

The diagram shows the MIT as a function of the temperature T and the
U
interaction U in units of the half-bandwidth 2W . The data sets Uc1 (T ) and
Uc2 (T ) were fitted using gnuplot with the following fitting functions

• f1 (x) = a1 + b1 xc1 f or Uc1

• f2 (x) = a2 + b2 x + c2 x2 f or Uc2

Fig. 22 shows the increase of both critical interactions (Uc1 (T ), Uc2 (T )) with
decreasing temperature and as well an approximated estimate for the position
of the critical endpoint, where the Uc1 (T ) and Uc2 (T ) lines merges. In the
phase-diagram the different phases of the Hubbard model are also shown
with a coexistence region broadening with decreasing temperature.

27
2.4 Phase-diagram of the 2D-square-lattice
Analogously we have also obtained the DMFT phase-diagram for the Hub-
bard model on a square lattice, shown in Fig. 23.

Figure 23: DMFT phase-diagram of the half-filled Hubbard model on a

square-lattice

The data sets for Uc1 (T ) and Uc2 (T ) were calculated for the same temper-
atures and similarly fitted with the program gnuplot. Due to the relatively
small variation w.r.t. U , it is not a priori obvious which function yields the
best fit for Uc1 (T ).
Several functions used to fit Uc1 (T ) are shown in Fig. 24, and listed in the
following.

• g1 (x) = a1 + b1 x

• g2 (x) = a2 + b2 xc2

• g3 (x) = a3 + b3 ec3 x
1
• g4 (x) = 0.043 − a4 (x−2.3)
1
• g5 (x) = 0.043 − a b
5 (x−2.3) 5

28
 1
• g6 (x) = 0.9 − a b
6 (x−2.3) 6

1
• g7 (x) = 0.9 − a7 (x−b c
7) 7

Figure 24: Uc1 (T ) of the square-lattice fitted by different functions.

The function g1 is a linear fit. The fits g2 and also g3 are similar to the linear
solution of g1 . Polynomials with order greater than one yield a non-physical
solution, because the slope of the functions changes the sign and thus these
functions are not shown in Fig. 24. The functions g4 , g5 , g6 and g7 are
constructed to get a closer fit. This works for the functions g4 , g5 and g6 .
But for g7 , which has a additional fitting parameter, the solution is similar to
the linear fit of g1 . The function which yields the best solution is not distinct.
The choice of the function, which has been used for DMFT phase-diagram,
is the linear fit because it is the solution of three different fits. The Functions
shown Fig. 23:
• h1 (x) = a1 + b1 xc1 f or Uc1
• h2 (x) = a2 + b2 ec2 x f or Uc2
Fig. 23 shows that with decreasing temperature, Uc1 (T ) decreases (contrary
to the Bethe case) while Uc2 (T ) increases. The corresponding critical end-
point and the phases of the Hubbard model on a square lattice are shown in
the same way as in Fig. 22.

29
3 Comparison and Discussion
The Bethe-lattice and the square-lattice DMFT phase-diagrams are shown
in a single figure (Fig. 25) for comparison.

Figure 25: Comparison of the DMFT coexistence region of the Hubbard

model on a Bethe-lattice and on a square-lattice, respectively. The red area
corresponds to the Bethe-lattice and the blue one to the square-lattice.

At first Fig. 25 shows that for all calculated temperatures the critical inter-
actions Uc1 (T ) and Uc2 (T ) of the square-lattice have lower values U than the
ones of the Bethe-lattice. Thus the square-lattice coexistence region (blue
area) is slightly shifted leftwards to the Bethe-lattice one (red area). Sec-
ondly, Fig. 25 shows that the width of both coexistence regions are nearly
equal. The behavior of Uc2 (T ) for both lattices are similar and thus also
the fitted functions. A quantitative difference between the phase-diagrams
is the behavior of the Uc1 (T ). For the Uc1,bethe (T ) the critical interaction
values increase with decreasing T , below the critical endpoint. On the other
side, for Uc1,square (T ) the critical interaction values decrease. This gives the
square-lattice phase-diagram a qualitatively different shape to the Bethe-
lattice phase-diagram.

30
4 Conclusion
In this Bachelor thesis we studied the difference of the MIT in the DMFT so-
lution of the Hubbard model on the Bethe-lattice and the 2D-square-lattice.
Starting from the corresponding DOSes, exploiting the w2dynamics-package
as impurity solver, DMFT calculations have been done.
A typical output of the calculations is the self-energy and the double occu-
pancy, which are used to estimate Uc1 (T ) and Uc2 (T ), separately for both
lattices to get a lead whether the differences in the physical MIT can effect
the low temperature behavior of the irreducible vertex divergences.

(a) Phase-diagram of the Mott-Hubbard- (b) First divergence line for the 2D
transition on the Bethe-lattice and the square lattice (solid line) and the Bethe
2D-square-lattice lattice (dashed line) at low temperature,
taken from [9]

Figure 26: Comparison of the irreducible vertex divergences in the Hubbard

model and the Mott-Hubbard transition, of the Bethe-lattice and the 2D-
square-lattice

Fig. 26 shows the phase-diagrams and the divergence lines next to each other
for a better comparison. This illustrates that for decreasing temperature the
divergence lines are getting closer to the shape of Uc1 (T ) lines of the phase-
diagram. So the overall slope of divergence lines might indeed be related to
the Uc1 (T ) lines of the MIT. To confirm this assumption more studies are
necessary, also because the divergence of the irreducible vertexes might else
appear in model without a MIT [10].
More over, it would be also interesting to determine, through is a minimiza-
tion of the free energy whether the first order transition line of MIT is also
effected by the lattice details.

31
References
[1] A. Toschi. “Local electronic correlations”. In: QFT-Lecture TU Wien,
LVA 138.062 (2016).
[2] M. P. Marder. Condensed Matter Physics, 2nd Edition. 2nd ed. Wiley,
2010. Chap. 17.5.
[3] J. Hubbard. “Electron Correlations in Narrow Energy Bands”. In: Pro-
ceedings of the Royal Society of London. Series A, Mathematical and
Physical Sciences 276.1365 (1963), pp. 238–257. url: www.jstor.org/
stable/2414761.
[4] E. Pavarini et al. “The LDA+DMFT approach to strongly correlated
materials”. In: Schriften des Forschungszentrums Jülich Reihe Modeling
and Simulation 1 (2011).
[5] M. W. Kinza. “Single Impurity Anderson Model and Dynamical Mean
Field Theory”. In: Dissertation - RWTH Aachen (2003). url: http://
publications.rwth-aachen.de/record/229234/files/4979.pdf.
[6] A. Georges et al. “Dynamical mean-field theory of strongly correlated
fermion systems and the limit of infinite dimensions”. In: Rev. Mod.
Phys. 68 (1996), pp. 13–125. doi: 10.1103/RevModPhys.68.13.
[7] A. Toschi. “Strong Electronic Correlation in Dynamical Mean Field
Theory and beyond”. In: Habilitationsschrift TU Wien (2011).
[8] G. Kotliar and D. Vollhardt. “Strongly Correlated Materials: Insights
From Dynamical Mean-Field Theory”. In: Physics Today 57 (2004),
pp. 53–59.
[9] T. Schäfer et al. “Nonperturbative landscape of the Mott-Hubbard
transition: Multiple divergence lines around the critical endpoint”. In:
Phys. Rev. B. 94 (2016). doi: 10.1103/PhysRevB.94.235108. arXiv:
1606.03393.
[10] P. Chalupa et al. “Divergences of the irreducible vertex functions in cor-
related metallic systems: Insights from the Anderson Impurity Model”.
In: ArXiv e-prints (2017). arXiv: 1712.04171.
[11] N. Blümer. “Mott-Hubbard Metal-Insulator Transition and Optical
Conductivity in High Dimensions”. In: Dissertation - Universität Mainz
(2003). url: http://www.komet337.physik.uni-mainz.de/Dateien/
bluemer_color.pdf.

32
[12] M. Wallerberger. “w2dynamics: Continuous time quantum Monte Carlo
calculations of one- and two-particle propagators”. In: Dissertation -
TU Wien (2016). url: http://repositum.tuwien.ac.at/obvutwhs/
content/titleinfo/1370376.
[13] M. Wallerberger et al. “w2dynamics: Local one- and two-particle quan-
tities from dynamical mean field theory”. In: ArXiv e-prints (2018).
arXiv: 1801.10209.

33