F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach
We derive an alternative representation for the f-electron spectral function of the Falicov-Kimball model from the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted to the real time axis, allowing us to go to arbitrarily low temperatures. The g...
Gespeichert in:
| Veröffentlicht in: | Condensed Matter Physics |
|---|---|
| Datum: | 2008 |
| Hauptverfasser: | , |
| Format: | Artikel |
| Sprache: | English |
| Veröffentlicht: |
Інститут фізики конденсованих систем НАН України
2008
|
| Online Zugang: | https://nasplib.isofts.kiev.ua/handle/123456789/119339 |
| Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Zitieren: | F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach / A.M. Shvaika, J.K. Freericks // Condensed Matter Physics. — 2008. — Т. 11, № 3(55). — С. 425-442. — Бібліогр.: 25 назв. — англ. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
nasplib_isofts_kiev_ua-123456789-119339 |
|---|---|
| record_format |
dspace |
| spelling |
Shvaika, A.M. Freericks, J.K. 2017-06-06T13:42:24Z 2017-06-06T13:42:24Z 2008 F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach / A.M. Shvaika, J.K. Freericks // Condensed Matter Physics. — 2008. — Т. 11, № 3(55). — С. 425-442. — Бібліогр.: 25 назв. — англ. 1607-324X PACS: 71.10.-w, 71.27.+a, 71.30.+h, 02.30.Rz DOI:10.5488/CMP.11.3.425 https://nasplib.isofts.kiev.ua/handle/123456789/119339 We derive an alternative representation for the f-electron spectral function of the Falicov-Kimball model from the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted to the real time axis, allowing us to go to arbitrarily low temperatures. The general formula for the retarded Green's function involves two determinants of continuous matrix operators that have the Toeplitz form. By employing the Wiener-Hopf sum equation approach and Szeg¨ o's theorem, we can derive exact analytic formulas for the large-time limits of the Green's function; we illustrate this for cases when the logarithm of characteristic function (which de nes the continuous Toeplitz matrix) does and does not wind around the origin. We show how accurate these asymptotic formulas are to the exact solutions found from extrapolating matrix calculations to the zero discretization size limit. Отримано нове представлення для спектральної функцiї f-електронiв моделi Фалiкова-Кiмбала, яке альтернативне оригiнальному представленню Брандта i Урбанека. У новому представленнi усi розрахунки виконуються тiльки на дiйснiй часовiй осi, що дозволяє розглядати як завгодно низькi температури. Загальний вираз для запiзнюючої функцiї Ґрiна включає два детермiнанти неперервних матричних операторiв зi структурою типу Теплiца. Застосовуючи дискретний пiдхiд Вiнера-Гопфа i теорему Сеґо, отримано точнi аналiтичнi формули для довгочасової поведiнки функцiй Ґрiна; розглянуто випадки, коли логарифм характеристичної функцiї (яка визначає неперервну матрицю Теплiца) робить i не робить виток навколо початку координат. Показано наскiльки точними є данi асимптотичнi вирази у порiвняннi з точними розв’язками, якi отримуються при екстраполяцiї прямих матричних розрахункiв до границi нульової дискретизацiї. It is with great pleasure that we honor Prof. Stasyuk with this contribution. This publication is based on work supported by Award No. UKP2{2697{LV{06 of the U.S. Civilian Research & Development Foundation (CRDF). J.K.F. also acknowledges support by the National Science Foundation under grant number DMR{0705266. en Інститут фізики конденсованих систем НАН України Condensed Matter Physics F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach Спектральна функцiя f-електронiв для моделi Фалiкова-Кiмбала i дискретний пiдхiд Вiнера-Гопфа Article published earlier |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| title |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach |
| spellingShingle |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach Shvaika, A.M. Freericks, J.K. |
| title_short |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach |
| title_full |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach |
| title_fullStr |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach |
| title_full_unstemmed |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach |
| title_sort |
f-electron spectral function of the falicov-kimball model and the wiener-hopf sum equation approach |
| author |
Shvaika, A.M. Freericks, J.K. |
| author_facet |
Shvaika, A.M. Freericks, J.K. |
| publishDate |
2008 |
| language |
English |
| container_title |
Condensed Matter Physics |
| publisher |
Інститут фізики конденсованих систем НАН України |
| format |
Article |
| title_alt |
Спектральна функцiя f-електронiв для моделi Фалiкова-Кiмбала i дискретний пiдхiд Вiнера-Гопфа |
| description |
We derive an alternative representation for the f-electron spectral function of the Falicov-Kimball model from
the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted
to the real time axis, allowing us to go to arbitrarily low temperatures. The general formula for the retarded
Green's function involves two determinants of continuous matrix operators that have the Toeplitz form. By employing
the Wiener-Hopf sum equation approach and Szeg¨ o's theorem, we can derive exact analytic formulas
for the large-time limits of the Green's function; we illustrate this for cases when the logarithm of characteristic
function (which de nes the continuous Toeplitz matrix) does and does not wind around the origin. We show
how accurate these asymptotic formulas are to the exact solutions found from extrapolating matrix calculations
to the zero discretization size limit.
Отримано нове представлення для спектральної функцiї f-електронiв моделi Фалiкова-Кiмбала, яке альтернативне оригiнальному представленню Брандта i Урбанека. У новому представленнi усi розрахунки виконуються тiльки на дiйснiй часовiй осi, що дозволяє розглядати як завгодно низькi температури. Загальний вираз для запiзнюючої функцiї Ґрiна включає два детермiнанти неперервних матричних операторiв зi структурою типу Теплiца. Застосовуючи дискретний пiдхiд Вiнера-Гопфа i теорему Сеґо, отримано точнi аналiтичнi формули для довгочасової поведiнки функцiй Ґрiна; розглянуто випадки, коли логарифм характеристичної функцiї (яка визначає неперервну матрицю Теплiца) робить i не робить виток навколо початку координат. Показано наскiльки точними є данi асимптотичнi вирази у порiвняннi з точними розв’язками, якi отримуються при екстраполяцiї прямих матричних розрахункiв до границi нульової дискретизацiї.
|
| issn |
1607-324X |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/119339 |
| citation_txt |
F-electron spectral function of the Falicov-Kimball model and the Wiener-Hopf sum equation approach / A.M. Shvaika, J.K. Freericks // Condensed Matter Physics. — 2008. — Т. 11, № 3(55). — С. 425-442. — Бібліогр.: 25 назв. — англ. |
| work_keys_str_mv |
AT shvaikaam felectronspectralfunctionofthefalicovkimballmodelandthewienerhopfsumequationapproach AT freericksjk felectronspectralfunctionofthefalicovkimballmodelandthewienerhopfsumequationapproach AT shvaikaam spektralʹnafunkciâfelektronivdlâmodelifalikovakimbalaidiskretniipidhidvineragopfa AT freericksjk spektralʹnafunkciâfelektronivdlâmodelifalikovakimbalaidiskretniipidhidvineragopfa |
| first_indexed |
2025-11-25T18:25:28Z |
| last_indexed |
2025-11-25T18:25:28Z |
| _version_ |
1850521375272337408 |
| fulltext |
Condensed Matter Physics 2008, Vol. 11, No 3(55), pp. 425–442
F-electron spectral function of the Falicov-Kimball
model and the Wiener-Hopf sum equation approach
A.M.Shvaika1, J.K.Freericks2
1 Institute for Condensed Matter Physics of the National Academy of Sciences of Ukraine, 1 Svientsitskii Str.,
79011 Lviv, Ukraine
2 Department of Physics, Georgetown University, 37th and O Sts. NW, Washington, DC 20057, USA
Received May 21, 2008, in final form June 20, 2008
We derive an alternative representation for the f -electron spectral function of the Falicov-Kimball model from
the original one proposed by Brandt and Urbanek. In the new representation, all calculations are restricted
to the real time axis, allowing us to go to arbitrarily low temperatures. The general formula for the retarded
Green’s function involves two determinants of continuous matrix operators that have the Toeplitz form. By em-
ploying the Wiener-Hopf sum equation approach and Szegö’s theorem, we can derive exact analytic formulas
for the large-time limits of the Green’s function; we illustrate this for cases when the logarithm of characteristic
function (which defines the continuous Toeplitz matrix) does and does not wind around the origin. We show
how accurate these asymptotic formulas are to the exact solutions found from extrapolating matrix calculations
to the zero discretization size limit.
Key words: F-electron spectral function, Falicov-Kimball model, Wiener-Hopf approach, dynamical
mean-field theory
PACS: 71.10.-w, 71.27.+a, 71.30.+h, 02.30.Rz
1. Introduction
The Falicov-Kimball model [1] has often been studied as one of the simplest models of strong
electron correlations. Indeed, the exact solution in the limit of large dimensions displays charge
density wave order, the Mott metal-insulator transition, and phase separation [2]. It is particularly
relevant to present a paper on this model in an issue dedicated to Prof. Stasyuk, as he has worked
on this model and closely related ones throughout his career; this began with his first article on the
Falicov-Kimball model [3] while his most relevant work to this contribution involves approximate
solutions for the f -electron spectral functions [4–6].
Our focus in this work is on the spectral function of the f -electrons in the exact solution
of Falicov-Kimball model with dynamical mean-field theory. This problem was first investigated
by Brandt and Urbanek [7] and by Janis [8]. The numerical approach of Brandt and Urbanek
was extended by Zlatić et al. [9] and Freericks, Turkowski and Zlatić [10–12] within the original
Brandt-Urbanek framework. An alternative approach was taken by employing the numerical renor-
malization group by Anders and Czycholl [13]. Some analytic work was completed by Liu [14] which
was related to early work on this problem [15] and proceeded from the perspective of the X-ray
edge problem [16,17]. Here, we develop a new approach which is numerically more tractable than
other approaches and allows us to develop asymptotically exact formulas for the Green’s function
in the time representation.
The single-site Hamiltonian of the Falicov-Kimball model [1] involves two types of electrons—
conduction electrons, denoted by the letter d and localized electrons, denoted by the letter f
Hloc = Efnf + Undnf − µ (nd + nf ) = H0(1− nf ) +H1nf , (1.1)
H0 = −µnd , (1.2)
H1 = Ef − µ+ (U − µ)nd , (1.3)
c© A.M.Shvaika, J.K.Freericks 425
A.M.Shvaika, J.K.Freericks
-i β
0 t
Figure 1. The Kadanoff-Baym-Keldysh contour. The contour starts at time 0, moves along the
real time axis to time t, moves back along the real axis to time 0, then extends down the
imaginary axis to time −iβ.
where nd = d†d and nf = f †f are the number operators for the d- and f electrons, respectively, U
is the on-site Coulomb interaction between the d and f electrons, and Ef is the site-energy of the
f state. The full Hamiltonian on the lattice includes a repeat of this local Hamiltonian for each
lattice site and a hopping of the conduction electrons between nearest-neighbor sites. The density
matrix of the single-impurity problem is equal to
ρ = e−βHlocTc exp
{
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t
′, t′′)d(t′′)
}
1
Z , (1.4)
where the time-ordering and integration are performed over the Kadanoff-Baym-Keldysh con-
tour [18,19] and β = 1/T is the inverse temperature; the Kadanoff-Baym-Keldysh contour is shown
in figure 1 – it starts at t = 0, runs out along the real axis to a maximal time tmax, then returns
along the real axis back to t = 0, and finally runs along the negative imaginary axis down to −iβ.
Here, the dynamical mean-field λc(t
′, t′′) and chemical potential µ are taken from the equilibrium
solution of the conduction electron problem via dynamical mean-field theory [2]. In other words,
the dynamical mean field satisfies
λc(t, t
′) = − i
π
∞∫
−∞
dωImλ(ω)e−iω(t−t′)[f(ω)−Θc(t, t
′)], (1.5)
where λ(ω) is the ordinary dynamical mean field on the real axis, extracted from the dynamical
mean-field theory solution of the model, f(ω) = 1/[1 + exp(βω)] is the Fermi-Dirac distribution
function, and Θc(t, t
′) is the Heaviside function on the contour, equal to 1 when t is ahead of t′
on the contour, equal to 0 when t is behind t′ on the contour and equal to 1/2 when t = t′ on
the contour (note that this definition of the λc field is i times the definition in [2]). The partition
function for the single-site problem contains two terms:
Z = Tr e−βHlocTc exp
{
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t′, t′′)d(t′′)
}
= Z0 + Z1 , (1.6)
which correspond to the different occupations of the f -particles, as follows:
Z0 =
[
1 + eβµ
]∏
m
(
1− λm
iωm + µ
)
, nf = 0; (1.7)
Z1 = eβ(µ−Ef )
[
1 + eβ(µ−U)
]∏
m
(
1− λm
iωm + µ− U
)
, nf = 1. (1.8)
Here, iωm = iπT (2m+ 1) is the fermionic Matsubara frequency and
λm =
β∫
0
dτeiωmτλ(τ) (1.9)
426
F-electron spectral function of the Falicov-Kimball model . . .
is the dynamical mean field evaluated at the mth Matsubara frequency [the λc field is a function
only of the difference of the two time arguments when both times lie on the imaginary part of the
Kadanoff-Baym-Keldysh contour and we use the notation λ(τ) = −iλc(−iτ, 0), which corresponds
to the conventional definition of the dynamical mean field for imaginary time].
2. Real-time Green’s functions
We consider the contour ordered Green’s function for the f -electrons (with times t and t′ both
lying on the contour)
Gc
f (t, t′) = −i
〈
Tcf(t)f †(t′)
〉
, (2.1)
where the time ordering is taken along the contour, the angle brackets denote taking a trace
(weighted by the density matrix for the equilibrium system), and the time dependence of the
fermionic operators is taken with respect to the local Hamiltonian Hloc. Depending on the location
of the time arguments on the contour, we obtain different real-time Green’s functions:
the greater Green’s function
G>
f (t− t′) = −i
〈
f(t)f †(t′)
〉
= −i
1
Z
∑
pq
e−βεp 〈p |f | q〉
〈
q
∣∣f †
∣∣ p
〉
ei(εp−εq)(t−t′); (2.2)
the lesser Green’s function
G<
f (t− t′) = i
〈
f †(t′)f(t)
〉
= i
1
Z
∑
pq
e−βεq
〈
q
∣∣f †
∣∣ p
〉
〈p |f | q〉 ei(εp−εq)(t−t′); (2.3)
the time-ordered Green’s function
Gt
f (t− t′) = −i
〈
Ttf(t)f †(t′)
〉
= Θ(t− t′)G>
f (t− t′) + Θ(t′ − t)G<
f (t− t′); (2.4)
and the antitime-ordered Green’s function
Gt̄
f (t− t′) = −i
〈
Tt̄f(t)f †(t′)
〉
= Θ(t′ − t)G>
f (t− t′) + Θ(t− t′)G<
f (t− t′). (2.5)
Here |p〉 and |q〉 denote many-body eigenstates of the local Hamiltonian with eigenvalues εp and
εq, respectively, and Θ(t) is the ordinary Heaviside function. We are also interested in the retarded
and advanced Green’s functions,
Gr
f (t− t′) = −iΘ(t− t′)
〈[
f(t), f †(t′)
]
+
〉
= Θ(t− t′)
[
G>
f (t− t′)−G<
f (t− t′)
]
, (2.6)
Ga
f (t− t′) = iΘ(t′ − t)
〈[
f(t), f †(t′)
]
+
〉
= Θ(t′ − t)
[
G<
f (t− t′)−G>
f (t− t′)
]
, (2.7)
which are easily constructed from taking combinations of the lesser and greater Green’s functions
(note the symbol [. . . ]+ denotes the anticommutator of the two operators in the brackets). In
equilibrium, the Green’s functions depend only on the time difference of the real-time variables
and the greater and lesser Green’s functions also satisfy
[
G>
f (t− t′)
]∗
= −G>
f (t′ − t),
[
G<
f (t− t′)
]∗
= −G<
f (t′ − t). (2.8)
It is often useful to represent quantities in terms of frequencies, instead of time. The Fourier
transforms of the greater and lesser Green’s functions (to frequency space) are equal to:
G>
f (ω) = −2πi[1− f(ω)]Af (ω), G<
f (ω) = 2πif(ω)Af (ω), (2.9)
where
Af (ω) =
1
Z
∑
pq
[e−βεp + e−βεq ] 〈p |f | q〉
〈
q
∣∣f †
∣∣ p
〉
δ(ω + εp − εq) (2.10)
427
A.M.Shvaika, J.K.Freericks
is the f -electron spectral density. The corresponding Fourier transforms of the retarded and time-
ordered Green’s functions are as follows:
Gr
f (ω) =
+∞∫
−∞
dω′ Af (ω′)
ω − ω′ + iδ
, ImGr
f (ω) = −πAf (ω), (2.11)
and
Gt
f (ω) = ReGr
f (ω) + i ImGr
f (ω) tanh
βω
2
. (2.12)
We use the Kadanoff-Baym-Keldysh approach to calculate the real-time Green’s function be-
cause there is no simple way of performing the analytical continuation of the Matsubara frequency
Green’s function to the real axis [7]. We are interested in the retarded Green’s function, so we
consider the greater Green’s function for positive time t > 0
G>
f (t) = −i
1
Z Tr
{
e−βHlocTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t
′, t′′)d(t′′)
]
f(t)f †(0)
}
(2.13)
and the lesser Green’s function for negative time t < 0 (or t̄ = −t > 0)
G<
f (t) = i
1
Z Tr
{
e−βHlocTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t
′, t′′)d(t′′)
]
f †(t̄)f(0)
}
; (2.14)
we use the fact that f †(0)f(t) can be replaced by f †(t̄)f(0) in the trace due to the time-translation
invariance that we have in equilibrium. The greater and lesser Green’s functions can be found for
other time values by using the relations in equation (2.8).
The first step in solving for the Green’s function is to write equations of motion for the f -
operators:
df(t)
dt
= −i [Und(t) +Ef − µ] f(t), (2.15)
df †(t̄)
dt̄
= i [Und(t̄) +Ef − µ] f †(t̄), (2.16)
and substitute their solutions into equations (2.13) and (2.14), which yields the following expres-
sions for the greater Green’s function (t > 0)
G>
f (t) =− i
1
Z Tr
{
e−βHlocTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t
′, t′′)d(t′′)
− i
∫
c
dt′Uc(t, t
′)nd(t′)− i(Ef − µ)t
]
f(0)f †(0)
}
and for the lesser Green’s function (t̄ = −t > 0)
G<
f (t) =i
1
Z Tr
{
e−βHlocTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λ(t′, t′′)d(t′′)
+ i
∫
c
dt′Uc(t̄, t
′)nd(t′) + i(Ef − µ)t̄
]
f †(0)f(0)
}
,
where Uc(t, t
′) is a new time-dependent field, which is nonzero only on the upper branch of the
Kadanoff-Baym-Keldysh contour and is equal to U for t′ ∈ [0, t] and is zero otherwise. It is because
this Uc field is not time-translation invariant that we need to use the Kadanoff-Baym-Keldysh
formalism for the analytic continuation. Now we take into account the fact that the f -particle
occupation number nf = f †f is a conserved quantity of the Hamiltonian; for the greater Green’s
428
F-electron spectral function of the Falicov-Kimball model . . .
function we have a projection onto states without f -particles (nf = 0) while for the lesser Green’s
function we have a projection onto states with one f -particle (nf = 1), which yields:
G>
f (t) = −i
1
Z e−i(Ef−µ)t
× Tr
{
eβµndTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t′, t′′)d(t′′)− i
∫
c
dt′Uc(t, t
′)nd(t′)
]}
(2.17)
and
G<
f (t) = i
1
Z eβ(µ−Ef )ei(Ef−µ)t̄
× Tr
{
eβ(µ−U)ndTc exp
[
−i
∫
c
dt′
∫
c
dt′′d†(t′)λc(t
′, t′′)d(t′′),+i
∫
c
dt′Uc(t̄, t
′)nd(t′)
]}
,
(2.18)
where the trace is now over the d-electrons only.
We next expand the contour-ordered-exponent into a series that we can ultimately sum exactly.
To begin, we need to recall Wick’s theorem, where the expectation value of any number of pairs
of fermionic operators can be expanded in terms of products of the two-operator contractions
�
d(t1)d†(t2) = igα(t1, t2), (2.19)
gα(t1, t2) = −i
〈
Tcd(t1)d†(t2)
〉
α
= ie−iεα(t1−t2) [f(εα)−Θc(t1, t2)]
when the Hamiltonian is quadratic in the fermionic operators. Here we have ε0 = −µ for sites with
nf = 0 (α = 0) and ε1 = U − µ for sites with nf = 1 (α = 1). The traces in equations (2.17) and
(2.18) can then be written in the following diagrammatic expression
(
1 + e−βεα
)
exp
{
+
1
2
+
1
3
+ . . .
}
, (2.20)
where the arrows denote the zero-order Green functions gα(t′, t′′) in equation (2.19) and the wavy
lines denote a generalized time-dependent hopping (λ-field)
λ̃t(t
′, t′′) = λ(t′, t′′) + Uc(t, t
′)δc(t
′, t′′) (2.21)
for the greater Green’s function and we have to replace
Uc(t, t
′)→ −Uc(t̄, t
′) (2.22)
for the lesser Green’s function.
Motivated by approaches that involve the integration over a coupling constant, we modify our
time-dependent field by introducing a dependence on some new parameter that we denote by x
Uc(t, t
′)→ Uc(t, t
′|x), (2.23)
with the following limiting behavior
Uc(t, t
′|0) = 0, Uc(t, t
′|1) = Uc(t, t
′). (2.24)
Next, we take the derivative of the diagrammatic series in equation (2.20) with respect to x and
find:
d
dx
{
. . .
}
= −
∫
c
dt′
dUc(t, t
′|x)
dx
G
(α)
U (t′, t′|x), (2.25)
429
A.M.Shvaika, J.K.Freericks
where we introduced a parameter-dependent Green’s function G
(α)
U (t′, t′′|x), which is the solution
of the following Dyson equation:
G
(α)
U (t′, t′′|x) =gα(t′, t′′) +
∫
c
dt1
∫
c
dt2 gα(t′, t1)λ̃t(t1, t2|x)G
(α)
U (t2, t
′′|x)
=gα(t′, t′′) +
∫
c
dt1
∫
c
dt2 gα(t′, t1)λc(t1, t2)G
(α)
U (t2, t
′′|x)
+
∫
c
dt1 gα(t′, t1)Uc(t, t1|x)G
(α)
U (t1, t
′′|x). (2.26)
This expression holds for the greater Green’s function (α = 0). For the lesser Green’s function, one
has to use the substitution in equation (2.22) and set α = 1. Finally, we find explicit expressions
for the greater
G>
f (t) = −iw0 exp
−i(Ef − µ)t−
1∫
0
dx
∫
c
dt1
dUc(t, t1|x)
dx
G
(0)
U (t1, t1|x)
(2.27)
and lesser
G<
f (t) = iw1 exp
i(Ef − µ)t̄+
1∫
0
dx
∫
c
dt1
dUc(t̄, t1|x)
dx
G
(1)
U (t1, t1|x)
(2.28)
Green functions, where
w0 = 〈1− nf 〉 =
Z0
Z0 + Z1
, w1 = 〈nf 〉 =
Z1
Z0 + Z1
(2.29)
are the average densities of sites without (w0) and with (w1) f -electrons.
We can rewrite equation (2.26) in the following two forms:
G
(α)
U (t′, t′′|x) = gaux
Uα (t′, t′′|x) +
∫
c
dt1
∫
c
dt2 g
aux
Uα (t′, t1|x)λc(t1, t2)G
(α)
U (t2, t
′′|x), (2.30)
or
G
(α)
U (t′, t′′|x) = Gα(t′, t′′) +
∫
c
dt1Gα(t′, t1)Uc(t, t1|x)G
(α)
U (t1, t
′′|x). (2.31)
The first form emphasizes the lambda field as an effective self-energy, while the second emphasizes
the U field as the effective self-energy. In the top equation, we introduced the auxiliary Green’s
function
gaux
Uα (t′, t′′|x) = gα(t′, t′′) +
∫
c
dt1 gα(t′, t1)Uc(t, t1|x)gaux
Uα (t1, t
′′|x) (2.32)
as first defined by Brandt and Urbanek [7], while we introduced the bare time-ordered Green’s
function for the effective medium
Gα(t′, t′′) = gα(t′, t′′) +
∫
c
dt1
∫
c
dt2 gα(t′, t1)λc(t1, t2)Gα(t2, t
′′), (2.33)
in the bottom equation. The bare time-ordered Green’s function for the effective medium can be
determined directly on the real axis via
Gα(t′, t′′) = − i
π
+∞∫
−∞
dω ImGα(ω)e−iω(t′−t′′) [f(ω)−Θc(t
′, t′′)] , (2.34)
with
Gα(ω) =
1
ω + iδ − εα − λ(ω + iδ)
(2.35)
430
F-electron spectral function of the Falicov-Kimball model . . .
being the retarded effective medium on the real frequency axis.
The first representation (2.30) was originally used by Brandt and Urbanek [7] and improved
by Zlatić et al. [9] and Freericks, Turkowski and Zlatić [10] (see for review Freericks and Zlatić [2])
and requires the calculation of continuous matrix operator determinants over the entire contour.
In what follows, we shall instead use the second representation in equation (2.31) because in
equations (2.27) and (2.31) all times in Gα
U (t1, t1|x) are only on the real axis (0 6 t1 6 t) [or
(0 6 t1 6 t̄)] and we can replace the contour integral in equation (2.31) by an integral over the
upper real time branch of the Kadanoff-Baym-Keldysh contour where the time-dependent field
Uc(t, t1|x) is nonzero
Gα
U (t′, t′′|x) = Gα(t′ − t′′) +
+∞∫
−∞
dt1Gα(t′ − t1)U(t, t1|x)Gα
U (t1, t
′′|x). (2.36)
Integrals over the entire Kadanoff-Baym-Keldysh contour survive only in the definition of Gα(t)
in equation (2.33), which can instead be calculated directly on the real axis by picking t and t′ on
the upper branch of the contour in equation (2.34)
Gα(t′ − t′′) = − i
π
+∞∫
−∞
dω ImGα(ω)e−iω(t′−t′′) [f(ω)−Θ(t′ − t′′)] . (2.37)
Also, for further discussions we will need its Fourier transform
Gα(ω) =
+∞∫
−∞
dt eiωtGα(t) = ReGα(ω) + i tanh
βω
2
ImGα(ω). (2.38)
We are now ready to summarize the final formulas for the Green’s function by evaluating all of
the terms in equation (2.31) and then performing the integration over the parameter x. We begin
with the formal matrix solution for Gα
U (t′, t′′|x)
Gα
U (t′, t′′|x) =
∫
c
dt̃
[
‖δc(t1, t2)−Gα(t1, t2)Uc(t, t2|x)‖−1
]
t′ t̃
Gα(t̃, t′′). (2.39)
Here δc(t, t
′) is the Dirac delta function along the contour defined by
∫
c dt′δc(t, t
′)f(t′) = f(t).
After substituting this result into equation (2.27) and then integrating over the parameter x, we
obtain for the greater Green’s function the following results:
G>
f (t) =− iw0 exp
{
− i(Ef − µ)t
−
1∫
0
dx
∫
c
dt′
∫
c
dt̃
[
‖δc(t1, t2)−G0(t1, t2)Uc(t, t2|x)‖−1
]
t′ t̃
G0(t̃, t′)
dUc(t, t
′|x)
dx
}
=− iw0 exp
−i(Ef − µ)t−
1∫
0
dxTrc
[
‖I−G0Uc(t|x)‖−1
G0
dUc(t|x)
dx
]
=− iw0 exp {−i(Ef − µ)t+ ln detc ‖I−G0Uc(t)‖}
=− iw0e−i(Ef−µ)t detc ‖δc(t1, t2)−G0(t1, t2)Uc(t, t2)‖ , (2.40)
where Trc denotes the trace of a continuous matrix operator expressed as a line integral over the
contour and detc is the determinant of the continuous matrix operator defined on the contour.
However, in our case, the time-dependent field Uc(t, t2) (and hence all nondiagonal components)
431
A.M.Shvaika, J.K.Freericks
are nonzero only on the upper real time branch of the contour for 0 6 t2 6 t, which allows us to
reduce the problem to the determinant of a continuous matrix operator over the finite real time
interval:
G>
f (t) = −iw0e−i(Ef−µ)t det[0,t] ‖δ(t1 − t2)−G0(t1 − t2)U‖ . (2.41)
It is the restriction to this finite real-time interval that makes the numerical calculations much
more efficient than other methods used previously. Note that all of the temperature dependence is
in the Green’s function G0, which is the time-ordered Green’s function for the effective medium of
the dynamical mean-field theory.
In a similar fashion, we can derive the expression for the lesser Green’s function (t̄ = −t > 0):
G<
f (t) = iw1ei(Ef−µ)t̄ det[0,t̄] ‖δ(t1 − t2) +G1(t1 − t2)U‖ . (2.42)
Here, the Green’s function G1(t1− t2) is defined in equation (2.37) and is the time-ordered Green’s
function for the effective medium when there is an f -electron on the impurity site.
When we perform numerical calculations to evaluate the continuous matrix determinants, we
first discretize the time contour with a discretization step ∆t, and then discretize the continuous
matrix operator into a discrete matrix, which can be manipulated by standard computational
libraries like LAPACK. The discretized representation for the Green’s functions is then
G>
f (t) = −iw0e−i(Ef−µ)t det[0,t] ‖δij −WiG0(ti − tj)U‖ (2.43)
and
G<
f (t) = iw1ei(Ef−µ)t̄ det[0,t̄] ‖δij +WiG1(ti − tj)U‖ , (2.44)
where Wi are the quadrature weights for the discrete set of times {ti}; for a uniform grid Wi =
W = ∆t, we have to calculate the determinant of a Toeplitz matrix. This will allow us to use
Szegö’s theorem and the Wiener-Hopf approach to find analytic expressions of the continuous
matrix operator determinant which are accurate for long times. Note, that while it appears that
these expressions are equally valid for all temperatures, the numerical solution, in the ∆t → 0
limit, becomes more difficult at lower temperatures.
The retarded Green’s function is found from the greater and lesser Green functions via
Gr
f (t) = Θ(t)
{
G>
f (t)−G<
f (t)
}
. (2.45)
Using the relation in equation (2.8), then yields
Gr
f (t) = −iΘ(t)e−i(Ef−µ)t
{
w0 det[0,t] ‖I−G0U‖+ w1
(
det[0,t] ‖I + G1U‖
)∗}
(2.46)
for real times, and
Gr
f (ω) = −i
+∞∫
0
dtei(ω+µ−Ef )t
{
w0 det[0,t] ‖I−G0U‖+ w1
(
det[0,t] ‖I + G1U‖
)∗}
(2.47)
for real frequencies.
3. The Wiener-Hopf sum equation approach and Szegö’s theorem
We closely follow the development of the Wiener-Hopf sum equation approach in McCoy and
Wu [20], which is based on the original integral equation development by Wiener and Hopf [21],
as summarized in Krein [22]. The Wiener-Hopf integral equation approach has been widely used
in the X-ray edge problem and in the Falicov-Kimball model [8] primarily at T = 0; we will show
that the Wiener-Hopf sum equation approach is very powerful for finite temperature calculations.
Before we can apply this technique to the problem of calculating the f -electron spectral function,
we should begin with some mathematical preliminaries (further details can be found in [20]).
432
F-electron spectral function of the Falicov-Kimball model . . .
Consider a sequence of numbers cn which satisfy
∑∞
n=−∞ |cn| < ∞ and the set of N + 1 linear
equations
N∑
m=0
cn−mx
(N)
m = yn (3.1)
for 0 6 n 6 N with cn and yn known. We will eventually be taking the limit where N →∞, so we
will require
∑∞
n=−∞ |yn| <∞ too. Note that the matrix cn−m that appears in equation (3.1) is a
Toeplitz matrix, because it depends only on the difference n−m and not on n and m separately.
Our job is to find the solution x
(N)
m to these equations. We define
x(N)
n = yn = 0 (3.2)
for n < 0 and n > N . In addition, we need to define two more sets of coefficients:
un =
{ ∑N
m=0 cN+n−mx
(N)
m for n > 0
0 for n 6 0
and
vn =
{ ∑N
m=0 c−n−mx
(N)
m for n > 0
0 for n 6 0
.
Next, we relax the condition on equation (3.1) to apply it for all n and we allow the sum over m
to extend from −∞ to ∞ (although all cases have just a finite number of terms because x
(N)
m is
potentially nonzero for only N + 1 terms). The result is
+∞∑
m=−∞
cn−mx
(N)
m = yn + Θ(n > N)un−N + Θ(n < 0)v−n . (3.3)
We consider a variable ξ = exp[iθ] which lies on the unit circle, so that |ξ| = 1. We multiply both
sides of equation (3.3) by ξn and sum over n to find
C(ξ)X(ξ) = Y (ξ) + U(ξ)ξN + V (ξ−1), (3.4)
with
C(ξ) =
+∞∑
n=−∞
cnξ
n, X(ξ) =
N∑
n=0
x(N)
n ξn, Y (ξ) =
N∑
n=0
ynξ
n,
U(ξ) =
∞∑
n=1
unξ
n, and V (ξ) =
∞∑
n=1
vnξ
n. (3.5)
It may appear that we have made the problem more complicated, because we have replaced a
problem with one unknown X(ξ) by a problem with three unknowns X(ξ), U(ξ) and V (ξ). But
it turns out that this more complicated problem can actually be solved by carefully analyzing its
analytic structure.
The crucial step of Wiener and Hopf is to find a unique factorization of C(ξ) into two factors,
one analytic inside the unit circle and continuous on the circle, and the other analytic outside
the unit circle and continuous on the circle. Before we show how to do this, we need some more
mathematical definitions. We call a function a + function if it can be expanded as a Laurent series
for |ξ| = 1 (
∑∞
n=0 anξ
n) and the coefficients satisfy
∑∞
n=0 |an| < ∞. In this case, the function is
analytic within the unit circle and continuous on the unit circle. Similarly, we call a function a −
function if it can be expanded in a Laurent series for |ξ| = 1 (
∑−1
n=−∞ anξ
n) and the coefficients also
satisfy
∑−1
n=−∞ |an| < ∞. This function is analytic outside the unit circle and continuous on the
unit circle; it vanishes as |ξ| → ∞. Note that any function defined on the unit circle via a Laurent
expansion f(ξ) =
∑∞
n=−∞ anξ
n with
∑∞
n=−∞ |an| < ∞ can obviously be uniquely decomposed
into the sum of a + and a − functions via f(ξ) = f+(ξ) + f−(ξ).
433
A.M.Shvaika, J.K.Freericks
The function C(ξ) is factorized in a so-called canonical factorization, where
C(ξ) = P (ξ)−1Q(ξ−1)−1, (3.6)
and both P (ξ) and Q(ξ) are nonzero for |ξ| 6 1; the factorization is made unique by requiring
Q(0) = 1. The functions P and Q can then be found in terms of + and − functions as
P (ξ) = eG+(ξ) and Q(ξ−1) = eG−(ξ), (3.7)
which automatically satisfies G−(∞) = 0 [and Q(0) = 1] and makes both P (ξ) and Q(ξ) nonvan-
ishing inside the unit circle. An explicit calculation of G±(ξ) is nontrivial, but it is easy to write
down a formal solution via
G+(ξ) = − [lnC(ξ)]+ and G−(ξ) = − [lnC(ξ)]− , (3.8)
where the ± subscripts denote a restriction of the Laurent series for the logarithm to its nonnegative
or negative power terms, respectively. Note that the condition for properly defining the + and −
functions requires the logarithm to be single-valued after ξ winds around the unit circle, or, in
other words, we require IndC(ξ) = 0, where
IndC(ξ) =
1
2iπ
[
lnC(e2iπ)− lnC(ei0)
]
. (3.9)
We will also need to consider situations where IndC(ξ) = −1. In this case, a new function C̄(ξ) =
−ξC(ξ) will have zero index, and the canonical factorization can be applied to C̄(ξ) instead.
We are now ready to solve equation (3.4) by substituting in the factorization from equation (3.6).
Next, multiply both sides of the equation by Q(ξ−1) to get
[P (ξ)]−1X(ξ) = Q(ξ−1)Y (ξ) +Q(ξ−1)U(ξ)ξN +Q(ξ−1)V (ξ−1). (3.10)
The term on the left hand side is a + function, because P (0) 6= 0, and [P (ξ)]−1 can be expanded in a
power series with just positive powers for |ξ| < 1 and X(ξ) is a + function. Similarly, Q(ξ−1)V (ξ−1)
is a − function. The other two terms are neither + nor − functions, but they can be uniquely
decomposed into + and − function pieces. We do this, and move all + functions to the left hand
side and all − function pieces to the right hand side. We are left with
[P (ξ)]−1X(ξ)−
[
Q(ξ−1)Y (ξ)
]
+
−
[
Q(ξ−1)U(ξ)ξN
]
+
=
[
Q(ξ−1)Y (ξ)
]
−
+
[
Q(ξ−1)U(ξ)ξN
]
−
+Q(ξ−1)V (ξ−1). (3.11)
The left hand side of equation (3.11) is an analytic function for |ξ| < 1 and is continuous on the
unit circle, the right hand side is an analytic function for |ξ| > 1 and is continuous on the unit
circle, and both functions agree on the unit circle, due to the equality, so we can define a function
that is analytic in the entire complex plane, and hence is an entire function. But this function
vanishes as |ξ| → ∞, and the only entire function that vanishes for large argument is the function
that is identically equal to 0. Hence we learn that
X(ξ) = P (ξ)
[
Q(ξ−1)Y (ξ)
]
+
+ P (ξ)
[
Q(ξ−1)U(ξ)ξN
]
+
(3.12)
and
V (ξ−1) = −
[
Q(ξ−1)
]−1
{[
Q(ξ−1)Y (ξ)
]
−
+
[
Q(ξ−1)U(ξ)ξN
]
−
}
. (3.13)
Using the fact that X(ξ−1)ξN is a + function [recall X(ξ) is an order N polynomial in ξ], we
can substitute ξ → ξ−1 in equation (3.11), multiply both sides of the equation by ξN , and then
separate into the + and − function pieces to create another vanishing entire function, and learn
that both
X(ξ−1)ξN = Q(ξ)
{[
P (ξ−1)Y (ξ−1)ξN
]
+
+
[
P (ξ−1)V (ξ)ξN
]
+
}
(3.14)
434
F-electron spectral function of the Falicov-Kimball model . . .
and
U(ξ−1) = −
[
P (ξ−1)
]−1
{[
P (ξ−1)Y (ξ−1)ξN
]
−
+
[
P (ξ−1)V (ξ)ξN
]
−
}
(3.15)
hold. These are all of the necessary relations we will need for determining the Toeplitz determinants.
Now we should take another mathematical interlude to state Szegö’s theorem. This theorem
determines how a Toeplitz determinant exponentially decays for large matrix size [23] and even
finds a constant prefactor for the exponential decay [24,25]. The proof can be carried out using
the Wiener-Hopf sum equation approach [20], but it would take too much space for us to repeat
it here, and it is rife with many technical mathematical manipulations that would take us away
from the physical phenomena we wish to discuss here. So we will instead simply state the result,
and then show how to evaluate asymptotic expressions for the f -electron Green’s function in real
time.
If we consider the determinant DN of the N ×N Toeplitz matrix defined by
DN =
∣∣∣∣∣∣∣∣∣∣∣∣∣
c0 c−1 c−2 . . . c−N+1
c1 c0 c−1 . . . c−N+2
c2 c1 c0 . . . c−N+3
...
...
...
...
...
cN−2 cN−3 cN−4 . . . c−1
cN−1 cN−2 cN−3 . . . c0
∣∣∣∣∣∣∣∣∣∣∣∣∣
,
then Szegö’s theorem says that
lim
N→∞
DN = exp
[
Ng0 +
∞∑
n=1
ngng−n
]
, (3.16)
with
gn =
1
2π
π∫
−π
dθe−inθ lnC(eiθ). (3.17)
The conditions that IndC(ξ) = 0, C(ξ) is continuous on the unit circle, and
∑∞
n=−∞ |cn| <∞ are
all required for the theorem to hold.
Now, with all of the mathematical preliminaries complete, we are ready to apply these tech-
niques to finding the Toeplitz determinants needed for the f -electron Green’s function. We need
to evaluate the two determinants in equations (2.43) and (2.44). We will explicitly show how to
calculate the determinant for the greater Green’s function. Modifications for the lesser Green’s
function are straightforward.
Case 1: no winding. We start by discretizing the time axis [0, t] with N + 1 points given by
tn = n∆t, with ∆t = t/N . Then we define coefficients cn to satisfy
cn = δn0 −∆tUG0(n∆t). (3.18)
Since G0(t) is bounded, and decays exponentially fast in time, it is easy to see that
∑∞
n=−∞ |cn| <
∞. We will assume for the moment that IndC(ξ) = 0 (which turns out to be true at half filling
when U . 0.866 on the hypercubic lattice, see below). Using the definition for cn, we find
C(eiθ) =
∞∑
n=−∞
cneinθ = 1−∆tU
∞∑
n=−∞
G0(n∆t)einθ. (3.19)
If we let n∆t = t′ and ω∆t = θ, we recognize that in the limit where ∆t→ 0, the sum in the right
hand side of the above equation approaches the Fourier transform of G0(t′), so we write
C(ω) = 1− U
∞∫
−∞
dt′G0(t′)eiωt′ = 1− UG0(ω), (3.20)
435
A.M.Shvaika, J.K.Freericks
where G0(ω) is defined by equation (2.38). The coefficients gn in equation (3.17) can then be
written as
g(t′) =
gn
∆t
=
1
2π
π/∆t∫
−π/∆t
dωe−in∆tω lnC(ei∆tω) =
1
2π
∞∫
−∞
dωe−iωt′ lnC(ω), (3.21)
where we again used t′ = n∆t. Now, defining DN to satisfy
D(t) = D(N∆t) = DN = det[0,t] ‖δij −∆tG0(ti − tj)U‖ , (3.22)
we find from equation (3.16) that
lim
t→∞
D(t) = exp
t
2π
∞∫
−∞
dω lnC(ω) +
∞∫
0
dt′t′g(t′)g(−t′)
, (3.23)
which is the asymptotically exact result in the limit of large t when the index of C(ξ) is equal
to zero. Since the coefficient of t in the exponent can be a complex number, the determinant can
oscillate while exponentially decaying at long times. This occurs when the system has undergone
the Mott transition and the f -electron spectral function displays a gap at low frequency.
However, we can improve upon the result in equation (3.23) to provide finite-time corrections.
Now we illustrate how to do this using the Wiener-Hopf sum equation approach. We start with our
original N + 1 simultaneous equations that we need to solve in equation (3.1), with the coefficients
cn defined in equation (3.18). Next, we choose yn = δn0. Cramers rule, applied to the first column
of the c matrix immediately tells us that
x
(N)
0 =
DN
DN+1
, (3.24)
which further implies that
DN =
(
∞∏
M=N
x
(M)
0
)
×DM→∞ . (3.25)
If we assume that x
(M)
0 = (1−∆tδM ) exp(−g0) (which we will verify below), then we find
D(t) = exp
[
−∆t
∞∑
M=N
δM
]
exp
t
2π
∞∫
−∞
dω lnC(ω) +
∞∫
0
dt′t′g(t′)g(−t′)
. (3.26)
This will provide the next order of corrections to the asymptotic limit of the determinant.
We need to find x
(M)
0 for large M . Our strategy is to set U(ξ) = 0 in equation (3.13) to
find V (ξ−1), then solve for U(ξ−1) with equation (3.15) and the approximate V (ξ), and finally
substitute U(ξ) into equation (3.12) to find X(ξ). The term x
(M)
0 follows from taking the ξ → 0
limit. Using the facts that Y (ξ) = 1 for our equation and Q(ξ−1)−1 is a − function, we immediately
find that equation (3.13) can be solved by
V (ξ−1) = −1 +
[
Q(ξ−1)
]−1
. (3.27)
Substituting into equation (3.15) (after replacing ξ → ξ−1), yields
U(ξ−1) = −
[
P (ξ−1)
]−1
[
P (ξ−1) {Q(ξ)}−1
ξN
]
−
. (3.28)
Now we should replace ξ → ξ−1 in the above equation. This requires care because we will convert
the − function into a + function, but with no constant term in the power series expansion. So we
write
U(ξ) = − [P (ξ)]−1
[
P (ξ)
{
Q(ξ−1)
}−1
ξ−N
]′
+
, (3.29)
436
F-electron spectral function of the Falicov-Kimball model . . .
with the prime indicating there is no constant term in the + function. This result can now be
substituted into equation (3.12) to yield
X(ξ) = P (ξ)− P (ξ)
[
Q(ξ−1) {P (ξ)}−1 ξN
[
P (ξ)
{
Q(ξ−1)
}−1
ξ−N
]′
+
]
+
, (3.30)
where we used the fact that [Q(ξ−1)]+ = 1. We will need to evaluate this in the limit ξ → 0.
Note that we have the following two identities for P and Q
P (ξ) = exp
[
−
∞∑
n=0
gnξ
n
]
and P (ω) = exp
−
∞∫
0
dtg(t)eiωt
, (3.31)
and
Q(ξ−1) = exp
[
−
−1∑
n=−∞
gnξ
n
]
and Q(−ω) = exp
−
−∆t∫
−∞
dtg(t)eiωt
. (3.32)
Using the definition for g(t) in equation (3.21), one can immediately verify that P (ω)Q(−ω) =
1/C(ω) in the limit as ∆t→ 0, as it should. Using the first relation in equation (3.31), immediately
shows us that P (ξ = 0) = exp[−g0]. Evaluating the other term in equation (3.30) requires more
work. In order to create the relevant + functions, we should first take the functions of ω, Fourier
transform them to functions of t and then reverse Fourier transform back, but include only the
positive time contributions. For example, we can write
[
P (ξ)
{
Q(ξ−1)
}−1
ξ−N
]′
+
=
1
2π
∞∫
∆t
dt′eiωt′
∞∫
−∞
dω′ P (ω′)
Q(−ω′)
e−iω′(t+t′). (3.33)
The limit of ∆t is to remind us that there should be no constant term in the expansion. Using
this strategy, we can immediately write down the final answer for x
(N)
0 . We obtain x
(N)
0 = (1 −
∆tδN ) exp[−g0] with
δN =
1
4π2
∞∫
−∞
dω
∞∫
∆t
dt′
∞∫
−∞
dω′Q(−ω)
P (ω)
P (ω′)
Q(−ω′)
ei(ω−ω′)(t+t′), (3.34)
where t = N∆t. Putting this all together finally yields the asymptotic expression for D(t)
D(t) = exp
t
2π
∞∫
−∞
dω lnC(ω) +
∞∫
0
dt′t′g(t′)g(−t′)−
∞∫
t
dt̄
∞∫
∆t
dt′h(t+ t′)h′(−t− t′)
, (3.35)
with
h(t) =
1
2π
∞∫
−∞
dω
P (ω)
Q(−ω)
e−iωt and h′(t) =
1
2π
∞∫
−∞
dω
Q(−ω)
P (ω)
e−iωt. (3.36)
Case 2: winding of –1. Szegö’s theorem fails when the index of C(ξ) is nonzero, and lnC(ξ)
winds around the origin. This occurs on the hypercubic lattice at half filling when U & 0.866, see
below. We can still derive asymptotic formulas for the Green’s function by solving an auxiliary
problem, which has the winding removed. We now show how this is done.
In all numerical cases we have examined, the index of C(ξ) satisfies IndC(ξ) = −1 when U is
large enough. Such a winding can be removed by considering C̄(ξ) = −ξC(ξ) which has index zero
(the minus sign is introduced for convenience, as will be clear below). The coefficients obviously
satisfy
c̄n = −cn−1. (3.37)
437
A.M.Shvaika, J.K.Freericks
We examine the auxiliary Toeplitz determinant
D̄N =
∣∣∣∣∣∣∣∣∣∣∣
c̄0 c̄−1 c̄−2 . . . c̄−N+1
c̄1 c̄0 c̄−1 . . . c̄−N+2
c̄2 c̄1 c̄0 . . . c̄−N+3
...
...
...
...
...
c̄N−1 c̄N−2 c̄N−3 . . . c̄0
∣∣∣∣∣∣∣∣∣∣∣
=
∣∣∣∣∣∣∣∣∣∣∣
−c−1 −c−2 −c−3 . . . −c−N
−c0 −c−1 −c−2 . . . −c−N+1
−c1 −c0 −c−1 . . . −c−N+2
...
...
...
...
...
−cN−2 −cN−3 −cN−4 . . . −c−1
∣∣∣∣∣∣∣∣∣∣∣
,
and consider the N + 1 simultaneous equations in equation (3.1) with the barred variables
N∑
m=0
c̄n−mx̄
(N)
m = ȳn (3.38)
with ȳn = δn0. Then Cramers rule, evaluated for the N + 1st column, tells us
x̄
(N)
N = DN−1/D̄N , (3.39)
where the two factors of (−1)N cancel. Szegö’s theorem can be applied to C̄(ξ), so we immediately
learn that
lim
N→∞
D̄N = exp
[
Nḡ0 +
∞∑
n=1
nḡnḡ−n
]
, (3.40)
with
ḡn =
1
2π
π∫
−π
dθe−inθ ln C̄(eiθ). (3.41)
We need to calculate x̄
(N)
N in order to find the determinant DN .
First note that equations (3.12–3.15) hold for the barred functions when we have a winding
around the origin. We will once again solve for V̄ (ξ−1) setting Ū(ξ) = 0, then substitute into
equation (3.14) to find X̄(ξ−1)ξN . Taking the limit ξ → 0 will then give us x̄
(N)
N . Just like before,
we find
V̄ (ξ) = −1 +
[
Q̄(ξ)
]−1
, (3.42)
which then gives
X̄(ξ−1)ξN = Q̄(ξ)
[{
Q̄(ξ)
}−1
P̄ (ξ−1)ξN
]
+
. (3.43)
We take the limit ξ → 0 by integrating the function on the right hand side over the unit circle and
dividing by 2π. Using θ = ∆tω, then yields
x̄
(N)
N =
∆t
2π
π/∆t∫
−π/∆t
dωeiωtP̄ (−ω)/Q̄(ω), (3.44)
where we used t = N∆t again. Note that the limits on the integration cannot be simply extended
to ±∞ as we did before. Since the function lnC(exp[iθ]) winds once around the origin, it displays a
discontinuity of −2iπ to its imaginary part as θ runs from −π to π. We remove this discontinuity by
adding a linear function that is equal to 0 at θ = −π and is equal to 2iπ at θ = π (the minus sign in
C̄ is needed to move the branch cut of the logarithm from the positive real axis to the negative real
axis). Since this linear function cannot be extended to infinity, we need to work with a finite value
of ∆t in the calculations (for a further discussion, look at figure 2 and the corresponding discussion
below). We simply take ∆t small enough, that we see the numerical results do not change, and the
system has approached its ∆t→ 0 limit. The final result for the determinant is then
D(t) = exp
t
2π
∞∫
−∞
dω ln C̄(ω) +
∞∫
0
dt′t′ḡ(t′)ḡ(−t′)
∆t
2π
π/∆t∫
−π/∆t
dω′eiω′tP̄ (−ω′)/Q̄(ω′). (3.45)
438
F-electron spectral function of the Falicov-Kimball model . . .
The limits on the first integral have been extended to ±∞, because it turns out that the linear
piece gives no contribution to the integral. The definitions of ḡ(t), P̄ (ω) and Q̄(ω) are obvious from
the above discussion.
It should be noted that, in contrast to the case without winding where the third term in the
exponent of equation (3.35) gives finite time corrections to the large t exponential asymptotics of
equation (3.23), in the case with winding, equation (3.45) has a time dependent prefactor with a
nontrivial ∆t → 0 and t → ∞ limit which cannot be considered as a vanishing correction to the
exponential decay of equation (3.23).
4. Numerical results
We begin our numerical discussion by describing the different behavior of C(ξ) when it has no
winding or when it winds with an index of −1 around the origin. In figure 2, we show a plot of
lnC(ξ) for two cases: (a) the case with no winding, and (b) the case with an index of −1. These
two cases are distinguished by the sign of the real part of C(ξ) at ξ = 0, where its imaginary part
changes sign. On the hypercubic lattice at half filling we have ReC(0) > 0 for U < 0.866 and there
is no winding and ReC(0) < 0 for U > 0.866 and now the index is equal to −1. Notice how the
imaginary part has a steep change in its value near ξ = 0 [panel (a)], which evolves into a discrete
jump by 2πi when U & 0.866 [panel (b)], the jump at the origin is compensated by the linear
shift which runs from the minimal to maximal values of the frequency used in the calculations
(|ω| = 10π here). While it is obvious that integration over ω can easily be extended to ω = ±∞
when there is no winding, one needs to carefully include contributions present at the given value
of ∆t when there is winding, and one cannot extend the integration limits in this case.
-20 -10 0 10 20
ω
-0.5
0
0.5
ln
C
(ω
)
Re ln C(ω)
Im ln C(ω)
(a)
-30 -20 -10 0 10 20 30
ω
-3
-2
-1
0
1
2
3
ln
C
(ω
)
Re ln C(ω)=Re ln C(ω)
Im ln C(ω)
Im ln C(ω)
(b)
Figure 2. (Color online) Plot of (a) ln C(ω) for the case with no winding [U = 0.5, IndC(ξ) = 0]
and (b) ln C(ω) and ln C̄(ω) for the case with winding [U = 1.5, ∆t = 0.1, IndC(ξ) = −1].
Next we compare our different asymptotic expansions to the results derived from a direct
numerical calculation with the matrix formalism from equation (2.43) when T = 0.1. This involves
a straightforward approach where we use three different discretization sizes (∆t = 0.1, 0.0666,
and 0.0333) which we extrapolate to ∆t → 0, and we compare those results to the calculations
with different approximations. In the case with no winding, we use both the asymptotic form for
large t in equation (3.23) and the more correct form for smaller t in equation (3.35). In figure 3,
we show the results for the greater Green’s function at half filling on the hypercubic lattice for
U = 0.5 [panel (a)] and U = 0.7 [panel (b)]. Note how the exponential form, that comes from
Szegö’s theorem, is quite accurate for large times, but becomes poor for small times, and how
the corrections arising from the Wiener-Hopf approach produce remarkable agreement with the
numerically exact results for essentially all t. The errors are the largest near t = 0, but even there,
they lie below the few percent level. This shows that the analytic formulas are quite accurate for
all t when there is no winding.
Next we examine what happens in the case with winding. We examine two values of U in
439
A.M.Shvaika, J.K.Freericks
0 10 20 30 40
t
-0.5
-0.4
-0.3
-0.2
-0.1
0
G
f> (t
)
matrix determinant
asymptotic
finite-time
Real part
Imaginary part
(a)
0 10 20 30 40 50
t
-0.5
-0.4
-0.3
-0.2
-0.1
0
G
f> (t
)
matrix determinant
asymptotic
finite time
Real part
Imaginary part
(b)
Figure 3. (Color online) Real and imaginary parts of the greater Green’s function for cases with
no winding (a) U = 0.5 and (b) U = 0.7 (T = 0.1). The solid (black) lines represent exact
results of the functional determinant in equation (2.43), the dot-dashed (red) line represents
the asymptotic result in equation (3.23), and the dashed (blue) line represents the expression
with finite-time corrections in equation (3.35). Note that we cannot really discern any difference
between the exact matrix results and the asymptotic Wiener-Hopf results all the way down to
t = 0.
figure 4: U = 1.5 which is near the critical interaction for the metal-insulator transition Uc =
√
2,
and U = 2 which is a small gap insulator. We show the results from the scaled matrix calculations
with a solid line (black) and the asymptotic Wiener-Hopf results with a (blue) dashed line. In this
case, the Green’s function shows a behavior that does not approach the asymptotic exponential
behavior for short times, and the analytic formula is less accurate for small times (with maximal
error of the order of 10%).
0 10 20 30 40
t
-0.5
-0.4
-0.3
-0.2
-0.1
0
G
f> (t
)
determinant
WH approximation
Imaginary part
Real part
(a)
0 10 20 30 40
t
-0.5
-0.4
-0.3
-0.2
-0.1
0
0.1
G
f> (t
)
determinant
WH approximation
Real part
Imaginary part (b)
Figure 4. (Color online) Real and imaginary parts of the greater Green’s function for the case
with winding [IndC(ξ) = −1] (a) U = 1.5 and (b) U = 2 (T = 0.1). The solid lines represent the
extrapolated matrix results for the functional determinant in equation (2.43) and the dashed
line represents the asymptotic expression with the finite-time corrections in equation (3.45).
Now we can specify the differences between the cases without winding and with winding.
For U < 0.866 there is no winding and our exact results of the functional determinants rapidly
approach the asymptotic result (with finite-time corrections) in equation (3.35) and we do not
observe a crossing of the zero axis at any temperature and at T → 0 the long-time behavior is
replaced by a power law. On the other hand, for 0.866 < U < Uc we always observe a crossing of the
axis at high temperatures, which originates from the time dependent prefactor in equation (3.45).
With decreasing temperature, the crossing point shifts to larger times producing some kind of
exponential decay for intermediate values of t and at T = 0 we observe a crossover to power law
behavior when the zero crossing is pushed out to infinity. For U > Uc the crossing of the axis is
440
F-electron spectral function of the Falicov-Kimball model . . .
observed at finite time values for any temperature and does not transform into power law at zero
temperature.
The Fourier transforms to the spectral formula produces the results essentially equivalent to
those already shown in [10–12], so we do not repeat them here.
5. Conclusions
In this work, we have shown a new representation for the f -electron Green’s function, which
allows for efficient numerical computation. By using the property that Uc(t, t
′) is nonzero only on
the interval [0, t], we restrict the determinant of the continuous matrix operators to a finite time
interval instead of over the three branches of the Kadanoff-Baym-Keldysh contour. This produces
huge savings in the computational effort, because the size of the matrices does not grow with
temperature and they are less than one half the size of the matrices used in previous numerical
calculations. As a result, we can examine much wider regions of parameter space.
In addition, we used the Wiener-Hopf sum equation approach, along with Szegö’s theorem,
to derive exact analytical expressions for the Toeplitz determinants required to find the Green’s
functions. While these expressions are asymptotically exact in the limit where t→∞, we find that
they have errors less than 10% all the way down to t = 0. This allows us to use the matrix results
for small times, and then append them by the asymptotic expressions for large times in order to
find accurate results for the Green’s function at all times. Then they are Fourier transformed to
real frequencies to determine the f -electron spectral function.
The f -electron spectral function displays the expected behavior as well. For U values smaller
than the critical U for the Mott transition (U =
√
2), the spectral function develops a sharp
peak at low T , which ultimately diverges as an inverse power law in ω as T → 0 [rigorously
speaking, our asymptotic formulas are not valid at T = 0, because the function lnC(ξ) becomes
discontinuous due to the jump in the Fermi factor, but the power law behavior can be extracted
via other techniques—our focus here was on the finite-T behavior]. When U is larger than the
critical U , a gap forms in the f -electron spectral function (rigorously speaking, it is a pseudogap
on the hypercubic lattice) but subgap states rapidly enter as a function of T for low T .
Acknowledgements
It is with great pleasure that we honor Prof. Stasyuk with this contribution. This publica-
tion is based on work supported by Award No. UKP2–2697–LV–06 of the U.S. Civilian Research
& Development Foundation (CRDF). J.K.F. also acknowledges support by the National Science
Foundation under grant number DMR–0705266.
References
1. Falicov L.M., Kimball J.C., Phys. Rev. Lett., 1969, 22, 997.
2. Freericks J.K., Zlatić V., Rev. Mod. Phys., 2003, 75, 1333.
3. Stasyuk I.V., Shvaika A.M., J. Phys. Studies, 1999, 3, 177.
4. Stasyuk I.V., Hera O.B., Condens. Matter Phys., 2003, 6, 127.
5. Stasyuk I.V., Hera O.B., Phys. Rev. B, 2005, 72, 045134.
6. Stasyuk I.V., Hera O.B., Eur. Phys. J. B, 2005, 48, 339.
7. Brandt U., Urbanek M.P., Z. Phys. B: Condens. Matter, 1992, 89, 297.
8. Janis V., Phys. Rev. B, 1994, 49, 1612; Physica B, 1996, 223–224, 616; Int. J. Mod. Phys. B, 1997
11, 433.
9. Zlatić V., Freericks J.K., Lemański R., Czycholl G, Philos. Mag. B, 2001, 81, 1443.
10. Freericks J.K., Turkowksi V.M., Zlatić V., Phys. Rev. B, 2005, 71, 115111.
11. Freericks J.K., Turkowksi V.M., Zlatić V., in Proceedings of the HPCMP Users Group Conference
2004, Williamsburg, VA, June 7–11, 2004 edited by R. E. Peterkin, IEEE Computer Society, Los
Alamitos, CA, 2004, p. 7.
12. Freericks J.K., Turkowksi V.M., Zlatić V., Physica B, 2005 359–361, 684.
441
A.M.Shvaika, J.K.Freericks
13. Anders F.B., Czycholl G., Phys. Rev. B, 2005, 71, 125101.
14. Liu Y.-L., Phys. Rev. B, 2005, 72, 045123.
15. Möller G., Ruckenstein A.E., Schmitt-Rink S., Phys. Rev. B, 1992, 46, 7427.
16. Mahan G.D., Phys. Rev., 1967, 163, 612.
17. Nozières P., De Dominicis C., Phys. Rev., 1969, 178, 1097.
18. Kadanoff L.P., Baym G., Quantum statistical mechanics, Benjamin, New York, 1962.
19. Keldysh L.V., Zh. Eksp. Teor. Fiz., 1964, 47, 1945 (in Russian) [Sov. Phys. JETP, 1965, 20, 1018].
20. McCoy B.M.,Wu T.T., The two-dimensional Ising model, Harvard University Press, Cambridge, MA,
1973.
21. Wiener N., Hopf E., Sitzber. Preuss. Akad. Wiss. Berlin, Kl. Math. Phys. Tech., 1931, 696.
22. Krein M.G., Uspekhi Mat. Nauk, 1958, 13, 3 (in Russian) [Am. Math. Soc. Transl., 1962, 22, 163].
23. Szegö G., Math. Z., 1920, 6, 167; Math. Z., 1921, 9, 167.
24. Szegö G., Commun. Sem. Math. Univ. Lund, suppl. ded. Marcel Reisz, 1952, 228.
25. Grenander U., Szegö G., Toeplitz forms and their applications, University of California Press, Berkeley
and Los Angeles, 1958.
Спектральна функцiя f -електронiв для моделi
Фалiкова-Кiмбала i дискретний пiдхiд Вiнера-Гопфа
А.М.Швайка1, Дж.К.Фрiрiкс2
1 Iнститут фiзики конденсованих систем НАН України, вул. Свєнцiцького 1, 79011 Львiв
2 Фiзичний факультет, Унiверситет Джорджтауну, Вашингтон, округ Колумбiя 20057, США
Отримано 21 травня 2008 р., в остаточному виглядi – 20 червня 2008 р.
Отримано нове представлення для спектральної функцiї f -електронiв моделi Фалiкова-Кiмбала, яке
альтернативне оригiнальному представленню Брандта i Урбанека. У новому представленнi усi роз-
рахунки виконуються тiльки на дiйснiй часовiй осi, що дозволяє розглядати як завгодно низькi тем-
ператури. Загальний вираз для запiзнюючої функцiї Ґрiна включає два детермiнанти неперервних
матричних операторiв зi структурою типу Теплiца. Застосовуючи дискретний пiдхiд Вiнера-Гопфа i
теорему Сеґо, отримано точнi аналiтичнi формули для довгочасової поведiнки функцiй Ґрiна; розгля-
нуто випадки, коли логарифм характеристичної функцiї (яка визначає неперервну матрицю Теплiца)
робить i не робить виток навколо початку координат. Показано наскiльки точними є данi асимпто-
тичнi вирази у порiвняннi з точними розв’язками, якi отримуються при екстраполяцiї прямих матри-
чних розрахункiв до границi нульової дискретизацiї.
Ключовi слова: спектральна функцiя f -електронiв, модель Фалiкова-Кiмбала, пiдхiд Вiнера-Гопфа,
теорiя динамiчного середнього поля
PACS: 71.10.-w, 71.27.+a, 71.30.+h, 02.30.Rz
442
|