Table Of ContentSpectral functions and time evolution from the Chebyshev recursion
F. Alexander Wolf,1 Jorge A. Justiniano,1 Ian P. McCulloch,2 and Ulrich Schollw¨ock1
1Department of Physics, Arnold Sommerfeld Center for Theoretical Physics,
LMU Munich, Theresienstrasse 37, 80333 Mu¨nchen, Germany
2Centre for Engineered Quantum Systems, School of Physical Sciences,
The University of Queensland, Brisbane, Queensland 4072, Australia
(Dated: April 28, 2015)
WelinklinearpredictionofChebyshevandFourierexpansionstoanalyticcontinuation. Wepush
the resolution in the Chebyshev-based computation of T = 0 many-body spectral functions to a
much higher precision by deriving a modified Chebyshev series expansion that allows to reduce
the expansion order by a factor ∼ 1. We show that in a certain limit the Chebyshev technique
5 6
becomes equivalent to computing spectral functions via time evolution and subsequent Fourier
1
transform. This introduces a novel recursive time evolution algorithm that instead of the group
0
operator e−iHt only involves the action of the generator H. For quantum impurity problems, we
2
introduceanadapteddiscretizationschemeforthebathspectralfunction. Wediscusstherelevance
r oftheseresultsformatrixproductstate(MPS)basedDMRG-typealgorithms,andtheirusewithin
p
dynamical mean-field theory (DMFT). We present strong evidence that the Chebyshev recursion
A
extractslessspectralinformationfromHthantimeevolutionalgorithmswhenfixingagivenamount
7 of created entanglement.
2
] I. INTRODUCTION terms, and by that reliably extracts much information
el about an underlying function from its local knowledge
- ExpandingthespectraldensityA(ω)ofanoperatorH {tn,gn}. In order for this to be meaningful, the underly-
r
t in the monomes ωn via the moments ingfunction,e.g.aGreen’sfunction,mustbecompatible
s
with (1).
. (cid:90)
at µmon = dωA(ω)ωn, In particular, we note that (1) can serve as an ansatz
n
m for analytic continuation of a zero-temperature Green’s
- is a tool that originates in the early days of quan- function
d tum mechanics.1 Computing these moments iteratively
n though is numerically unstable2,3 and one replaced ex- G(t)=−i(cid:104)ψ0|e−i(H−E0)t|ψ0(cid:105), (2)
o
c pansions in ωn by expansions in such polynomials pn(ω) where |ψ0(cid:105) is a single-particle excitation of the ground-
[ ofdegreenthatcanbestablycomputed.4,5 Aprominent state |E (cid:105) of H, for example the creation of a fermion
2 exampleforpn(ω)areChebyshevpolynomials,whoseas- |ψ0(cid:105) =0c†|E0(cid:105). Note that in the case of fermions,
sociated three-term recursion is stable as it does not ad- (2) describes only the t > 0 contribution (then usu-
v
6 mit a so-called minimal2 solution. ally more precisely denoted G>(t)) of the full fermionic
1 After the development of stable recursions the next Green’sfunction. G(t)isanalyticeverywhereinthecom-
2 stepinthemid1990swastheintroductionofkernelsthat plex plane except for t → i∞ and thereby allows for
7 damptheerroneousGibbs oscillations oftruncatedpoly- an analytic continuation of G(t) from a local descrip-
0 nomial expansions of discontinuous functions,6–8 which tion {t ,G(t )} to the domain [t ,∞). This analytic
. n n 0
1 lead to the kernel polynomial approximation. It deals continuation is highly different from the ill-conditioned
0 with redefined series expansions that represent the con- problem of continuing the frequency-space represented
5 volutionoftheexpandedfunctionwithabroadening ker- Green’s function from a domain in the complex plane
1 nel, like a Gaussian or Lorentzian. This technique has (e.g. the imaginary-frequency axis or a parallel of the
:
v been reviewed in Ref. 1 and more recently in Ref. 9 from real-frequencyaxis)tothereal-frequencyaxis,wherethe
i a numerical linear algebra perspective. frequency-space Green’s function has poles.
X
Inthispaper,wedroptheideaofsuchbroadening ker- In the context of Green’s functions, linear prediction
ar nelsinfrequencyspaceortheequivalentdamping orwin- has for the first time been used to extrapolate the time
dowing kernels in the associated Fourier or Chebyshev evolution of the spin structure factor in the one dimen-
expansions. Instead, we employ the fundamentally dif- sional Heisenberg model.11,12 While for the spin-1 model
ferent technique of linear prediction.10 Linear prediction it was clear that the ansatz (1) is justified as the time
is a linear recursive reformulation (Appendix C) of the evolution is dominated by a small number of magnons
non-linear problem to fit the surrogate function whose excitation energies correspond directly to the fre-
g(t)=(cid:88)α eiωit, α ,ω ∈C, t∈R, (1) quencies ωi in (1),11 this was not the case for the spin-12
i i i model.12 In the latter, spinons dominate which lead to
i an(infinitely)highnumberofpolesonthereal-frequency
to given numerical data {t ,g }. Due to linearity, linear axis, and the direct correspondence of pole energies and
n n
prediction is able to treat superpositions of hundreds of frequencies ω in (1) is lost. Still the ansatz works12 in
i
2
anapproximatesensebyextractingeffectivefrequencies. Analogously, any integrable function f(ω)| ,
ω∈[−a,a]
For the computation of spectral functions, the use of where a∈R, can be expanded in a Fourier series 2 2
linear prediction for the time evolution of Green’s func-
∞
tionsprovidesahighlyattractivealternativeapproachto f(ω)= 1 (cid:88) eiωtnf(t ), (7a)
the usual damping or windowing in real-time or broad- 2aπ n
n=−∞
ening in frequency space: an approach that enhances
(cid:90) a/2 n
resolution in frequency space. Up to now, it is not en- f(t )= dωe−iωtnf(ω), t = , (7b)
tirelyclearinwhichcasesthisiscontrolled. Ontheother n −a/2 n a
hand, the approach of damping the truncated series ex-
which represents a Fourier transform for a→∞.
pansioncannotbeconsideredcontrolled,too: Althougha
broadenedfunctionf (ω),whichisforGaussianbroaden-
η
ing given by fη(ω)= √1 (cid:82) dω(cid:48)e−(ω(cid:48)−ω)2/2η2f(ω(cid:48)), con- B. Expansion of a spectral function
η 2π
vergesuniformlytotheunderlyingoriginalfunctionf(ω)
forη →0,extractionofinformation(deconvolution)from Consider now the expansion of the spectral function
fη(ω)aboutf(ω)isuncontrolledasitcorrespondstothe A(ω) of a Hamilton operator H with respect to a refer-
problem of analytic continuation from a domain in the ence energy E and a state |ψ (cid:105) as in (2)
ref 0
complex plane to the real axis.
A(ω)=(cid:104)ψ |δ(ω−(H −E ))|ψ (cid:105). (8)
Recently, Ref. 13 suggested to extrapolate the Cheby- 0 ref 0
shevexpansionofaspectralfunctionusinglinearpredic-
ThespectralfunctionisrelatedtotheGreen’sfunctionof
tion,albeitonlyjustifiedbytheempiricalsuccess. Inthe (2) via its Fourier transform: A(ω)=−1ImG(ω+i0+).
remainder of this introduction, we place these results in π
The coefficients of the Fourier expansion can be com-
the context of the preceding discussion, and by that put (cid:80)
puted by inserting an identity of eigenstates |E (cid:105)(cid:104)E |
this approach on more firm grounds. i i i
in the integral over the delta function δ(ω−(H−E ))
ref
(cid:90) a/2
f(t )= dωe−iωtnA(ω)=(cid:104)ψ |ψ(t )(cid:105), (9a)
n 0 n
A. Chebyshev and Fourier transformation basics
−a/2
n
The Chebyshev polynomials of the first kind |ψ(tn)(cid:105)=e−i(H−Eref)tn|ψ0(cid:105), tn = a. (9b)
In order for (9) to hold true, a must be chosen large
Tn(x)=cos(narccos(x)) (3) enough for that the support of A(ω) is contained in
[−a,a]. A sufficient condition for that is spec(H −
2 2
can be generated by the recursion E ) ⊂ [−a,a], which is possible as we consider opera-
ref 2 2
tors H with bounded spectra. Eq. (9b) makes it obvious
Tn(x)=2xTn−1(x)−Tn−2(x), T1 =x, T0 =1, (4) that a has the meaning of an inverse time step.
To compute the coefficients for the Chebyshev expan-
whichisnumericallystableif|x|≤1. Chebyshevpolyno- sion, we need to consider a spectral function whose sup-
mialsareorthonormalwithrespecttotheweightedinner portiscontainedin[−1,1]. Forthis,introducearescaled
product and shifted version of H with appropriately chosen con-
stants a and b
(cid:90) 1
dxwn(x)Tm(x)Tn(x)=δnm, (5a) H = H −Eref +b, x= ω +b, (10)
−1 a,b a a
2−δ
w (x)= √ n0 . (5b) where a can again be considered an “inverse time step”
n
π 1−x2 and H is dimensionless. Note that in Ref. 14, the defini-
tion of b differed from the one here by a factor a. Then
Anyintegrablefunctionf(x)| canbeexpandedin
x∈[−1,1]
T (x) A (x)=(cid:104)ψ |δ(x−H )|ψ (cid:105) (11)
n a,b 0 a,b 0
∞ yields the original spectral function via A(ω)= 1A(cid:0)ω +
f(x)= (cid:88)wn(x)µnTn(x), (6a) b(cid:1), where we omitted to specify the indices a,ab, asa in
n=0 most of the rest of this paper. The Chebyshev moments
(cid:90) 1 for A(x) can be computed analogously to the Fourier
µ = dxT (x)f(x), (6b)
n n coefficients (9)
−1
(cid:90) 1
where the definition of the so-called Chebyshev moments µ = dxA(x)T (x)=(cid:104)ψ |ψ (cid:105), (12a)
n n 0 n
µ viathenon-weighted innerproduct(6b)followswhen −1
n
(cid:82)1 |ψ (cid:105)=T (H)|ψ (cid:105). (12b)
applying dxT (x)... to both sides of (6a). n n 0
−1 m
3
Inserting the recursive definition (4) of T (H) in the II. SPECTRAL FUNCTIONS IN THE
n
definition (12b) of |ψ (cid:105), one obtains a practical calcula- THERMODYNAMIC LIMIT
n
tionschemeforthepowerseriesexpansionofT (H),and
n
by that for the Chebyshev states |ψn(cid:105) in (12): A spectral function for a system of finite size L has
a finite-size peak structure due to an agglomeration of
|ψ (cid:105)=2H|ψ (cid:105)−|ψ (cid:105), |ψ (cid:105)=H|ψ (cid:105). (13) eigenvalues that is not present in the thermodynamic
n n−1 n−2 1 0
limit. In a weakly interacting system, this agglomer-
ation happens around the positions of the eigenvalues
of the corresponding noninteracting (single-particle) sys-
tem. This argument gives us the best, though still very
C. Analytic continuation
rough, estimate W /L for the spacing of finite-size
single
peaks,whereW isthesingle-particlebandwidth. At
single
Acomparisonof (9b)and(12b)clarifywhylinearpre- amuchsmallerspacingthanthat,spectralfunctionshave
diction is an equally justified approach for a Chebyshev an underlying delta peak structure, as is obvious from
as for a Fourier expansion. definition (8), which can be rewritten as
Rewriting the “evolution operators” that appear in (cid:88)
A(ω)= W δ(ω−(E −E )), (15)
(9b) and (12b) as i i ref
i
exp(−inHa,b=0) and cos(narccos(Ha,b)) (14) with weights Wi = |(cid:104)ψ0|Ei(cid:105)|2. The delta peak structure
merges to a (section-wise) smooth function only in the
thermodynamic limit.
makesitclearthatwedealwithanalytic functionsofnif
Expandingthespectralfunctionofafinite-sizesystem
we consider n as a continuous complex variable. By (9a)
in orthogonal polynomials is a very efficient way to not
and (12a), this makes the Fourier and the Chebyshev
resolveeitherfinite-sizepeaksorthedeltapeakstructure,
coefficients f(t ) and µ analytic functions of n, too.
n n buttoextractonlythesmoothfunctionofthethermody-
In addition, Re exp(−inH ) = cos(nH ) shows that
a,b a,b namiclimit,ase.g.discussedinRef.14. Itisthisfunction
the real part of the functional dependence of the time
ofthethermodynamiclimitthatweareinterestedin,and
evolution operator on n is the same as for the “Cheby-
for which we start our discussion.
shev evolution operator”, when neglecting a redefinition
ofoscillationsfrequencies. Asthisredefinitioncanbeac-
counted for by the fitting procedure, the particular form
A. Discontinuity of spectral functions
of the surrogate function g(t) in (1) is equally suited to
analytically continue both types of expansions. A fun-
The state |ψ (cid:105) and the energy E in (15) are gener-
damental theorem from complex analysis then tells us 0 ref
allyassociatedtothegroundstateofacertainsymmetry
that if linear prediction provides us with a function g(t)
sector N of H, which for fermions is typically a particle
that locally agrees with f(t ) or µ , we know that this
n n number. The reference energy for |ψ (cid:105) = c†|E (cid:105) then
function globally agrees with f(t ) or µ . Of course, in 0 0
n n is the Fermi energy, which is the ground state energy
practicetheseargumentsaretobetakenwithcare,aswe
E = EN−1 of the contiguous symmetry sector of |ψ (cid:105)
will never numerically find a function g(t ) that agrees ref 0 0
n (or E = EN+1 for a hole excitation). The weights
exactly with the local data {g ,t }. ref 0
n n W =|(cid:104)ψ |E (cid:105)|2 inthespectralfunction(15)canbenon-
i 0 i
zeroonlyforeigenstates|E (cid:105)andeigenvaluesE fromthe
i i
sector N. The particular meaning of E as a ground
ref
state energy then implies that even if the global spectral
D. Outline of the paper
function
(cid:88)
A (ω)= δ(ω−(E −E )) (16)
We will first study the convergence properties of the global i ref
Chebyshev expansion of discontinuous (spectral) func- i
tionsinthethermodynamiclimit. Thisallowstoderivea issmooth,theweightsW generallyintroduceadisconti-
i
newschemeforaChebyshevseriesdefintionthatleadsto nuity at ω =0 (we use the term global here, as |ψ (cid:105) usu-
0
exponential convergence and allows to reduce expansion ally is a local excitation associated with a certain quan-
orders in practical calculations by a factor ∼ 1 (Sec. II). tum number).
6
We then apply these results to the computation of spec-
tralfunctionsforfinitesystems(Sec.III),anddiscussthe
relevance for matrix product state (MPS) based compu- B. Convergence of Chebyshev series expansions
tations(Sec.IV).Afterthat,wedescribetheapproximate
equivalenceoftheChebyshevrecursiontotimeevolution The convergence of the Chebyshev moments µ → 0
n
and show how this leads to a novel time evolution algo- of a function f(x) in the limit n → ∞ can be charac-
rithm(Sec.V).Finally, weconcludethepaper(Sec.VI). terized by the degree of differentiability of f(x), similar
4
to a Fourier expansion.15 Let k denote the highest inte-
ger for which the kth derivative of f(x) is integrable: If
f(x) is smooth (k = ∞), the envelope of µn converges 1 )1
exponentially to zero with respect to n; if f(x) is a step (
) >
function (k = 1), the envelope converges algebraically ( A
> v
with 1; and if f(x) is a delta function, the envelope re- A
n v
mains constant. In general, the order of convergence is 0
0 2 4
at least 1 . Although in Ref. 15, this is stated for mo- /v
nk
ments computed with the weighted inner product (5b),
it also holds for moments computed using (6b) (see Ap- 0
E -E 0 E -E
min ref max ref
pendix A). In practice, we are not interested in the limit
/v
n→∞, but rather in intermediate values of n: but also
here, the degree of differentiability of A(ω) helps us to 101
10-1 >
learn something about the convergence of µ . s n
n nt >
A>C(oωn)siadsershaowtnypiincatlhedistcoopntpinanueolusofspFeigc.tr1a.l fIutnsctcioorn- ome 100 nn
responding Chebyshev moments µ>n are computed by v m10-1 10-20 10 20
numerically integrating (6b) and shown in the bottom e
h
panel of Fig. 1 as blue circles. The blue line in the in- s
y
set shows the envelope of µ>n, which evidently decreases heb10-2
algebraically to zero. C
NownotethatcontinuityofA>(ω)atω =0caneasily 10-3
0 5 10 15 20 25
be restored by defining
n/a (1/v)
A(cid:101)>(ω)=A>(ω)−A>(0). (17)
FIG. 1. (Color online) Top: Typical example of a discon-
The green crosses (lines) in the bottom panel of Fig. 1 tinuous spectral function due to the restriction to a given
showthattheChebyshevmomentsµ(cid:101)>n ofA(cid:101)>(ω)converge symmetry sector. In this case, this is the particle contribu-
exponentially for the values of n considered in the plot, tion of the spectral function of the single-impurity Anderson
i.e.,qualitatively differentlythanµ>. Thisisobservedal- model(SIAM)withsemi-ellipticbathdensityofstatesofhalf-
n bandwidth2v andinteractionU/v=4,16 takenfromRef.14.
thoughA(cid:101)>(ω)isnotsmooth,butonlyoncedifferentiable The spectral function is given by (8) using |ψ (cid:105) = c†|E (cid:105),
0 0
(kink in first derivative at ω =0).
where|E (cid:105)denotesthehalf-filledgroundstateandE isthe
0 ref
WhiletheconstructionofA(cid:101)>(ω)iscompletelygeneral, Fermi energy E0. Here, only the shape of the scalar function
for the particular case of a fermionic spectral function, is of importance, therefore we postpone the model definition
another way of constructing a continuous function from to (30). The same spectral function is obtained for the local
A>(ω) has been favored: In the appendix of Ref. 17, it density of states of the first site for spinless fermions hop-
ping on a semi-infinite chain with tunneling v and and an
was mentioned that the Chebyshev expansion of the full
interaction of U/v = 4 that acts only on the first site. Bot-
spectral function
tom: ComparisonofconvergenceoftheChebyshevmoments
A(ω)=A>(ω)+A<(−ω), (18) of A>(ω) with its redefinitions A(cid:101)>(ω) and A(ω), giving rise
to moments µ>, µ> and µ , respectively. The full spectral
≷ ≷ ≷ n (cid:101)n n
A (ω)=(cid:104)ψ0 |δ(ω−(H −E0))|ψ0 (cid:105), function(18)forthisexampleisA(ω)=A>(ω)+A>(−ω)as
here, A<(ω) = A>(ω) due to particle-hole symmetry. All of
obtainedbysummingoverparticle(>)andhole(<)con-
this is for the setup b = 0 using a rescaling of a = 100v in
tributions, is much better suited for a Chebyshev expan-
(10).
sion than A≷(ω), as it lacks the discontinuity. In Ref. 13
it was then pointed out that the full A(ω) is smooth and
therefore,Chebyshevmomentsshoulddecreaseexponen-
tially, which would allow to use linear prediction. In C. Comparison of setups b=0 and b(cid:39)−1
general, it is not true that A(ω) is smooth, due to the
possibilityofvanHove singularities,asappeare.g.forthe In Ref. 14, we pointed out that the choice b = 0 in
U = 0 case of the spectral function of the single impu- (10)iscomputationallymuchlessefficientthanthechoice
rity Anderson model (SIAM) (see Appendix B, Fig. 11). b(cid:39)−1(called“b(cid:39)−a”inRef.14). Whereasconstruct-
Still, A(ω) is likely to be smooth, and for the present ex- ing a Chebyshev expansion of the full spectral function
ample, it is. The bottom panel of Fig. 1 therefore shows A(ω) requires choosing b = 0, this is not the case for
thattheChebyshevmomentsµn forA(ω)decreaseatthe A(cid:101)(ω). For A(cid:101)(ω), we can therefore use the exponential
same exponential rate as the moments µ(cid:101)>n of A(cid:101)>(ω). rate of convergence to quantify the amount of spectral
The statements about the qualitatively different con- information that the Chebyshev recursion extracts from
vergence behaviors of A>(ω), A(cid:101)>(ω) and A(ω) are con- H in the setups b = 0 and b (cid:39) −1, and by that under-
firmed for further typical examples in Appendix B. stand the observations of Ref. 14 quantitatively.
5
150 150 inthesetupb=−0.995canbeunderstoodbylookingat
n=0
the natural stretching of the frequency scale of Cheby-
n=50
100 100 n=100 shev polynomials close to the boundaries of [−1,1]. Ex-
n=150
n=200 pressing the integral (19) by substituting x=cosθ
x) 50 x) 50 n=250
( (
n n n=300
T T
x) 0 x) 0 (cid:90) 0
>(a,bA >(a,bA µ>n =− dθA>a,b(cosθ)cos(nθ)sinθ, (20)
a -50 a -50 π
-100 -100 one arrives at a convolution with the regularly oscillat-
bb==00 bb==--00..999955 ing cos(nθ). Consider now the interval of width 0.05 on
-150 -150 [−1,1], which corresponds to the (single-particle) sup-
0 0.02 0.04 0.06 -1 -0.98 -0.96 -0.94
x x port of A> (x) in the example of Fig. 2. By comput-
101 ing the inat,ebgral widths under the map x = cosθ, one
> b=-0.995
ments 100 nn>> bb==0-0.995 l[−ea0r.n9s95th,−at0p.9la5c],inagstrheesuslutpspfoorrtbin=th−e0“.9b9o5u,nidnacrryearseegsiotnh”e
mo10-1 n resolution by a factor ∼ 6.4 compared to placing it in
the “center region” [0,0.05], as results for b = 0. These
v
he10-2 effectsarewellknownboundaryeffectsoftheChebyshev
ys polynomialsthatareexploitedalsointhesolutionofdif-
eb10-3 ferential equations.15
h
C The bottom panel of Fig. 2 shows the Chebyshev mo-
10-4 mentsobtainedintheb=−0.995setup,forA>(ω)(blue
0 2 4 6 8 10
n/a (1/v) circles) and for A(cid:101)>(ω) (green pluses). As mentioned be-
fore, for this setup, no Chebyshev expansion of the full
FIG. 2. (Color online) Top: Integrand for computation of spectral function A(ω) is possible. Instead, we compare
moments for the spectral function shown in Fig. 1. in the the b = −0.995 results to the b = 0 results, depicted as
two setups b = 0 (left) and b = −0.995 (right). Bottom: red crosses. It is evident that the Chebyshev expansion
Comparison of convergence of moments computed with the
in the b = −0.995 setup converges much faster than the
integrands shown in the top panels. In all of that, a=100v.
one in the b = 0 setup. After n = 1000 iterations, the
magnitude differs by more than 100.
This difference directly appears in the error of the
The key observation to make is that the integral
Chebyshevseries,asstatedbythefollowinggeneralrule:
(cid:90) 1 theorderoftheerrorεofaChebyshev(orFourier)series
µ>n = dxA>a,b(x)Tn(x) (19) representation of a function that is truncated at n = N
−1 can be estimated by (see Ref. 15, Chap. 2.12)
extracts a highly different amount of information about
thestructureofA>(ω)dependingonhowaandbin(10) ε=O(µ ) if µ converges exponentially, (21a)
N n
are chosen when generating A> (x). ε=O(Nµ ) if µ converges algebraically. (21b)
a,b N n
Throughout the whole paper, we keep a = 100v fixed
toguaranteethenumericalstabilityoftheChebyshevre-
cursionforthetypicalsystemsizesofaroundL≥80that
arelargeenoughtodisplay“thermodynamiclimitbehav- D. Linear prediction for the Chebyshev expansion
ior”. Ifwechoseasmaller,wecouldonlystablycompute
“small” systems or we would have to resort to the tech-
The main motivation for studying the convergence of
nique of energy truncation, which is strongly prone to
different Chebyshev expansions in the previous sections
errors.17 Furthermore, in the MPS context, it is impor-
lies in the possibility to extrapolate exponentially de-
tant to compare only computations in which a is kept
creasingsequenceswithlinearprediction. Asdiscussedin
constant: constant a means constant effective hopping
theintroduction,thelatterallowsanextremelyhighgain
energies v inH,andbythataconstantamountofentan-
a in resolution, if its application is justified. For details on
glementproductioninasingleiterationstepof (13). The
linear prediction, see Appendix C.
parameterb,bycontrast,canbechosenfreelywithoutaf-
In what follows, we compare the known approach of
fecting the numerical stability, and in principle, without
using linear prediction for the Chebyshev expansion of
affectingentanglementproductioninMPScomputations.
A(ω), as suggested in Ref. 13, with the new approach of
The top panels of Fig. 2 show the convolution of
A> (x) with Chebyshev polynomials T (x) of different extrapolating the Chebyshev expansion of A(cid:101)>(ω).
a,b n
degree n for the two setups b=0 and b=−0.995(cid:39)−1. We first compute the Chebyshev moments of the step
Thehighlyincreasedoscillationfrequencythatisevident function that has the discontinuity of A>(ω) at ω = 0,
6
which transforms to x=−b for A>(x), as 100
1.2 input b=0 full
b=0 full b=-0.995
(cid:90) 1 b=-0.995 )|
µsntep = dxTn(x) (22) 1 (ut10-1
−b p
1(cid:16)cos[(n+1)arccosx] )0.8 Ain
= 2 cos[(nn−+11)arccosx](cid:17)(cid:12)1 vA(0.6 ()-onst10-2
− n−1 (cid:12)(cid:12)−b. 0.4 v|Arec10-3
0.2
The Chebyshev moments of A(cid:101)>(x) are then given by
0 10-4
0 1 2 3 4 5 0 1 2 3 4 5
µ> =µ>−A>(0)µstep (23) /v /v
(cid:101)n n n
)
and are accessible by linear prediction, as they decrease -1(v 10-1 bb==0-0 f.u9l9l5
exponentially. )| b=0
The core problem in this new approach is that the (n10-2
value A>(0) of the spectral function is, in general, un- )-Ai10-3
known prior to linear prediction. But it fulfills the fol- (
c
lowingself-consistencyproblem, whichcanbeiteratively Are10-4
solved: ChoosingastartvalueA>(0)forA>(0),wecom- x|
puteµ ,extrapolatethesequence0uptoconvergence,and ma 10-5
(cid:101)n 200 400 600 800 1000120014001600
thenusetheextrapolatedsequencetoreconstructA>(ω), N
which provides us with a new value A>(0). We repeat
1
the procedure until the new and the old version A>(0) FIG. 3. (Color online) Top: Input spectral function and re-
i
and A> (0) agree. This procedure is found to converge constructedspectralfunctionsusinglinearpredictionforN =
i+1
stably and quickly for all examples studied (see also Ap- 200 computed Chebyshev moments. We compare our pro-
pendix B). posal with the original proposal,13 where for the b=0 setup
the expansion of the full spectral function was extrapolated.
Fig.3comparestheapproachofreconstructingthefull
Left: In the first case, we use the Chebyshev expansion of
spectral funciton A(ω) from A>(ω), using linear predic-
tion for the expansion of A(ω) in the b = 0 setup, with A(cid:101)>(ω)intheb=−0.995setup(reddashedlines),andinthe
second, we use the Chebyshev expansion of the full spectral
thenewapproachofusinglinearpredictionofA(cid:101)(ω)inthe function A(ω) in the b = 0 setup (solid green line). Right:
b=−0.995setup. WetakethefunctionofFig.1asinput Errorofthesefunctions|A (ω)−A (ω)|forbothse-
reconst input
function that shall be reconstructed. In the top panels tups. Bottom: The error max |A (ω)−A (ω)|
ω≥0 reconst input
of Fig. 3, we compare both setups for N = 200 com- versusnumberN ofcomputedChebyshevmoments. Here,we
puted moments that are then extrapolated to N (cid:29)1000 alsoshowresultsforusingtheChebyshevexpansionofA(cid:101)>(ω)
until they converge to a value of 10−6. We choose this in the b = 0 setup (light green crossed line). This is differ-
comparatively small number of computed moments, as ent from using the full spectral function A(ω) in the b = 0
in MPS algorithms the number of moments that can be setup. Lower error levels only occur for much higher higher
computed in a controlled way is strongly limited.14 expansion orders than shown in the panel.
The upper left panel of Fig. 3 shows that already for
N = 200, our approach (dashed red line) allows a very
Only at very high expansion orders, that in practice
good reconstruction of the input function. In the up-
canoftennotbereachedbyMPScomputations,theorig-
per right panel we show the error of this reconstruction,
inal approach allows to reach smaller error levels for the
which becomes maximal at the second peak of the input
function and is of order 10−2, i.e., a relative error of a presently studied generic example. More examples are
studied in Appendix B.
few per cent. The situation is very different for the ex-
trapolation scheme of the full A(ω) that uses the b = 0
setup. For N =200 computed moments, large errors are
observed in both top panels of Fig. 3. III. SPECTRAL FUNCTIONS FOR FINITE
In the bottom panel of Fig. 3, we plot the maximal SYSTEMS
error, defined as max |A (ω)−A (ω)|, ver-
ω≥0 reconst input
sus different values of the number of computed moments Let us now study the case of finite systems, where a
N. Anordersofmagnitudereductionoftheerrorisseen discretizedrepresentationofthespectralfunctionisused
upon using our over the previous approach. If one com- for reconstruction. The general previous arguments are
parestheexpansionorderN forwhichanerrorof5·10−3 still valid, but several technical details have to be taken
is reached (N ∼250 in the b(cid:39)−1 setup, and N =1200 into account. In particular, we suggest a new discretiza-
in the b=0 setup), one recovers the factor ∼6 that has tion scheme suited for reconstruction with Chebyshev
been derived in the previous section. expansions. Such a discretization scheme can be used
7
for problems that allow to manipulate the discretization As we already know ω from (25), we only have to
L+1
of the spectral function. This is e.g. the case for im- fixtheintermediateintervalboundaries{ω ,...,ω }.
L/2+1 L
purity models, for which the discretization of the input We define
bath spectral function determines the discretization of
(cid:0) (cid:1)
the spectral function. Still, the following discussion is ω =a cos(θ +l∆θ)−b , l=1,...,L/2
L/2+l L/2
also relevant to e.g. finite lattice models for which the ∆θ = 2(θ −θ ),
discretization is physically constrained. L L+1 L/2
θ =arccosb,
To construct a discrete representation of a continu- L/2
ous function A(ω), we employ the scheme that is used θ =arccos(b+ ωL+1). (28)
L+1 a
todiscretizethehybridization function ofimpuritymod-
els in the numerical renormalization group.18 This pro- Using these definitions, a discrete representation
ceeds as follows. For L given discretization intervals A (ω) of A(ω) (in the sense that A (ω) → A(ω)
discr discr
[ω ,ω ], l = 1,...,L, we compute discrete weights V2 for L→∞) is given by
l l+1 l
and eigenvalue positions (cid:15) by
l
A (ω)=(cid:104)ψ |δ(ω−H)|ψ (cid:105),
discr 0 0
V2 =(cid:90) ωl+1dωA(ω), (24) Hll(cid:48) =(cid:15)lδll(cid:48), l,l(cid:48) =1,...,L
l
1ωl (cid:90) ωl+1 |ψ0(cid:105)l =Vl, (29)
(cid:15) = dωωA(ω).
l V2 where H ∈ RL×L and |ψ (cid:105) ∈ RL, and the parameters (cid:15)
l ωl 0 l
and V are given in (24). This is consistent with defini-
l
The first line associates a weight and the second line tion (8) if we realize that this is a single-particle Hamil-
a representative energy with an interval of energies tonian for a particle that is in either of (cid:15) energy states
l
[ωl,ωl+1]. Fortheenergy,onecoulde.g.take9 thesimple with probability V2. The reference energy would be the
average 12(ωl + ωl+1). Eq. (24), by contrast, produces ground state energly of the vacuum Eref = 0. To obtain
an average using the weighting function 1 A(ω), which the step function behavior of A>(ω), we project out the
V2
attributes more weight to peaks of A(ω). l positive energy contributions from the initial state |ψ0(cid:105).
Wechoosetheleftboundaryofthefirstintervalω and In Fig. 4, we show the reconstruction of spectral func-
1
therightboundaryofthelastintervalω suchthatthe tionsbasedonthelinearpredictionofthemomentscom-
L+1
distance ωL+1−ω1 is minimized but14 putedforA(cid:101)>discr(ω)usingtheoperator-valuedChebyshev
expansionpresentedinSec.IBforthe“Hamiltonian”de-
(cid:90) ωL+1 (cid:90) ∞ fined in (29). This is analogous to the top left panel of
dωA(ω)≥0.999 dωA(ω), (25)
Fig. 3, which treated the thermodynamic limit.
ω1 −∞
For the finite-size system, the specific choice of dis-
where the integrand is non-negative, which guarantees cretization is important and we compare the linear and
that [ω1,ωL+1] contains almost the complete support of the cosine discretization in the top panels of Fig. 4 for
A>(ω), but minimizes finite-size effects. The intermedi- the expansion order N = 200 and a system size L = 80.
ate values of the discretization intervals {ω2,...,ωL} can From the large error at ω = 0 for the linear discretiza-
be chosen using a logarithmic discretization, as done in tion(reddashedline)seenintherighttoppanelofFig.4,
NRG.18 Butifanunbiased resolutioniswanted,oneusu- which was not present in the thermodynamic limit (red
ally chooses a linear discretization13,14 dashed line in right top panel of Fig. 3), we conclude
that the linear discretization starts resolving finite-size
ω =(l−1)∆ω+ω , l=2,...,L
l 1 features close to ω = 0 already for N = 200. The lower
∆ω = 1(ω −ω ). (26) panels then show how the error behaves as a function of
L L+1 1
the number of computed moments for different system
As Chebyshev polynomials do not show an unbiased sizes L. While the cosine discretization follows the error
energy resolution as they oscillate much quicker at the of the thermodynamic limit quite closely for low values
boundaries of [-1,1] than in the center, the linear dis- of N and lattice sizes of L ≥ 80, it almost saturates in
cretizationwillfirstresolvethefinitesize(discrete)struc- a plateau for higher expansion orders, and only starts
ture close to the boundaries of [−1,1]. We suggest to increasing slightly for very high expansion orders. For
adapt the discretization to account for the cosine map- thelineardiscretization,neithertheclosecorrespondence
ping (20) of the energy scale that is responsible for this with the thermodynamic limit is observed, nor does the
phenomenon. error only moderately depend on the expansion order:
Let us study the case of even L (for odd L, see Ap- Instead, the error increases exponentially for high values
pendix D) and assume without loss of generality that we of N, as then, finite-size features are inhomogeneously
want as many intervals {ωl,ωl+1} on the positive half- resolved. Both features make it difficult to determine
axis as on the negative half-axis, which implies the value of N for which the computation of Chebyshev
momentsshouldbestoppedinordertoobtainaminimal
ω =0. (27)
L/2 error.
8
Asanexample,wecomputethespectralfunctionofthe
1.2 ilcnionpsut )| 10-1 lcions saisngalecoimmmpuornitybeAncnhdmerasrokn13m,14o,d16e,l22(,2S3IAaMnd),iswhhiigchhlyserrevlees-
1 (ut vant as it is at the core of dynamical mean-field theory
p
()0.8 )-Ain10-2 (DTMhFeTH)a.2m4–i2lt7onian of the SIAM is given as,
A
v0.6 (nst
0.4 Areco10-3 HSIAM =Himp+Hbath+Hhyb, (30)
v| H =U(cid:0)n − 1(cid:1)(cid:0)n − 1(cid:1),
0.2 imp 0↑ 2 0↓ 2
0 0 1 2 3 4 5 10-4 0 1 2 3 4 5 Hbath =(cid:88)Lb (cid:88)(cid:15)lc†lσclσ,
/v /v
l=1 σ
1) 1) H =(cid:88)Lb (cid:88)(cid:16)V c† c +H.c.(cid:17).
-v -v hyb l 0σ lσ
)|( 10-1 )|( 10-1 l=1 σ
(put (put By a unitary transform effected by Lanczos tridiagonal-
n n
Ai Ai
)- )- ization,thiscanbemappedontheso-calledchain geom-
(nst10-2 (nst10-2 etry. Butasthisleadstohigherentanglement,wesimply
eco L=40 eco L=40 orderbathstatesbytheirpotentialenergywhichdirectly
x|Ar LL==81020 x|Ar LL==81020 gives a one-dimensional array that can be treated with
ma 10-3 ccooss L=160 ma 10-3 lliinn L=160 MPS.23 We solve the model for the semi-elliptic bath
L= L=
density of states
102 103 102 103
N N
(cid:114)
1 1 (cid:16)ω(cid:17)2
− ImΛ(ω)= 4− , (31)
FIG.4. (Coloronline)Reconstructionofthespectralfunction π 2vπ v
of Fig. 1 represented by the discrete Hamiltonian (29). Top
left: Input function and reconstructed functions using L = whichisdiscretizedaccordingtotheprocedurediscussed
80, N = 200, a = 100v, b = −0.995. We compare the linear in Sec. III, and then yields the parameters (cid:15) and V . It
l l
discretization (26) with the cosine (28) discretization. Top is important to realize that here, we discretize the bath
right: Differenceof input andreconstructed functionsofthe hybridization function whereas in Sec. III, we discretized
topleftpanel. Bottomleft: Theerrormaxω≥0|Areconst(ω)− the spectral function. While Sec. III did this to illustrate
A (ω)|versusnumberN ofcomputedChebyshevmoments
input the effect of discretization for a toy model for which the
usingacosinediscretization. Bottomright: Errorforlinear
spectral function was known from the beginning, in the
discretization.
present case, a true many-body compuation is involved.
Inthepresentcase,therelevantdiscretizationparameter
is the bath size L = L−1, and no longer the system
IV. IMPLICATIONS FOR MPS b
REPRESENTATIONS size.14
We compute the spectral function (18) of the impu-
rity Green’s function, where the initial states are single-
Whatistherelevanceofthepreviousresultsformatrix
product state (MPS) based computations of the spec- particle excitations of the ground state: |ψ0>σ(cid:105)=c†0σ|E0(cid:105)
tralfunctionforagivenmatrixproductoperator(MPO) and |ψ0<σ(cid:105) = c0σ|E0(cid:105). As we consider the particle-hole
H?19RepeatedMPOoperationsonMPScreateentangle- and spin-symmetric case of (30), we only need to com-
ment, which eventually makes manipulating and storing pute one Chebyshev recursion; to be precise: |ψ0(cid:105) =
MPS computationally very costly. Manipulations, such c†0↑|E0(cid:105). We compare our results with dynamic DMRG
as applying H to states |t (cid:105) in the recursion (13), or results from Ref. 16, which are believed to be highly re-
n
performing subsequent time evolution steps e−iH∆t, can liable. In particular, we compare computations in the
therefore only be carried out up to a certain recursion formerly suggested setup13,17 that uses the Chebyshev
order n or time t, before hitting an exponential wall in recursionforb=0in(10)andreconstructsthefull spec-
computation cost. For time evolution algorithms, this tral function A(ω) using linear prediction,13 and the one
has long been known,20,21 but this also limits computa- suggested here, that uses b = −0.995 and reconstructs
tions using the Chebyshev recursion.14 theshifted spectralfunctionA(cid:101)(ω)usinglinearprediction.
In the following, we show that the method intro-
duced in the previous sections outperforms the previous InthetopleftpanelofFig.5,weshowcomputationsof
approach:13 it extracts more spectral information from thespectralfunctionoftheSIAMforL=80forN =260
H when creating the same amount of entanglement or, in the b=−0.995 setup, and N =900 in the b=0 setup
which is equivalent up to technical details of the algo- and compare it with the result of Ref. 16. We choose
rithm, using the same computation time. thesetwoexpansionorders,astheyleadtoacomparable
9
0.07 makesgeneralcomparisonsforthespeedupdifficult. But
Raas (2004) b=-0.995 N=260
the lower right panel of Fig. 5 still shows that with only
1.4 b=-0.995 N=260 0.06 b=0 full N=900
b=0 full N=900 )| a few exceptions, the solid (b=−0.995) lines are always
1.2 (as0.05 clearly below the dashed (b = 0) lines. The logarithmic
vA()00..681 ()-AonstRa00..0043 abscissa therefore indicates a high speedup.
0.4 Arec0.02 V. COMPARISON TO TIME EVOLUTION
v|
0.2 0.01
It is interesting to compare the efficiency of the
0 0 available MPS algorithms to extract spectral informa-
0 1 2 3 4 5 -1 0 1 2 3 4 5
/v /v tion from H. Candidates are, aside from the dy-
0.06 0.06 namicDMRG,28 whichisbelievedtobecomputationally
-1v) LL==4800 -1v) LL==81020 highly costly, time evolution and recursive algorithms.
()|( 0.05 LL==112600 ()|( 0.05 Tpohleynloamtteiarlsa1r7ea,nidntphaerLtiacunlcazro,seaxlgpoarnitshiomns.29i,n30CLhaenbcyzsohseivs
as as
Ra0.04 Ra0.04 numericalinstableasthebasisthatitspansloosesitsor-
A A
)- )- thogonality for high numbers of iterations.31 This seems
(nst0.03 (nst0.03 to to disqualify Lanczos as a high-performing candidate.
o o Therefore the main question is whether the Chebyshev
ec ec
Ar0.02 Ar0.02 recursion can more efficiently extract spectral informa-
x| x|
a a tion from H than time evolution algorithms.
m m
0.01 0.01 To answer this question, in the following, we exploit
102 103 101 102 103
the fact that for b = 0 in the limit a → ∞, the Cheby-
N cpu time (min)
shev expansion becomes a Fourier expansion, and the
FIG. 5. (Color online) Comparison of MPS computed spec- Chebyshev states directly describe the time evolved sys-
tral functions in the two setups studied in the previous sec- tem. This is different from the procedure of computing
tions,withdatabyRaaset al.16. Solidlinesrefertothenew the time dependence of a Green’s function via a Fourier
methodthatusesA(cid:101)(ω),dashedlinesrefertothemethodthat transformofitsspectralfunction.17,32–34Forcomparison,
uses the full spectral function A(ω). Top left: For L = 80 we summarize the latter technique in Appendix E.
and two exemplary expansion orders. Top right: Errors of The approximate equivalence of the Chebyshev recur-
comparisonintopleftpanel. Bottom left: Plotofthemax-
sion to time evolution can be used as a novel time evo-
imum error versus expansion order N. Bottom right: Plot
lution algorithm. This is interesting as for long-range
of the maximum error versus computer time.
interacting Hamiltonians H, the MPO representation of
e−iHt isnotavailable,oronlyapproximately.35 Although
it is possible to use so-called Krylov algorithms for such
maximumerror,asshowninthetoprightpanelofFig.5. problems, this requires some programming effort, and is
In the b = −0.995 setup, this maximum error is slightly in general believed to be numerically rather inefficient
smaller. Around ω = 0, by contrast, the error in the as compared ot other time evolution algorithms. Long-
b = −0.995 setup is much smaller. If we compare the rangeinteractingproblemsappeare.g.ifmappingatwo-
computation time that is needed to reach this precision dimensionalsystemonaone-dimensionalchain,orinthe
(max|Areconst(ω)−Ainput| (cid:39) 0.015/v), we find that the solution of a SIAM using a star geometry23 as in (30).
b = −0.995 setup required ∼145min whereas the b = 0
setup required ∼434min. If one makes this comparison
for a slightly larger error (max|Areconst(ω) − Ainput| (cid:39) A. Statement of approximate equivalence
0.025/v), realized for expansion order N = 120 for the
b = −0.995 setup, and for expansion order N = 200 for
The time evolution of a state |ψ (cid:105)
the b = 0 setup, the comparison in computation times 0
reads ∼12min versus ∼160min.
|ψ(t)(cid:105)=exp(−iHt)|ψ (cid:105), (32)
0
When studying the convergence of the maximum er-
ror with respect to expansion order N in the lower left can be approximately linked to the sequence of Cheby-
panel of Fig. 5, we see that this is, after a sharp de- shev vectors generated by starting from |ψ (cid:105) as follows.
0
crease for low expansion orders, not monotonously de- Choose a reference energy E in (10) that is charac-
ref
creasing. The previously mentioned choices, N = 260 teristic for the initial state of the time evolution and the
in the b = −0.995 setup and N = 900 in the b = 0 Chebyshev recursion. When computing the time evo-
setup, bothcorrespondtoaminimumintheoscillations, lution of the Green’s function with |ψ (cid:105) = c†|E (cid:105), one
0 0
as seen when inspecting the green solid (dashed) lines chooses E =E , if |ψ (cid:105) is not an eigenstate, we choose
ref 0 0
for the b=−0.995 (b=0) setup. The non-monotonicity E =(cid:104)ψ |H|ψ (cid:105).
ref 0 0
10
Then define |φ(t)(cid:105) = exp(iE t)|ψ(t)(cid:105) and H = (H − We can now compute the time evolution
ref
E )/a as in (10) in the b = 0 setup. Here a has the
ref
meaning of an inverse time step of unit energy. With |φcos(tn(cid:48))(cid:105)=|φcheb(tn(cid:48))(cid:105)+(cid:15)(tn(cid:48))|ψ0(cid:105),
these definitions, (32) reads as |φ (t )(cid:105)=T (H)|ψ (cid:105), n(cid:48) =0,4,8,... (40)
cheb n(cid:48) n(cid:48) 0
|φ(t)(cid:105)=exp(−iaHt)|ψ0(cid:105) via the recursion (13)
(cid:0) (cid:1)
= cos(aHt)−isin(aHt) |ψ (cid:105)
0 |φ (t )(cid:105)=2H|φ (t )(cid:105)−|φ (t )(cid:105),
cheb n cheb n−1 cheb n−2
≡|φ (t)(cid:105)−i|φ (t)(cid:105). (33)
cos sin |φ (t )(cid:105)=H|ψ (cid:105), n=0,1,2,.... (41)
cheb 1 0
Let us discretize time by defining t = n, then
n a
|φcos(tn)(cid:105)=cos(nH)|ψ0(cid:105), (34a) 1 ReiG>(t) 101 ReiG>(t)
|φsin(tn)(cid:105)=sin(nH)|ψ0(cid:105). (34b) n n n 100 |ReiG>(t)- 4n|
T|ψh0Wi(cid:105)seuissnionnwgotawparoensctsuitrbosliecoonwmtitphhuatttehotenhlsyetaainncvtdioaolrnvdeosrfetchcouesr(asncioHtnio)|nfψoo0rf(cid:105)tHohne. e evolutio 0.05 e evolutio1100--21
cosine function, as shown in Appendix F1. m m
Let us instead consider the action of the Chebyshev rt-ti g-ti10-3
o-0.5 n
polynomials h o
s l10-4
Tn(H)=cos(narccos(H)) (35) -1 10-5
0 0.5 1 1.5 2 2.5 3 0 100 200 300 400
t=n/a (1/v) t=n/a (1/v)
on |ψ (cid:105). This action approximately reproduces the ac-
0
tionoftheplanecosinefunction,ifweconsiderevery4th
FIG.6. (Coloronline)Left: Timeevolutionofaparticlecre-
iteration, i.e. introduce the new index n(cid:48) =4n, n∈N:
atedonthefirstsiteofachainoflengthL=100withhopping
Tn(cid:48)(H)|ψ0(cid:105)=cos(cid:16)n(cid:48)(cid:16)π2 −H(cid:17)(cid:17)|ψ0(cid:105)+(cid:15)(n(cid:48))|ψ0(cid:105) (36) vtio=n o1ftthhaetGisre|eψn0’(cid:105)s =funcc†0t|ivoanc(cid:105)i.G>W(et)c=omepiEa0rte(cid:104)ψth0|eψt(itm)(cid:105)e(sehvoowlun-
as blue crosses) with the Chebyshev moments µ = (cid:104)ψ |ψ (cid:105)
(cid:16) (cid:17) n 0 n
=cos n(cid:48)H |ψ (cid:105)+(cid:15)(n(cid:48))|ψ (cid:105). n(cid:48) =0,4,8,... (shownasgreendots)obtainedwhencomputingtherecursion
0 0
|ψ (cid:105) = 2H/a|ψ (cid:105)−|ψ −2(cid:105). Hopping amplitudes v are
n n−1 n l
obtained from the discretization of the spectral function of
In the first line, we used the Taylor expansion
arccos(H)= π −H+ 1H3+···, that leads to the error Fig. 1. A qualitatively equivalent behavior is obtained for a
2 6 chain with homogeneous hopping v = v. Right: Long-time
function(cid:15)(n(cid:48))(AppendixF5),andinthesecondline,we l
evolution of the Green’s function (blue crosses), and differ-
usedn(cid:48)π =2πn,n∈N,whichobviouslydropsoutofthe
2 enceoftheGreen’sfunctionandtheChebyshevmoments(red
argument of the cosine (also see Appendix F2).
dots). The horizontal dashed line marks the the prefactor of
The error (cid:15)(n(cid:48)) is bounded by (F16) theerrorestimate(37),whichiscomputedasσ3/a2 =8·10−4.
t
|(cid:15)(t )|= if t<t
n(cid:48) t err
err
a2 B. Numerical examples
t = (37)
err σ3
1. Single-particle computation for SIAM
where σ is the spectral width of the initial state |ψ (cid:105) ≡
0
|ψ (cid:105) around E ,
0 ref
Fig. 6 shows the numerically exact time evolution of a
σ = max |E −E |, (38) single particle created on the first site |ψ (cid:105) = c†|E (cid:105) of
k ref 0 0 0
|Ek(cid:105)∈|ψ0(cid:105) a chain of 100 lattice sites. The spectral width therefore
is σ = 2v. The left panel of Fig. 6 plots the Cheby-
where “|E (cid:105) ∈ |ψ (cid:105)” refers to the decomposition of the
k 0
shev moments µ =(cid:104)ψ |t (cid:105) obtained with a=100v and
initial state in eigenstates |E (cid:105) of H n 0 n
k
thetimeevolutionofthecorrespondingGreen’sfunction
|ψ0(cid:105)=(cid:88)ck|Ek(cid:105). (39) iG>(t)=eiE0t(cid:104)ψ0|ψ(t)(cid:105) for the time step a4. Every forth
Chebyshev moment agrees with a value of the Green’s
k
function. In the right panel, we show the long time be-
The“spectralwidth”σ isusuallysmallcomparedtorea- havior of G>(t) and the difference G>(4n/a)−µ . The
4n
sonably high values of the inverse time step a. If one differenceisseentobeclearlybelowtheconservativeup-
is unsure of whether a was chosen large enough, one re- perbound(37),itremainsoftheorderofσ3/a2 =8·10−4
runs a calculation with a higher value of a and checks up to very high times that correspond to 100 hopping
convergence. processes (t=100/v).
