The large system asymptotics of persistent currents in mesoscopic quantum rings
A. Gendiar1,2 , R. Krcmar1 , and M. Weyrauch3

arXiv:0902.4385v3 [cond-mat.str-el] 11 May 2009

1

Institute of Electrical Engineering, Slovak Academy of Sciences, SK-841 04, Bratislava, Slovakia
2
Institute for Theoretical Physics C, RWTH University Aachen, D-52056 Aachen, Germany
3
Physikalisch-Technische Bundesanstalt, D-38116 Braunschweig, Germany
(Dated: May 29, 2018)

We consider a one-dimensional mesoscopic quantum ring filled with spinless electrons and threaded
by a magnetic flux, which carries a persistent current at zero temperature. The interplay of Coulomb
interactions and a single on-site impurity yields a non-trivial dependence of the persistent current on
the size of the ring. We determine numerically the asymptotic power law for systems up to 32 000
sites for various impurity strengths and compare with predictions from Bethe Ansatz solutions
combined with Bosonization. The numerical results are obtained using an improved functional
renormalization group (fRG) method. We apply the density matrix renormalization group (DMRG)
and exact diagonalization methods to benchmark the fRG calculations. We use DMRG to study
the persistent current at low electron concentrations in order to extend the validity of our results
to quasi-continuous systems. We briefly comment on the quality of calculated fRG ground state
energies by comparison with exact DMRG data.
PACS numbers: 71.10.Pm, 73.23.Ra, 73.63.-b, 05.10.Cc

I.

INTRODUCTION

Quantization of the magnetic flux was first observed
in superconducting cylinders [1, 2]. In normal electronic
systems the flux quantum φ0 = hc/e determines the period of the persistent currents, which can be observed in
mesoscopic one-dimensional disordered metallic rings [3].
Persistent currents in metallic rings are induced only if
the circumference of the ring is not larger than the electron phase coherence length and/or the electron’s meanfree path. Such conditions can be easily obtained for
sufficiently low temperatures at which inelastic scattering from phonons is suppressed. Normally, this happens
at temperatures below 1 K and ring sizes of about 1 µm.
Impurities (disorder) in such metallic rings act as elastic
scatterers, and it is not surprising that persistent currents
can flow despite the non-zero resistance of the rings.
Interest in persistent current phenomena remains
strong both experimentally [4, 5, 6, 7, 8, 9, 10, 11] and
theoretically [12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25]. The main reason for this is the fact that there
still is a discrepancy between experiment and theory. Experiments [4, 8] show persistent currents which are two
orders of magnitude larger than theoretical predictions
for rings filled with non-interacting electrons [13, 19, 26].
If electron-electron interactions are included, persistent
currents increase for repulsive as well as attractive electron interactions [14, 16, 27], however, the calculated currents were still about five times smaller than the experimental data. Recently, Bary-Soroker et al. [24] provided
an explanation for the magnitude of the persistent currents if magnetic impurities and the attractive interactions are considered.
Considerable theoretical effort has been invested into
the question whether the interplay between strong
electron-electron interactions, electron density, and disorder strength can explain the large persistent currents

observed experimentally. Disorder usually gives rise to
elastic scattering and consequently should reduce the
persistent current. Bouzerar et al. [15] claim that both
electron-electron interactions and disorder decrease the
persistent current, whereas the self-consistent HartreeFock calculations by Cohen et al. [18] suggest that the
current may be enhanced by a moderate disorder due to
screening effects. Other theoretical studies conclude that
electron-electron interaction could reduce or enhance the
current depending on the type of disorder present in the
system [14, 28, 29]. Abraham and Berkovits [14] conjectured that the increase would be rather small and onedimensional spinless models at half-filling cannot explain
the results of single ring experiments.
Among the theoretical methods used to study persistent current phenomena are Hartree-Fock approximations [17, 18], configuration interaction and quantum
Monte Carlo simulations [30, 31], Bosonization [16], conformal field theory [19, 20], the functional renormalization group (fRG) [21], and the density matrix renormalization group (DMRG) [21, 28, 29].
In the present paper we study the interplay between
electron-electron interactions and a single impurity in a
one-dimensional system. For one-dimensional interacting systems, the presence of a single impurity is known
to affect their physical properties [32, 33, 34, 35]. We
restrict our investigations to the Luttinger liquid phase
which is characterized by a power law decay of correlation
functions in a sufficiently large system [36]. Our aim is to
investigate the asymptotic behavior of the persistent current I as a function of the ring length L, from which the
power law exponent may be obtained I ∝ L−αB −1 . The
exponent αB = K −1 − 1 is related to the Luttinger liquid
parameter K [34, 35, 36] and depends on the electronelectron interaction, but it does not depend on the impurity strength if L → ∞.
It is known that the asymptotics of the currents is

2

II.

MODEL FOR SPINLESS FERMIONS

We consider a quantum ring of interacting spinless electrons at zero temperature. The ring is pierced by a magnetic flux φ and contains one single impurity. This system
is modelled by the tight-binding Hamiltonian
L

L

c† cℓ+1 e−iφ/L
ℓ

H = −t
ℓ=1

nℓ nℓ+1 +V n1 ,

+ h.c. +U

-1.985
-1.990

ε2

-2.000

ε0

-1.985

-2.000

ε1

V = 0.1

ε4
ε3

-1.990
-1.995

V=0

ε4
ε3

-1.995

ε0 , ε1 , ε2 , ε3 , and ε4

reached for smaller system sizes if sufficiently strong impurities are considered [16, 37, 38]. Nevertheless, it is
necessary to investigate relatively large systems in order to reach the asymptotics. We therefore need suitable many-body techniques which enable such calculations. We have chosen the functional renormalization
group [39, 40, 41, 42, 43], which proved to be a rather
useful tool in a previous study [21]. However, for numerical reasons the asymptotic regime could not really
be reached in that work. In the present paper we improve the fRG technology such that system sizes up to
about 32000 lattice sites can be investigated. In order
to benchmark the fRG results we use DMRG calculations [44, 45, 46, 47].
We also study the asymptotic power law decay of the
persistent current at low electron concentrations in order
to clarify differences between discrete lattice models (investigated in this paper) and continuous systems studied
by quantum Monte Carlo and configuration interaction
methods [30, 31]. For this reason, we apply DMRG to
our model away from half-filling (n = 1 ), in particular,
2
we study the case of quarter-filling (n = 1 ) as well as
4
1
1
n = 8 and n = 16 .
The paper is organized as follows. In Sec. 2 we define the model for our calculations and briefly review the
theory of the persistent current. We introduce the fRG
method in Sec. 3 and related Appendices. Section 4 is
devoted to the DMRG method with a short discussion of
its limitations. Results are presented in Sec. 5 gathering
the numerical data obtained by fRG, DMRG, and exact
diagonalization. We calculate the persistent current and
extract the effective exponent αeff related to the power
law decrease with respect to the ring size L. Finally we
show DMRG results for systems with low-filling. The
results are summarized in Sec. 6.

ε2
ε0

-1.985

ε4

-1.990

ε1

ε3

ε2
-1.995 ε
1
ε0
-2.000
-2π

V=1

-π

0

φ

2π

FIG. 1: (Color online) The five lowest-laying single-particle
energy levels for a non-interacting system as functions of the
magnetic flux φ for a ring of size L = 128. The three panels
correspond to V = 0, V = 0.1, and V = 1.0, respectively.

with strength V is placed on the lattice site ℓ = 1 along
the ring (ℓ = 1, 2, ..., L). Periodic boundary conditions
are imposed by L + 1 ≡ 1.
Without electron-electron interaction (U = 0) and impurity (V = 0) the single particle energy levels ǫm calculated from the Hamiltonian Eq. (1) depend quadratically on the magnetic flux φ. The persistent current
Im = −vm /L at energy level m with the wave vector
k = 2π (φ + m) /L is calculated from the electron velocity vm ,
vm (φ) = 2π

∂ǫm (φ)
∂ǫm(φ)
=L
.
∂k(φ)
∂φ

(2)

Summation over all occupied energy levels then yields
the total persistent current at temperature T = 0,

ℓ=1

(1)
written in terms of the electron creation and annihilation operators c† and cℓ as well as the electron density
ℓ
operator nℓ = c† cℓ .
ℓ
The Peierls factor e−iφ/L with φ = 2πΦ/Φ0 describes
the influence of the magnetic field in terms of the magnetic flux Φ [48]. The flux quantum Φ0 = hc/e is set to
one in the following. The hopping amplitude t between
neighboring sites is also set to one in order to define the
energy scale. The nearest-neighbor Coulomb interaction
is denoted by U . An (external) on-site impurity potential

π

mmax

I(φ) = −
m=0

∂ǫm (φ)
.
∂φ

(3)

The upper panel in Fig. 1 shows the parabolic energy
band structure (V = 0, U = 0). Consequently, the persistent current is proportional to the magnetic flux φ and
inversely proportional to the length L of the ring. For an
odd number of electrons in the ring one finds
I(φ) = −

vF
φ,
πL

−π ≤ φ < π ,

(4)

3
whereas for an even filling the current is given by

 φ+π,
vF
I(φ) = −
×
πL  φ − π ,

-63.85
-63.90

−π ≤ φ < 0,
(5)

2

T (kF ) =

4t
,
+V2

4t2

(6)

∂E0 (φ)
,
∂φ

I ∝ L−1−αB .

(8)

The exponent αB is a function of the electron-electron interaction (regardless of the impurity strength V ). Only
at half-filling, the exponent αB can be obtained analytically from a Bethe Ansatz solution [21, 36, 50]
2
U
arccos −
π
2

-63.45

E4

-63.50

E3

-63.55
-63.60

E2

V=1

E1
E0

-63.65
-63.25
-63.30
-63.35

V=2

-63.45
-2π

-π

0

φ

π

2π

FIG. 2: (Color online) The five lowest-laying many-body energy levels obtained by DMRG as a function of the magnetic
flux φ for a ring of size L = 128 at half-filling and electronelectron interaction U = 1. Upper panel: V = 0, center panel:
V = 1, lower panel: V = 2.

(7)

where E0 (φ) is the many-body ground-state energy at
zero temperature as plotted in Fig. 2. In general the
system is a Luttinger liquid. However, at half-filling it is
in the Luttinger liquid phase only for |U | ≤ 2. At U = 2,
the system exhibits a charge-density wave instability, and
phase separation characterizes the system for U < −2.
In this paper we focus on the Luttinger liquid phase
for which Bosonization techniques and additional approximations predict that asymptotically (i.e. for large L)
the persistent current decays algebraically with increasing system size,

αB =

-64.00

-63.40

which holds at half-filling in the non-interacting
limit [49]. With this relation, the strengths of various
impurities V used in this paper may be compared.
If the electron-electron interaction U is switched on,
the many-body energy levels are shown in Fig. 2 for
U = 1 and various impurity strengths V . The persistent currents are obtained from the relation
I(φ) = −

E0 , E1 , E2 , E3 , and E4

0≤φ<π

with vF being the Fermi velocity.
An impurity in the ring rounds off the cusps seen in
this function as shown in the central and lower panels of
Fig. 1 for a weak and an intermediate impurity, respectively. For strong impurities and odd filling, one finds
I ∝ −vF sin(φ)/πL [16]. The coupling strength V of an
on-site impurity can be related to the physically more relevant transmission coefficient at the Fermi wave vector
kF

V=0

-63.95

−1.

(9)

Throughout this paper we often consider the case U =
1
1 with the corresponding αB = 3 . We determine αB
from the Hamiltonian Eq. (1) using three different manybody techniques: fRG, DMRG, and exact diagonalization
(ED).

III.

FUNCTIONAL RG SCHEME FOR
INTERACTING FERMIONS

A scheme for one-dimensional Fermionic systems based
on the functional renormalization group, has been developed by Meden, Metzner, Sch¨nhammer, and collaborao
tors [51] based on seminal work by Wetterich [40] and
Morris [41]. It has been applied to study transport properties of quantum wires [51, 52, 53] and rings [21].
The fRG scheme of the present work is similar to the
one applied in the references cited above, but uses a
somewhat different truncation procedure. Details of our
method have been presented recently [54]. Here we review the essential steps relevant for a system of spinless
interacting Fermions.
The effective average action Γk [40] for interacting
Fermions evolves according to the fRG equation [40, 55]
∂
1
∂Rk
(2)
Γk [φ∗ , φ] = − Tr [Γk [φ∗ , φ] + Rk ]−1
∂k
2
∂k

. (10)

The system of L spinless electrons is described by an
L-component vector of Grassmann variables φ(τ ) =
(φ1 (τ ), . . . , φL (τ )), which describes the evolution of the
interacting Fermionic system in imaginary time τ . The
functional derivative of Γk with respect to φ∗ and φ is de(2)
noted by Γk . The trace in Eq. (10) is to be performed
over the Fermionic states j = 1, . . . , L. The regulator

4
Rk is introduced in order to suppress thermal and quantum fluctuations at energy or momentum scales k larger
than any physical scale relevant for our problem. With
decreasing k, the regulator gradually “switches on” such
fluctuations until they are fully included at k = 0, i.e.
at k = 0 the regulator Rk vanishes. As a regulator we
choose Rk = Ckθ(k 2 − ω 2 ) with C a large constant and
it satisfies all requirements for a useful regulator as discussed in more detail in Ref. [54]. It has an additional
advantage that the integration over the frequencies ω can
be done analytically.
In order to solve the functional differential equation (10) a particular functional form for Γk in terms
of the Grassmann variables φ∗ , φ is assumed
L

β

∂
φα (τ ) + Uk (φ∗ (τ ), φ(τ )) ,
∂τ
0
α=1
(11)
where the ‘effective potential’ Uk does not depend on
derivatives of the Grassmann variables with respect to
the imaginary time τ , i.e. it is represented as a Grassmann polynomial in terms of φ∗ and φ
Γk [φ∗ , φ] =

φ∗ (τ )
α

dτ

L

ak,jj φ∗ φj + ak,j,j+1 φ∗ φj+1
j
j

Uk (φ∗ , φ) = ak,0 +
j=1

+ ak,j,j−1 φ∗ φj−1
j
L

φ∗ φj φ∗ φj+1
j
j+1

+ Uk

(12)

j=1

with the understanding that φ0 = φL and φL+1 = φ1 .
Of course, in general this polynomial should have higher
order terms, but those are neglected in the hope they are
small.
For half-filling, the density-density interaction strength
Uk renormalizes according to
Uk =

U
1+

U
2π

k−

2+k2
√
4+k2

(13)

as was suggested in Ref. [50]. This result can be easily derived within the fRG scheme applied here by using
analogous approximations for the two-particle vertex as
in Ref. [50]. (The interaction parameter U is given in the
Hamiltonian Eq. (1)).
Inserting this Ansatz for the effective average action on
both sides of the flow equation, performing the integration over the frequencies ω, and comparing terms with
the same Grassmann structure, it is straightforward to
obtain the flow equations for the coefficients ak,0 and
ak,jj ′
a′
k,0 =
a′
k,jj

1
2π

Uk
=
2π

+

eiλ0 ln det
λ=±k
+

1 −1
g (iλ) ,
iλ k

eiλ0 gk,(j+s)(j+s) (iλ),
λ=±k s=±1

a′
k,j(j±1) = −

Uk
2π

+

eiλ0 gk,j(j±1) (iλ)

(14)

λ=±k
−1

with gk (iλ) = (ak + iλ1) and ak = (ak,jj ′ ). The con+
vergence factor eiλ0 is only needed at the beginning of
the flow from infinity down to a large constant k0 .
At k = ∞ the effective average action Γk is given by
the classical action S [40, 55], which is determined by the
Hamiltonian Eq. (1). Therefore, the initial conditions at
k = ∞ for the solution of the flow equations (14) are
directly obtained from Eq. (1),
a∞,0 = 0,
a∞,jj = −V δj1 + µδjj ,
a∞,j(j±1) = e

∓iφ/L

(15)

,

with the understanding that aL,L+1 = aL,1 and a1,0 =
a1,L . For convenience we have added a chemical potential
term −µnℓ to the Hamiltonian. Our initial conditions
differ from those of Ref. [21]: in the present work the
dependence on the magnetic flux φ resides in the initial
conditions for the off-diagonal ajj ′ only, while in Ref. [21]
this dependence partly resides in the free propagator, corresponding to a different definition of the kinetic energy.
Moreover, here we also include the chemical potential
into the ajj ′ and not in the free propagator.
In practice we must start the renormalization flow at
a finite k = k0 . The initial flow from k = ∞ to k = k0
is obtained from an analytical solution of the Eqs. (14).
For convenience this calculation is briefly reviewed in Appendix A with the result
V
L U
+
−µ ,
2
2 2
= (µ − U )δjj − V δj1 ,

ak0 ,0 =
ak0 ,jj

ak0 ,j(j±1) = e

∓iφ/L

(16)

.

At half-filling the chemical potential is given by µ = U
because of particle-hole symmetry. From the ground
state energy E0 (φ) = a0 (k = 0) + µ n , the persistent
current is calculated using Eq. (3).
For given L, Eq. (14) represents a (large) set of 3L − 1
non-linear coupled complex-valued ordinary differential
equations. The calculation of the right-hand side of these
equations requires the inversion of a (potentially) large
cyclic tridiagonal matrix as well as the computation of
the logarithm of the determinant of such a matrix at
each step of the numerical integration of the differential
equations. To accomplish this in an efficient manner is
described in some detail in Appendix B.
IV.

DMRG

The density matrix renormalization group is a numerical technique for the diagonalization of the very large
matrices typically encountered in quantum many-body

5
calculations. The technique is described in detail in
Refs [44, 45, 46, 47]. We use DMRG for the calculation of
the ground-state energy E0 (φ). In our case, the matrices
to be diagonalized are complex-valued due to the Peierls
factor e−iφ/L entering Eq. (1). Treating such complexvalued systems by DMRG does not lead to numerical
instabilities, apart from a few rare cases which occur at
very low electron fillings. Here, however, the numerical
diagonalization routines for complex matrices experience
convergence difficulties for very large system sizes. Also,
the superblock diagonalization routines seldom require
more than 104 cycles for convergence as compared to the
standard ∼ 102 cycles typically necessary for real-valued
matrices. The memory requirements increase by a factor
of 1.8 and the calculation time raises by a factor of 2.1
compared to a standard real-valued DMRG.
We kept the DMRG truncation error as small as
ε < 10−9 for system sizes L ≥ 128, for smaller systems
ε ≈ 10−15 could be achieved. The number of states kept
are left to vary such that the above truncation error condition could be satisfied. It is well known that the efficiency and accuracy of DMRG decreases substantially if
periodic boundary conditions are imposed. We checked
the accuracy of our results comparing with data obtained
from the exact diagonalization of the Hamiltonian for systems with up to L = 24 sites at various fillings.
We would like to point out that the calculations for
strong impurities V > 10 and for large rings with L > 100
require an extremely careful treatment because differences between the ground state energies for different magnetic flux values are extremely small |E0 (φ = 0)−E0 (φ =
π)|/|E0 (φ = 0)| < 10−6 . For this reason, both the
number of states kept and the number of the finite-size
method sweeps must be sufficiently large.

V.

the currents
∞

I(φ) =

Ik sin(kφ)

(17)

k=1

can be calculated without numerical differentiation using
an integration by parts
Ik =

2k
π

π

dφE0 (φ) cos(kφ).

(18)

0

Nevertheless, small numerical errors in the calculation of
the ground state energy may significantly affect the calculation of the power law exponent αB (cf. Eq. (8)). The
Fourier analysis for strong impurity strengths in Eq. (17)
is well justified due to the sine-like shape of the current,
and the higher coefficients Ik , k ≥ 2 decay to zero with
increasing L faster than I1 . However, if the impurity
V ≪ 1, the first Fourier coefficient I1 may not suffice
to characterize the decay of the persistent current accurately.
For system sizes L ≤ 24 and for impurities V < 100,
ED and DMRG yield essentially exact ground state energies in mutual agreement. The ground state energy
calculated from fRG differs from the exact results of ED
and DMRG at large interactions, as shown in Fig. 3 for
interactions within the Luttinger liquid regime. This deviation signals the drastic approximations involved in the
fRG method used here, in particular, the discarding of
higher order vertex functions. However, the shape of the
ground state energy as a function of the magnetic flux
can be reproduced quite well by this method as can be
seen from Fig. 4. Therefore, we can use fRG calculations
in order to obtain the persistent currents for systems as
large as L = 3 · 104. For very strong impurities (V ≫ 10)
DMRG and fRG may run into numerical difficulties and
we analyzed such cases by ED as is discussed in more
detail below.

RESULTS
A.

In this section we present and analyze numerical data
obtained by three different methods: fRG, DMRG, and
ED. We study the dependence of the persistent currents
on the ring size L, the impurity strength V , the magnetic
flux φ, and on the electron concentration n. Varying the
electron-electron interaction U within the Luttinger liquid regime −2 ≤ U ≤ 2 does not change physical properties qualitatively as we shall see below. However, the
weaker the electron-electron interaction the larger is the
system size L needed to reach the asymptotic power law
decay of the persistent current [21, 30].
The ground state energies E0 (φ) are calculated for 0 ≤
φ ≤ π and extended to the first Brillouin zone −π ≤ φ ≤
π using the reflection symmetry of the energy around the
origin φ = 0 (cf. Figs. 1 and 2). The persistent current in
Eq. (7) obtained by DMRG and ED is calculated using
numerical differentiation. The Fourier coefficients Ik of

fRG and DMRG analysis at half-filling

Figure 4 shows persistent currents as functions of the
1
flux calculated at half-filling (n = 2 ) using DMRG and
fRG for various ring lengths L and for an intermediate
impurity strength V = 2 (corresponding to transmission
1
T = 2 ). The data indicate that for small and intermediate ring lengths L, the DMRG and fRG calculations
agree rather well. DMRG calculations are hardly feasible for L > 256, and only fRG data are available in order
to study the large L limit.
Bosonization theory including additional approximations predict a power law decay of the persistent current at large L, cf. Eq. (8). In Fig. 5 we compare this
asymptotic decay (dashed lines) with the numerically determined first Fourier coefficient of the persistent current
LI1 in a logarithmic plot. From Bosonization one expects
that LI1 decays asymptotically as L−αB . The open circles represent the fRG data and the filled triangles the

6
4

-5
10

8

16

32

64

128

256

512

1024 2048 4096 8192 16384 32768

V = 0.1

0

V = 0.5
V = 1.0

-15

10

DMRG
ED
fRG

-20

-1

V = 2.0

L I1

E0 (φ = 0)

-10

10

V = 10.0

-2

-25
-30
-2.0

-1.5

-1.0

-0.5

0.0

0.5

1.0

1.5

2.0

10

V = 100.0

-3

U

10

1

10

2

10

3

10

4

L
FIG. 3: (Color online) Ground state energies at φ = 0 calculated by DMRG (filled triangles), ED (open circles), and fRG
(crosses) as functions of U within the Luttinger liquid regime
for L = 16 and V = 1.
0.8

FIG. 5: (Color online) The first Fourier coefficient I1 versus
the length of the ring L for U = 1 and various V . The open
circles and the filled triangles, respectively, correspond to the
fRG and DMRG data. The dashed straight lines show the
1
asymptotics αB = 3 . The upper x-axis shows the exact number of the sites L at which the calculations were performed.

0.4

LI

0.2
0.0

L=
4
L=
8
L = 16
L = 32
L = 64
L = 128
L = 256

-0.2
-0.4
-0.6
-0.8
-π

-π/2

0

φ

π/2

π

FIG. 4: (Color online) Persistent currents calculated by fRG
(lines) and DMRG (symbols ×) at half-filling for U = 1,
V = 2, and various ring lengths L. The currents with the
largest amplitude correspond to L = 4 and then the amplitudes decrease with increasing L.

DMRG results. Again we see the good agreement between both methods, and it is obvious that the asymptotic decay does not depend on the impurity strength (if
at least V > 0.1), which ranges from rather weak V = 0.1
(T = 0.9975) to very strong V = 100 (T = 0.0004). The
asymptotic regime is reached at smaller L for stronger
impurities.
Notice a tiny deviation of the fRG data from the
dashed lines with the common power-law exponent αB
at large L. This is due to the approximations involved in
the fRG method and is known from studies of quantum
wires [22]. However, there are also numerical limitations:
for large systems the differences between currents calculated for different magnetic fluxes are so tiny that they

αeff = − ∂ log10 (L I1) / ∂ log10 (L)

0.6
V = 10.0

0.4

V = 2.0

0.3
V = 100.0
V = 1.0

0.2

V = 0.5

0.1
V = 0.1

0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1 / log10 (L )
FIG. 6: (Color online) The effective exponent αeff versus the
inverse of log10 (L) for the interaction U = 1 and various impurities V . The open and filled symbols correspond to fRG
and DMRG data, respectively.

can not be resolved by the floating point data representation of the computer, and the current cannot be calculated reliably. This numerical limitation also prevents us
to calculate even larger systems.
In order to quantify the deviation of the numerical data
from the expected power law, we define an effective exponent
αeff = −

∂ log10 (LI1 )
,
∂ log10 (L)

(19)

which is shown in Fig. 6 as a function of 1/ log10 (L).

7
1.8

0.50

1.6

V L | I1 |

1.4
1.2

V= 1
V= 2
V= 4
V = 10
V = 16
V = 20
V = 100
V = 10000

0.45
0.40

αeff

V=1
V=2
V=4
V = 10
V = 16
V = 20
V = 100
V = 10000

0.35
0.30

1.0

0.25
0.20

0.8
4

6

8

10

12

14

16 18

0.15
0.0

0.2

0.4

0.6

FIG. 7: (Color online) The persistent current versus L for a
variety of V ’s obtained by ED at U = 1 and half-filling.

0.8

1.0

1.2

1.4

1.6

1 / log10 (L )

L

FIG. 8: (Color online) Dependence of the effective exponent
on 1/ log10 (L) for the same parameters as in Fig. 7.

V
1 2 4 10

From our data analysis it emerges that the effective
exponents for larger L are also non-negligibly influenced
by the numerical limitations discussed above, which are
seen in the figure as small irregularities in the approach
of the numerical data to saturation. Notice that we are
dealing with extremely tiny effects observed in the exponent αeff . Thus, extremely accurate calculations of
the ground state energy is the key to the analysis of the
effective exponent. There is another feature of the effective exponent worth mentioning: it shows a characteristic
minimum at lengths 8 < L < 100 within the range of V .
∼ ∼
The effective exponent determined from the DMRG
data (full symbols) is systematically below the fRG results for all but the very small systems. Since the DMRG
data are expected to be accurate, one may conjecture
that αeff approaches αB = 1 for L → ∞, if DMRG cal3
culations for such large systems would be feasible.

2

10

3

10

4

10

5

10

6

60

∆α = (αB - αeff ) / αB [%]

The data labeled by the open and filled symbols are obtained from the fRG and DMRG results shown in Fig. 5
using a numerical derivative. If the impurity is weak
(V = 0.1), the effective exponent αeff deviates significantly from αB = 1 (the horizontal dashed line). It is ex3
pected that substantially larger systems would need to be
considered to approach the asymptotics in this case. For
V ≥ 0.5, the numerically determined exponents tend to
saturate at αeff ≈ 0.35 which is about 5% larger than the
1
expected αB = 3 . Our result is in good agreement with
the leading U -behavior, where the correction to the exponent for the quantum wire is known to be U = 0.318 [22].
π
As already stated, we attribute this deviation to the approximations involved in the fRG. Similar results were
obtained in a study on open chains [53], where the same
power law decay is extracted from the decay of the spectral weight near the impurity site.

10

40
20
0
L= 6
L= 8
L = 10
L = 12
L = 18

-20
-40
-60

0

1

2

3

4

5

6

log10 (V )
FIG. 9: (Color online) Relative difference between the effective and asymptotic exponents as a function of the impurity
strength V for several sizes L.

B.

Non-monotonous behavior of the effective
exponent

Another interesting feature can be seen in Fig. 6.
Both fRG and DMRG data show that the result for the
strongest impurity V = 100 appears to be out of sequence
if compared to the other V ’s. To elucidate this further,
we performed a series of calculations using ED at smaller
sizes L. Figure 7 shows the first Fourier coefficient of the
persistent current I1 as a function of L on a logarithmic
scale for on-site impurities in the range 1 ≤ V ≤ 104 . In
order to gather the data, the Fourier coefficient is scaled
with V L in Fig. 7. Furthermore, we plot the absolute
value of the Fourier coefficient in order to remove the

8

1.1
1.0

10

0

0.8
10

0.7

U=0
U=1
U=2
U=0
U=1
U=2
U=0
U=1
U=2

0.6
0.5
0.4

(n = 1/2)
(n = 1/2)
(n = 1/2)
(n = 1/4)
(n = 1/4)
(n = 1/4)
(n = 1/8)
(n = 1/8)
(n = 1/8)

10

6

8

10 12

16

V = 2,
V = 2,
V = 2,
V = 2,
V = 10,
V = 10,
V = 10,
V = 10,

-2

10

4

-1

L I1

| LI | (φ = π / 2)

0.9

n = 1/2
n = 1/4
n = 1/8
n = 1/16
n = 1/2
n = 1/4
n = 1/8
n = 1/16

-3

4

8

16

L
FIG. 10: (Color online) Persistent currents as a function of
the system size L for magnetic flux φ = π , various electron
2
concentrations n, and interactions U . Plotted is the absolute
value of the persistent current in order to remove the sign
change between even and odd electron fillings.

32

64

128

256

L

20 24

FIG. 11: (Color online) Persistent currents versus system size
L for U = 1 at various electron fillings n. The open and
full symbols correspond to the impurity V = 2 and V = 10,
respectively. The slopes of the dashed lines coincide with the
asymptotic exponents αB .
0.6

C.

DMRG study at low fillings

Studies of the persistent current at low fillings try to
address the question of differences between discrete lattice models and continuous models of one-dimensional
systems. Such questions have been addressed in a recent comprehensive numerical study [30, 31]. Discrete
models mimic continuous electron gas models if the electron concentration n is such that the cosine dispersion of
the discrete system can reasonably well approximate the
parabolic dispersion of the electron gas. Moreover, lat-

αeff

0.4

V=
V=
V=
V=

2.0
2.0
2.0
2.0

n = 1/2
n = 1/4
n = 1/8
n = 1/16

V = 10.0
V = 10.0
V = 10.0
V = 10.0

n = 1/2
n = 1/4
n = 1/8
n = 1/16

0.3
0.2
0.1
0.0
0.6
0.5
0.4

αeff

even-odd effect discussed in Section II. The straight dotdashed line corresponds to αB = 1 .
3
Figure 8 shows corresponding calculations of the effective exponent αeff for the same set of impurity strengths.
The dotted lines only serve as guides to the eye and
should not be taken as extrapolations. Again, it is obvious that the data for the cases V ≥ 102 appear to be
out of sequence as in Fig. 6. In Fig. 9 we plot the relative difference ∆α = (αB − αeff )/αB in order to quantify
the dependence of the effective exponent on the impurity
strength V for several system sizes L. A minimum of
∆α appears at V ≈ 4 for U = 1. We, therefore, identify
three regions in the figure (for reachable system sizes):
in the first region, 0 < V < 4, we observe the stan∼
dard behavior, i.e., the stronger the impurity, the faster
the convergence of αeff to αB . In the the second region
which starts at around V ≈ 4 and ends around V = 103 ,
this behavior is reversed. In the third region V > 103
the effective exponent αeff saturates and does not change
any more with increasing impurity strength.

0.5

0.3
0.2
0.1
0.0
0.0

0.2

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1 / log10 (L )
FIG. 12: (Color online) The effective exponents versus the
inverse of log10 (L) for the data displayed in Fig. 11. The horizontal dashed lines correspond to the asymptotic exponents
αeff (n) with n as listed in the legend.

tice models violate Galilean invariance [56]. In a Galilean
invariant model without impurity (V = 0) the persistent
current should not depend on the interaction U .
First we check whether the interaction dependence of
the current really disappears at lower fillings: Figure 10
shows persistent currents as a function of system size
L for U = 0, U = 1, and U = 2 at various fillings
1
1
n = 1 , 4 , and 8 . The currents at half-filling (full sym2

9
bols) strongly depend on the electron-electron interaction, while the electron gas is better approximated at
quarter-filling (shaded symbols). A negligibly small de1
pendence on U is found at n = 8 (open symbols) where
Galilean invariance is almost satisfied.
Figure 11 shows the decay of the persistent current
as a function of the ring size L for electron concentra1 1
1
tions n = 1 , 4 , 8 , and 16 on a logarithmic scale for
2
V = 2 (open symbols) and V = 10 (full symbols). The
asymptotic values αB (n) expected from the Bosonization
analysis [36, 50, 57] are indicated by the dashed straight
lines corresponding to exponents αB (n = 1 ) = 1 ,
2
3
1
1
αB ( 4 ) = 0.18383, αB ( 1 ) = 0.08825, αB ( 16 ) = 0.042940.
8
On the scale of the figure the data seemingly reach the
asymptotics. A more detailed picture of how the effective
power-law exponents αeff approach the expected asymptotic values αB is shown in Fig. 12. The upper and lower
panels show the effective exponent calculated for V = 2
and V = 10, respectively, at various fillings n. It is obvious that the αeff (n) do not converge to αB (n), and
substantially larger system sizes L would need to considered in order to see convergence. From our data we,
therefore, cannot support the suggestion [30, 31] that the
asymptotic regime can be reached at a small number of
electrons nL in continuous models.

VI.

SUMMARY

We studied the asymptotic power law decay of persistent currents using three different numerical methods
(fRG, DMRG, ED). To this end we improved the fRG
method so that systems with up to 3 · 104 sites could be
treated. We used DMRG and ED in order to benchmark
the fRG results and confirmed a sufficiently good agreement between fRG, DMRG, and ED calculations. We
confirmed that the fRG is a suitable method for extracting the asymptotic Luttinger power law exponent αB .
We found that αB which describes the decay of the
persistent current is about 5% overestimated by the fRG
results compared to expectations from Bosonization and
additional approximation. This deviation we attribute
to the approximation of the fRG method, and the effective exponent is known to be correct to the leading order
in the electron-electron interaction from the calculations
of the same exponent for quantum wires. The accurate
DMRG analysis cannot, unfortunately, be extended to
sufficiently large systems. However, there are indications
that the expected asymptotic values could be reached, if
larger systems could be calculated. It is found that the
effective exponent for a fixed impurity strength V shows
a typical minimum for lengths 8 < L < 100 at half-filling.
∼ ∼
Generally, it is known that the stronger the impurity,
the faster the asymptotic regime is reached. However,
we identified three regions for U = 1: in the first region,
0 < V < 4, we observed the standard behavior. (The
∼
case V = 4 corresponds to T = 1 .) In the the second
5
region 4 < V < 103 this behavior is reversed. In the third
∼ ∼

region V > 103 the dependence of the effective exponent
αeff on the impurity strength V saturates.
Since discrete lattice models are not Galilean invariant,
1
we decreased the electron filling down to n = 16 such
that the discrete Hamiltonian mimics a continuous electron gas with a quadratic dispersion law. From our data
we see that even at low fillings, i.e. for quasi-continuous
models, we need large systems in order to reach the expected power law exponents.

Acknowledgments

We thank S. Andergassen, V. Meden, and
U. Schollw¨ck for stimulating discussions. This work is
o
partially supported by the Slovak Agency for Science
and Research grant APVV-51-003505, APVV-VVCE0058-07, QUTE, and VEGA grant No. 1/0633/09 (A.G.
and R.K.). A.G. also acknowledges support from the
Alexander von Humboldt foundation. A.G. and R. K.
thank PTB for hospitality and support.

APPENDIX A: SOLUTION OF THE FRG
EQUATIONS FOR LARGE k

Here we show how the initial conditions given in
Eqs. (16) are determined. We start from the equations
(14) for the coefficients ak for large k
a′
k,jj =

2U sin(k 0+ )
,
π
k

a′
k,j(j±1) = 0,

(A1)

noting that Uk = U for large k (0+ denotes a positive
infinitesimal increment). This equation is easily solved
using the initial conditions at k = ∞ given in Eq. (15)
with the result
ak,jj = −U

1−

ak,j(j±1) = 0.

2
Si(k 0+ ) δjj − V δj1 + µδjj ,
π
(A2)

Here, Si(z) is the sine integral function. The impurity
with strength V is placed at site j = 1. This result can
now be used in order to solve the equation for ak,0 for
large k,
Lµ − V sin(k 0+ ) LU sin(k 0+ )
−
π
k
π
k

2
Si(k 0+ )
π
(A3)
where the expansion ln(1 + x) ≈ x for x < 1 has been
used. Integration of Eq. (A3) yields

a′ =
k,0

ak,0 =

L
V
+
2
2

1−

U
−µ .
2

where the limit 0+ → 0 has been performed.

(A4)

10
APPENDIX B: NUMERICAL SOLUTION OF THE
FLOW EQUATIONS FOR THE SELF ENERGIES
AND THE GROUND STATE ENERGY

The solution of the flow equations for the self-energies
and ground state energy (14) requires the calculation of
the inverse and determinant of (potentially) large cyclic
tridiagonal matrices at each step of the integration of
a set of ordinary differential equations. Of course, this
could be achieved with standard library routines. However, due to the special form of the Eqs. (14) we only need
to determine the cyclic tridiagonal part of these matrices. Therefore, in the following we will develop methods adapted to our special needs. As a consequence we
achieve a significant speed-up of the numerical calculation as well as considerable reduction of the computer
memory requirements. Such improvements enable us to
treat large system sizes.
1.

Inversion of cyclic tridiagonal matrices

For matrices M which can be represented as the sum
of a tridiagonal matrix T and a direct product of two
vectors u and v,
M = T + u ⊗ v,

(B1)

the inverse can be obtained using the Sherman-Morrison
formula [58]
M−1 = T−1 −

z⊗w
1+v·z

For the inversion of the tridiagonal part T, we use
a method inspired by Andergassen et al. [50], which is
based on Ref. [59]. However, in our case the tridiagonal
matrix is not complex symmetric, therefore, we have to
use an LU decomposition of T instead of the LDU decomposition employed in Ref. [50]. Details are described
in the following subsection.

2.

Implementation

According to Eq. (B2) the inversion of a general cyclic
tridiagonal matrix requires two essential steps: (i) the inversion of a tridiagonal matrix and (ii) the determination
of the vectors z and w.
(i) A tridiagonal matrix can be represented as a prod˜˜
˜
uct of two matrices T = LU = UL, where U, U are
˜ are lower diagonal matrices.
upper diagonal and L, L
From the LU decomposition
T = LU =

d1
 c2 d2

..
..

.
.


c

n−1


dn−1
cn dn







(B7)


1 u1
1 u2
.. ..
. .
1 un−1
1







with
d1 = α1 ,

(B2)

di = αi − ui−1 ci ,

ui = bi /di

(B8)

one easily finds for 1 < i < j < n the relations

with
z = T−1 u,

w = (T−1 )T v.

A general cyclic tridiagonal matrix

a1 b 1
0
···
c1
 c2 a2 b 2
0


.
..
..
..
.
.
.
.
.
M= 0

 .
.
 .
cn−1 an−1 bn−1
bn · · · 0
cn
an

(B3)

(T−1 )i,j = −ui (T−1 )i+1,j .

(B9)

˜˜
Similarly, from the UL decomposition of T,









(B4)

can be written in the form (B1) with the tridiagonal part
T given by


α1 b1
0
···
0
 c2 α2 b2
0 



. 
..
..
..
. 
.
.
.
. 
T= 0
(B5)

 .

 .
.
cn−1 αn−1 bn−1 
0 ··· 0
cn
αn

with α1 = 2a1 , αn = an + c1 bn /a1 and αi = ai for
1 < i < n. The vectors u and v are defined by
u = (−a1 , 0, · · · , 0, bn ) ,

(T−1 )n,n = 1/dn ,

v = (1, 0, · · · , 0, −c1 /a1 ) .
(B6)

˜˜
T = UL =

1 u1
˜
1 u2
˜


.. ..

. .


1 un−1
˜
1

(B10)

 d
˜1
˜

 c2 d2




..
..


.
.




˜
cn−1 dn−1
˜n
cn d

with
˜
dn = αn

˜
di = αi − ui ci+1
˜

˜
ui = bi /di+1 ,
˜

(B11)

one obtains for 1 < i < j < n the relations
˜
(T−1 )1,1 = 1/d1 ,

(T−1 )i,j = −˜j−1 (T−1 )i,j−1 . (B12)
u

Combining Eqs. (B9) and (B12) yields a recursion relation for the diagonal elements of T−1
(T−1 )i+1,i+1 =

ui −1
˜
(T )i,i
ui

(B13)

11
with the initial condition given in Eq. (B12). From the
diagonal elements one then generates the upper diagonal
as well as the upper corner element using the relations
˜˜
given above. Alternatively, the UL decomposition yields
the relation
ci
(T−1 )i,j = − (T−1 )i−1,j for 1 < j < i < n, (B14)
˜
di
which is used to generate the lower diagonal and corner
element, respectively.
(ii) The vector z is obtained from Eq. (B3), Tz =
LUz = u, as follows: first we solve recursively Lz′ = u
for z′
′
z1 = u1 /d1 ,

ci

bi +
i=1

n

n

n

(−1)n+1

+ Tr

i=1

i=1

ai −bi−1 ci
1
0

with b0 = bn . For large systems, the products in this
formula can easily underflow or overflow. This must be
carefully controlled appropriately.

′
′
zi = (ui −ci zi−1 )di for 1 < i ≤ n, (B15)

and then obtain z recursively from Uz = z′
′
zn = zn ,

′
zi = zi − ui zi+1

for 1 ≤ i < n.

(B16)

In a similar way one determines w from TT w = v.
3.

Determinant of a cyclic tridiagonal matrix

The determinant of a cyclic tridiagonal matrix is calculated using a formula given in Ref. [60]


a1 b 1
0
···
c1
 c2 a2 b 2
0 



. 
..
..
..
. =
 0
.
.
.
. 
det 
(B17)

 .
 .
.
cn−1 an−1 bn−1 
bn · · · 0
cn
an

REFERENCES

[1] B.S. Deaver Jr. and W. M. Fairbank. Phys. Rev. Lett.,
7:43, 1961.
[2] N. Byers and C. N. Yang. Phys. Rev. Lett., 7:46, 1961.
[3] M. B¨ttiker, Y. Imry, and R. Landauer. Phys. Lett.,
u
96A:365, 1983.
[4] L.P. L´vy, G. Dolan, J. Dunsmuir, and H. Bouchiat.
e
Phys. Rev. Lett., 64:2074, 1990.
[5] V. Chandrasekhar, R. A. Webb, M. J. Brady, M. B.
Ketchen, W. J. Gallagher W J, and A. Kleinsasser. Phys.
Rev. Lett., 67:3578, 1991.
[6] D. Mailly, C. Chapelier, and A. Benoit. Phys. Rev. Lett.,
70:2020, 1993.
[7] B. Reulet, M. Ramin, H. Bouchiat, and D. Mailly D.
Phys. Rev. Lett., 75:124, 1995.
[8] E. M. Q. Jariwala, P. Mohanty, M. B. Ketchen, and R. A.
Webb. Phys. Rev. Lett., 86:1594, 2001.
[9] W. Rabaout, L. Saminadayar, D. Mailly, K. Hasselbach,
A. Benoˆ and B. Etienne. Phys. Rev. Lett., 86:3124,
it,
2001.
[10] N. A. J. M. Kleemans et al. Phys. Rev. Lett., 99:146808,
2007.
[11] H. Bouchiat. Physics, 1:7, 2008.
[12] R. Landauer and M. B¨ttiker. Phys. Rev. Lett., 54:2049,
u
1985.
[13] H. F. Cheung, Y. Gefen, E. K. Riedel, and W. H. Shih.

Phys. Rev. B, 37:6050, 1988.
[14] M. Abraham and R. Berkovits. Phys. Rev. Lett., 70:1509,
1993.
[15] G. Bouzerar, D. Poiblanc, and G. Montambaux. Phys.
Rev. B, 49:8258, 1994.
[16] A. O. Gogoglin and Prokof’ev. Phys. Rev. B, 50:4921,
1994.
[17] H. Kato and D. Yoshioka. Phys. Rev. B, 50:4943, 1994.
[18] A. Cohen, K. Richter, and R. Berkovits. Phys. Rev. B,
57:6223, 1998.
[19] M. Henkel and D. Karevski. Eur. Phys. J. B, 5:787, 1998.
[20] S. Jaimungal, M. H. S. Amin, and G. Rose. Int. J. Mod.
Phys. B, 13:3171, 1999.
[21] V. Meden and U. Schollw¨ck. Phys. Rev. B, 67:035106,
o
2003.
[22] T. Enss, V. Meden, S. Andergassen, X. Barnab´e
Th´riault, W. Metzner, and K. Sch¨nhammer. Phys.
e
o
Rev. B, 71:155401, 2005.
[23] S. Friederich and V. Meden. Phys. Rev. B, 77:195122,
2008.
[24] H. Bary-Soroker, O. Entin-Wohlmann, and Y. Imry.
Phys. Rev. Lett., 101:057001, 2008.
[25] J. Heinrichs. J. Phys.: Condens. Matter, 20:345232,
2008.
[26] B. L. Altshuler, Y. Gefen, and Y. Imry Y. Phys. Rev.

12
Lett., 66:88, 1991.
[27] V. Ambegaokar and U. Eckern. Europhys. Lett., 13:733,
1990.
[28] E. Gambetti-C´sare, D. Weinmann, R. A. Jalabert, and
e
Ph. Brune. Eurphys. Lett., 60:120, 2005.
[29] E. Gambetti. Phys. Rev. B, 72:165338, 2005.
[30] R. Krˇm´r, A. Gendiar, M. Moˇko, R. N´meth,
c a
s
e
P. Vagner, and L. Mitas. Physica E, 40:1507, 2008.
[31] R. N´meth, M. Moˇko, R. Krˇm´r, A. Gendiar, M. Ine
s
c a
dlekofer, and L. Mitas. arXiv:0902.2225.
[32] A. Luther and I. Peschel. Phys. Rev. B, 9:2911, 1974.
[33] W. Apel and T. M. Rice. Phys. Rev. B, 26:7063, 1982.
[34] C. L. Kane and M. P. A. Fisher. Phys. Rev. Lett., 68:1220,
1992.
[35] C. L. Kane and M. P. A. Fisher. Phys. Rev. B, 46:15233,
1992.
[36] F. D. M. Haldane. Phys. Rev. Lett., 45:1358, 1980.
[37] K. A. Matveev, D. Yue, and L. I. Glazmann. Phys. Rev.
Lett., 71:3351, 1993.
[38] D. Yue, L. I. Glazmann, and K. A. Matveev. Phys. Rev.
B, 49:1966, 1994.
[39] J. Polchinski. Nucl. Phys. B, 231:269, 1984.
[40] C. Wetterich. Phys. Lett. B, 301:90, 1993.
[41] T. R. Morris. Int. J. Mod. Phys. A, 9:2411, 1994.
[42] M. Salmhofer. Commun. Math. Phys., 194:249, 1998.
[43] M. Salmhofer and C. Honerkamp. Prog. Theor. Phys.,
105:1, 2001.
[44] S. R. White. Phys. Rev. Lett., 69:2863, 1992.
[45] S. R. White. Phys. Rev. B, 48:10345, 1993.
[46] I. Peschel, X. Wang, M. Kaulke, and K. Hallberg (Eds.).
Lecture Notes in Physics 528 Density-Matrix Renormal-

[47]
[48]
[49]
[50]

[51]
[52]
[53]

[54]
[55]
[56]
[57]
[58]

[59]
[60]

ization, A New Numerical Method in Physics. Springer,
Berlin, 1999.
U. Schollw¨ck. Rev. Mod. Phys., 77:259, 2005.
o
R. E. Peierls. Z. Phys. B: Condens. Matter, 80:763, 1933.
V. Meden, P. Schmitteckert, and N. Shannon. Phys. Rev.
B, 57:8878, 1998.
S. Andergassen, T. Enss, V. Meden, W. Metzner,
U.Schollw¨ck, and K. Sch¨nhammer. Phys. Rev. B,
o
o
70:075102, 2004.
V. Meden, W. Metzner, U. Schollw¨ck, and
o
K. Sch¨nhammer. Phys. Rev. B, 65:045318, 2002.
o
V. Meden, W. Metzner, U. Schollw¨ck, and
o
K. Sch¨nhammer.
o
J. Low Temp. Phys., 126:1147,
2002.
S. Andergassen, T. Enss, V. Meden, W. Metzner,
U. Schollw¨ck, and K. Sch¨nhammer. Phys. Rev. B,
o
o
70:075102, 2004.
M. Weyrauch and D. Sibold. Phys. Rev. B, 77:125309,
2008.
J. Berges, N. Tetradis, and C. Wetterich. Phys. Rep.,
363:223–386, 2002.
M. Schick. Phys. Rev., 166:144, 1967.
S. Qin, M. Fabrizio, L. Yu, M. Oshikawa, and I. Affleck.
Phys. Rev. B, 56:9766, 1997.
S.A. Teukolsky, W. T. Vetterling, and B.P.Flannery. Numerical Recipes in Fortran 77. Cambrige University
Press, Cambridge, MA, 1986.
G. Meurant. SIAM J. Matrix Anal. Appl., 13:707, 1992.
L. G. Molinari. arXiv:0712.0681v3.