arXiv:cond-mat/0403517v2 [cond-mat.str-el] 31 Aug 2004

Functional renormalization group for Luttinger liquids
with impurities
S. Andergassen,1 T. Enss,1 V. Meden,2 W. Metzner,1
U. Schollw¨ck,3 and K. Sch¨nhammer2
o
o
1 Max-Planck-Institut
2 Institut
3 Institut

f¨r Festk¨rperforschung, D-70569 Stuttgart, Germany
u
o

f¨r Theoretische Physik, Universit¨t G¨ttingen, D-37077 G¨ttingen, Germany
u
a o
o

f¨r Theoretische Physik, Technische Hochschule Aachen, D-52056 Aachen, Germany
u
(Dated: March 19, 2004)

We improve the recently developed functional renormalization group (fRG) for
impurities and boundaries in Luttinger liquids by including renormalization of the
two-particle interaction, in addition to renormalization of the impurity potential.
Explicit flow equations are derived for spinless lattice fermions with nearest neighbor interaction at zero temperature, and a fast algorithm for solving these equations
for very large systems is presented. We compute spectral properties of single-particle
excitations, and the oscillations in the density profile induced by impurities or boundaries for chains with up to 106 lattice sites. The expected asymptotic power-laws at
low energy or long distance are fully captured by the fRG. Results on the relevant
energy scales and crossover phenomena at intermediate scales are also obtained. A
comparison with numerical density matrix renormalization results for systems with
up to 1000 sites shows that the fRG with the inclusion of vertex renormalization is
remarkably accurate even for intermediate interaction strengths.
PACS: 71.10.Pm, 72.10.-d, 73.21.Hb

I.

INTRODUCTION

Low dimensional electron systems exhibit a rich variety of surprising effects which
are due to the cooperative interplay of impurities and interactions. In one dimension
even clean metallic systems are always strongly affected by interactions: at low energy
scales physical properties obey anomalous power-laws, known as Luttinger liquid behavior, which is very different from the conventional Fermi liquid behavior describing most
higher dimensional metals.1,2 In Luttinger liquids with repulsive interactions already a single static impurity is known to affect the low-energy properties drastically.3,4,5,6,7,8,9 The
impurity potential in a repulsive Luttinger liquid becomes dressed by long-range oscillations which suppress the spectral weight for single-particle excitations near the impurity
and also the conductance through the impurity down to zero in the low-energy limit.
The asymptotic low-energy properties of Luttinger liquids with a single impurity are
rather well understood. Universal power-laws and scaling functions have been obtained by

2
bosonization, conformal field theory, and exact solutions for the low-energy asymptotics
in special integrable cases.10 What remains to be developed, however, is a many-body
method for microscopic models of interacting fermions with impurities, which does not
only capture correctly the universal low-energy asymptotics, but allows one to compute
observables on all energy scales, providing thus also non-universal properties, and in
particular an answer to the important question at what scale the ultimate asymptotics
sets in. That scale can indeed be surprisingly low, and the properties above it very
different from the asymptotic behavior. Some of the non-universal properties can be
computed numerically by the density matrix renormalization group (DMRG),11 but this
method is limited to lattice systems with about 1000 sites, and only a restricted set of
observables can be evaluated with affordable computational effort.
In the last few years it has been realized that the functional renormalization group
(fRG) is a source of powerful new computation tools for interacting Fermi systems, especially for low-dimensional systems with competing instabilities and entangled infrared
singularities. The starting point of this approach is an exact hierarchy of differential flow
equations for the Green or vertex functions of the system, which is obtained by taking
derivatives with respect to an infrared cutoff Λ.12 Approximations are then constructed
by truncating the hierarchy and parametrizing the vertex functions with a manageable
set of variables or functions.13
A relatively simple fRG approximation for impurities in spinless Luttinger liquids
has been developed recently by some of us.14 The scheme starts from an fRG hierarchy for one-particle irreducible vertex functions, as first derived in a field theoretical
context by Wetterich15 and Morris,16 and for interacting Fermi systems by Salmhofer
and Honerkamp.17 A Matsubara frequency cutoff is used as the flow parameter. The hierarchy is then truncated already at first order, such that the 2-particle vertex remains
unrenormalized, and the flow of the self-energy, which describes the renormalized impurity
potential, is determined by the bare 2-particle interaction. No simplified parametrization
of the self-energy is necessary such that the full spatial dependence of the renormalized
impurity potential can be obtained for very large lattice systems. In spite of the striking
simplicity of this scheme it was shown that the effects of a single static impurity in a
spinless Luttinger liquid are fully captured qualitatively, and in the weak coupling limit
also quantitatively.14 In particular, one obtains that impurity potentials in repulsive Luttinger liquids become effectively stronger at lower energy scales, and act ultimately as a
weak link between two otherwise separate wires, as predicted by Kane and Fisher.7 The
fRG also correctly yields the universal low energy power-laws with exponents that do
not depend on the bare impurity strength. In addition, it was shown that the asymptotic behavior holds typically only at very low energy scales and for very large systems,
except for very strong bare impurities.14 The fRG approach to impurities in Luttinger
liquids was originally developed and tested for spectral densities of single-particle excitations, but has been applied very recently also to transport problems, such as persistent

3
currents in mesoscopic rings18,19 and the conductance of interacting wires connected to
non-interacting leads.19,20 The full power of the fRG with its ability to deal naturally also
with complex crossover phenomena emerges most convincingly in the multi-scale problem
posed by the transport through a resonant double barrier.21
In the present work we further develop the fRG approach for impurities in Luttinger
liquids by including 2-particle vertex renormalization, that is we go one step further in
the hierarchy of flow equations than previously.22 For spinless fermions this extension
does not matter qualitatively, as the lowest order is already qualitatively correct, but the
quantitative accuracy of the results improves considerably in particular at intermediate
1
interaction strengths. For spin- 2 systems, which we will treat in a subsequent work,23
vertex renormalization is necessary to take into account that backscattering of particles
with opposite spins at opposite Fermi points scales to zero in the low energy limit. A
crucial point is to devise an efficient parametrization of the vertex by a managable number
of variables. Here we focus on spinless lattice fermions with nearest neighbor interaction
as a prototype model. The bulk model is supplemented by site or hopping impurities.
We also analyze the influence of boundaries, which can be viewed as infinite barriers or
infinitesimal weak links. We choose to parametrize the vertex by a renormalized nearest
neighbor interaction, which allows us to capture various advantageous features: the low
energy flow of the vertex at kF in the pure system is obtained correctly to second order
in the renormalized couplings; the non-universal contributions at finite energy scales are
correct to second order in the bare interaction; the algorithm for the flow of the self-energy
remains as fast as in the absence of vertex renormalization, such that one can easily deal
with up to 106 lattice sites. We compute spectral properties of single-particle excitations
near an impurity and the oscillations in the density profile induced by an impurity, or
by a boundary. The accuracy of the calculation is checked by comparing with DMRG
results for systems with up to 1000 lattice sites and with exact results for the asymptotic
behavior, which can be obtained from the Bethe ansatz and bosonization.
The article is structured as follows. In Sec. II we introduce the microscopic model
and various types of impurities. The fRG formalism is developed in Sec. III, and worked
out explicitly for the spinless fermion model with nearest neighbor interaction. Sec. IV is
dedicated to results for the renormalized impurity potential, spectral properties, and the
density profile. We finally conclude in Sec. V with an outline of promising extensions of
the present work.
II.

MICROSCOPIC MODEL

We consider the spinless fermion lattice model with nearest neighbor interaction and
various types of impurity potentials. The Hamiltonian has the form
H = H0 + HI + Himp

(1)

4
where the kinetic energy
c† cj + c† cj+1
j+1
j

H0 = −t

(2)

j

is given by nearest neighbor hopping with an amplitude t, and
HI = U

nj nj+1

(3)

j

is a nearest neighbor interaction of strength U. We use standard second quantization
notation, where c† and cj are creation and annihilation operators on site j, respectively,
j
†
and nj = cj cj is the local density operator. The impurity is represented by
V j ′ j c† ′ cj
j

Himp =

(4)

j,j ′

where Vj ′j is a static potential. For ”site impurities”
Vj ′j = Vj δjj ′

(5)

this potential is local. For the important special case of a single site impurity
Vj = V δjj0

(6)

the potential acts only on one site j0 . We also consider ”hopping impurities” described
by the non-local potential
Vj ′j = Vjj ′ = −tj,j+1 δj ′ ,j+1 .

(7)

For the special case of a single hopping impurity,
tj,j+1 = (t′ − t) δjj0 ,

(8)

the hopping amplitude t is replaced by t′ on the bond linking the sites j0 and j0 + 1. In
the following we will set the bulk hopping amplitude t equal to one, that is all energies
are given in units of t.
The clean spinless fermion model H0 + HI is exactly solvable via the Bethe ansatz.24
The system is a Luttinger liquid for all particle densities n and any interaction strength
except at half-filling for |U| > 2. For U > 2 a charge density wave with wave vector π
forms, for U < −2 the system undergoes phase separation. The Luttinger liquid parameter
Kρ , which determines all the critical exponents of the liquid, can be computed exactly
from the Bethe ansatz solution.25 At half-filling Kρ is related to U by the simple explicit
formula
U
2
−1
(9)
Kρ = arccos −
π
2

5
for |U| ≤ 2.
We will compute physical properties for the above models on finite chains with open
boundary conditions. The fRG equations can be solved easily for chains with up to 106
sites, the DMRG can be applied for up to 103 sites. To avoid oscillations emerging from
the ends of the chain we will sometimes attach non-interacting semi-infinite leads to the
finite interacting chain, with a smoothly decaying interaction at the contact.
III.

FUNCTIONAL RENORMALIZATION GROUP

This section is dedicated to a detailed presentation of the fRG equations for onedimensional Fermi systems with impurities, including vertex renormalization, at zero
temperature. We use the one-particle irreducible version of the functional RG.15,16 The
starting point is an exact hierarchy of flow equations for the irreducible vertex functions,
which is obtained by cutting off the infrared part of the free propagator on a scale Λ,
and differentiating the generating functional for the vertex functions with respect to this
scale.17
A.

Cutoff and flow equations

The cutoff can be imposed in many different ways. The only requirement is that
the infrared singularities must be regularized such that the flow equations allow for a
regular perturbation expansion in powers of the renormalized two-particle vertex. For our
purposes a sharp Matsubara frequency cutoff is the most efficient choice. A momentum
cutoff is less suitable here since the impurity spoils momentum conservation. The cutoff is
imposed by excluding modes with frequencies below scale Λ from the functional integral
representation of the system, or equivalently, by introducing a regularized bare propagator
GΛ (iω) = Θ(|ω| − Λ) G0 (iω) .
0

(10)

Here G0 is the bare propagator of the pure system, involving neither interactions nor
impurities. Instead of the sharp cutoff imposed by the step function Θ one may also
choose a smooth cutoff function, but the sharp cutoff has the advantage that it reduces
the number of integration variables on the right hand side of the flow equations. Note
that we will frequently write expressions which are well defined only if the sharp cutoff
is viewed as a limit of increasingly sharp smooth cutoff functions. The suppression of
frequencies below scale Λ affects all Green and vertex functions of the interacting system,
which become thus functions of Λ. The original system is recovered in the limit Λ → 0.
The first two equations in the exact hierarchy of flow equations for the one-particle
irreducible m-particle vertex functions are depicted graphically in Figs. 1 and 2. The first
equation yields the flow of the self-energy, which is related to the interacting propagator

6
by the usual Dyson equation
GΛ = (GΛ )−1 − ΣΛ
0

−1

.

(11)

Here and below GΛ , ΣΛ etc. are operators, which do not refer to any particular singleparticle basis, unless we write matrix indices explicitly. The right hand side of the flow
equation involves the 2-particle vertex ΓΛ , and the so-called single scale propagator S Λ
defined as
S Λ = GΛ ∂Λ (GΛ )−1 GΛ = −
0

1
∂GΛ
1
0
.
Λ Λ
1 − G0 Σ ∂Λ 1 − ΣΛ GΛ
0

(12)

Note that S Λ has support only on a single frequency scale: |ω| = Λ. The precise flow
equation for the self-energy (sketched graphically in Fig. 1) reads
∂ Λ ′
1
Σ (1 , 1) = −
∂Λ
β

+

eiω2 0 S Λ (2, 2′ ) ΓΛ (1′ , 2′ ; 1, 2) ,

(13)

2,2′

where β is the inverse temperature. The numbers 1, 2, etc. are a shorthand for Matsubara
frequencies and labels for single-particle states such as site and spin indices. Note that
′
′
ω1 = ω1 and ω2 = ω2 due to Matsubara frequency conservation. The exponential factor
in the above equation is irrelevant at any finite Λ, but is necessary to define the initial
conditions of the flow at Λ = Λ0 → ∞ (see below).
The right hand side of the flow equation for the 2-particle vertex ΓΛ , shown in Fig.
2, involves ΓΛ itself, but also the 3-particle vertex ΓΛ . The latter could be computed
3
from a flow equation involving the 4-particle vertex ΓΛ , and so on. To avoid this un4
managable proliferation of vertex functions we make our first approximation: we neglect
the contribution of the 3-particle vertex to the flow of ΓΛ . The coupled system of flow
equations for ΓΛ and ΣΛ is then closed. The contribution of ΓΛ to ΓΛ is small as long
3
as ΓΛ is sufficiently small, because ΓΛ is initially (at Λ0 ) zero and is generated only from
3
terms of third order in ΓΛ .17 A comparison of the RG results to exact DMRG results and
exact scaling properties shows that the truncation error is often surprisingly small even
for rather large interactions.14 The explicit form of the truncated flow equation for the
2-particle vertex reads
1
∂ Λ ′ ′
Γ (1 , 2 ; 1, 2) =
∂Λ
β
×

GΛ (3, 3′) S Λ (4, 4′)
3,3′ 4,4′

ΓΛ (1′ , 2′ ; 3, 4) ΓΛ(3′ , 4′ ; 1, 2)
− ΓΛ (1′ , 4′ ; 1, 3) ΓΛ(3′ , 2′ ; 4, 2) − (3 ↔ 4, 3′ ↔ 4′ )
+ ΓΛ (2′ , 4′ ; 1, 3) ΓΛ(3′ , 1′ ; 4, 2) + (3 ↔ 4, 3′ ↔ 4′ )

The various contributions are shown diagrammatically in Fig. 3.

(14)

7
For a sharp frequency cutoff the frequency sums on the right hand side of the flow equations can be carried out analytically in the zero temperature limit, where the Matsubara
sums become integrals. At this point one has to deal with products of delta functions
δ(|ω| − Λ) and expressions involving step functions Θ(|ω| − Λ). These at first sight ambiguous expressions are well defined and unique if the sharp cutoff is implemented as a
limit of increasingly sharp broadened cutoff functions Θǫ , with the broadening parameter ǫ tending to zero. The expressions can then be conveniently evaluated by using the
following relation,16 valid for arbitrary continuous functions f :
1

δǫ (x − Λ) f [Θǫ (x − Λ)] → δ(x − Λ)

f (t) dt ,

(15)

0

where δǫ = Θ′ǫ . Note that the functional form of Θǫ for finite ǫ does not affect the result
in the limit ǫ → 0.
Instead of writing down the frequency-integrated flow equations in full generality, we
first implement our second approximation: the frequency dependent flow of the renormalized 2-particle vertex ΓΛ is replaced by its value at vanishing (external) frequencies, such
that ΓΛ remains frequency independent. As a consequence, also the self-energy becomes
frequency independent. Since the bare interaction is frequency independent, neglecting
the frequency dependence leads to errors only at second order (in the interaction strength)
for the self-energy, and at third order for the vertex function at zero frequency. Besides
the quantitative errors we miss qualitative properties related to the frequency dependence
of the self-energy, such as the suppression of the one-particle spectral weight in the bulk
of a pure Luttinger liquid. On the other hand, a comparison with exact numerical results
and asymptotic analytical results shows that the impurity effects are not qualitatively
affected by the frequency dependence of Σ, at least for weak interactions.
Carrying out the frequency integration in the flow equation for the self-energy, and
inserting a frequency independent 2-particle vertex, one obtains
1
∂ Λ
+
˜
(16)
eiω0 GΛ ′ (iω) ΓΛ′ ,2′ ;1,2 ,
Σ1′ ,1 = −
1
2,2
∂Λ
2π ω=±Λ ′
2,2

where the lower indices 1, 2, etc. label single-particle states (not frequencies) and
˜
GΛ (iω) = G−1 (iω) − ΣΛ
0

−1

.

(17)

˜
Note that GΛ has no jump at |ω| = Λ, in contrast to GΛ . The frequency integrated flow
equation for the 2-particle vertex, evaluated at vanishing external frequencies, has the
form
∂ Λ
Γ′ ′
=
∂Λ 1 ,2 ;1,2
1 ˜Λ
1
˜ 4,4
G ′ (iω) GΛ ′ (−iω) ΓΛ′ ,2′ ;3,4 ΓΛ′ ,4′ ;1,2
1
3
2π ω=±Λ ′ ′
2 3,3
3,3

4,4

˜
˜
+ GΛ ′ (iω) GΛ ′ (iω) −ΓΛ′ ,4′ ;1,3 ΓΛ′ ,2′ ;4,2 + ΓΛ′ ,4′ ;1,3 ΓΛ′ ,1′ ;4,2
3,3
4,4
1
3
2
3

(18)

8
The flow is determined uniquely by the differential flow equations and the initial conditions at Λ = ∞. For Λ = ∞ the vertex functions are given by the bare interactions of the
system. In particular, the flow of the 2-particle vertex starts from the antisymmetrized
bare 2-particle interaction while m-particle vertices of higher order vanish at Λ = ∞, in
the absence of m-body interactions with m > 2. The self-energy at Λ = ∞ is given by
the bare one-particle potential, that is by those one-particle terms which are not included
already in G0 . We choose to include the translationally invariant bulk kinetic energy and
the chemical potential in G0 , while the impurity (site or hopping) potential V is assigned
to the initial condition of the self-energy. In a numerical solution the flow starts at some
large finite initial cutoff Λ0 . Here one has to take into account that, due to the slow decay
of the right hand side of the flow equation for ΣΛ at large Λ, the integration of the flow
from Λ = ∞ to Λ = Λ0 yields a contribution which does not vanish in the limit Λ0 → ∞,
˜
but rather tends to a finite constant. Since GΛ ′ (iω) → δ2,2′ /(iω) for |ω| = Λ → ∞, this
2,2
constant is easily determined as
−

1
lim
2π Λ0 →∞

Λ0

+

eiω0

dΛ
∞

ω=±Λ 2,2′

δ2,2′
1
I1′ ,2′ ;1,2 =
iω
2

I1′ ,2;1,2 ,

(19)

2

where I1′ ,2′ ;1,2 is the bare antisymmetrized interaction. In summary, the initial conditions
for the self-energy and the 2-particle vertex at Λ = Λ0 → ∞ are
ΣΛ0 ′ = V1,1′ +
1,1

1
2

I1′ ,2;1,2

(20)

2

ΓΛ′0 ′ ;1,2 = I1′ ,2′ ;1,2 ,
1 ,2

(21)
+

where V1,1′ is the bare impurity potential. For the flow at Λ < Λ0 the factor eiω0 in Eq.
(16) can be discarded.
B.

Parametrization of 2-particle vertex

For a finite lattice system with L sites the flow of the 2-particle vertex ΓΛ′ ,2′ ;1,2 involves
1
O(L3 ) independent flowing variables, if translation invariance is assumed, and O(L4 )
variables, if the influence of the impurity on the flow of the 2-particle vertex is taken into
account. For a treatment of large systems it is therefore necessary to reduce the number
of variables by a suitable approximate parametrization of the vertex. In the low energy
limit (small Λ) the flow is dominated by a very small number of variables anyway, the
others being irrelevant according to standard RG arguments.2 In particular, the frequency
dependence of the vertex, discarded already above, is irrelevant for the flow of ΓΛ at small
Λ. For larger Λ one can use perturbation theory as a guide for a simple but efficient
parametrization of ΓΛ .
˜
We neglect the feedback of the self-energy into the flow of ΓΛ and replace GΛ by G0 in
Eq. (18). The renormalization of ΓΛ includes thus only bulk, not impurity, contributions,

9
such that ΓΛ remains translation invariant. While this turns out to be sufficient for
capturing the effects of isolated impurities in otherwise pure systems, it is known that
impurity contributions to vertex renormalization become important in macroscopically
disordered systems.26 For a more concrete treatment of the vertex renormalization, we
have to focus on a specific model. In this article we restrict ourselves to the spinless
1
fermion model with nearest neighbor interaction, leaving an extension to spin- 2 systems
with more general interactions for subsequent work.23
For spinless fermions the 2-particle vertex and the self-energy are fully characterized
by either site or momentum variables. In the low energy limit, the flow of the vertex is
dominated by contributions with momenta close to the Fermi points, such that the right
hand side of the flow equation is determined by momentum components of the vertex
′
′
ΓΛ′ ,k′ ;k1 ,k2 with k1 , k2 , k1, k2 = ±kF . Due to the antisymmetry of the vertex, there is only
k1 2
one such component which is non-zero:
g Λ = ΓΛF ,−kF ;kF ,−kF .
k

(22)

In the low energy limit the momentum dependence of the vertex away from ±kF is irrelevant. There are therefore many possible choices for the functional form of ΓΛ′ ,k′ ;k1 ,k2 ,
k1 2
which all lead to the correct low energy asymptotics. For a model with a bare nearest
neighbor interaction U, a natural and efficient choice is to parametrize the flowing vertex
simply by a renormalized nearest neighbor interaction U Λ , which leads to a real space
vertex of the form
′
′
′
′
(23)
ΓΛ1 ,j2 ;j1 ,j2 = UjΛ ,j2 (δj1 ,j1 δj2 ,j2 − δj1 ,j2 δj2 ,j1 )
j′ ′
1
with UjΛ ,j2 = U Λ (δj1 ,j2−1 + δj1 ,j2 +1 ) . This yields the following structure in momentum
1
space:
(2π)
′
′
ΓΛ1 ,k2 ;k1 ,k2 = 2U Λ [cos(k1 − k1 ) − cos(k2 − k1 )] δk1 +k2 ,k′ +k′ ,
(24)
k′ ′
1

2

where the Kronecker delta implements momentum conservation (modulo 2π). The flowing
coupling constant U Λ is linked to the value of the vertex at the Fermi points by the relation
g Λ = 2U Λ [1 − cos(2kF )] .

(25)

The flow equation for g Λ becomes
1
∂g Λ
=
∂Λ
2π ω=±Λ

dp
( PP + PH + PH′ )
2π

(26)

with the particle-particle and particle-hole contributions
PP =

1 0
Gp (iω) G0 (−iω) ΓΛF ,−kF ;p,−p ΓΛ
−p
k
p,−p;kF ,−kF
2

PH = −[G0 (iω)]2 ΓΛF ,p;kF ,p ΓΛ F ;p,−kF
p
k
p,−k
PH′ = G0 F (iω) G0 F (iω) ΓΛ F ,p+kF ;kF ,p−kF ΓΛ F ,kF ;p+kF ,−kF ,
p−k
p+k
−k
p−k

(27)

10
where ΓΛ on the right hand side of the flow equation is given by Eq. (24). Using Eq. (25)
to replace ∂Λ g Λ by ∂Λ U Λ on the left hand side of Eq. (26), one obtains a flow equation
for U Λ of the simple form
∂Λ U Λ = h(Λ) (U Λ )2 .
(28)
The function h(Λ) depends only on the cutoff Λ and the Fermi momentum kF . An explicit
formula for h(Λ) can be obtained by carrying out the momentum integral in Eq. (26) using
the residue theorem, as sketched in Appendix A. For finite systems the momentum integral
should be replaced by a discrete momentum sum; however, this leads to sizable corrections
only for very small systems. The flow equation (28) can be integrated to
UΛ =

U
,
1 − U H(Λ)

(29)

where H(Λ) is the primitive function of h(Λ) with H(Λ) → 0 for Λ → ∞. Integrating
h(Λ) one obtains
H(Λ) = −

Λ − iµ0
1
(4 − µ2 )Λ2 − 2iµ0 (2 − µ2 )Λ + µ4 − 6µ2 + 8 iµ0
Λ
0
0
0
0
+ Re
sinh−1
−
2π π
2
2
2 (4 − µ2 ) Λ2 − 2iµ0 Λ + 4 − µ2
0
0

µ4
0
+
tanh−1
2 3/2
2(4 − µ0 )

4 + µ2 + iµ0 Λ
0
(4 + µ2 + iµ0 Λ)2 + 4(Λ − 2iµ0 )2
0

,

(30)

where µ0 = −2 cos kF , while sinh−1 and tanh−1 denote the main branch of the inverse
of the complex functions sinh and tanh, respectively. At half-filling, corresponding to
kF = π/2, the function h(Λ) is particularly simple (see Appendix A), such that also U Λ
has a relatively simple form:
U

UΛ =
1+ Λ−

2+Λ2
√
4+Λ2

.

(31)

U/(2π)

In Fig. 4 we show results for the renormalized nearest neighbor interaction U Λ as obtained
from the flow equation at various densities n, for a bare interaction U = 1. While
the renormalization does not follow any simple rule at intermediate scales Λ, all curves
saturate at a finite value U ∗ in the limit Λ → 0, corresponding to a finite g ∗ , as expected
for a Luttinger liquid fixed point.2
The above simplified flow equation yields not only the correct low energy asymptotics
to second order in the renormalized vertex, but contains also all non-universal second
order corrections to the vertex at ±kF at higher energy scales. Hence the resulting fixed
point coupling g ∗ , from which the Luttinger liquid parameter Kρ can be computed, is
obtained correctly to second order in U. The calculation of Kρ , presented in Appendix
B, illustrates how Luttinger liquid parameters can be related to the fixed point couplings
obtained from the fRG. A comparison of the RG result for Kρ with exact results from the
Bethe ansatz solution of the spinless fermion model shows that the vertex renormalization
scheme described above is not only very simple, but also surprisingly accurate.

11
Parametrizing ΓΛ by a renormalized nearest neighbor interaction has the enormous
advantage that the self-energy, as determined by the flow equation (16), is a tridiagonal
matrix in real space, that is only the matrix elements ΣΛ and ΣΛ
j,j
j,j±1 are non-zero. InsertΛ
ing Γ from Eq. (23) into (16), one obtains the following simple coupled flow equations
for the diagonal and off-diagonal matrix elements:
UΛ
∂ Λ
Σj,j = −
∂Λ
2π
∂ Λ
Σ
=
∂Λ j,j±1

UΛ
2π

˜
GΛ
j+r,j+r (iω)
ω=±Λ r=±1

˜
GΛ (iω) .
j,j±1

(32)

ω=±Λ

˜
Note that the self-energy enters also the right hand side of these equations, via GΛ =
(G−1 − ΣΛ )−1 . Since ΣΛ and G−1 are both tridiagonal in real space, the matrix inversion
0
0
˜
required to compute the diagonal and first off-diagonal elements of GΛ from ΣΛ can be
carried out very efficiently, and thus for very large systems. In Appendix C we describe an
algorithm for the numerical solution of the flow equation for Σ which scales only linearly
with the size of the system.
C.

Observables

In the next section we will present results for spectral properties of single-particle
excitations near an impurity or boundary, and for the density profile. Here we describe
how the relevant observables are computed from the solution of the flow equations.
1.

Single-particle excitations

Integrating the flow equation for the self-energy ΣΛ down to Λ = 0 yields the physical
(cutoff-free) self-energy Σ, and the single-particle propagator G = (G−1 − Σ)−1 . From
0
the latter spectral properties of single-particle excitations can be extracted. We focus on
local spectral properties, which are described by the local spectral function
1
ρj (ω) = − Im Gjj (ω + i0+ ) ,
π

(33)

where Gjj (ω + i0+ ) is the local propagator, analytically continued to the real frequency
axis from above.
In our approximation the self-energy is frequency-independent and can therefore be
viewed as an effective single-particle potential. The propagator G is thus the Green
function of an effective single-particle Hamiltonian. In real space representation this
Hamiltonian is given by the tridiagonal matrix heff = h0 + Σ, where the matrix elements
of h0 are the hopping amplitudes in H0 , Eq. (2). For a lattice with L sites this matrix has L

12
(including possible multiplicities) eigenvalues ǫλ and an orthonormal set of corresponding
eigenvectors ψλ . For the spectral function ρj (ω) one thus obtains a sum of δ-functions
ρj (ω) =

wλj δ(ω − ξλ ) ,

(34)

λ

where ξλ = ǫλ − µ, and the spectral weight wλj is the squared modulus of the amplitude
of ψλ on site j. For large L the level spacing between neighboring eigenvalues is usually
of order L−1 , except for one or a few levels outside the band edges which correspond to
bound states.
Due to even-odd effects etc. the spectral weight wλj generally varies quickly from
one eigenvalue to the next one. A smooth function of ω which suppresses these usually
irrelevant finite-size details can be obtained by averaging over neighboring eigenvalues.
In addition it is useful to divide by the level spacing between eigenvalues to obtain the
conventional local density of states, which we denote by Dj (ω).
2.

Density profile

Boundaries or impurities induce a density profile with long-range Friedel oscillations,
which are expected to decay with a power-law at long distances.28 The expectation value
of the local density nj could be computed from the local one-particle propagator Gjj , if
G was known exactly. However, the approximate flow equations for Σ can be expected to
describe the asymptotic behavior of G correctly only at long distances between creation
and annihilation operator in time and/or space, while in the local density operator time
and space variables coincide. In the standard RG terminology nj is a composite operator,
which has to be renormalized separately.29
To derive a flow equation for nj , we follow the usual procedure for the renormalization
of correlation functions involving composite operators: we add a term φj nj with a small
field φj to the Hamiltonian, and take derivatives with respect to φj in the flow equations.
The local density is given by
∂Ω(φj )
,
(35)
nj =
∂φj φj =0
where Ω(φj ) is the grand canonical potential of the system in the presence of the field φj .
Note that we use the same symbol nj for the density operator and its expectation value.
In the presence of a cutoff Λ the grand canonical potential obeys the exact flow-equation
∂ Λ
1
Ω =
∂Λ
β

+

tr eiω0 [∂Λ GΛ (iω)−1 ] [GΛ (iω) − GΛ (iω)]
0
0

(36)

ω

in a short-hand matrix notation. This flow-equation follows from the functional flow
equation for the generating functional for vertex functions17 and the relation between
the grand canonical potential and the zero-particle vertex. At zero temperature the

13
Matsubara frequency sum becomes an integral which, for the sharp frequency cutoff (10),
can be carried out analytically. This yields
1
∂ Λ
+
tr eiω0 log[1 − G0 (iω) ΣΛ (iω)]
Ω =
∂Λ
2π ω=±Λ

.

(37)

We attribute φj to the ”interaction” part of the Hamiltonian, not to H0 , such that G0
remains independent of φj . The self-energy is modified via the additional local and
frequency-independent contribution φj δjj ′ to its initial value ΣΛ0′ at scale Λ0 .
jj
The density profile can be obtained from the above equations and the flow equation
for ΣΛ by computing the shift of Ω generated by a small finite perturbation φj , that is by
numerical differentiation. Alternatively one may carry out the φj -derivative analytically
in the flow equations, which yields a flow equation for the density in terms of the density
response vertex. Taking the φj derivative in (37) yields
∂ Λ
1
+
Λ
˜
tr eiω0 GΛ (iω) Rj (iω)
nj = −
∂Λ
2π ω=±Λ

(38)

with the density response vertex
Λ
Rj (iω) =

∂ΣΛ (iω)
∂φj

,

(39)

φj =0

˜
where the propagator GΛ is computed as previously, that is in the absence of φj . We
compute the self-energy ΣΛ in the presence of φj within the same approximation as previously. It is thus determined from the flow equation (16) with a frequency-independent
two-particle vertex ΓΛ . Taking a derivative of that equation with respect to φj at φj = 0
yields the flow equation for the response vertex
1
∂ Λ
Rj;1′ ,1 = −
∂Λ
2π

Λ
˜
˜
GΛ ′ (iω) Rj;3′,3 GΛ ′ (iω) ΓΛ′ ,2′ ;1,2 .
1
3,2
2,3

(40)

ω=±Λ 2,2′ 3,3′

Λ
Note that Rj is frequency-independent in our approximation and that there is no contribution from the φj -derivative of ΓΛ , since we neglect self-energy contributions in the flow
of the two-particle vertex.
For spinless fermions with a (renormalized) nearest neighbor interaction, the matrix
Λ
Λ
Λ
Rj is tridiagonal, that is only the components Rj;l,l and Rj;l,l±1 are non-zero, and their
flow is given by

∂ Λ
UΛ
Rj;l,l = −
∂Λ
2π
UΛ
∂ Λ
Rj;l,l±1 =
∂Λ
2π

Λ
˜′
˜
GΛ ′ (iω) Rj;l′,l′ +r′ GΛ+r′ ,l+r (iω)
l
l+r,l
ω=±Λ l′

r=±1 r ′ =0,±1
Λ
˜
˜′
GΛ ′ (iω) Rj;l′,l′ +r′ GΛ+r′ ,l±1 (iω) .
l,l
l

ω=±Λ

l′

r ′ =0,±1

(41)

14
The right hand sides involve only one unrestricted lattice summation, and can be computed in O(L) time (see Appendix C). The initial condition for the response vertex is
Λ0
1
Rj;l,l′ = δjl δll′ . The initial condition for the density is nΛ0 = 2 , for any filling, due to
j
the slow convergence of the flow equation (38) at large frequencies, which yields a finite
contribution to the integrated flow from Λ = ∞ to Λ0 for arbitrarily large finite Λ0 , as in
the case of the self-energy discussed in more detail above.
To avoid the interference of Friedel oscillations emerging from the impurity or one
boundary with those coming from the (other) boundaries of our systems we suppress the
influence of the latter by coupling the finite chain to semi-infinite non-interacting leads,
with a smooth decay of the interaction at the contacts. The presence of leads modifies
the self-energy in the interacting region only via a boundary term. In particular, a lead
coupled to the first site of the interacting region via nearest neighbor hopping (with
amplitude t = 1) yields the additional contribution
Σlead (iω) =
j,j ′

iω + µ0
2

1−

1−

4
(iω + µ0 )2

δ1j δjj ′

(42)

to the first diagonal element of self-energy matrix in real space.30 Note that this term
is frequency dependent and independent of Λ. A lead coupled to the last site yields an
analogous contribution to the last diagonal element of the self-energy.
IV.

RESULTS

We now present and discuss explicit results for the self-energy, spectral properties
near the impurity or boundary, and the density profile, as obtained from the fRG. A
comparison with exact DMRG results is made for the spectral weight at the Fermi level
and for the density profile. Some fRG results for the self-energy and spectral properties
of the spinless fermion model have been published already previously,14 but only for
half-filling and without any vertex renormalization. Here we present results also away
from half-filling and demonstrate the quantitative improvement of accuracy obtained by
including vertex renormalization. In addition we present for the first time fRG results
for the density profile. Using a faster algorithm (see Appendix C) than previously we are
now able to solve the flow equations for much larger systems.
The typical shape of the self-energy in the vicinity of an impurity can be seen in Fig. 5,
where we plot the diagonal elements Σj,j and the off-diagonal elements Σj,j+1 near a site
impurity of strength V = 1.5 added to the spinless fermion model with interaction strength
U = 1 at quarter-filling. Recall that the self-energy is tridiagonal in real space and
frequency independent within our treatment. The diagonal elements can be interpreted as
a local effective potential, the off-diagonal elements as a non-local effective potential which
renormalizes the hopping amplitudes. At long distances from the impurity both Σj,j and

15
Σj,j+1 tend to a constant. The former describes just a bulk shift of the chemical potential,
the latter a bulk renormalization of the hopping amplitude toward larger values. The
oscillations around the bulk shifts are generated by the impurity. The wave number of the
oscillations is 2kF = π/2, where kF is the Fermi wave vector of the bulk system at quarterfilling. The amplitude of the oscillations of Σ decays slower than the inverse distance
from the impurity at intermediate length scales, but approaches a decay proportional to
1/|j −j0 | for |j −j0 | → ∞. This can be seen most clearly by plotting an effective exponent
βj for the decay of the oscillations, defined as the negative logarithmic derivative of the
oscillation amplitude with respect to the distance |j − j0 |. In Fig. 6 we show the effective
exponent resulting from the oscillations of Σj,j as a function of the distance from a site
impurity, for U = 1 and half-filling. The impurity is situated at the center of a long chain
with L = 218 + 1 sites. To avoid interferences with oscillations from the boundaries we
have attached semi-infinite non-interacting leads to the ends of the interacting chain, as
described in the last paragraph of III.C.2. Only for relatively large impurity strengths
the asymptotic regime corresponding to βj = 1 is reached before finite size effects set
in. For small V one can see that βj increases from values below one, but the asymptotic
long-distance behavior is cut off by the boundaries of the interacting region. For very
small V (for example V = 0.01 in Fig. 6) we observe a plateau in βj for intermediate
distances from the impurity site. In this regime βj is close to Kρ which can be understood
by analytically solving our flow equations for small V .14
The long-range 2kF -oscillations of the self-energy lead to a marked suppression of the
spectral weight for single-particle excitations at the Fermi level, that is at ω = 0. In Fig.
7 we show the local density of states Dj (ω) on the site next to a site impurity of strength
V = 1.5 for the spinless fermion model at half-filling. The result for the interacting system
at U = 1 is compared to the non-interacting case. Even-odd effects have been eliminated
by averaging over neighboring eigenvalues (see Sec. III.C.1). δ-peaks outside the band
edges corresponding to bound states are not plotted. The interaction leads to a global
broadening of the band, which is due to an enhancement of the bulk hopping amplitude,
and also to a strong suppression of Dj (ω) at low frequencies which is not present in the
non-interacting system. For a finite system (here L = 1025) the spectral weight at the
Fermi level remains finite, but tends to zero with increasing system size. In Fig. 8 we
show results for the density of states choosing the same parameters as in Fig. 7, but now
for densities away from half-filling: n = 1/4 and n = 3/4. In addition to the dip near
ω = 0 a second singularity appears at a finite frequency. This effect is due to the fact that
a long-range potential with a wave number 2kF does not only strongly scatter states with
momenta near kF , but also those with momenta close to π − kF . Indeed the singularity is
situated at ω = ǫπ−kF −µ, where ǫk is the renormalized (bulk) dispersion. In the half-filled
case only one singularity is seen simply because π − kF = kF for kF = π/2.
The spectral weight at the Fermi level is expected to vanish asymptotically as a power−1
law L−αB with increasing system size, where αB = Kρ − 1 is the boundary exponent

16
describing the power-law suppression of the density of states at the boundary of a semiinfinite chain.7 That exponent depends only on the bulk parameters of the model, not
on the impurity strength. For the spinless fermion model it can be computed exactly
from the Bethe ansatz solution.25 We therefore analyze the asymptotic behavior of the
spectral weight at the Fermi level by defining an effective exponent α(L) as the negative
logarithmic derivative of the spectral weight with respect to the system size, such that
α(L) tends to a (positive) constant in case of a power-law suppression. In Fig. 9 we
show results for α(L) as obtained from the fRG for the spinless fermion model at quarterfilling with up to about 106 sites, for a weak (U = 0.5) and an intermediate (U = 1.5)
interaction parameter. The spectral weight has been computed either at a boundary,
or near a hopping impurity of strength t′ = 0.5. Results obtained from the RG without
(upper panel) and with (lower panel) vertex renormalization are compared to exact DMRG
results (for up to 512 sites) and the exact boundary exponents αB , plotted as horizontal
lines. The fRG results follow a power-law for large L, with the same asymptotic exponent
for the boundary and impurity case, confirming thus the expected universality. However,
the asymptotic regime is reached only for fairly large systems, even for the intermediate
interaction strength U = 1.5. The comparison with the exact DMRG results and exact
exponents shows that the fRG is also quantitatively rather accurate, and that the inclusion
of vertex renormalization leads to a substantial improvement at intermediate coupling
strength. Results for the effective exponent α in the case of a site impurity are shown in
Fig. 10, at quarter-filling and for an interaction strength U = 1. The comparison of the
different curves obtained for different impurity strengths confirms once again the expected
asymptotic universality, and also how the asymptotic regime shifts rapidly toward larger
systems as the bare impurity strength decreases.
We finally discuss results for the density profile nj . Boundaries and impurities induce
Friedel oscillations of the local density with a wave vector 2kF . In a non-interacting
system these oscillations decay proportionally to the inverse distance from the boundary
or impurity. In an interacting Luttinger liquid the Friedel oscillations are expected to
decay as |j − j0 |−Kρ at long distances |j − j0 |. For a very weak impurity one expects a
slower decay proportional to |j − j0 |1−2Kρ at intermediate distances, and a crossover to
the asymptotic power-law with exponent Kρ at very long distances.28 At intermediate
distances the response of the density to a weak impurity can be treated in linear response
theory, such that the density modulation is determined by the density-density response
function at 2kF , which leads to the power-law decay with exponent 2Kρ − 1. In Fig. 11 we
show fRG and DMRG results for the density profile nj for a spinless fermion chain with
128 sites and interaction strength U = 1 at half-filling. The Friedel oscillations emerge
from both boundaries and interfere in the center of the chain. The accuracy of the RG
results is excellent for all j. For incommensurate filling factors the density profile looks
more complicated. This can be seen in Fig. 12, where functional RG results are shown for
the density modulation |nj − n| near the boundary of a system with an average density

17
n = 0.393 and 8192 sites. For long distances from the boundary the oscillation amplitude
has a well-defined envelope which fits to a power-law as a function of j. We now analyze
the long-distance behavior of the amplitudes more closely for the half-filled case, and
compare to exact results for the asymptotic exponents. In Fig. 13 we show fRG results
for the amplitude of density oscillations emerging from an open boundary, for a very long
spinless fermion chain with 219 + 1 sites and various interaction strengths U at half-filling.
The other end of the chain (opposite to the open boundary) is smoothly connected to a
non-interacting lead. In a log-log plot (upper panel of Fig. 13) the amplitude follows a
straight line for almost all j, corresponding to a power-law dependence. Deviations from
a perfect power-law can be seen more neatly by plotting the effective exponent αj , defined
as the negative logarithmic derivative of the amplitude with respect to j (see the lower
panel of Fig. 13). The effective exponent is almost constant except at very short distances
or when j approaches the opposite end of the interacting chain, which is not surprising.
From a comparison with the exact exponent (horizontal lines in the figure) one can assess
the quantitative accuracy of the fRG results. Effective exponents describing the decay
of Friedel oscillations generated by site impurities of various strengths are shown in Fig.
14, for a half-filled spinless fermion chain with 218 + 1 sites and interaction U = 1. Both
ends of the interacting chains are coupled to non-interacting leads to suppress oscillations
coming from the boundaries. For strong impurities the results are close to the boundary
result (cf. Fig. 13), as expected. For weaker impurities the oscillations decay more slowly,
that is with a smaller exponent, and approach the boundary behavior only asymptotically
at large distances (beyond the range of our chain for V < 1). For very weak impurities
(V = 0.01 in Fig. 14) the oscillation amplitude follows a power-law corresponding to the
linear response behavior with exponent 2Kρ − 1 at intermediate distances. We finally
present some results for the effective exponent of the density oscillation decay in the case
of an attractive interaction U = −1, see Fig. 15. In that case the effective impurity
strength should scale to zero at low energies and long distances.7 Indeed, for weak and
moderate bare impurity potentials the effective exponent in Fig. 15 approaches the linear
density response exponent 2Kρ −1. Only for very strong impurities the density oscillations
decay with the smaller boundary exponent over several length scales.
V.

CONCLUSION

In summary, we have shown that the fRG provides an ideal tool for computing the
intriguing properties of Luttinger liquids with static impurities. The method yields ”ab
initio” results for microscopic model systems at all energy scales from the Fermi energy
to the ultimate low energy limit. We have demonstrated the power of the method by
computing spectral properties of single-particle excitations as well as the oscillations in
the density profile induced by impurities or boundaries for a spinless fermion model with
nearest neighbor interaction. With the inclusion of vertex renormalization, in addition to

18
the renormalization of the effective impurity potential, our results agree remarkably well
with exact asymptotic results and numerical DMRG data even for intermediate interaction
strength.
There is a broad range of interesting further applications and extensions of the fRG
for impurities in Luttinger liquids beyond the scope of the present article:
Spin- 1 fermions: The inclusion of spin degrees of freedom requires a parametrization
2
of the two-particle vertex with several coupling constants. We have already derived flow
equations for the (extended) Hubbard model, with an effective vertex parametrized by
local and nearest neighbor interactions.23
Feedback of impurities on vertex: For isolated impurities the influence of impurities on
the vertex renormalization is irrelevant for the asymptotic low energy or long distance
behavior, although it may contribute quantitatively at intermediate scales. In the case of
disordered systems with a finite impurity density the influence of the latter on the twoparticle vertex is crucial and must be taken into account.26 In principle this is possible by
computing the vertex flow with full propagators, which contain the renormalized impurity
potential via the self-energy.
Finite temperature: The fRG approach can be extended without major complications
to finite temperature. This is particularly useful for studying the temperature dependence
of transport properties.30
Transport: In the absence of inelastic processes (ImΣ = 0) the conductance of the
interacting wire can be computed from the one-particle Green function in the presence of
leads.30 Several results have already been presented in short articles.19,20,21
Inelastic processes: Inelastic processes appear at second order in the interaction and
can be included in the flow equations by inserting the second order vertex into the flow
equation for the self-energy without neglecting its frequency dependence. This procedure
would also capture the anomalous dimension of the bulk system, which is missing in the
present work.

Acknowledgements:
We are grateful to Manfred Salmhofer for useful discussions. This work benefitted from the
workshop ”Realistic theories of correlated electron materials” in fall 2002 at ITP-UCSB.
V.M. thanks the Bundesministerium f¨ r Bildung und Forschung for financial support.
u
APPENDIX A: EVALUATION OF VERTEX FLOW

Here we sketch some details concerning the explicit evaluation of the flow equations
for the 2-particle vertex. Inserting the momentum structure of ΓΛ , (24), into the flow

19
equation (26) and replacing g Λ by U Λ on the left hand side yields
(U Λ )2
∂U Λ
=
∂Λ
2π sin2 kF

2π
ω=±Λ

0

dp
f (p, ω)
2π

(A1)

where
f (p, ω) =

2 sin2 kF sin2 p
(cos kF − cos p)2
[cos(2kF ) − cos p]2
. (A2)
−
+
0
0
0
0
0
(iω − ξp )(−iω − ξ−p )
(iω − ξp )2
(iω − ξp−kF )(iω − ξp+kF )

0
Here ξk = −2 cos k − µ0 with µ0 = −2 cos kF is the bare dispersion relation relative to the
bare Fermi level. Since f (p, ω) can be written as a rational function of cos p and sin p, the
p-integral can be carried out analytically using the substitution z = eip and the residue
theorem. After a lengthy but straightforward calculation one obtains the following result
for the coefficient h(Λ) in (28):

h(Λ) = −
×

1
i
− Re (µ0 + iΛ)
2π
2

1−

4
(µ0 + iΛ)2

3iµ4 − 10µ3 Λ − 12iµ2 (Λ2 + 1) + 6Λ3 µ0 + 18Λµ0 + 6iΛ2 + iΛ4
0
0
0
π(2µ0 + iΛ)(4 − µ2 + Λ2 − 2iΛµ0)2
0

. (A3)

At half-filling (µ0 = 0), this reduces to
h(Λ) = −

1
2π

1−Λ

6 + Λ2
(4 + Λ2 )3/2

.

(A4)

APPENDIX B: CALCULATION OF Kρ

In this appendix we show how the Luttinger liquid parameter Kρ , which determines the
critical exponents of Luttinger liquids, can be computed from the fixed point couplings as
obtained from the RG. A comparison of the RG result for Kρ with the exact Bethe ansatz
result for the bulk model (without impurity) serves also as a check for the accuracy of
our vertex renormalization. A relation between the fixed point couplings and Kρ can be
established via the exact solution of the fixed point Hamiltonian of Luttinger liquids, the
Luttinger model.
For spinless fermions, Kρ is determined by the Luttinger model parameters g and vF
as
1 − g/(2πvF )
Kρ =
(B1)
1 + g/(2πvF )
where g is the interaction between left and right movers and vF the effective Fermi velocity
of the model, that is the slope of the (linear) dispersion relation, with a possible shift
due to interactions between particles moving in the same direction (g4 -coupling) already
included.2 We therefore need to extract g and vF from the RG flow in the limit Λ → 0.

20
In order to obtain Kρ correctly to order U 2 , it is sufficient to obtain vF correctly to linear
order in U.
The Luttinger model interaction g and the fixed point coupling g ∗ = ΓΛ→0 F ;kF ,−kF
kF ,−k
from the RG are not simply identical, in contrast to what one might naively expect. To
find the true relation between g and g ∗ , one has to take into account that the forward
scattering limit of the dynamical 2-particle vertex is generally not unique (in the absence
of cutoffs), and depends on whether momentum or frequency transfers tend to zero first.
This ambiguity is well-known in Fermi liquid theory, where it leads to the distinction
between quasi-particle interactions and scattering amplitudes,27 but is equally present in
Luttinger liquids, for the same reason in all cases: the ambiguity of the small momentum,
small frequency limit of particle-hole propagators contributing to the vertex function. In
the dynamical limit, where the momentum transfer q vanishes first, the singular particlehole propagators do not contribute. In Fermi liquids this limit yields the quasi-particle
interaction. In the opposite static limit the frequency transfer ν vanishes first and particlehole propagators yield a finite contribution. In the presence of an infrared cutoff Λ > 0
the forward scattering limit of the vertex function is unique, since the ambiguity in the
particle-hole propagator is due to the infrared pole of the single-particle propagator. Hence
ΓΛF ,−kF ;kF ,−kF is well defined. However, ΓΛF ,−kF ;kF ,−kF and also its limit for Λ → 0 depend
k
k
on the choice for the cutoff function. For a momentum cutoff, which excludes states with
excitation energies below Λ around the Fermi points, particle-hole excitations with small
momentum transfers q are impossible. Hence particle-hole propagators with infinitesimal
q do not contribute to the vertex at any Λ > 0, such that ΓΛF ,−kF ;kF ,−kF converges to the
k
dynamical forward scattering limit, which is simply given by the bare coupling constant g
in the Luttinger model. For a frequency cutoff the particle-hole propagators with vanishing
momentum and frequency transfer yield a finite contribution at Λ > 0, which tends to
the static limit for Λ → 0. This can be seen directly by integrating ω=±Λ dp [G0 (iω)]2
p
Λ
over Λ from infinity to zero. Hence the vertex ΓkF ,−kF ;kF ,−kF obtained from our frequency
cutoff RG tends to the static forward scattering limit.
For the Luttinger model, the static forward scattering limit of the vertex can be obtained from the dynamical effective interaction between left and right movers D(q, iν),
which is defined as the sum of particle-hole chains
D(q, iν) = g + g Π0 (q, iν) g Π0 (q, iν) g + · · · =
−
+
where
Π0 (q, iν) = ±
±

g
1−

g 2 Π0 (q, iν) Π0 (q, iν)
−
+

1
q
2π iν ∓ vF q

(B2)

(B3)

is the bare particle-hole bubble for right (+) and left (-) movers. Note that only odd
powers of g contribute to the effective interaction between left and right movers. This
effective interaction appears naturally in the exact solution of the Luttinger model via

21
Ward identities.31 For the static limit one obtains
lim D(q, 0) =

q→0

g
1 − [g/(2πvF )]2

(B4)

which we identify with our fixed point coupling g ∗ as obtained from the RG with frequency
cutoff. Inverting this relation between g and g ∗ we obtain
g=

2πvF
g∗

−πvF +

(πvF )2 + (g ∗ )2

(B5)

For spinless fermions the difference between g and g ∗ appears only at third order in the
coupling, but for models with spin the distinction becomes important already at second
order.
The Fermi velocity vF can be computed from the (frequency independent) self-energy
in momentum space as
0
vF = vF + ∂k Σk |kF ,
(B6)
0
where vF = ∂k ǫk |kF is the bare Fermi velocity. The self-energy is computed from the flow
equation (32), which can be rewritten in momentum space as

∂ Λ
UΛ
Σk = −
∂Λ
π

ω=±Λ

dp 1 − cos(k − p)
2π iω − ξp − ΣΛ
p

(B7)

where ξp = ǫp − µ. The chemical potential µ has to be fixed by the final condition
ξkF + ΣkF = 0, where kF = πn depends only on the density, not the interaction. From
the tridiagonal structure of Σ in real space, but also from the above expression it follows
that ΣΛ has the form ΣΛ = aΛ + bΛ cos k. The functional flow equation for ΣΛ yields a
k
k
k
coupled set of ordinary flow equations for the coefficients aΛ and bΛ , with initial conditions
aΛ0 = U and bΛ0 = 0. The momentum integrals can be evaluated analytically via the
residue theorem, such that the remaining set of two coupled differential equations (with
U Λ as input) can be easily solved numerically. The result for vF is correct at least to
first order in U, but not necessarily to second order, since our simplified parametrization
captures the 2-particle vertex correctly to second order only at the Fermi points.
Inserting g and vF into the Luttinger model formula (B1) we can now compute Kρ as a
function of U and density for the microscopic spinless fermion model. In Fig. 16 we show
results for Kρ (U) for various fixed densities, which are compared to the exact result from
the Bethe ansatz solution25 in the inset. The RG results are indeed correct to order U 2
for small U, as expected, and they are surprisingly accurate also for larger values of U.
APPENDIX C: NUMERICAL SOLUTION OF Σ FLOW

For Λ < Λ0 < ∞ the flow equation for the self-energy (16) can be written as
1
∂ Λ
Σ1′ ,1 = −
∂Λ
2π

˜
ΓΛ′ ,2′ ;1,2 · 2 Re GΛ ′ (iΛ) .
1
2,2
2,2′

(C1)

22
In order to compute its right-hand side, one needs to invert the tridiagonal matrix
T = G−1 (iΛ) − ΣΛ ,
0

(C2)

where T is complex symmetric (not hermitean) with diagonal elements ai := iΛ + µ − ΣΛ ,
i,i
Λ
i = 1, . . . , L, and first off-diagonal elements bi := t − Σi,i+1 , i = 1, . . . , L − 1. Note that
Im(ai ) = Λ > 0 such that T is non-singular and its inverse well-defined.
˜
The inverse GΛ (iΛ) = T −1 is not tridiagonal but a full matrix which can be computed
by standard methods in O(L2 ) time. However, for an interaction that does not extend
˜
beyond nearest neighbors on the lattice, only the tridiagonal part of G is required, which
can be computed in O(L) time, such that much larger lattices can be treated. We shall
first explain how this is done and then present the resulting algorithm that can directly
be incorporated into a computer program.
Under certain assumptions (see below), a matrix can be uniquely factorized into a
lower unit triangular matrix L, a diagonal matrix D, and an upper unit triangular matrix
U (“LDU factorisation”): T = L · D · U.32 For a tridiagonal matrix T the unit triangular
matrices L and U are in fact unit bidiagonal: their matrix elements are one on the
diagonal, and only the first off-diagonal is nonzero. Since our T is symmetric we have
L = U T . Thus we obtain a factorization of the form




+
+
1 U1
D1
1

  1 U+

 +
+



U1 1
D2
2
+T + +



T =U D U =
.. 
+
+


D3
U2 1
.
1




..
.. ..
..
.
. .
.
where the label “+” distinguishes this factorization from another one used below. The
+
prescription to compute the elements Di and Ui+ is well known and can be found for
example in Ref.32 . Starting in the upper left corner one proceeds to increasing row and
column numbers until one arrives at the lower right corner of T :
+
D1 := a1 ,

+
Ui+ := bi /Di ,

+
Di+1 := ai+1 − bi Ui+

(i = 1, . . . , L − 1) .

(C3)

+
This works well since in our case Im(Di ) ≥ Λ > 0, such that one never divides by zero.
˜
To compute the inverse G = T −1 , one could directly calculate (U + )−1 (D + )−1 (U +T )−1 .
It is however easier and more accurate to find the inverse by solving the linear system of
˜
equations T G = 1, where 1 is the identity matrix, by “backsubstitution”. To be specific,
˜
˜
consider the ith column vector G·,i of G:

˜
˜
ei = T G·,i = U +T (D + U + G·,i ) = U +T gi ,

˜
U + G·,i = (D + )−1 gi

(C4)

where ei is the ith unit vector. The first step is to solve the linear system U +T gi = ei
˜
˜
for gi , and the second step to solve U + G·,i = (D + )−1 gi for G·,i . To solve a tridiagonal

23
˜
linear system for one vector takes O(L) time, so solving for the full inverse matrix G takes
O(L2 ) time.
˜
Now we shall derive an algorithm to compute the elements of gi and G·,i . Begin with
the last column i = L: U +T gL = eL can be solved from the first to the last row and
˜
gives gL = eL . Next U + G·,L = (D + )−1 eL can be solved starting from the last row,
+
˜
˜
˜
GL,L = 1/DL . From there one can work upwards by backsubstitution, Gj,L = −Uj+ Gj+1,L
(j = 1, . . . , L − 1). For the other columns i < L, one cannot take the shortcut and has
˜
to solve both linear systems for gi and G·,i . But it is now important to realize that for
˜
˜
any column vector G·,i+1, if we somehow know the diagonal element Gi+1,i+1 , the next
element above the diagonal is
˜
˜
Gi,i+1 = −Ui+ Gi+1,i+1

(i = 1, . . . , L − 1) .

(C5)

˜
Thus we have a prescription how to go up one row in G. Together with the symmetry
˜
˜
˜
of G, that is Gi,i+1 = Gi+1,i , which follows from the symmetry of T , we get the first offdiagonal element one column to the left without solving the two linear systems in (C4).
So it is possible to compute directly the tridiagonal part of the inverse. However, there is
another algorithm which is much more accurate for near-singular matrices at the end of
˜
the RG flow: the double factorization.33 It does not rely on the symmetry of G but uses
the complementary “UDL” factorization
T = U − D − L− = U − D − U −T ,

(C6)

where the matrix elements are obtained as
−
DL := aL ,

−
Ui− := bi /Di+1 ,

−
Di := ai − bi Ui−

(i = L − 1, . . . , 1) .

(C7)

We proceed as for the LDU factorization above and get
−
˜
G1,1 = 1/D1
˜
˜
Gi,i+1 = −Ui− · Gi,i .

(C8)
(C9)

We can combine equations (C5) and (C9) to relate consecutive diagonal elements:
+
−
˜
˜
˜
˜
Gi+1,i+1 = −Gi,i+1 /Ui+ = Gi,i · Ui− /Ui+ = Gi,i · Di /Di+1 .

(C10)

Thus we start with (C8) and use the U − D − L− decomposition to go one matrix element
to the right in the inverse matrix, from the diagonal to the first off-diagonal (C9), while
the L+ D + U + decomposition allows to go down by one, back to the next diagonal element
(C5). There is no need to compute the full inverse matrix.
One can implement the algorithm without knowing the derivation by using equations
(C3) and (C7) - (C10). One can further eliminate the U’s using equations (C3) and
(C7) and implement the algorithm such that only the input vectors ai , bi and the output

24
˜
˜
vectors Gi,i , Gi,i+1 enter the temporary storage. This double factorization is numerically
accurate to more than 10 significant digits (using double precision) even for large lattices
(106 sites) and almost singular matrices with |ai | ∼ 10−15 which appear at the end of the
flow for half filling.
The right hand side of the flow equations (41) for the density response vertex RΛ can
be computed in O(L) time using the fact that the upper triangular part of the inverse
˜
of a tridiagonal matrix (here G) is the upper triangular part of the outer product of two
33,34
vectors.

1

T. Giamarchi, Quantum Physics in One Dimension (Oxford University Press, New York,
2003).

2

For a review on Luttinger liquids, see J. Voit, Rep. Prog. Phys. 58, 977 (1995).

3

A. Luther and I. Peschel, Phys. Rev. B 9, 2911 (1974).

4

D. C. Mattis, J. Math. Phys. 15, 609 (1974).

5

W. Apel and T.M. Rice, Phys. Rev. B 26, 7063 (1982).

6

T. Giamarchi and H.J. Schulz, Phys. Rev. B 37, 325 (1988).

7

C.L. Kane and M.P.A. Fisher, Phys. Rev. B 46, 15233 (1992).

8

S. Eggert and I. Affleck, Phys. Rev. B 46, 10866 (1992).

9

D. Yue, L.I. Glazman, and K.A. Matveev, Phys. Rev. B 49, 1966 (1994).

10

For a review on isolated impurities in Luttinger liquids, see Ref. [1], chapter 10.

11

For a recent comprehensive review of the DMRG, see U. Schollw¨ck, to appear in Rev. Mod.
o
Phys.

12

M. Salmhofer, Renormalization (Springer, Berlin, 1998).

13

First computational applications of the fRG to Fermi systems were performed for the 2D
Hubbard model, see D. Zanchi and H.J. Schulz, Phys. Rev. B 61, 13609 (2000); C.J. Halboth
and W. Metzner, ibid 61, 7364 (2000); C. Honerkamp, M. Salmhofer, N. Furukawa, and T.
M. Rice, ibid. 63, 035109 (2001).

14

V. Meden, W. Metzner, U. Schollw¨ck, and K. Sch¨nhammer, Phys. Rev. B 65, 045318
o
o
(2002); J. Low Temp. Phys. 126, 1147 (2002).

15

C. Wetterich, Phys. Lett. B 301, 90 (1993).

16

T.R. Morris, Int. J. Mod. Phys. A 9, 2411 (1994).

17

For generic interacting Fermi systems, a concise derivation of the full hierarchy of flow equations is presented by M. Salmhofer and C. Honerkamp, Prog. Theor. Phys. 105, 1 (2001).

18

V. Meden and U. Schollw¨ck, Phys. Rev. B 67, 035106 (2003).
o

19

V. Meden and U. Schollw¨ck, Phys. Rev. B 67, 193303 (2003).
o

20

V. Meden, S. Andergassen, W. Metzner, U. Schollw¨ck, and K. Sch¨nhammer, Europhys.
o
o
Lett. 64, 769 (2003).

25
21

V. Meden, T. Enss, S. Andergassen, W. Metzner, and K. Sch¨nhammer, cond-mat/0403655
o
(unpublished).

22

Vertex renormalization has already been applied in two letters on transport properties, Ref.
[20,21], but without any derivation of the method.

23

S. Andergassen, T. Enss, V. Meden, W. Metzner, U. Schollw¨ck, and K. Sch¨nhammer
o
o
(unpublished).

24

C.N. Yang and C.P. Yang, Phys. Rev. 150, 321 (1966); ibid 327 (1966).

25

F.D.M. Haldane, Phys. Rev. Lett. 45, 1358 (1980).

26

See, for example, T. Giamarchi and H. Maurey, in Correlated Fermions and Transport in
Mesoscopic Systems, edited by T. Martin, G. Montambaux, and J. Tran Thanh Van (Editions
Frontieres, Paris, 1997).

27

J.W. Negele and Orland, Quantum Many-Particle Systems (Addison-Wesley, Reading, 1987).

28

R. Egger and H. Grabert, Phys. Rev. Lett. 75, 3505 (1995).

29

See, for example, J. Zinn-Justin, Quantum Field Theory and Critical Phenomena (Oxford
University Press, New York, 1996).

30

The derivation will be given in an upcoming paper containing a detailed discussion of transport properties, see T. Enss, S. Andergassen, V. Meden, W. Metzner, U. Schollw¨ck, and K.
o
Sch¨nhammer (unpublished).
o

31

I.E. Dzyaloshinskii and A.I. Larkin, Zh. Eksp. Teor. Fiz. 65, 411 (1973) [Sov. Phys. JETP
38, 202 (1974); W. Metzner and C. Di Castro, Phys. Rev. B 47, 16107 (1993).

32

See, for example, W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery, Numerical
Recipes in Fortran 77 (Cambridge Univ. Press, 1986).

33

For a pedagogical review on the inverse of symmetric tridiagonal matrices, see G. Meurant,
SIAM J. Matrix Anal. Appl. 13, 707 (1992).

34

For more details, see T. Enss, Ph.D. thesis, University of Stuttgart, 2004 (unpublished).

26

Ë£

£

¦£

 £

FIG. 1: Flow equation for the selfenergy ΣΛ .


£

 £

 £

Ë £  £
£

Ë£
+

 £
¿

FIG. 2: Flow equation for the 2-particle vertex ΓΛ .

27
1

2


£

1’

2’

 £
1

1’

1

1’

2

2’

2

2’

1

2’

1

2’

2

1’

2

1’

FIG. 3: Flow equation for the 2-particle vertex at 1-loop level with particle-particle and particlehole channels written explicitly.

1.2
n = 1/2
1.1
n = 3/8
UΛ

1
0.9

n = 1/4

0.8
0.7
0.001

n = 1/8
0.01

0.1

1
Λ

10

100

1000

FIG. 4: Flow of the renormalized nearest neighbor interaction U Λ for the spinless fermion model,
for U = 1 and various densities n.

28
0.5

0.4

Σj,j

0.3
-0.12
Σj,j+1

-0.16
-0.2
-75

-50

-25
j − j0

0

25

FIG. 5: Self-energy near a site impurity of strength V = 1.5 for the spinless fermion model at
quarter-filling and interaction strength U = 1; the impurity is situated at the center of a chain
with L = 1025 sites.

exponent

1

0.9

0.8

0.7
101

102

103
j − j0

104

105

FIG. 6: Effective exponent for the decay of oscillations of Σj,j as a function of the distance
from a site impurity of strengths V = 0.01, 0.1, 0.3, 1, 10 (from bottom to top), for the spinless
fermion model at half-filling and interaction strength U = 1; the impurity is situated at the
center of a chain with L = 218 + 1 sites.

29
0.4
U =0
U =1

Dj0 −1

0.3

0.2

0.1

0
-3

-2

-1

0
ω

1

2

3

FIG. 7: Local density of states on the site next to a site impurity of strength V = 1.5 for spinless
fermions at half-filling and U = 1; the impurity is situated at the center of a chain with L = 1025
sites; the non-interacting case U = 0 is shown for comparison.

0.5
n = 1/4
n = 3/4

Dj0 −1

0.4
0.3
0.2
0.1
0
-4

-3

-2

-1

0
ω

1

2

3

4

FIG. 8: Local density of states on the site next to a site impurity as in Fig. 7 (same parameters),
but now for densities n = 1/4 and n = 3/4.

30
0.4
0.3

α

0.2
0.1
0
-0.1
102

103

104
L

105

106

102

103

104
L

105

106

0.4
0.3

α

0.2
0.1
0
-0.1

FIG. 9: Logarithmic derivative of the spectral weight at the Fermi level near a boundary (solid
lines) or hopping impurity (dashed lines) as a function of system size L, for spinless fermions
at quarter-filling and interaction strength U = 0.5 (circles) or U = 1.5 (squares); upper panel:
without vertex renormalization, lower panel: with vertex renormalization; the open symbols
are fRG, the filled symbols DMRG results; the horizontal lines represent the exact boundary
exponents for U = 0.5 and U = 1.5. In the boundary case (solid lines) the spectral weight has
been taken on the first site of a homogeneous chain, in the impurity case (dashed lines) on one
of the two sites next to a hopping impurity t′ = 0.5 in the center of the chain.

31

0.2

α

0.1
0
boundary
V =3
V =1
V = 0.3
V = 0.1

-0.1
-0.2
102

103

104
L

105

106

FIG. 10: Logarithmic derivative of the spectral weight at the Fermi level near a boundary (solid
line) or near a site impurity (dashed lines) as a function of system size L, for spinless fermions
at quarter-filling and interaction strength U = 1. In the boundary case the spectral weight has
been taken on the first site of a homogeneous chain, in the impurity case on the site next to a
site impurity of strength V in the center of the chain. The horizontal line represents the exact
boundary exponent for U = 1.

0.7
DMRG
RG

nj

0.6

0.5

0.4
0

20

40

60

80

100

120

j

FIG. 11: Density profile nj for a spinless fermion chain with 128 sites and interaction strength
U = 1 at half-filling; fRG results are compared to exact DMRG results.

32
0.1

|nj − n|

0.01

0.001

0.0001
10

100

1000
j

FIG. 12: Density modulation |nj − n| as a function of the distance from a boundary, for spinless
fermions with interaction strength U = 1 and (average) density n = 0.393 on a chain with 8192
sites; the broken straight line is a power-law fit to the envelope of the oscillation amplitudes.

33
10−1
10−2

amplitude

10−3
10−4
10−5
10−6
10−7
10−8 1
10

U
U
U
U
102

= 0.1
= 0.5
= 1.0
= 1.5
103

104

105

104

105

j

1

exponent

0.9
0.8
0.7
0.6
0.5 1
10

102

103
j

FIG. 13: Amplitude (envelope) of oscillations of the density profile nj induced by a boundary
as a function of the distance from the boundary, for spinless fermions with various interaction
strengths U at half-filling; the interacting chain with 219 + 1 sites is coupled to a semi-infinite
non-interacting lead at one end (opposite to the boundary); upper panel: log-log plot of the
amplitude, lower panel: effective exponents for the decay, and the exact asymptotic exponents
as horizontal lines.

34

exponent

0.7

0.6

0.5

0.4 1
10

102

103
j − j0

104

105

FIG. 14: Effective exponent for the decay of density oscillations as a function of the distance
from a site impurity of strengths V = 0.01, 0.1, 0.3, 1, 3, 10 (from bottom to top); the impurity
is situated at the center of a spinless fermion chain with 218 + 1 sites and interaction strength
U = 1 at half-filling; the interacting chain is coupled to semi-infinite non-interacting leads at
both ends.

exponent

1.8

1.6

1.4

1.2 1
10

102

103
j − j0

104

105

FIG. 15: Effective exponent for the decay of density oscillations as a function of the distance
from a site impurity of strengths V = 0.1, 1, 10, 100, 1000 (from top to bottom) for the same
chain as in Fig. 14 but now with an attractive interaction, U = −1.

35
1.6
0.05

∆Kρ

1.4
0

Kρ

1.2
-0.05

1

-1

0

1

2

U

0.8
0.6
-1

-0.5

0

0.5
U

1

1.5

2

FIG. 16: Luttinger liquid parameter Kρ as a function of U at various densities (as in Fig. 4) for
the spinless fermion model; the inset shows the difference between the RG result and the exact
Bethe ansatz result for Kρ .