Papers in Physics, vol. 9, art. 090005 (2017)

Received: 31 March 2017, Accepted: 06 June 2017
Edited by: D. Domı́nguez
Reviewed by: A. Feiguin, Northeastern University, Boston, United States.
Licence: Creative Commons Attribution 3.0
DOI: http://dx.doi.org/10.4279/PIP.090005

www.papersinphysics.org

ISSN 1852-4249

An efficient impurity-solver for the dynamical mean field theory
algorithm

Y. Núñez Fernández,1∗ K. Hallberg1

One of the most reliable and widely used methods to calculate electronic structure of
strongly correlated models is the Dynamical Mean Field Theory (DMFT) developed over
two decades ago. It is a non-perturbative algorithm which, in its simplest version, takes
into account strong local interactions by mapping the original lattice model on to a single
impurity model. This model has to be solved using some many-body technique. Several
methods have been used, the most reliable and promising of which is the Density Matrix
Renormalization technique. In this paper, we present an optimized implementation of
this method based on using the star geometry and correction-vector algorithms to solve
the related impurity Hamiltonian and obtain dynamical properties on the real frequency
axis. We show results for the half-filled and doped one-band Hubbard models on a square
lattice.

I. Introduction

Materials with strongly correlated electrons have
attracted researchers in the last decades. The
fact that most of them show interesting emergent
phenomena like superconductivity, ferroelectric-
ity, magnetism, metal-insulator transitions, among
other properties, has triggered a great deal of re-
search.

The presence of strongly interacting local or-
bitals that causes strong interactions among elec-
trons makes these materials very difficult to treat
theoretically. Very successful methods to calculate
electronic structure of weakly correlated materials,
such as the Density Functional Theory (DFT) [1],
lead to wrong results when used in some of these
systems. The DFT-based local density approxima-

∗E-mail: yurielnf@gmail.com

1 Centro Atómico Bariloche and Instituto Balseiro, CNEA,
CONICET, Avda. E. Bustillo 9500, 8400 San Carlos de
Bariloche, Ŕıo Negro, Argentina

tion (LDA) [2] and its generalizations are unable
to describe accurately the strong electron correla-
tions. Also, other analytical methods based on per-
turbations are no longer valid in this case so other
methods had to be envisaged and developed.

More than two decades ago, the Dynamical Mean
Field Theory (DMFT) was developed to study
these materials. This method and its successive
improvements [3–8] have been successful in incor-
porating the electronic correlations and more reli-
able calculations were done. The combination of
the DMFT with LDA allowed for band structure
calculations of a large variety of correlated ma-
terials (for reviews, see Refs. [9, 10]), where the
DMFT accounts more reliably for the local corre-
lations [11, 12].

The DMFT relies on the mapping of the cor-
related lattice onto an interacting impurity for
which the fermionic environment has to be deter-
mined self-consistently until convergence of the lo-
cal Green’s function and the local self-energy is
reached. This approach is exact for the infinitely

090005-1



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

coordinated system (infinite dimensions), the non-
interacting model and in the atomic limit. There-
fore, the possibility to obtain reliable DMFT solu-
tions of lattice Hamiltonians relies directly on the
ability to solve (complex) quantum impurity mod-
els.

Since the development of the DMFT, several
quantum impurity solvers were proposed and used
successfully; among these, we can mention the it-
erated perturbation theory (IPT) [13, 14], exact di-
agonalization (ED) [15], the Hirsch-Fye quantum
Monte Carlo (HFQMC) [16], the continuous time
quantum Monte Carlo (CTQMC) [17–20], non-
crossing approximations (NCA) [21], and the nu-
merical renormalization group (NRG) [22, 23]. All
of these methods imply certain approximations.
For a more detailed description, see [24].

Some years ago, we proposed the Density Ma-
trix Renormalization Group (DMRG) as a reliable
impurity-solver [25–27] which allows to surmount
some of the problems existing in other solvers, giv-
ing, for example, the possibility of calculating dy-
namical properties directly on the real frequency
axis. Other related methods followed, such as in
[28, 29]. This way, more accurate results can be ob-
tained than, for example, using algorithms based
on Monte Carlo techniques. The scope of this pa-
per is to detail the implementation of this method
and to show recent applications and potential uses.

II. DMFT in the square lattice

We will consider the Hubbard model on a square
lattice:

H = t
∑
〈ij〉σ

c
†
iσcjσ + U

∑
i

ni↑ni↓ −µ
∑
i

ni, (1)

where ciσ

(
c
†
iσ

)
annihilates (creates) an electron

with spin σ =↑,↓ at site i, niσ = c
†
iσciσ is the den-

sity operator, ni = ni↓ + ni↑, U is the Coulomb
repulsion, µ is the chemical potential, and 〈ij〉 rep-
resents nearest neighbor sites.

Changing to the Bloch basis d
†
k, the non-

interacting part becomes:

H0 =
∑
k,σ

t(k)d
†
kσdkσ, (2)

with t(k) = 2t (cos kx + cos ky) − µ. The Green’s
function for (1) is hence given by:

G(k,ω) = [ω − t(k) − Σ(k,ω)]−1 , (3)

where Σ(k,ω) is the self-energy.
The DMFT makes a local approximation of

Σ(k,ω), that is, Σ(k,ω) ≈ Σ(ω). This locality of
the magnitudes allows us to map the lattice prob-
lem onto an auxiliar impurity problem that has the
same local magnitudes G(ω) and Σ(ω). The im-
purity is coupled to a non-interacting bath, which
should be determined iteratively. The Hamiltonian
can be written:

Himp = Hloc + Hb, (4)

where Hloc is the local part of (1)

Hloc = −µn0 + Un0↑n0↓, (5)

and the non-interacting part Hb representing the
bath is:

Hb =
∑
iσ

λib
†
iσbiσ +

∑
iσ

vi

[
b
†
iσc0σ + H.c.

]
, (6)

where b
†
iσ represents the creation operator for the

bath-site i and spin σ, label “0” corresponds to the
interacting site.

The algorithm is summarized as:

(i) Start with Σ(ω) = 0.

(ii) Calculate the Green’s function for the local in-
teracting lattice site:

G(ω) =
1

N

∑
k

G(k,ω) (7)

=
1

N

∑
k

[ω − t(k) − Σ(ω)]−1 .

(iii) Calculate the hybridization

Γ(ω) = ω + µ− Σ(ω) − [G(ω)]−1 . (8)

(iv) Find a Hamiltonian representation Himp with
hybridization Γd(ω) to approximate Γ(ω). The
hybridization Γd(z) is characterized by the pa-
rameters vi and λi of Himp through:

090005-2



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

Γd(ω) =
∑
i

v2i
ω −λi

. (9)

(v) Calculate the Green’s function Gimp(ω) at
the impurity of the Hamiltonian Himp using
DMRG.

(vi) Obtain the self-energy

Σ(ω) = ω + µ− [Gimp(ω)]
−1 − Γd(ω). (10)

Return to (ii) until convergence.

At step (iv) we should find the parameters vi and
λi by fitting the calculated hybridization Γ(ω) us-
ing expression (9). At half-filling, because of the
electron-hole symmetry, we have Γ(ω) = Γ(−ω)
and hence λ−i = −λi, and v−i = vi, where the
bath index i goes from −p to p, and it does not
include i = 0 for an even number of bath sites 2p.

Almost all of the computational time is spent
at step (v), where the dynamics of a single impu-
rity Anderson model (SIAM) (see Fig. 1) is cal-
culated. We use the correction-vector for DMRG
following [30]. The one-dimensional representation
of the problem (needed for a DMRG calculation) is
as showed in Fig. 1, except that for the spin degree
of freedom we duplicate the graph, generating two
identical chains, one for each spin. Moreover, it
should be noticed that this is not a local or short-
range 1D Hamiltonian (usually called chain geom-
etry, where the DMRG is supposed to work very
well). However, we refer to [31, 32] where strong
evidence of better performance of the DMRG for
this kind of geometry (star geometry) compared to
chain geometry is presented.

The correction-vector for DMRG consists of tar-
geting not only the ground state |E0〉 of the system
but also the correction-vector |Vi〉 associated to the
frequency ωi (and its neighborhood), that is:

(ωi + iη −Himp −E0) |Vi〉 = c
†
0 |E0〉 , (11)

where a Lorentzian broadering η is introduced to
deal with the poles of a finite-length SIAM. For
a better matching between the ω windows (with
width approximately η), we target the correction
vectors of the extremes of the window. Once the
DMRG is converged, the Green’s function is eval-
uated for a finer mesh (around 0.2 of the original

Figure 1: Schematic representation of the impurity
problem for the DMFT. The circles (square) repre-
sent the non-interacting (interacting) sites, and the
lines correspond to the hoppings. Top: star geome-
try drawn in two ways. Bottom: 1D representation
as used for DMRG calculations.

window) [30]. In this way, a suitable renormalized
representation of the operators is obtained to cal-
culate the properties of the excitations around ωi,
particularly the Green’s function.

In what follows, we present results for a paradig-
matic correlated model using the method described
above.

III. Results

We have used this method to calculate the density
of states (DOS) of the Hamiltonian (Eq. 1) on a
square lattice with unit of energy t = 0.25, for sev-
eral dopings, given by the chemical potential. We
consider a discarded weight of 10−11 in the DMRG
procedure for which a maximum of around m = 128
states were kept, even for the largest systems (50
sites). For these large systems, the ground state
takes around 20 minutes to converge and each fre-
quency window, between 5 and 20 minutes. This is
an indication of the good efficiency of the method.

The metal-insulator Mott’s transition at half-
filling is showed in Fig. 2. The transition occurs
between U = 3 and U = 4. In Fig. 3, we observe
that the metallic character of the bands remains
robust under doping for a given value of the in-
teraction, showing a weight transfer between the
bands due to the correlations. The metallic char-
acter is also seen in the variation of the filling with
µ.

090005-3



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

 0

 0.1

 0.2

 0.3

 0.4

 0.5

 0.6

 0.7

−4 −3 −2 −1  0  1  2  3  4

−
1

/π
 I

m
[G

(ω
+

 i
η

)]

U=1
U=2
U=3
U=4

−30

−25

−20

−15

−10

−5

 0

−4 −3 −2 −1  0  1  2  3  4

Im
[Σ

(ω
+

 ι
η

)]

ω

U=1
U=2
U=3
U=4

Figure 2: Top: Density of states for U = 1, 2, 3, 4 at
half-filling. We use a bath with 30-50 sites per spin
and a Lorentzian broadening η = 0.12. The Fermi
energy is located at ω = 0. Bottom: Imaginary
part of the self-energy.

Figure 4 shows our results for a larger value of the
interaction U, for which we find a regime of dop-
ing having an insulating character. However, for a
large enough doping (obtained for a large negative
value of the chemical potential), the systems turn
metallic and acquire a large density of states at
the Fermi energy. While the system is insulating,
changing the chemical potential only results in a
rigid shift of the density of states. The small finite
values of the DOS at the Fermi energy for the insu-
lating cases are due to the Lorentzian broadening
η, see Eq. (11).

IV. Conclusions

We have presented here an efficient algorithm to
calculate dynamical properties of correlated sys-

 0

 0.05

 0.1

 0.15

 0.2

 0.25

 0.3

 0.35

 0.4

 0.45

 0.5

−3 −2 −1  0  1  2  3  4

−
1

/π
 I

m
[G

(ω
+

 i
η

)]

ω

µ=−0.0
µ=−0.5
µ=−1.0

 0.8

 0.85

 0.9

 0.95

 1

−1 −0.8 −0.6 −0.4 −0.2  0

fi
ll

in
g

µ

Figure 3: Top: Density of states for U = 3, same
parameters as in Fig. 2, and several chemical po-
tentials (µ = 0 corresponds to the half-filled case).
Bottom: Filling vs chemical potential showing a
metallic behavior.

tems such as the electronic structure for any dop-
ing. It is based on the Dynamical Mean Field
theory method where we use the Density Matrix
Renormalization Group (DMRG) as the impurity
solver. By using the star geometry for the hy-
bridization function (which reduces the entangle-
ment enhancing the performance of the DMRG for
larger bath sizes) together with the correction vec-
tor technique(which accurately calculates the dy-
namical response functions within the DMRG) we
were able to obtain reliable real axis response func-
tions, in particular, the density of states, for any
doping, for the Hubbard model on a square lattice.
This improvement will allow for the calculation of
dynamical properties on the real energy axis for
complex and more realistic correlated systems.

090005-4



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

 0

 0.1

 0.2

 0.3

 0.4

 0.5

−3 −2 −1  0  1  2  3  4  5

−
1

/π
 I

m
[G

(ω
+

 i
η

)]

ω

µ=−0.0
µ=−0.5
µ=−1.0
µ=−1.5

 0.8

 0.85

 0.9

 0.95

 1

−1.6 −1.4 −1.2 −1 −0.8 −0.6 −0.4 −0.2  0

fi
ll

in
g

µ

Figure 4: Top: Density of states for U = 4, same
parameters as in Fig. 2, and several chemical po-
tentials (µ = 0 corresponds to the half-filled case).
Bottom: Filling vs chemical potential showing the
transition from a metal to an insulator.

Acknowledgements - We thank Daniel Garćıa for
useful discussions.

[1] P Hohenberg, W Kohn, Inhomogeneous elec-
tron gas, Phys. Rev. 136, B864 (1964).

[2] R O Jones, O Gunnarsson, The den-
sity functional formalism, its applications and
prospects, Rev. Mod. Phys. 61, 689 (1989).

[3] G Kotliar, D Vollhardt, Strongly correlated
materials: Insights from dynamical mean-field
theory, Physics Today 57, 53 (2004).

[4] A Georges, G Kotliar, W Krauth, M J Rozen-
berg, Dynamical mean-field theory of strongly
correlated fermion systems and the limit of in-

finite dimensions, Rev. Mod. Phys. 68, 13
(1996).

[5] G Kotliar, S Y Savrasov, G Pálsson, G Biroli,
Cellular dynamical mean field approach to
strongly correlated systems, Phys. Rev. Lett.
87, 186401 (2001).

[6] T Maier, M Jarrell, T Pruschke, M H Hettler,
Quantum cluster theories, Rev. Mod. Phys.
77, 1027 (2005).

[7] M H Hettler, A N Tahvildar-Zadeh, M Jar-
rell, T Pruschke, H R Krishnamurthy, Nonlo-
cal dynamical correlations of strongly interact-
ing electron systems, Phys. Rev. B 58, R7475
(1998).

[8] D Sénéchal, D Perez, M Pioro-Ladrière, Spec-
tral weight of the Hubbard model through clus-
ter perturbation theory, Phys. Rev. Lett. 84,
522 (2000).

[9] M Imada, T Miyake, Electronic structure cal-
culation by first principles for strongly corre-
lated electron systems, J. Phys. Soc. Jpn. 79,
112001 (2010).

[10] K Held, Electronic structure calculations using
dynamical mean-field theory, Adv. in Phys. 56,
829 (2007).

[11] V I Anisimov, A I Poteryaev, M A Korotin,
A O Anokhin, G Kotliar, First-principles cal-
culations of the electronic structure and spec-
tra of strongly correlated systems: dynamical
mean-field theory, J. Phys. Condens. Mat. 9,
7359 (1997).

[12] A I Lichtenstein, M I Katsnelson, Ab initio
calculations of quasiparticle band structure in
correlated systems: LDA++ approach, Phys.
Rev. B 57, 6884 (1998).

[13] A Georges, G Kotliar, Hubbard model in in-
finite dimensions, Phys. Rev. B 45, 6479
(1992).

[14] M J Rozenberg, G Kotliar, X Y Zhang, Mott-
Hubbard transition in infinite dimensions. II,
Phys. Rev. B 49, 10181 (1994).

090005-5



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

[15] M Caffarel, W Krauth, Exact diagonalization
approach to correlated fermions in infinite di-
mensions: Mott transition and superconduc-
tivity, Phys. Rev. Lett. 72, 1545 (1994).

[16] J E Hirsch, R M Fye, Monte Carlo method
for magnetic impurities in metals, Phys. Rev.
Lett. 56, 2521 (1986).

[17] A N Rubtsov, V V Savkin, A I Lichten-
stein, Continuous-time quantum Monte Carlo
method for fermions, Phys. Rev. Lett. 72,
035122 (2005).

[18] P Werner, A Comanac, L de Medici, M Troyer,
A J Millis, Continuous-Time Solver for Quan-
tum Impurity Models, Phys. Rev. Lett. 97,
076405 (2006).

[19] H Park, K Haule, G Kotliar, Cluster Dynami-
cal Mean Field Theory of the Mott Transition,
Phys. Rev. Lett. 101, 186403 (2008).

[20] E Gull, A J Millis, A I Lichtenstein, A N
Rubtsov, M Troyer, P Werner, Continuous-
time Monte Carlo methods for quantum impu-
rity models, Rev. Mod. Phys. 83, 349 (2011).

[21] T Pruschke, D L Cox, M Jarrell, Hubbard
model at infinite dimensions: Thermodynamic
and transport properties, Phys. Rev. Lett. 47,
3553 (1993).

[22] K G Wilson, The renormalization group: Crit-
ical phenomena and the Kondo problem, Rev.
Mod. Phys. 47, 773 (1975).

[23] R Bulla, Zero temperature metal-insulator
transition in the infinite-dimensional hubbard
model, Phys. Rev. Lett. 83, 136 (1999);

R Bulla, A C Hewson, T Pruschke, Numeri-
cal renormalization group calculations for the
self-energy of the impurity Anderson model, J.
Phys. Condens. Mat. 10, 8365 (1998).

[24] K Hallberg, D J Garćıa, P Cornaglia, J Fa-
cio, Y Núñez-Fernández, State-of-the-art tech-
niques for calculating spectral functions in
models for correlated materials, EPL 112,
17001 (2015).

[25] D J Garćıa, K Hallberg, M J Rozenberg, Dy-
namical mean field theory with the density ma-
trix renormalization group, Phys. Rev. Lett.
93, 246403 (2004).

[26] D J Garćıa, E Miranda, K Hallberg, M
J Rozenberg, Mott transition in the Hub-
bard model away from particle-hole symmetry,
Phys. Rev. B 75, 121102 (2007);

E Miranda, D J Garćıa, K Hallberg, M J
Rozenberg, The metal-insulator transition in
the paramagnetic Hubbard Model, Physica B:
Cond. Mat. 403, 1465 (2008);

D J Garćıa, E Miranda, K Hallberg, M J
Rozenberg, Metal-insulator transition in cor-
related systems: A new numerical approach,
Physica B: Cond. Mat. 398, 407 (2007);

S Nishimoto, F Gebhard, E Jeckelmann, Dy-
namical density-matrix renormalization group
for the Mott-Hubbard insulator in high dimen-
sions, J. Phys. Condens. Mat. 16, 7063 (2004);

M Karski, C Raas, G Uhrig, Electron spectra
close to a metal-to-insulator transition, Phys.
Rev. B 72, 113110 (2005);

M Karski, C Raas, G Uhrig, Single-particle
dynamics in the vicinity of the Mott-Hubbard
metal-to-insulator transition, Phys. Rev. B 75,
075116 (2008);

C Raas, P Grete, G Uhrig, Emergent Col-
lective Modes and Kinks in Electronic Disper-
sions, Phys. Rev. Lett. 102, 076406 (2009).

[27] Y Núñez Fernández, D Garćıa, K Hallberg,
The two orbital Hubbard model in a square lat-
tice: a DMFT + DMRG approach, J. Phys.:
Conf. Ser. 568, 042009 (2014).

[28] M Ganahl et al, Efficient DMFT impurity
solver using real-time dynamics with matrix
product states, Phys. Rev. B 92, 155132
(2015).

[29] F Wolf, J Justiniano, I McCulloch, U
Schollwöck, Spectral functions and time evolu-
tion from the Chebyshev recursion, Phys. Rev.
B 91, 115144 (2015).

[30] T D Kühner, S R White, Dynamical correla-
tion functions using the density matrix renor-
malization group, Phys. Rev. B 60, 335 (1999).

090005-6



Papers in Physics, vol. 9, art. 090005 (2017) / Y. Núñez Fernández et al.

[31] A Holzner, A Weichselbaum, J von Delft, Ma-
trix product state approach for a two-lead mul-
tilevel Anderson impurity model, Phys. Rev.
B 81, 125126 (2010).

[32] F Alexander Wolf, I McCulloch, U Schollwöck,
Solving nonequilibrium dynamical mean-field
theory using matrix product states, Phys. Rev.
B 90, 235131 (2014).

090005-7