Chebyshev matrix product state approach for spectral functions

1 MB·

by  Andreas Holzner

Chebyshev matrix product state approach for spectral functions
Andreas Holzner,1 Andreas Weichselbaum,1 Ian P. McCulloch,2 Ulrich Schollw¨ock,1 and Jan von Delft1
1Physics Department, Arnold Sommerfeld Center for Theoretical Physics,
and Center for NanoScience, Ludwig-Maximilians-Universit¨at Mu¨nchen, D-80333 Mu¨nchen, Germany
2School of Physical Sciences, University of Queensland, Brisbane, Queensland 4072, Australia
(Date: February 1, 2011)
We show that recursively generated Chebyshev expansions offer numerically efficient represen-
tations for calculating zero-temperature spectral functions of one-dimensional lattice models using
matrix product state (MPS) methods. The main features of this Chebychev matrix product state
1 (CheMPS)approachare: (i)itachievesuniformresolutionoverthespectralfunction’sentirespectral
1 width;(ii) it can exploit thefact that thelatter can bemuchsmaller thanthemodel’s many-body
0 bandwidth; (iii) it offers a well-controlled broadening scheme that allows finite-size effects to be
2 either resolved or smeared out, as desired; (iv) it is based on using MPS tools to recursively calcu-
lateasuccessionofChebychevvectors|t i,(v)whoseentanglement entropies werefoundtoremain
n n
bounded with increasing recursion order n for all cases analyzed here; (vi) it distributes the total
a
J entanglement entropy that accumulates with increasing n over the set of Chebyshev vectors |tni,
which need not be combined into a single vector. In this way, the growth in entanglement entropy
1
thatusuallylimitsdensitymatrixrenormalization group(DMRG)approachesispackagedintocon-
3
veniently manageable units. We present zero-temperature CheMPS results for the structure factor
] of spin-21 antiferromagnetic Heisenberg chains and perform a detailed finite-size analysis. Making
l comparisons to three benchmark methods, we find that CheMPS (1) yields results comparable in
e
- qualitytothoseofcorrectionvectorDMRG,atdramaticallyreducednumericalcost;(2)agreeswell
r with Bethe Ansatz results for an infinite system, within the limitations expected for numerics on
t
s finitesystems;(3)canalsobeappliedinthetimedomain,whereithaspotentialtoserveasaviable
t. alternative to time-dependent DMRG (in particular at finite temperatures). Finally, we present a
a detailed error analysis of CheMPS for thecase of the noninteracting resonant level model.
m
- PACSnumbers: 02.70.-c,75.10.Pq,75.40.Mg,78.20.Bh
d
n
o I. INTRODUCTION to achieve an accuracy beyond that of the expected L-
c dependent errors, and having convenient control of this
[
Consider a one-dimensionallattice model amenable to accuracy can significantly reduce numerical costs.

1 treatment by the density matrix renormalization group In this paper, we show that Chebyshev expansions
v (DMRG),1–4 with Hamiltonian Hˆ, ground state 0 and offer numerically efficient representations for calculat-
5 ground state energy E . This paper is concerne|diwith ing spectralfunctions using matrix product state (MPS)
9 0
8 zero-temperature spectral functions of the form methods,4,6–10 with numerical costs that compare fa-
5 vorably to those of other established DMRG-based ap-
. BC(ω)= 0 ˆδ(ω Hˆ +E0) ˆ0 , (1) proaches. In particular, the Chebychev MPS approach
1 A h |B − C| i
presented here, to be called CheMPS, allows the above-
0 whichrepresentstheFouriertransform dteiωtGBC(t)of
1 2π mentioned control of accuracy and resolution to be im-
the correlator
1 R ported into the DMRG/MPS arena.
v: GBC(t)= 0 ˆ(t)ˆ(0)0 . (2) The historically first approach for calculating spec-
h |B C | i tral functions with DMRG is the continued-fraction
i
X One possible framework for calculating such spectral expansion.11Whilethismethodrequiresonlymodestnu-
r functions isto expandthemintermsofChebychevpoly- merical resources, it is limited to low frequencies and it
a
nomials, as advocated in Ref. 5. Such a Chebyshev ex- is difficult to produce reliable results with it in the case
pansionofferspreciseandconvenientcontroloftheaccu- ofcontinua(however,algorithmicimprovementswerere-
racyandresolutionwithwhichaspectralfunctionistobe ported recently12). At present, the most accurate, but
computed. This is very useful, particularly when broad- alsomosttime-consumingapproachesare: (i)thecorrec-
ening the spectral function of a length-L system, which tion vector (CV) method,13–15 and (ii) time-dependent
exhibits finite-size subpeaks with spacing ω 1/L, in DMRG (tDMRG),7,9,16–18 in particular when combined
L
order to mimick that of an infinite system: i∼f the lat- withlinearpredictiontechniques.19–22 Sinceanynewap-
ter has structures (e.g. sharp or diverging peaks) which proachmustmeasureuptotheirstandards,letusbriefly
are not yet properly resolvedat the scale ω , the broad- summarize their key ideas, advantages and drawbacks.
L
enedversionofthefinite-sizespectralfunctioninevitably (i) To calculate BC(ω) using the CV approach, it is
A
bears L-dependent errors in the vicinity of these struc- expressed as
tures. Hence, when calculating the finite-size version of
thesestructuresforthelength-Lsystem,thereisnoneed BC(ω)= 0 ˆ , (3a)
ω
A h |B|Ci
2
in terms of the so-called correction vector (a)
1 1
lim Im ˆ0 . (3b)
|Ciω ≡−η→0π ω Hˆ +E +iη C| i
(cid:20) − 0 (cid:21)
Thecorrectionvectorcanbecalculated(forfinite broad- (b)
eningparameterη)usingeitherconventionalDMRG13–15
orvariationalmatrix productstate(MPS) methods.23 A
majoradvantageofthis approachisthatarbitrarilyhigh
spectral resolution can be achieved by reducing η and
sampling enough frequency points. However, this comes
at considerable numerical costs: first, a separate calcu- Fig.1. (a)Sketchofaspectralfunctionwhosespectralwidth
lation is required for every choice of ω (though in doing WA is much smaller than the many-body bandwidth W. Be-
fore making a Chebychev expansion, we rescale the interval
so,resultsfor ’s frompreviousfrequencies canbe in-
corporated); a|nCdiωsecond, the calculation of |Ciω involves ωint∈erv[a0l,Wω∗′],∈w[i−thWe′ff,eWct′i]v,eshboawndnwiindt(hb)W,∗w=ith2rWeAsc,aolendtohathlfe-
aconnodpiteiroanteodr,inevveerrsmioonrpersoobltehme stmhaatlliesrnηumis.ericallypoorly bandwidth W′=1− 21ǫt and a safety factor ǫt ≃0.025.
(ii) An alterative possibility is to use tDMRG to cal-
culate the time-domain correlatorGBC(t), Fourier trans-
needforcalculatingasinglestatewithsuchhighaccuracy
forming to the frequency domain only at the very end.
andinsteadallowsnumericalresourcestobe focusseddi-
To this end, one expresses
rectly on the calculation of the relevant expectation val-
GBC(t)=eiE0t 0 ˆ ˜ (4a) ues. This can be achieved by representing the spectral
h |B|Cit function via a Chebychev expansion,5,24,25 whose coef-
in terms of the time-evolved state ficients, the so-called Chebyshev moments, can be cal-
culated recursively using MPS tools. Below, we briefly
˜ e−iHˆtˆ0 (4b)
t summarizethestructureandmainfeaturesofsuchanex-
|Ci ≡ C| i
pansion, thereby providing both an introduction and an
and uses tDMRG to calculate the latter. Two attractive
overview of the material developed in detail in the main
features of this strategy are: first, it builds on an exten-
part of this paper.
sive body of algorithmic knowledge for efficiently calcu-
lating time-evolution;7,16,17 and second, a simple linear- The Chebychev polynomials Tn(x) form an orthonor-
prediction scheme19–22 can be used to extrapolate the amraelvseertyowfeplloslytundoimediamlsaotnhetmheatiinctaelrlyv,a2l6–x28∈a[n−d1a,r1e].wTidheelyy
time-dependence calculated for short and intermediate
used for function expansions since they have very favor-
time scales to longer times, thereby improving the qual-
able convergence properties. As will be described in de-
ity of results at low frequency at hardly any additional
tail below, the spectral function can be represented ap-
numerical cost. However, obtaining reliable results over
proximately by a so-called Chebychev expansion, which
a sufficiently large time interval can, in itself, be nu-
becomes exact for N , of the following form:
merically very expensive, since the time-evolution of the →∞
many-bodystate ˜ isaccompaniedbya stronggrowth
in entanglement e|nCtirtopy. This unavoidably also implies BC(ω)= 2W′/W∗ g µ +2N−1g µ T (ω′) . (5)
a growth of tDMRG truncation errors. AN π√1 ω′2 " 0 0 n n n #
Notethatinbothoftheschemesoutlinedabove,signif- − nX=1
icant (often heroic) amounts of numerical resources are HeretheChebyshev moments µ = 0 ˆt areobtained
n n
dorevo˜tedfotrogciavlecnulat,tiansgaaccsuinragtleelystaatse,po|Cssiiωblefo;rthgeivoenverω- from the Chebychev vectors |tni =hTn|B(H|ˆ′)iCˆ|0i, and the
t g are known damping factors that influence broaden-
laps|Coir expectation values of interest, namely 0 ˆ ω inng effects. The primes indicate that the Hamiltonian
h |B|Ci
for 0 ˆ ˜t, are only calculated at the end, in a single, Hˆ and frequency ω were expressed in terms of rescaled
finahl s|tBe|pC,iafter or ˜ have been fully determined. and shifted versions, Hˆ′ and ω′, in such a manner that
ω t
|Ci |Ci
Actually, these states are calculated so accurately that an interval ω [0,W ] that contains the entire spectral
∗
theywouldhavebeenequallysuitableforcalculatingany weight is map∈ped onto a rescaled band ω′ [ W′,W′]
other quantity (correlator or matrix element) involving of halfwidth W′ <1. ∈ −
that state. In a sense, DMRG is asked to work harder This representation has several useful features:
than necessary: it is usedto calculate a single state with (i)Itresolvestheintervalω [0,W ]withauniform res-
∗
∈
“general-purposeaccuracy”,whereas the accurate calcu- olution of (W /N).
∗
O
lation of a particular expectation value involving that (ii)Therangeoffrequenciesoverwhichthespectralfunc-
state would have been sufficient. tionhas nonzero weight,sayW (to be calledits spectral
A
The main motivation for the present work is to at- width) is often significantly smaller than the many-body
tempttoreducethiscalculationaloverheadbyemploying bandwidth of the Hamiltonian, say W, as depicted in
a representation of the spectral function that avoids the Fig 1. By choosing the effective bandwidth W to be of
∗
3
order W instead of W, huge gains in resolutionare pos- and 8, respectively). In Sec. V we perform an exten-
A
sible. sive error analysis of the CheMPS approach using the
(iii)Awell-controlledbroadeningscheme,encodedinthe quadratic resonantlevel model, and discuss some salient
damping factors g , is availablethat allowsfinite-size ef- features of density matrix eigenspectra in Sec. VI. Sec-
n
fects to be either resolved or smeared out, as desired. tionVII summarizesourmainconclusions,andSec.VIII
(iv) The Chebyshev vectors t are calculated using a presents a brief outlook towards possible future appli-
n
| i
(numerically stable) recursion scheme, which exploits cations,involvingtime dependence orfinite-temperature
Chebychev recurrence relations to calculate t from correlators. An Appendix gives a detailed account of
n
H′ t and t (seeEq.(30)). Thus,theex|pecitation CheMPS results for the resonant level model used for
n−1 n−2
| i | i
valuesfromwhichthespectralfunctionisconstructedare the error analysis of Sec. V.
built up in a series of recursive steps (see Eq. (7) below)
instead of being calculated at the end in one final step.
(v) The bond entropy of successive Chebyshev vectors II. CHEBYSHEV EXPANSION OF ABC(ω)
t isfoundempiricallytoremain bounded withincreas-
n
| i
ingrecursionnumbern,thusthecomplexityofthesevec-
A. Chebyshev basics
tors remains managable up to arbitrarily large n.
(vi)Finally,andfromthe perspectiveofnumericalcosts,
Letusstartbybrieflysummarizingthosepropertiesof
most importantly: CheMPS efficiently copes with the
Chebyshev polynomials that will be needed below. We
growth in bond entropy with increasing iteration num-
follow the notation of Ref. 5, which gives an excellent
ber that usually limits DMRG approaches. It does so by
general discussion of Chebyshev expansion techniques
distributing this entropy over all t , thereby packaging
n
| i (though without mentioning possible DMRG/MPS ap-
it into managable units (see (v)). In particular, when
plications).
constructing and using the states t , one never needs
n
| i Chebyshevpolynomialsofthefirstkind,T (x), hence-
to know more than three at a time (and after use may n
forth simply called Chebyshev polynomials, are defined
delete them frommemory). Hence, it is not necessaryto
by the recurrence relations
combineallinformationcontainedinall t intoasingle
n
| i
MPS.
T (x)=2xT (x) T (x),
Let us constrast this with the CV or tDMRG ap- n+1 n − n−1 (8)
proaches: imagine expanding the correction vector or T0(x)=1, T1(x)=x.
time-evolvedstateintermsoftheChebyshevvectors t ,
n
i.e. expressing them as linear combinations of the for|mi They also satisfy the useful relation (for n n′)
≥
N−1 N−1 T (x)=2T (x)T (x) T (x). (9)
n+n′ n n′ n−n′
Cn t , ˜ C˜n t , (6) −
|Ciω ≃ ω| ni |Cit ≃ t | ni
n=0 n=0 Two useful explicit representations are:
X X
respectively. (The coefficients Cωn and C˜tn are related Tn(x)=cos[narccos(x)]=cosh[narccosh(x)]. (10)
by Fourier transformation.) Now, the CV or tDMRG
approaches in effect attempt to accurately represent the On the interval I = [ 1,1] the Chebyshev polynomials
entire linear combination using a single MPS. This en- constitute an orthogon−al system of polynomials (over a
devourisnumericallyverycostly,sincetheentanglement weight function (π√1 x2)−1), in terms of which any
entropyofthislinearcombinationgrowsrapidlywithN. −
piecewise smooth and continuous function f(x) can
The Chebychev approach avoids this problem by taking |x∈I
be expanded. In fact, the T (x) are optimally suited for
n
expectation values before performing the sum on n:
thispurpose,sincetheyhavetheuniqueproperty(setting
them apart from other systems of orthogonal polynomi-
N−1 N−1
0 ˆ ˆ Cnµ , 0 ˆ ˜ C˜nµ . (7) als) that on I their values are confined to Tn(x) 1,
h |B|Ciω ≃ ω n h |B|Cit ≃ t n withallextremalvaluesequalto1or 1. Th|isisev|i≤dent
nX=0 nX=0 from the first equality in Eq. (10); t−he second equality
Thus, the Chebychev expansion very conveniently orga- implies that for x / I, T (x) grows rapidly with in-
n
nizes the calculation into many separate, and hence nu- creasing x. These∈prop|erties |are illustrated in Fig. 2.
merically less costly, packages or subunits. | |
Our paper is organized as follows. We introduce the There are several ways of constructing Chebyshev ap-
ChebyshevexpansionforspectralfunctionsinSec.IIand proximations for f(x) (see Weisse et al.,5, Section
discussitsimplementationusingMPSincludinganewal- II.A). The Chebyshev|xe∈xIpansion that is practical for
gorithmforperformingaprojectioninenergyinSec.III. present purposes has the form
In Sec. IV we present CheMPS results for the structure
factor of a spin-1 Heisenberg chains, perform a detailed ∞
2 1
analysisoffinite-sizeeffects(seeFig.5),andcompareour f(x)= µ0+2 µnTn(x) , (11)
results to CV, Bethe Ansatz and tDMRG (see Figs. 4, 6 π√1−x2 " nX=1 #
4
(a) mostly employ Jackson damping, given by
1
(N n+1)cos πn +sin πn cot π
gJ = − N+1 N+1 N+1. (15)
n N +1
(x)n 0
T This is usually the best choice,since it guaranteesanin-
tegratederrorof (1)forf (x). Whenusedtoapproxi-
O N N
mateaδ-functionδ(x x¯)sittingatx¯ I,Jacksondamp-
-1 ingyieldsanearlyGa−ussianpeakofw∈idth√1 x¯2π/N.
-1 -0.5 0 0.5 1 −
On one occasion we will also employ Lorentz damping,
x
(b) 102 gL = sinh λ 1− Nn , (16)
T0 T3 T6 n,λ sinhλ
T T T (cid:2) (cid:0) (cid:1)(cid:3)
1 4 7
T T T where λ is a real parameter. Lorentz damping preserves
 2 5 8
(x)n101 analytical properties (causality) of Green’s function and
T broadens a δ-function δ(x x¯) into a peak whose shape,

−
forthechoiceλ=4usedhere(followingRef.5),isnearly
Lorentzian, of width √1 x¯2λ/N.
100 −
Tosummarize: the order-N Chebychev-reconstruction
-1.5 -1 -0.5 0 0.5 1 1.5
f (x) with Jackson or Lorentzian damping with λ = 4
N
x
yieldsaresultthatisveryclosetothebroadenedfunction
Fig.2. (Coloronline)Chebyshevpolynomialsofthefirstkind, 1
Tarne(xlo),cafotrednwuipthtoin8t.h(ea)inAtellrvzaerloIs=and[−e1x,t1r]e,maandofaelvleerxytTrenm(xa)l fNX(x)= dx¯KηXN′ ,x¯(x−x¯)f(x¯), (17)
Z
values equal 1 or −1. (b) Chebyshev polynomials |T (x)|for −1
n
x∈[−1.5,1.5]. The|Tn>0(x)|growrapidlywhen|x|increases (X =J,L) with broadening kernels and widths given by
beyond 1.
e−x2/(2(η′2) π
where the Chebychev moments µn are given by KηJ′(x)= √2πη′ , ηN′ ,x¯ = 1−x¯2 N , (18a)
1 η′/π p 4
KL(x)= , η′ = 1 x¯2 , (18b)
µ = dxf(x)T (x). (12) η′ x2+η′2 N,x¯ − N
n n
p
Z
−1 respectively. Thus,f (x)resolvestheshapeoff(x)with
N
a resolution of (1/N).
An approximate representation of order N is obtained O
For purposes of illustration, Fig. 3(a) shows three
for f(x) if only the first N terms (i.e. n N 1) are re-
≤ − Chebyshev reconstructions of a δ-function at x¯ = 0:
tained. However,suchatruncationingeneralintroduces
without damping, giving Gibbs oscillations; with Jack-
artificial oscillations, of period 1/N, called Gibbs os-
≃ son damping, yielding a near-Gaussian peak; and with
cillations. These can be smoothened by employing cer-
Lorentz damping, yielding a near-Lorentzian peak. Fig-
tainbroadeningkernels,whichineffectrearrangethe in-
ure3(b)showsaJackson-dampedChebyshevreconstruc-
finite series (11) before truncation. This leads to a re-
tionofacombofGaussianpeaks KJ (x x¯ ),whose
constructed expansion of the form α η¯α′ − α
widths η¯′ are all equal. It illustrates how increasing N
α P
N−1 reducestheamountofbroadeninguntiltheoriginalpeak
1
fN(x)= g0µ0+2 gnµnTn(x) , (13) form is recovered for sufficiently large N. It also shows
π√1−x2 " nX=1 # that the broadened peak widths depend on the peak po-
sitions, reflecting the fact that convolving a Gaussian of
which(for properlychosenkernels)convergesuniformly:
widthη¯′ withanear-Gaussianofwidthη′ (Eq.(17))
α N,x¯α
produces a near-Gaussianof width
N→∞
max f(x) f (x) 0. (14)
N
−1<x<1| − | −→
η′ η¯′2+η′2 . (19)
α ≃ α N,x¯α
The reconstructed series (13) contains the same Cheby-
q
shevmomentsµ asEq.(12), buttheyaremultipliedby ToevaluatetheT (x)thatoccurinEq.(13),weusethe
n n
damping factors g , real numbers whose form is charac- first equality of Eq. (10). Although numerically more ef-
n
teristic of the chosen kernel. Several choices have been ficientmethodsexistforthispurpose,5 theirusebecomes
proposed, which damp out Gibbs oscillations in some- advisableonlyforexpansionordersmuchlargerthanthe
what different ways (see Ref. 5 for details). We will N < (103) that we will need in this work.
∼O
5
(a) slightlysmallerthantherequisite[ 1,1]bychoosingthe
8 gn 1 L, λ=4 undampeδdJ rseasfecatyledfachtaolfr-5baonfdǫwtidt0h.0W25′.tToobbees−mexapllliecritt,hwane d1e,fiwnieth a
J ≃
6 0.5 KJ(x) ω W
δL ω′ = W′ , a= ∗ , (21a)
f(x)N 4 0 0 2n5 50 KL(x) N=50 Hˆ′ = Haˆ−−E0 W′ , 2W′ (21b)
a −
2
where Hˆ′ has ground state energy E′ = W′. Then we
0 −
0 express the spectral function (1) as
-0.4 -0.2 0 0.2 0.4 BC(ω)= 1 0 ˆδ(ω′ Hˆ′) ˆ0 , (22)
x A ah |B − C| i
N=1000 N=100 N=50
(b)
20 (with ω′ = ω′(ω) and Hˆ′ given by Eqs. (21)), which by
construction has no weight for ω′ / [ W′,W′].
∈ −
One possible choice for W is to equate it to the width
15 ∗
ofthe many-bodyspectrum of H, givenby W =E
max
−
E . When using DMRG, E is usually already known
x) 0 0
Jf(N 10 fbreomfoucnaldc,uela.gti.n,gbtyhecaglrcouulantdinsgta29teth|0eigorfoHun,dansdtaEtemoafx caHn
−
(reduced DMRG accuracy relative to usual groundstate
5
calculations is sufficient, since only E is of interest
max
here.)
0 AdisadvantageofthechoiceW =W isthatthemany-
∗
-1 -0.5 0 0.5 1
body bandwidth W typically is large (it scales with sys-
x
tem size), whereas optimal spectral resolution requires
W tobeassmallaspossible: sinceanN-thorderCheby-
Fig. 3. Three Chebyshev reconstructions of δ(x), with N = ∗
shevexpansionyieldsaresolutionof (1/N)ontheinter-
50: the undamped case (gn = 1) yields Gibbs oscillations O
val [ 1,1], its resolution on the original interval [0,W ]
(central peak has height δ50(0) = 16.23); Jackson damping − ∗
(δJ) mimicks a Gaussian peak KJ(x) of width π/N; Lorentz will be (W∗/N), which evidently becomes better the
O
damping(δL) forλ=4 mimicksaLorentzian peak KL(x)of smaller W∗. If ˆand ˆare single-particle operators, the
width λ/N. Inset: Jackson and Lorentz damping factors, gJ spectral widthBW ofC BC(ω) is independent of system
n A
and gnL,λ=4, respectively, plotted for N = 50. Both decrease size and hence much smAaller than the many-body band-
monotonically from 1 to 0, but in somewhat different ways. width W. In this case, it is advisable to chooseW to be
∗
(b) Jackson-damped reconstruction of a comb of normalized
of similar order (though still larger) than W . We will
Gaussians (dashed line), all of width η¯′ =0.02, for three val- A
chooseW =2W ,whichis typically W,asillustrated
uesofN (solidlines). Thex-dependenceofthepeakheightsis ∗ A ≪
in Fig. 1.
givenby[2π(η¯′2+η′2 )]−1/2 (dash-dottedline),seeEq.(19).
N,x
C. Chebychev expansion in frequency domain
B. Rescaling of ω and Hˆ
To expand the δ-function in Eq. (22) in Chebyshev
To construct a Chebyshev expansion of the spectral
polynomials, we use fˆ(x) = δ(x Hˆ′) with x = ω′
function BC(ω)ofEq.(1),weneedtorescaleandshift5 −
A in Eq. (12), and obtain from Eq. (13) a reconstructed
the Hamiltonian Hˆ Hˆ′ and the frequency ω ω′ in
7→ 7→ Chebyshev operator expansion of the form:
suchawaythatthespectralrangeof (ω),i.e.theinter-
A
val[0,WA]withinwhichithasnonzeroweight,ismapped 1 N−1
into theinterval[ 1,1]. Rescaled,dimensionlessenergies δ (ω′ Hˆ′)= g +2 g T (Hˆ′)T (ω′) .
and frequencies w−ill always carry primes. As safeguards N − π√1 ω′2 " 0 n n n #
− nX=1
against “leakage” beyond [ 1,1] due to numerical inac- (23)
curacies, we choose the line−ar map (see Fig. 1) Insertingthis intoEq.(22)for BC(ω)yields the Cheby-
A
shev expansion (5), with Chebyshev moments given by
ω [0,W ] ω′ [ W′,W′], W′ =1 1ǫ , (20)
∈ ∗ 7→ ∈ − − 2 t µ = 0 ˆT (Hˆ′)ˆ0 . (24)
n n
which entails two precautionary measures: first, the ω- h |B C| i
intervalistakentobelargerthantherequisite[0,W ]by Thusµ isagroundstateexpectationvalueofann-thor-
A n
choosingtheeffectivebandwidth W tobelargerthanthe der polynomial in Hˆ′, whose constructionmight a priori
∗
spectral width W ; second, the ω′-interval is taken to be appear to become increasingly daunting as n increases.
A
6
Fortunately, this challenge can be dealt with recursively, A. Recurrence fitting
by expressing the moments as
To initialize the Chebyshev expansion, we calculate
µn = 0 ˆtn , tn =Tn(Hˆ′)ˆ0 , (25) groundstate 0 andgroundstate energy E of Hˆ, make
h |B| i | i C| i 0
| i
aspecificchoiceforW andW′,andconstructHˆ′accord-
and calculating the Chebyshev vectors t by exploiting ∗
| ni ing to Eq.(21b). Then comes the maintask, namely the
the Chebyshev recurrence relations (8). The details of
recursive calculation of the moments µ . This is done
this recursive scheme will be discussed in Section III. n
starting from
t = ˆ0 , t =Hˆ′ t , (29)
0 1 0
D. Chebychev expansion in time domain | i C| i | i | i
andusingtherecurrencerelation(obtainedfromEq.(8))
The Chebyshev expansion can also be employed for
t =2Hˆ′ t t . (30)
studying time evolution in general, and the correlator | ni | n−1i−| n−2i
GBC(t) in particular. To this end, we express the time-
Eq.(30)canbeimplementedusingtheso-calledcompres-
evolution operator as sionorfittingprocedure32(see4,Sec.4.5.2fordetails). It
finds an MPS representation for t , at minimal loss of
1 n
| i
Uˆ(t)=e−iHˆt = dω′e−i[a(ω′+W′)+E0]tδ(ω′ Hˆ′), (26) informationforgivenMPSdimensionm,byvariationally
− minimizing the fitting error
Z
−1
2
and insert Eq. (23) (without damping, gn = 1) into the ∆fit = |tni− 2Hˆ′|tn−1i−|tn−2i . (31)
latter. This yields30,31 (cid:13) (cid:16) (cid:17)(cid:13)
Wewillcallthis(cid:13)procedurerecurrencefitting.(cid:13)Inpractice,
(cid:13) (cid:13)
N−1 the variational minimization proceeds via a sequence of
UˆN(t)=e−i(E0+aW′)t c0(t)+2 Tn(Hˆ′)cn(t) , (27a) fitting sweeps back and forth along the chain. These
" #
n=1 are continued until the state being optimized becomes
X
1 stationary, in the sense that the overlap
e−iatω′T (ω′)
c (t)= n dω′ =( i)nJ (at). (27b)
n −Z1 π√1−ω′2 − n ∆c = 1− thtn|t′nti′ (32)
(cid:12) k| nikk| nik(cid:12)
HereJ (at)istheBesselfunctionofthefirstkindoforder (cid:12) (cid:12)
n between the states t (cid:12) and t′ before a(cid:12)nd after one fit-
n. It decays very rapidly with n once n > at. Hence, | n(cid:12)i | ni (cid:12)
ting sweep, drops below a specified fitting convergence
an expansion of given order N gives an essentially exact
threshold (typically in the range 10−6 to 10−8). The
representation of Uˆ(t) for times up to t < N, while
max a maximum expansion order for which t is obtained us-
cN−1(t) provides an estimate of the error. ∼ ing recurrence fitting will be denoted|bnyiN .
Inserting Eqs. (27) into Eqs. (4) for GBC(t) we find max
The MPS dimension m needed to achieve accurate re-
currence fitting turns out to be surprisingly small (see
N−1
GBC(t)=e−iaW′t µ J (at)+2 ( i)nµ J (at) , Sec.Vforadetailedanalysis). Forexample,m=32suf-
N " 0 0 − n n # ficedfortheantiferromagneticHeisenbergchainoflength
n=1
X (28) L=100 discussed in Section IV. The reason for this re-
where the Chebyshev moments µ are again given by markableandeminentlyusefulfeatureliesinthefactthat
n
Eq. (25). Thus the Chebyshev expansions of GBC(t) and theChebychevrecurrencerelations(30)containonlytwo
BC(ω) aregovernedby the same set ofmoments µ , as terms on the right-hand side, whose addition requires
n
iAs to be expected for functions linked by Fourier trans- only modest computational effort. In contrast, CV or
formation. tDMRG typically require much larger m, since they at-
tempt to represent the sum of many states, see Eq. (6),
in terms of a single MPS.
III. MPS EVALUATION OF THE CHEBYSHEV For the special but common case that ˆ= ˆ†, Eq. (9)
B C
MOMENTS µ yields a relation between different moments,
n
µ =2 t t µ . (33)
n+n′ n n′ n−n′
We now present a recursive scheme for calculat- h | i−
ing the Chebyshev moments µn. The manipulations This can be used to effectively double the order of the
described below were implemented using MPS-based expansion to 2N without calculating any additional
max
methods,4,6–10whichareveryconvenientforconstructing Chebyshev vectors, by setting n′ =n 1 or n:
the states of interest, while matrix-product operators10 −
(MPOs) simplify the implementation of the shift- and µ˜2n−1 =2 tn tn−1 µ1,
h | i− (34)
rescaling transformation Eq. (21b) of the Hamiltonian. µ˜ =2 t t µ .
2n n n 0
h | i−
7
We use tildes to distinguish µ˜ -moments calculated in choiceofW . WehavefoundthecombinationW =2W
n ∗ ∗ A
thismannerfromtheµ -momentsobtainedviaEq.(25). and ε = 1.0 to work well (but other choices, involving,
n P
Although they should nominally be identical, in numer- e.g. smaller W and larger ε would be possible, too.)
∗ P
ical practice µ˜ -moments are less accurate (by up to a In the following, we describe the procedure just out-
n
factor of 5 in Fig. 9(c) below), since they depend on linedinmoredetailforasinglesite,usingstandardMPS
two Chebyshev vectors, whereas µ -moments depend on nomenclature. Let the effective local Hilbert space for
n
only one. Our Chebyshev reconstructions thus gener- thissitebespannedbythe left,localandrightbasisvec-
ally employ the µ -moments, and unless stated other- tors l , σ and r , and expand the Chebyshev vector
n
| i | i | i
wise, µ˜ -moments are used only for results requiring ψ = t in this basis:
n n
| i | i
N n<2N .
max max
≤
ψ = A[σ] l σ r . (35)
| i lr | i| i| i
lσr
B. Energy truncation X
To construct a projection operator P that projects out
We have argued above that in order to optimize spec- the high energy components for this site, ψ P ψ ,
| i 7→ | i
tral resolution, it may be desirable to choose the effec- one may proceed as follows:
tbiavnedbwainddtwhiWdth. WIf∗thtoisbiessdmonalel,erhothwaenvetrh,eitfuilslmesasennyt-ibaoldtyo spaFnirstl,bσuirld aaKnrdylcoavlcsuulbastepatcheeomfdaitmrixeneslieomnednKtswoifthHˆin′
{| i| i| i}
include an additional energy truncation step into the re- within it (no truncation necessary):
cursion procedure, to ensure that each t remains free
n
from“high-energy”components, i.e. Hˆ′-|eigienstateswith ˜i =(Hˆ′)i−1 ψ , i=1...dK, (36a)
| i | i
eigenenergiesE′ >1,whichfalloutsidetherange[ 1,1] ˜i i orthonormalize via Gram-Schmidt,
thatisadmissabkleforargumentsofChebyshevpolyn−omi- | i7→| i (36b)
als. If W < W, numerical noise causes the state t to
∗ | ni (Hˆ′ ) = iHˆ′ j , Hˆ′K CdK×dK. (36c)
containsuchhigh-energycontributionsinspiteofthepre- K ij h | K| i ∈
cautionarymeasuresdescribedafterEq.(20),becausethe
application of Hˆ′ to tn−1 in Eq. (30) entails a DMRG Next,fully diagonalizeHˆK′ toobtainalleigenenergiesε′α
| i and eigenvectors e
truncation step, which is not performed in the eigenba- α
| i
sis of Hˆ′. If such high-energy components were fed into
subsequentrecursionsteps,thenorms t ofsuccessive dK
Chebyshev vectors would diverge rapikdlnyik(as would the Uˆ†HˆK′ Uˆ = |eαiε′αheα|. (37)
resulting moments µn), because this effectively amounts αX=1
to evaluating Chebyshev polynomials Tn(x) for x > 1, Then construct the projection operator
| |
where T 1 (see Fig. 2(b)).
n
| |≫
Asaconsequence,afterobtaininganewstate t from
n P = e e (38)
| i α α
Eq.(30),wetaketheprecautionarymeasureofprojecting 1− | ih |
out any high-energy components that it might contain, α:εX′α≥εP
before proceeding to the next |tn+1i. This can be done for a certain energy threshold εP and apply it:
by performing several energy truncation sweeps. During
an energy truncation sweep, we focus on one site at a ψ P ψ . (39)
time, performanenergytruncationina localKrylovba- | i7→ | i
sisconstructedforthatsite,andthenmoveontothenext Performingthisprocedureonceforeverysiteofthechain
site. Shifting the current site is accomplished by stan- constitutes a truncation sweep. The state obtained after
dard MPS means, without any truncation, as a DMRG severaltruncationsweeps,say t ,isstrippedfromthe
n tr
| i
truncation would counteract the energy truncation. (As unwantedhigh-energycomponentsof t ,aswellaspos-
n
| i
a consequence, an energy truncation in terms of two-site sible within a Krylov approximation. After fitting and
sweeps has not been implemented.) truncation have been completed, the resulting (unnor-
The truncation must take place in the energy eigen- malized)state t isrenamed t , usedforcalculating
n tr n
basis of the Hamiltonian Hˆ′. Of course, its complete µ , and fed int|o tihe next recursi|onistep.
n
eigenbasis is not accessible, thus we build a Krylov sub- To quantify the effects of energy truncation, we con-
space of dimension dK within the effective Hilbert space sidertwomeasuresofhowmuch tn changesduringtrun-
ateverysite. Alternatively,energytruncationcanalsobe cation. First, foragiventruncat|ioni sweep,we define the
performed in the bond representation ψ = Blr lk rk . average truncatedweight per site (averagedoverallsites)
| i | i| i
In this Krylov subspace, the effective Hamiltonian Hˆ′ by
K
of dimensiond can be fully diagonalizedandso we can
K
construct a projection operator to project out all eigen- 1
states with energy bigger than some energy truncation Ntsrweep =vL |hekα|ψi|2, (40)
threshold εP. Thechoiceofthisthresholddependsonthe u Xk α:εX′α≥εP
u
t
8
parameter recommended value description task
W∗ 2WA (or W) effectivebandwidth with (or without) energy truncation rescaling of H
ǫt 0.025 safety offset in rescaled half-bandwidth: W′=1− 21ǫt
m MPS dimension
recurrencefitting
∆c 10−6...10−8 fitting convergencethreshold
d 30 Krylov subspace dimension
K
n 10 numberof sweeps energy truncation
S
ε 1.0 energy truncation threshold (in rescaled units)
P
N dependson system size order of expansion, broadening
spectral reconstruction
g gJ choice of dampingfactors
n n
TABLEI.List ofCheMPSparametersthatcontrolvariousalgorithmic tasks,influencingtheirnumericalcostsandthequality
ofresults. Sincethetasks“recurrencefitting”and“energy truncation”arecarried outat everyrecursion step,theimportance
of the corresponding parameters is self-evident. However, W∗ and ǫt, which determine the rescaled Hamiltonian Hˆ′, turn out
to havea high impact on theresults, too, as thequality of theenergy truncation sweeps strongly dependson Hˆ′. For ε =1,
P
thechoiceof takingW∗ tobetwicethespectral bandwidthWA (orequaltothemany-bodybandwidthW)was foundtowork
well with (or without) energy truncation, respectively. N and g do not affect the calculated moments of the expansion, but
n
control the broadening of thereconstructed spectral function.
where ek are the vectors constituting the projector of IV. RESULTS: HEISENBERG
| αi
Eq. (38) at site k. Second, we define the truncation- ANTIFERROMAGNET
induced state change by
To illustrate the capabilities and power of the pro-
posed CheMPS approach, this section presents results
∆ = t t 2 . (41) for the spin structure factor of a one-dimensional spin-
tr n tr n
k| i −| ik 1 Heisenberg antiferromagnet (HAFM) and compares
2
them againstresults obtained from CV and tDMRG ap-
proaches.
Itmeasureschangesinthestateduetotheintendedtrun-
cationofhighenergyweight,butalsoduetounavoidable
numerical errors. In our experience, neither of the trun-
cationmeasuresNsweep and∆ showclearsignsofdecay A. Spin structure factor
tr tr
whenincreasingthenumberoftruncationsweeps,sayn
S
(see Fig. 10(c) below). This reflects the fact that energy We study the spin-1 HAFM for a lattice of length L
2
trunctation has the status of a precautionary measure,
not a variational procedure, and implies that there is no L−1
dynamic criterionwhen to stop truncation sweeping. As HˆHAFM =J Sˆj Sˆj+1, (42)
·
a consequence, one has to analyze how the accuracy of j=1
X
theresultsdependsonn andoptimizethelatteraccord-
S
ingly. This will be described in Sec. VB below. where Sˆ denotes the spin operator at site j. We
j
choose J = 1 as unit of energy throughout this section.
The numerical costs for energy truncation are as
This model exhibits SU(2) symmetry, which has been
follows: The cost for the steps in Eqs. (36) are
exploited33 in our calculations , accordingly all MPS di-
(d2 m3d2D ), where m is the MPS dimension, d the
O K H mensions noted for the HAFM are to be understood as
size of the local site basis and D the matrix product
H number of SU(2) (representative) states being kept. To
operator dimension of Hˆ′. The diagonalization of Hˆ′K
accountfor the open boundary conditions we define spin
is of (d3 ) where d is theoretically bounded by m2d.
O K K wave operators as:
In our experience, the purpose of the energy truncation,
which is solely to eliminate high-energy contributions, is
L
well accomplished already for a relatively small Krylov Sˆ = 2 sin jkπ Sˆ . (43)
subspace dimension of d =30 m2d. k L+1 L+1 j
K ≪ r j=1
X
AnoverviewofalltheparametersrelevantforCheMPS
is given in Table I. Where applicable, it also lists the The spin structure factor (spectral function) we are in-
values that we found to be optimal. A detailed error terested in is given by
analysis, tracing the effects of various choices for these
parameters, will be presented in Sec. V. S(k,ω)= S†k·Sk(ω). (44)
A
9
It is known from exact solutions34–38 that the dominant
part of the spin structure factor stems from two-spinon ChL: N=111
2
contributions, bounded from below and above by ChL: N=222
CV: η=0.1
π k CV: η=0.05
ω = sink and ω =π sin . (45) 1.5
1 2
2 | | (cid:12) 2(cid:12) ω) L=100, λ=4.0, a=3.19
(cid:12) (cid:12) 2,
Mdivoerregoevears,for an infinite system S(k,ω)(cid:12)(cid:12) is kn(cid:12)(cid:12)own37,38 to πS(k=/ 1 mCWV=∗1=060.30,, εmt=C0h.=02352
S(k,ω) [ω ω1]−21 ln[1/(ω ω1)],for k =π, (46a) 0.5
∼ − − 6
S(π,ω) ω−1 ln(1/pω), (46b)
∼
asω approachespthelowerthresholdω fromabove. This 0
1
divergence reflects the tendency towards staggered spin 0.4 0.5 0.6 0.7
order of the ground state of the Heisenberg antiferro- ω/π
magnet. It poses a severe challenge for numerics, which
always deals with systems of finite size, and hence will Fig. 4. (Color online) Comparison of CheMPS vs. correc-
never yield a true divergence. Instead, the divergence tion vector (CV) calculations of S(k = π/2,ω) for a Heisen-
will be cut off at ω ω 1/L, yielding a peak of finite berg chain: Lines show Chebyshev results reconstructed for
1
− ≃ N = 118 and 236 using Lorentz damping, Eq. (16), with
height,
λ = 4.0; symbols show CV results, obtained using broad-
1 ening parameters of η = 0.1 and 0.05. We expect and in-
maxS(k,ω) [LlnL]2 , for k =π, (47a)
∼ 6 deed find good agreement between lines and symbols, since
1
maxS(π,ω) L[lnL]2 . (47b) both Lorentz damping and the correction vector method in
∼
effect broaden the spectral function by Lorentzians, whose
Thus, the best one can hope to achieve with numerics is widths we equated by choosing aλ/N = η (with a = 3.19),
to capture the nature of the divergence as ω approaches see Eq. (48), and also Fig. 3. Since mCh ≪ mCV, the nu-
ω before it is cut off by finite size, or the scaling of the merical cost of obtaining an entire curve via Chebychev is
1
peak height with system size. dramatically cheaper than calculating a single point via CV,
as discussed in thetext.
Eq.(45)givesagoodguideforchoosingW . Wefound
∗
the choices W = 6.3 2π and ǫ = 0.025 to work well
∗ t
≃
for all k and have used them for all figures (4 to 6) of
to be very accurate,thoughalsocomputationally expen-
this section. As consistency checks, we verified that the
sive. The CV method involves a broadening parame-
resulting S(k,ω) is essentially independent of W , and
∗ ter η and broadens δ-functions into Lorentzian peaks of
that it agrees with a calculation that included the full
width η. This can be mimicked with CheMPS by us-
many-body bandwidth (W =W).
∗ ing Lorentz damping (with λ=4.0), since this also pro-
To haveanaccuratestartingpoint for allcalculations,
duces Lorentzian broadening, representing a δ-function
wethroughoutusedagroundstateobtainedbystandard δ(ω′ ω¯′) by a near-Lorentzian peak, albeit with a
DMRG with MPS dimension m = 512. From expansion −
frequency-dependent width, η′ = √1 ω¯′2λ/N (see
order n = 1 onwards, it turned out to be sufficient to N,ω¯′ −
Sec. IIA and Fig. 3). To compare CheMPS results with
represent all Chebyshev vectors t using a surprisingly
| ni CVresultsatgivenη,wethusidentifyη =aη′ ,where
small MPS dimension of m = 32, or m = 64 for some N,ω¯′
a is the scaling factor from Eqs. (21) and ω¯′ is taken to
resultsinvolvingverylargeiterationnumber,asindicated
betherescaledandshiftedversionofthefrequencyω
in every figure. (In retrospect, this implies that for the max
atwhichthepeakreachesitsmaximum. Thuswesetthe
groundstate,too,amuchsmallermwouldhavesufficed.)
expansion order used for reconstruction to
We have verified that the structure factor S(k,ω) is well
converged w.r.t. m nevertheless. Detailed evidence for
4a
this claim will be presented below. However, already N = 1 (ω /a W′)2. (48)
max
η − −
at this stage it is worth remarking that the ability of
p
CheMPS to get good results with comparatively small m-
Figure 4 shows such a comparison for the structure fac-
values is perhaps the single most striking conclusion of
tor S(π/2,ω) of a L = 100 Heisenberg chain. We
our work. This will be discussed in detail below.
used two choices of η that are large enough to avoid
finite-size effects, namely η = 0.1 and 0.05, and set
ω = π/2 (cf. ω of Eq. (45)). We used MPS dimen-
max 1
B. Comparison to CV sions of m = 1000 or m = 32 for CV or CheMPS
CV Ch
calculations,respectively. (Ourchoiceform aimedfor
CV
WebeginourdiscussionofCheMPSresultsbycompar- achieving highly accurate CV results; for η = 0.05 this
ing them to those of CV calculations, which are known required m =1000, but for η =0.1, a slightly smaller
CV
10
value for m would have sufficed.) We find excellent fact is important, since it implies that the structure fac-
CV
agreement between the two approaches without adjust- tor of a finite-size system is exhausted almost fully by
ing any free parameter,since N is fixedby Eq.(48). For the set of dominant subpeaks, with very small intrinsic
example, for η =0.05, N =255, the relative error is less widths.
than 3 % for all ω.
Since this level of agreementis obtained using m We have checked that there are (L) dominant sub-
m , we conclude that CheMPS with Lorentz damCphin≪g peaks within the spectral bandwidtOh of S(k,ω). Corre-
CV
gives results whose accuracy is comparable to those of spondingly, the average spacing between dominant sub-
CV, at dramatically reduced numerical cost. Indeed, for peaks,to be calledthe finite-size energy scale ωL,ispro-
η = 0.05 the calculation of the entire CheMPS spectral portionalto 1 (Fig. 5(c,d)). The weightofeachsubpeak
L
function was 25 times faster than that of a single CV decreases similarly, ensuring that the total weight in a
data point. given frequency interval converges as L . The in-
verse subpeak spacing ~/ω corresponds→to∞the Heisen-
L
berg time, i.e. the time within which a spin wave packet
propagates the length of the system.
C. Finite-size effects
Figures 5(e,f) illustrate two slightly different broaden-
Letusnowanalysetheroleoffinitesystemsize. Tothis
ing strategies. In Fig. 5(e), L is increased for fixed N:
end, it is of course important to understand broadening
the distinct subpeaks increasingly overlap,resulting in a
effects in detail. The fact that CheMPS offers simple
smooth spectral function once ω drops below η . In
L N
and systematic control of broadening via the choice of
Fig.5(f)optimalbroadeningis used(η justlargerthan
N
the expansion order (and damping factors), as will be
ω : now no subpeaks are visible, and the L-evolution of
L
illustrated below, is very convenient and may regarded
the main peak is revealed with better resolution.
as one of its main advantages.
Figure 5(a) shows CheMPS results for the spin struc- In both Figs. 5(e) and 5(f), the peak height shows no
ture factors S(k,ω) of four different momenta k, calcu- indications of converging with increasing L. (The same
lated for L = 100 using Jackson damping. They were is true for the data ofFig. 5(a).) This reflects the intrin-
reconstructedusing the largestexpansionorder,say NL, sicdivergenceofthepeakheightexpectedfromEq.(46).
that does not yet resolve finite-size effects, a choice that Figures 6(a) and 6(b) contain a quantitative analysis of
will be called optimal broadening. Each curve shows a this divergence, for S(π,ω) and S(π/2,ω), respectively.
dominant peak, and we are interested in finding its in- The shape of the divergencies for an infinite system are
trinsic shape S∞(k,ω) in the continuum limit of an in- shown by the thick solid lines, representing exact Bethe
finitely longchain(L ). Thus, the followinggeneral Ansatz results from Ref. 38. Thin dashed lines show re-
→∞
question arises: under what conditions will a spectrum sults from tDMRG from Ref. 22 for L = 100, and thin
calculatedforfinitesystemsizeLandreconstructedwith solid lines CheMPS results for several system sizes be-
finiteexpansionorderN,saySL(ω),correctlyreproduce tween L = 50 and 300. For CheMPS spectral recon-
N
thedesiredcontinuumspectrumS∞(ω)? Thegeneralan- struction, we determined the expansion order N that
300
swer,ofcourse,isthattheoptimallybroadenedspectrum ensures optimal broadening for L = 300, and used a
should have converged as function of L, i.e. the shape fixed ratio of N/L = N /300 for all curves (namely
300
of SL (ω) should not change upon increasing L. How- 0.42 or 0.67 for Figs. 4(a,b), respectively). CheMPS (for
NL
ever,foraspectrumwithanintrinsicdivergence,suchas L = 300) and tDMRG reproduce the peak’s tail and
Eq.(46), the peak’sheightwillneversaturatewithL; at flank well, but clearly and expectedly are unable to pro-
bestone canhope to observeL-convergenceofthe shape duce a true divergence at the lower threshold frequency.
ofits tail,andthe properscalingofits height(Eq.(47)). Nevertheless, the insets show that the manner in which
To illustrate the nature of finite-size effects and the the CheMPSpeak heights increase with L is indeed con-
role of N in revealing or hiding them, Fig. 5(b) shows sistent with Eq. (46). (For the limited range of avail-
S(π,ω)forL=100andseveralvaluesofN,bothsmaller ablesystemsizes,however,areliabledistinctionbetween
and larger than NL. As N is increased and the effective L[ln(L)]1/2, [Lln(L)]1/2 or L behavior is not possible.)
broadening η (W /N) decreases, the main peak
N ∗
≃ O
of the initially very broad and smooth spectral function It is also possible to determine the lower threshold
becomes sharper. Optimal broadening in Fig. 5(b) cor- frequency ω rather accurately from the CheMPS re-
1
respondstoN 70,beyondwhichadditional“wiggles” sults by doing an 1/L extrapolation. We illustrate this
L
≃
emerge. These develop, with beautifully uniform resolu- in Fig. 6(b) by extrapolating the frequencies at which
tion, into dominant subpeaks as N is increased further. S(π/2,ω) = 0.1 (triangles). Since the data exhibit a
The discrete subpeaks reflect the quantized energies of slight curvature when plotted against 1/L (see lower in-
spin-waveexcitationsinafinitesystem. Withsufficiently set of Fig. 6(b)), they were fitted using a second order
high resolution (N = 999 in Fig. 5(b)) numerous addi- polynomial in 1/L. Extrapolating the fit to 1/L = 0
tional minor subpeaks emerge, but their weight is very yields ω =0.496π (marked by a square), in good agree-
1
small compared to that of the dominant subpeaks. This ment with the prediction ω =π/2 from Eq. (45).
1