Table Of ContentGeneralized Moment Method for Gap Estimation
and Quantum Monte Carlo Level Spectroscopy
Hidemaro Suwa1,2 and Synge Todo2,3
1Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA
2Department of Physics, University of Tokyo, Tokyo 113-0033, Japan
3Institute for Solid State Physics, University of Tokyo, Kashiwa 277-8581, Japan
5 (Dated: June2, 2015)
1
We formulate a convergent sequence for the energy gap estimation in the worldline quantum
0
Monte Carlo method. The ambiguity left in the conventional gap calculation for quantum systems
2
is eliminated. Our estimation will be unbiased in the low-temperature limit and also the error
y bar is reliably estimated. The level spectroscopy from quantum Monte Carlo data is developed as
a an application of the unbiased gap estimation. From the spectral analysis, we precisely determine
M
the Kosterlitz-Thouless quantum phase-transition point of the spin-Peierls model. It is established
that thequantumphonon with a finite frequencyis essential tothe critical theory governed by the
0 antiadiabatic limit, i.e., the k=1 SU(2) Wess-Zumino-Witten model.
3
]
The excitation gap is one of the most fundamental necessarilyintroducessomebiasinthefitting procedure.
h
c physical quantities in quantum systems. The Haldane Our purpose in the present paper is to establish a versa-
e phase and the Z topological phase are characterized by tile and unbiased gap-estimation procedure free from an
2
m
the topologically protected gap[1]. Recently the exis- ambiguous fitting.
- tence of gapful/gapless quantum spin-liquid phases has In the meanwhile, the recently advanced technology
t
a been discussed in frustrated spin systems[2]. Not only hasallowedforquantumsimulatorsthatcanrealizeideal
t in the gapfulbut also criticalphases,the system-size de- quantummany-bodysystems[8]. Inparticular,thequan-
s
. pendence of the excitation gap is useful for the analysis tum phonon effect of trapped ions has caught a great
t
a ofthe quantumphase transition[3]. Particularly,the en- deal of attention, which provides rich physics and engi-
m ergy gap ∆ in the conformal quantum phases scales as neering,e.g.,spinfrustration[9],long-rangespininterac-
- ∆ xv/L, apart from possible logarithmic correction, tion[10],phononsuperfluids[11],quantumgates[12],etc.
d ∝
where L is the system size, x is the scaling dimension, As anapplication of the presentgap-estimationmethod,
n
and v is the velocity appearing in the conformal field wewillelucidatethequantumphasetransitionofthespin
o
theory[4]. An unbiased gap calculation thus allows for system coupled with quantum phonons, which is called
c
[ extracting the universal properties of the critical phases the spin-Peierlstransition[13–19]andis accessibleinthe
from finite-size data. ion system[20]. The level spectroscopy[3] from Monte
2
Carlo data is developed to overcome the difficulty of the
v
Thegapestimationforlargesystemsisnottrivial. For [Kosterlitz-Thouless(KT)]transitionthatmakesthecon-
7
4 small systems, it is possible to calculate the gap by the ventional approaches ineffective. We establish that the
8 exactdiagonalizationmethod. Thereachablesystemsize quantumphononeffectisessentialtothespin-Peierlssys-
0 is, however, strongly limited because of the explosion of tem and its critical phenomena.
2. required memory size and computation time. The den- The spectral information of a quantum system de-
0 sity matrix renormalization group (DMRG) method[5] scribed by a Hamiltonian H is encoded in the imaginary
4 works well for many one-dimensional systems, but it be- time (dynamical) correlation function:
1 comes less effective in gapless or degenerated phases.
v: In the meanwhile, the quantum Monte Carlo (QMC) C(τ) = Oˆ(τ)Oˆ† = 1tr eτHOˆe−τHOˆ†e−βH
methodbasedontheworldlinerepresentationisapower- h i Z
i
X ful method for various strongly correlated systems with- 1 h i
r out dimensional restriction[6]. In previous QMC cal- = Z bℓ,ℓ′e−τ(Eℓ−Eℓ′)e−βEℓ′
a ℓ,ℓ′
culations[7], the gap is extracted by the fitting of the X
correlation function; the tail of the exponential function b e−τ(Eℓ−E0) (β ), (1)
ℓ
→ →∞
is estimated as a fitting parameter [see Eq.(1) below]. ℓ≥1
X
Hereweencounteratrade-offbetweenthesystematicer-
ror and the statistical error. The lower the temperature whereOˆ isa chosenoperator,Z is the partitionfunction
is in a QMC simulation, the smaller the systematic er- and β = 1/T is the inverse temperature. Also, ℓ
{| i}
rorbecomesbutthelargerthestatisticalerrordoes. Itis is the complete orthogonal set of eigenstates, Eℓ is the
becausethecorrelationinlongimaginarytimehasanex- associated eigenenergy, bℓ,ℓ′ = ℓ′ Oˆ ℓ 2, and bℓ = bℓ,0,
|h | | i|
ponentially smallabsolutevalue. Since in practicewe do where the ground state is 0 , the first excited state is
| i
not have a prior knowledge for the optimal temperature 1 , and so are the higher excited states, respectively.
| i
and the range of imaginary time where the correlation We assume Oˆ 0 =0, 0Oˆ 0 = 0, and E >E (ℓ 1),
ℓ 0
| i6 h | | i ≥
function follows the asymptotic form, a choice of data i.e., the groundstate is certainlyexcited by the operator
2
and the gapis finite. The latter is the case for finite-size O(β−(2n+2)∆−(2n+3)) with coefficients x . It will be
ℓ,ℓ′ n,k
systems even if the ground state is degenerated in the dominated by the smallest gap ∆ in β and
1
→ ∞
thermodynamic limit. n . The coefficients x satisfy the following equa-
n,k
Here,letusconsiderthemomentoftheimaginarytime tio→ns:∞ n x k2m = δ (0 m n), where δ
k=0 n,k m,n ≤ ≤ mn
correlation function[21]: is the Kronecker delta. They are exactly solved for by
P
the formula of the inverse of Vandermonde’s matrix[24].
Ik = k1!Z0∞τkC(τ)dτ =Xℓ≥1 ∆bkℓℓ+1 (k ≥0), (2) Weℓ,ℓ′cgaℓ,nℓ′ω1a2(lnso−1)∆shℓ−o,ℓ(w′2n−(1−)/1Z)n+−1OP(βnk−=20nk∆2−ℓx,nℓ(′2,knC+˜1()ω)kb)y u=s-
ing the same x . Then, a sequence of the higher-order
n,k
where ∆ℓ ≡∆ℓ,0, and ∆ℓ,ℓ′ =Eℓ−Eℓ′. The higher mo- Pestimators is derived as
ment will be dominated by the contribution from the
fi∆rst<exc∆itati<on gap. aTshIekn ∼web1c/a∆nk1s+e1e baecuasuesfuel∆li1mi<t, ∆ˆ =ω n k2x C˜(ω ) n x C˜(ω ) (6)
(Ik2/Ik+m)31/m →···∆1 (k → ∞) ∀m ∈ N. However, (n,β) 1vuu−Xk=0 n,k k ,Xk=0 n,k k
we cannotuse the momentdirectly in finite-temperature t
with x = 1/ n (k + j)(k j). Importantly,
simulations. It is because the correlation function is pe- n,k j=0,j6=k −
riodic for bosons or anti-periodic for fermions and the ∆ˆ I /I (β ), and the systematic
(n,β) → 2(n−Q1) 2n → ∞
moment is not well defined. Then the Fourier series can error is expressed as
p
be exploited instead. Let us think of the bosonic case
because we will investigate spin excitation. The Fourier ∆ˆ(n,β) 1 bℓ ∆1 2n−1 ∆1 2n
1+ +O .(7)
cfroemqupeonnceyntωof=t2hπejc/oβrr(ejlatiZon) ifsunexctpiroenssaedt aasMatsubara ∆1 → 2Xℓ>1"b1 (cid:18)∆ℓ(cid:19) (cid:18)∆ℓ(cid:19) !#
j
∈ Note that I can be achieved only for even-numbered
k
β k since C˜(ω ) is real. As examples, (x ,x ,x ) =
C˜(ωj)= C(τ)eiτωjdτ (3) (1, 1, 1 ) fjor n = 2[23], and (x ,x2,0,x2,1,x2,2) =
Z1 Z0∆gℓ2,ℓ′∆+ℓ,ωℓ′2 (ωj 6=0) (t−e4n3−1d6,o3w214n1,2−th6e10b,i3a16s0)offotrheng=ap3.esWtimeahta3o,v0re(a63n),a1alnydt3i,cs2ahlolyw3,n3wrthite-
= Xℓ,ℓ′ ℓ,ℓ′ j following remarkable property:
Z1{EℓX6=Eℓ′ ∆gℓℓ,,ℓℓ′′ +EℓX=Eℓ′bℓ,ℓ′e−βEℓβ} (ωj =0), nl→im∞βl→im∞∆ˆ(n,β) =βl→im∞nl→im∞∆ˆ(n,β) =∆1. (8)
wheregℓ,ℓ′ =bℓ,ℓ′(e−βEℓ′ e−βEℓ). Inmanysimulations, Thatis,thesetwolimitsareinterchangeable(seetheSup-
the so-called second mom−ent[22] is used as the lowest- plementalMaterial[25]fordetails). Thisimportantprop-
erty makes our gapestimation greatlyrobust. Note that
order gap estimator:
the present approach works also in the stochastic series
C˜(ω ) I expansion QMC method by the time generation[26].
∆ˆ =ω 1 0 (β ). (4) Our generalized gap estimator is applicable to any
(1,β) 1sC˜(ω0) C˜(ω1) → rI2 →∞ quantum system. As an example, we will show an appli-
−
cation with the level spectroscopy to the following one-
Interestingly,thisestimatorwillbetheratioofthezeroth
dimensional S =1/2 spin-Peierls model:
andthesecondmomentinthelow-temperaturelimit. We
havetotakenoticeofthesystematicerrorcarefully. The
ωλ
error remains even in β →∞ as H = r (1+r 2 (ar+a†r))Sr+1·Sr+ r ωa†rar, (9)
∆ˆ 1 b ∆ ∆ 2 X X
(1,β) 1+ ℓ 1 +O 1 , (5) where ω is a dispersionless phonon frequency, λ is a
∆1 → 2 ℓ>1"b1 ∆ℓ (cid:18)∆ℓ(cid:19) !# spin-phonon coupling constant, Sr is the spin-21 opera-
X tor, a anda† arethe annihilationandcreationoperator
r r
which is typically a few percent of ∆1[23]. This cor- of the soft-core bosons (phonons) at site r, respectively.
rection hampers proper identification of the universality This spin-Peierlsmodelhas beeninvestigatedinthe adi-
class in the level spectroscopy analysis as we will see be- abatic limit (ω 0)[13], from the antiadiabatic limit
low. (ω )[14–18→, 27, 28], and in its crossover[19]. The
Our main idea in the present paper is to construct rele→van∞ce of the present model to real materials, such as
a sequence of gap estimators that converges to a ra- CuGeO [29], has been discussed[30]. The model is ex-
3
tio of higher-order moments in the low-temperature pected to exhibit a KT-type quantum phase transition
limit. We will consider an estimator that has a cor- between the Tomonaga-Luttinger (TL) liquid phase and
rection of O((∆1/∆ℓ)2n−1). Let us expand Eq.(3) in the dimerphaseatafinite spin-phononcoupling[15,19],
powers of (1/β∆ℓ,ℓ′) and make a linear combination of which is absent when either spin or phonon is classically
Fourier components so that the lowest orders cancel; treated[13]. Therealizationofthequantumphasetransi-
(−1)n nk=0xn,kC˜(ωk) = ℓ,ℓ′gℓ,ℓ′ω12n∆−ℓ,ℓ(′2n+1)/Z + tionwasrecentlyproposedinthetrappedionsystem[20].
P P
3
of the diagonal operators, our simulation is free from an
0.02 occupation-number cutoff of the soft-core bosons. The
Fourier components of the correlation function(3) are
0.01
directly calculated during the simulation. The worm-
scattering probability is optimized in rejection (bounce)
0
×10-2 1/n ratebybreakingthedetailedbalance[37]. Theboundary
n=1
ε -0.01 2 n=2 condition was periodic in the space and time directions.
0 n=3 More than 225( 3.4 107) Monte Carlo samples were
-0.02 -2 inflecnti=o∞n takenintotalaf≃ter218×( 2.6 105)thermalizationsteps.
≃ ×
-4 optimal fit The error bar of the gap estimates is calculated by the
-0.03 0 0.5 1
(τ - τ )/τ jackknife analysis[38].
2 1 2
First, the convergence of our gap estimate was tested
0 0.05 0.1 0.15 0.2 0.25 0.3
for L = 4, ω = 4, λ = 1/2, where L is the system
T
size. We set ω here is fairly larger than the actual spin
gap because this condition is satisfied for large systems
FIG. 1. Triplet-gap estimation (relative) error of our es-
in the relevant spin-phonon coupling region. The bo-
timators(6) (n = 1,2,3,∞), the inflection-point value of
son occupation number cutoff D was set to 4 only in
−dlogC(τ)/dτ,andtheoptimalfit(definedinthemaintext)
this test for comparing with the diagonalization result.
for the spin-Peierls model(9) with L = 4, ω = 4, λ = 1/2,
D = 4, which has ∆ ≃ 1.111388. The inset shows the n Fig.1 shows the calculated triplet-gap estimation errors,
1
dependence of the gap estimate (circles) for 1≤ n ≤ 10 and where Oˆ = Szeiπr is used in the dynamical corre-
r r
theτ dependenceofthelogC(τ)linear-fitresult(diamonds) lation function. We compared the gap estimators(6) to
1
for τ1 ≤ τ ≤ τ2 with τ2 = β/4, calculated from 220(∼ 106) thepreviousPapproach[7]wherethefirstgapisestimated
Monte Carlo steps at β=6. as dlogC(τ)/dτ from the asymptotic form(1). The
−
derivative will show a plateau at the gap value in an ap-
propriate τ region. When β(=1/T) is not large enough,
however, the plateau is indistinct. Then the inflection
Itisdifficultto preciselylocatethe transitionpointby point could be used, but it is hard to estimate in prac-
the conventional analyses. The huge Hilbert space with tice (here we calculatedit by longerQMC simulationfor
the soft-core bosons hinders the sufficient-size calcula- comparison). Asananotherpracticalandreasonablegap
tion by the diagonalization method[16]. The effective estimation,wetestalinearfitoflogC(τ)forτ τ τ ,
1 2
≤ ≤
spin-model approachby the perturbation[14] orthe uni- wherewefixτ =β/4. TheinsetofFig.1showsthefeasi-
2
tary transformation[17] does not take into account all bleconvergenceofthegapestimateinnandthedifficulty
marginal terms, e.g., the 4-spin and 6-spin interactions of finding appropriate τ for the linear fit. The function
1
examinedinRef.[31]. AsfortheDMRGmethod,itneeds logC(τ)ispoorlyfittedtoalinearformatsmallτ ,while
1
an additional symmetry breaking term, which blurs the ithaslargerstatisticalerroratlargeτ . Thenthegaper-
1
phase transition point, in the degenerated phase[27]. In rorresulting fromthe linear regressiontakes a minimum
addition, the method has difficulty in precise calculation value at optimal τ∗, which we call “optimal fit.” Even
1
of the relevant quantities, such as the central charge, though it seems reasonable, the optimal fit underesti-
around an essential singularity[32]. Also, the previous mates the gapat β =6,8and overestimatesit at β =12
QMC approach[28] suffers from the exponential diver- asshowninthemainpanelofthefigure. Meanwhile,the
gence of the correlation length and the logarithmic cor- second-moment estimator(n = 1) has a non-negligible
rection around the KT transition point. These difficul- bias even in T 0 as expected from Eq.(5). The esti-
→
ties mentioned above can be overcome by employing the matewithlargeenoughn(we callitthe n= estimate
∞
levelspectroscopymethod[3,33]combinedwithourpre- hereafter),ontheotherhand,exponentiallyconvergesto
cise gap estimation. In the TL liquid phase, the both of the exact value as the temperature decreases[25]. The
triplet and singlet excitation are gapless in the thermo- bias convergence is much faster than that of the inflec-
dynamic limit, but the lowest excited state of finite-size tion point (one of the best estimates from the fitting ap-
systems is the triplet because of the logarithmic correc- proach). Moreover, the higher-order estimator provides
tion[34]. Inthedimerphase,ontheotherhand,the first a reliable error bar, while the optimal fit significantly
excitedstateisthesingletforfinite-sizesystems. Itforms underestimates it[25]. Therefore, our approach is more
the degenerated ground states eventually in the thermo- preciseandstraightforwardthanthefittingapproach. In
dynamiclimit. Thustheexcitationgapsofthetripletand the present study, we have used a simple recipe to opti-
singletexcitationintersectat a spin-phononcoupling for mize n and β, minimizing both the systematic and the
finite-sizesystems. Thetransitionpointcanbeefficiently statistical error[25].
extrapolated from the gap-crossingpoints[3]. The scaling of the gap-crossing point for the spin-
We used the continuous-time worldline representa- Peierls model between the triplet and singlet excitation
tion[26, 35] and the worm (directed-loop) algorithm[36] isshowninFig.2. Forthesingletexcitationgap,weused
in the QMC method. Thanks to the exponential form Oˆ = S S eiπr. The bare excitation phonon gap
r r · r+1
P
4
0.55
0.3
S=1, n=1
0.28
0.5
0.26 S=1, n=∞
λ(L)c 0.24 ∆ x(L) S=0, n=1
0.22 0.45
0.2 0.1 0.2 0.3
0.18 λc(∞) = 0.2245(17) λ S=0, n=∞
0.4
0 0.0004 0.0008 0 0.0005 0.001
1/L2 1/L2
FIG. 2. Convergence of the gap-crossing point (circles)
between the triplet and the singlet excitation for L = FIG.3. System-sizedependenceofthescalingdimensioncor-
36,40,48,64,togetherwiththecrossingpointofthespinsus- respondingtothetripletorthesingletexcitationatthetran-
ceptibility (diamonds) between χ (L)/L and χ (L/2)/(L/2). sitionpoint(λ=0.2245),calculatedfromthesecond-moment
s s
The spin-phonon coupling dependence of the gaps is shown (n=1) or the n=∞ gap estimate.
in the inset for each L. The dashed line is the fitting curve
with λ (∞) fixed, which results in large χ2/dof ≈ 5.0. The
c
statistical errors are smaller than thesymbol size.
where the effective spin model is the frustrated J -J
1 2
chain[3]. Ourresultstronglyindicatesthatthequantum
phononeffectisrelevanttothespin-Peierlssysteminthe
was set to ω = 1/4 for the comparison with the previ-
sense that it necessarily triggers the universal KT phase
ous result[28]. The transition point λ = 0.2245(17) in
c
transition.
thethermodynamiclimitwasextrapolatedwithoutloga-
rithmic correction, which is much more precise than the In conclusion, we have presented the generalized mo-
previousestimate,0.176<λ <0.23[28]inournotation. ment method for the gap estimation. The advantages of
c
Also the spin susceptibility χ = β Sz(τ)Sz eiπrdτ our method over the previous approaches are as follows:
could be used for finding thes tran0sitiornh proint0(iFig.2). theunbiasedestimation[Eq.(8)],theabsenceofambigu-
Nevertheless, the gap-crossing poRintPprovides the much ousprocedure,thefasterconvergencew.r.t.thetempera-
morereliableextrapolationwith1/L2 correctionfromir- ture,andthereliableerror-barestimation. Weemphasize
relevantfields[3],whilethesusceptibilityislikelytohave thatourapproachisgenerallyapplicabletoanyquantum
some more complicated corrections. system. The QMClevelspectroscopywasdemonstrated,
forthefirsttime,fortheKTtransitioninthespin-Peierls
We have also calculated the velocity, the central
model. This spectral analysis will likely work in various
charge, and the scaling dimensions at the transition
systemsincludingmostconformalphases. Weelucidated
point, fixing λ = 0.2245. The velocity v = 1.485(8)
thatthequantumphononeffectisrelevanttothecritical
was calculated from the scaling form v(L) = ∆ /k =
v + a/L2 + b/L4 + o(1/L4), where ∆ is thek1trip1let theory of the spin-phonon system, which is expected to
k1 be universal in many kinds of one-dimensional systems,
gap at k = 2π/L, a and b are non-universal constants.
1
e.g., (spinless) fermion-phonon systems, by virtue of the
The central charge c = 0.987(13) was obtained from the
well-established transformations. The clarified quantum
finite-sizecorrection[34],E (L)=E πvc/6L+o(1/L).
0 0
− phase transitionandthe criticalitywouldbe directly ob-
Thescalingdimensioncorrespondingtothetripletorsin-
served in the quantum simulator[20].
glet excitation was calculated from the relation x(L) =
L∆ /2πv,where∆ isthelowest(tripletorsinglet)exci- The authors are grateful to Anders W. Sandvik,
π π
tation gapat k =π. As shownin Fig.3, the n= esti- Thomas C. Lang, and Hiroshi Ueda for the valuable dis-
∞
matesconvergedtox =0.502(3)andx =0.499(3) cussion. The most simulations in the presentpaper were
S=1 S=0
without logarithmic correction as expected only at the done by using the facility of the Supercomputer Center,
transition point[3]. Hence we conclude that this transi- Institute for Solid State Physics, University of Tokyo.
tionpointisdescribedbythek =1SU(2)Wess-Zumino- We also used the computational resources at the Insti-
Witten model[39] with c=1 and x=1/2. Onthe other tute for Information Management and Communication,
hand, the second-moment estimates (n = 1) failed to Kyoto University and Computing and Communications
approach1/2 as seen in Fig.3. This identification of the Center, Kyushu University through the HPCI System
criticaltheoryclearlydemonstratestheimportanceofthe Research Project (No.hp140204, hp140162). The sim-
higher-order estimator. The present study non-trivially ulation code has been developed based on the ALPS
clarified that the critical theory at the transition point library[40]. The authors acknowledge the support by
of the spin-Peierls model with a finite phonon frequency KAKENHI (No.23540438, 26400384) from JSPS and
coincides with that in the antiadiabatic limit (ω ) CMSIinSPIREfromMEXT,Japan. HSissupportedby
→ ∞
5
the JSPSPostdoctoralFellowships forResearchAbroad.
[1] M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82, 3045 [29] M. Hase, I. Terasaki, and K. Uchinokura, Phys. Rev.
(2010). Lett. 70, 3651 (1993).
[2] L. Balents, Nature464, 199 (2010). [30] G. S. Uhrig and H. J. Schulz, Phys. Rev. B 54, R9624
[3] K. Nomura, J. Phys. A:Math. Gen. 28, 5451 (1995). (1996).
[4] J. L. Cardy, Scaling and Renormalization in Sta- [31] Y. Tang and A. W. Sandvik, Phys. Rev. Lett. 107,
tistical Physics (Cambridge Universilty Press, 1996); 157201 (2011).
P. Di Francesco, P. Mathieu, and D. Senechal, Confor- [32] P. Chen, Z.-l. Xue, I. P. McCulloch, M.-C. Chung,
mal field theory (Springer, New York,1997). M. Cazalilla, and S.-K. Yip, J. Stat. Mech. , P10007
[5] S. R.White, Phys. Rev.Lett. 69, 2863 (1992). (2013).
[6] N.KawashimaandK.Harada,J.Phys.Soc.Jpn.73,1379 [33] Y.-C. Tzeng, Phys. Rev.B 86, 024403 (2012).
(2004); A.W.Sandvik,AIPConf.Proc.1297,135(2010). [34] I. Affleck, D. Gepner, H. J. Schulz, and T. Ziman, J.
[7] S. Yamamoto, Phys. Rev. Lett. 75, 3348 (1995); Z. Y. Phys. A: Math. Gen. 22, 511 (1989); R. R. P. Singh,
Meng, T. C. Lang, S. Wessel, F. F. Assaad, and A. Mu- M. E. Fisher, and R. Shankar, Phys. Rev. B 39, 2562
ramatsu, Nature464, 847 (2010). (1989); T. Giamarchi and H. J. Schulz, ibid. 39, 4620
[8] R. Blatt and C. F. Roos, Nature Physics 8, 277 (2012); (1989); S.Eggert, ibid. 54, R9612 (1996).
I.M.Georgescu,S.Ashhab, andF.Nori,Rev.Mod.Phys. [35] B. B. Beard and U. J. Wiese, Phys. Rev. Lett. 77, 5130
86, 153 (2014). (1996).
[9] A. Bermudez, J. Almeida, F. Schmidt-Kaler, A. Retzker, [36] N. V. Prokof’ev, B. V. Svistunov, and I. S. Tupitsyn,
and M. B. Plenio, Phys. Rev.Lett. 107, 207209 (2011). Sov. Phys. JETP 87, 310 (1998); O. F. Syljuasen and
[10] J. W. Britton, B. C. Sawyer, A. C. Keith, C.-C. A. W.Sandvik,Phys. Rev.E 66, 046701 (2002).
JosephWang,J.K.Freericks,H.Uys,M.J.Biercuk, and [37] H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603
J. J. Bollinger, Nature484, 489 (2012). (2010).
[11] D. Porras and J. I. Cirac, Phys. Rev. Lett. 93, 263602 [38] B.A.Berg,Markov Chain Monte CarloSimulations and
(2004). Their Statistical Analysis (World Scientific Publishing,
[12] S.-L.Zhu,C.Monroe, andL.-M.Duan,Phys.Rev.Lett. 2004).
97, 050505 (2006). [39] E. Witten,Commun. Math. Phys.92, 455 (1984).
[13] M. C. Cross and D. S. Fisher, Phys. Rev. B 19, 402 [40] B. Bauer et al.,J. Stat Mech. , P05001 (2011).
(1979).
[14] K.KubokiandH.Fukuyama,J.Phys.Soc.Jpn.56,3126
(1987).
[15] L.G.Caron andS.Moukouri,Phys.Rev.Lett.76,4050
(1996).
[16] G. Wellein, H. Fehske, and A. P. Kampf, Phys. Rev.
Lett. 81, 3956 (1998).
[17] A. Weiße, G. Wellein, and H. Fehske, Phys. Rev. B 60,
6566 (1999).
[18] A.Weiße,G.Hager,A.R.Bishop, andH.Fehske,Phys.
Rev.B 74, 214426 (2006).
[19] R. Citro, E. Orignac, and T. Giamarchi, Phys. Rev. B
72, 024434 (2005).
[20] A. Bermudez and M. B. Plenio, Phys. Rev. Lett. 109,
010501 (2012).
[21] S.Todo, Phys.Rev. B 74, 104415 (2006).
[22] F.Cooper,B.Freedman, andD.Preston,Nucl.Phys.B
210[FS6], 210 (1982).
[23] S.TodoandK.Kato,Phys.Rev.Lett.87,047203(2001).
[24] V.Neagoe, IEEE Signal Proc. Let. 3, 119 (1996).
[25] See Supplemental Material attached below for the
derivation of the systematic error in the case of dis-
crete/continuum spectrum, the error-bar comparison to
thefittingapproach,andtherecipeoftheerroroptimiza-
tion.
[26] A. W. Sandvik, R. R. P. Singh, and D. K. Campbell,
Phys. Rev.B 56, 14510 (1997).
[27] C.J.Pearson,W.Barford, andR.J.Bursill,Phys.Rev.
B 82, 144408 (2010).
[28] A.W.SandvikandD.K.Campbell,Phys.Rev.Lett.83,
195 (1999).
6
Supplemental Material: Generalized Moment Method for Gap Estimation
and Quantum Monte Carlo Level Spectroscopy
I. SYSTEMATIC ERROR OF HIGHER-ORDER GAP ESTIMATE
We providehere the detailed formofthe bias of the gapestimator (6) andits limiting formin n and β .
→∞ →∞
The gap estimator is rewritten as
n
1
v gℓ,ℓ′∆ℓ,ℓ′ ∆2 +ω2
∆ˆ =uu EℓX6=Eℓ′ jY=1 ℓ,ℓ′ j . (S1)
(n,β) u n n
uuuuuEℓX6=Eℓ′ ∆gℓℓ,,ℓℓ′′ jY=1∆2ℓ,ℓ′1+ωj2 +EℓX=Eℓ′bℓ,ℓ′βe−βEℓjY=1ω1j2
t
In the equation, as the reminder of the definitions, ∆ℓ,ℓ′ = Eℓ Eℓ′, bℓ,ℓ′ = ℓ′ Oˆ ℓ 2, gℓ,ℓ′ = bℓ,ℓ′(e−βEℓ′ e−βEℓ),
ω =2πj/β (j Z). Here we used the useful equations: − |h | | i| −
j
∈
n n
x 1
n,k =( 1)nω2n , (S2)
∆2+ω2 − 1 ∆2+ω2
k=0 k k=0 k
X Y
n k2x n 1
n,k =( 1)n−1ω2(n−1) , (S3)
∆2+ω2 − 1 ∆2+ω2
k=0 k k=1 k
X Y
where x = 1/ n (k+j)(k j) as defined in the main text. These equations are derived from the imposed
n,k j=0,j6=k −
condition on xn,k to cancel the lowest orders of (1/β∆ℓ,ℓ′) as explained in the main text:
Q
n n
x k2m =δ (0 m n) x σ ( k)=( 1)nδ (0 m n) (S4)
n,k m,n n,k n,m m,n
≤ ≤ ⇐⇒ − − ≤ ≤
k=0 k=0
X X
n
= k2x σ ( k)=( 1)n−1δ (0 m n 1), (S5)
n,k n,m m,n−1
⇒ − − ≤ ≤ −
k=1
X
whereσ ( k)=σ (0,1,22,32,...,(k 1)2,(k+1)2,...,n2),andσ istheelementarysymmetricpolynomialoforder
n,m m m
− −
m, e.g., σ (a ,a ,a )=1, σ (a ,a ,a )=a +a +a , σ (a ,a ,a )=a a +a a +a a , σ (a ,a ,a )=a a a .
0 1 2 3 1 1 2 3 1 2 3 2 1 2 3 1 2 2 3 3 1 3 1 2 3 1 2 3
The systematic error of the higher-order gap estimate is expressed as
∆ˆ 1+R (β) 1 b ∆ 2n−1 ∆ 2n
(n,β) n ℓ 1 1
= 1+ +O (β ), (S6)
∆1 s1+Fn(β)+Dn(β) → 2 ℓ>1"b1 (cid:18)∆ℓ(cid:19) (cid:18)∆ℓ(cid:19) !# →∞
X
where
gℓ,ℓ′ ∆1 −1
Rn(β)= h(n,β,∆ℓ,ℓ′), (S7)
Eℓ6=Eℓ′,X(ℓ,ℓ′)6=(1,0)(cid:18)g1,0(cid:19)(cid:18)∆ℓ,ℓ′(cid:19)
gℓ,ℓ′ ∆1
Fn(β)= h(n,β,∆ℓ,ℓ′), (S8)
Eℓ6=Eℓ′,X(ℓ,ℓ′)6=(1,0)(cid:18)g1,0(cid:19)(cid:18)∆ℓ,ℓ′(cid:19)
D (β)= bℓ,ℓ′ β∆ e−βEℓ n 1+ ∆21 , (S9)
n g 1 ω2
EℓX=Eℓ′(cid:18) 1,0(cid:19) jY=1 j !
n ∆2+ω2
1 j
h(n,β,∆ℓ,ℓ′)= ∆2 +ω2, (S10)
j=1 ℓ,ℓ′ j
Y
the summations for Fn and Rn are taken over Eℓ = Eℓ′ except (ℓ,ℓ′) = (1,0) and (0,1), Dn term comes from
6
x C˜(0), and b = b . We showed, in the main text, the limiting form: lim ∆ˆ = I /I , where I
n,0 ℓ ℓ,0 β→∞ (n,β) 2(n−1) 2n k
is the moment defined as Eq. (2). Then, lim lim ∆ˆ =∆ since (I /I )1/m ∆ (k ) m N.
n→∞ β→∞ (n,β) 1 k k+m →p 1 →∞ ∀ ∈
7
Next, let us consider the limiting form in n at a finite temperature. We use the product expansion form of
the hyperbolic function: sinh(πz)=πz ∞j=1(1→+z∞2/j2) with z =β∆ℓ,ℓ′/2π in Eq. (S6). The finite-n corrections are
expressed as
Q
n
sinh(πz) 1 1
(1+z2/j2) exp z2 +
∼ πz − n+1 2(n+1)(n+2)
jY=1 (cid:26) (cid:18) (cid:19)(cid:27) (S11)
sinh(πz) 1 1 z4
1 z2 + + .
∼ πz − n+1 2(n+1)(n+2) 2(n+1)2
(cid:20) (cid:18) (cid:19) (cid:21)
Here the asymptotic expansion of the Riemann zeta function was used;
n 1 π2 1 1
= . (S12)
k2 6 − n+1 − 2(n+1)(n+2) −···
k=1
X
The limiting form is then expressed as
∆ˆ 1+R (β) 1
(n,β) ∞
= +O , (S13)
∆1 s1+G∞(β) (cid:18)n(cid:19)
where
R∞(β)= bℓ,ℓ′ ∆1 −2e−β2∆ex, (S14)
Eℓ6=Eℓ′,X(ℓ,ℓ′)6=(1,0)(cid:18) b1 (cid:19)(cid:18)∆ℓ,ℓ′(cid:19)
G∞(β)= bℓ,ℓ′ e−β2∆ex, (S15)
b
(ℓ,ℓ′X)6=(1,0)(cid:18) 1 (cid:19)
and ∆ex =Eℓ+Eℓ′ E1 E0. Note that the correction terms coming from Fn and Dn are included in G∞. Then,
lim lim ∆ˆ − =−∆ because ∆ > 0. As a result, the finite-n corrections are O(1/n). It is attainable to
β→∞ n→∞ (n,β) 1 ex
confirm the convergence as shown in the inset of Fig. 1 in the main text. Furthermore, the n = estimate shows
∞
the exponentialconvergencew.r.t. the temperature,whichwasindeedobservedforthe testcaseinthe maintext(see
Fig. 1).
II. ASYMPTOTIC BEHAVIOR OF GAP ESTIMATE FOR CONTINUUM SPECTRUM
The systematic error (S6) was formulated for a discrete spectrum. This is the case in finite-size quantum systems
with a reasonable basis. We will show a generalizationto a continuum spectrum that will be achievedin the thermo-
dynamic limit. In practice, the crossoverfromthe discrete to continuum casewill be observedin fairly large systems.
First, let us generalize the formulation of the correlation functions. They are expressed as
∞
C(τ)= dǫS(ǫ)e−τǫ, (S16)
Z−∞
β 1 ∞ A(ǫ)
C˜(ω )= dτC(τ)eiτωj = dǫ , (S17)
j
2π ǫ iω
Z0 Z−∞ − j
where A(ǫ) = 2πS(ǫ)(1 ζe−βǫ) is the spectral function with ζ = 1 (bosons or spins) or ζ = 1 (fermions). The
− −
formulation of a discrete spectrum is recoveredby setting
1
S(ǫ)= bℓ,ℓ′e−βEℓ′δ(ǫ (Eℓ Eℓ′)). (S18)
Z − −
ℓ,ℓ′
X
Then the gap estimator (6) is generally expressed as
∞ n 1
dǫA(ǫ)ǫ
v ǫ2+ω2
∆ˆ =uuZ−∞ jY=1 j . (S19)
(n,β) uu ∞ A(ǫ) n 1
u dǫ
uuZ−∞ ǫ j=1ǫ2+ωj2
u Y
t
8
TABLE I. Asymptotic behavior of the gap estimate for typical continuum spectrum. In the spectral function row, θ(x)
is theunit step function; it takes 0 (x<0) or 1 (x≥0).
A+(ǫ) (ǫ−∆1)aθ(ǫ−∆1) e−c/(ǫ−∆1)aθ(ǫ−∆1) ǫa (a>0) e−c/ǫa (a,c>0)
(a>−1) (a,c>0)
∆ˆ(∞,β) ∆1[1+O(1/β)] ∆1[1+O(1/βa+11)] ∼1/β ∼1/βa+11
∆ˆ(n,∞) ∆1[1+O(1/n)] ∆1[1+O(1/na+11)] 0 ∼1/na1
Let us focus on the spin case here. The integrand with O = Sz or SzSz (i and j are site indices) in Eq. (S19) is
i i j
symmetric at ǫ = 0, so considering only for ǫ 0 suffices. Let us then write A(ǫ) = A (ǫ) (ǫ 0). When there
+
≥ ≥
is a delta peak at the first gap (b > 0) and a finite gap between the first gap (∆ ) and the second gap (∆ ), as
1 1 2
∆ ∆ > 0, the systematic error (S6) is expressed in a similar way converting the discrete summation to the
2 1
−
corresponding integral. The asymptotic form becomes
∆ˆ
(∞,β) =1+O e−β2(∆2−∆1) , (S20)
∆
1
(cid:16) (cid:17)
∆ˆ ∆ 2n−1
(n,∞) 1
=1+O . (S21)
∆ ∆
1 (cid:18) 2(cid:19) !
Nextletusthinkofthecasewherethesystemhasacontinuumspectrumabovethelowestgap(∆ ). Theasymptotic
1
behaviors of the gap estimate for some typical spectral function (with ∆ >0 or ∆ =0) are summed up in Table I.
1 1
The parameter constraint for each case comes from the bounded correlation function: C(τ)< .
∞
(i) A (ǫ) (ǫ ∆ )aθ(ǫ ∆ ) (a> 1)
+ 1 1
∼ − − −
When the spectrum is in a power of ǫ ∆ , the systematic errorwill be asymptotically O(1/n) or O(1/β); that is,
1
−
Q(a,2) 1
∆ˆ =∆ 1+O , (S22)
(∞,β) ∼sQ(a,0) 1 β
(cid:20) (cid:18) (cid:19)(cid:21)
T(a,2n 1) 1
∆ˆ − =∆ 1+O , (S23)
(n,∞) ∼sT(a,2n+1) 1 n
(cid:20) (cid:18) (cid:19)(cid:21)
where
Q(a,m)≡Z∆ǫ1M dǫ(ǫsi−nh∆1β)2aǫǫm ∼Z∆ǫ1M dǫ2(ǫ−∆1)aǫme−β2ǫ =2a+2β∆a+m11e−β∆21 Z0qM dqqa(cid:18)1+ β2∆q1(cid:19)me−q
= β∆a+m11e−β∆21 Q(cid:16)0+(cid:17)Qβ1 (β ≫1/∆1), (S24)
(cid:18) (cid:19)
ǫM (ǫ ∆ )a TM ∆a+1−m ta TM ∆a+1−m ta
T(a,m) dǫ − 1 = dt 1 dt 1
≡ ǫm ma+1 (1+ t )m ∼ ma+1 et(1 t2 )
Z∆1 Z0 m Z0 − 2m
∆a+1−m T
= 1 T + 1 (m 1). (S25)
ma+1 0 m ≫
(cid:18) (cid:19)
In Eq. (S24) q =β(ǫ ∆ )/2, in Eq. (S25) ǫ=(1+ t )∆ , and ǫ ,q ,t ,Q ,Q ,T ,T are constants.
− 1 m 1 M M M 0 1 0 1
(ii) A+(ǫ) e−c/(ǫ−∆1)aθ(ǫ ∆1) (a,c>0)
∼ −
Even when the spectral function vanishes exponentially at ǫ = ∆ > 0, the gap estimator will be asymptotically
1
9
unbiased. The convergence becomes slower than the power-law case:
U(a,c,2) 1
∆ˆ =∆ 1+O , (S26)
(∞,β) ∼sU(a,c,0) 1" βa+11!#
V(a,c,2n 1) 1
∆ˆ − =∆ 1+O , (S27)
(n,∞) ∼sV(a,c,2n+1) 1(cid:20) (cid:18)na+11(cid:19)(cid:21)
where
U(a,c,m)≡Z∆ǫ1M dǫe−si(nǫh−∆c1β)2aǫǫm =:Z∆ǫ1M dǫe−um(ǫ) ≈su′m′2(πǫ∗m)e−um(ǫ∗m) ∼∆−1me−mβ−a+11 (β ≫1/∆1), (S28)
(cid:16) (cid:17)
V(a,c,m)≡Z∆ǫ1M dǫe−(ǫ−∆c1)aǫ−m =:Z∆ǫ1M dǫe−vm(ǫ) ≈svm′′2(πǫ∗m)e−vm(ǫ∗m) ∼∆−1mm−2(aa++21)e−ma+a1 (m≫1),
(S29)
u′′(ǫ) and v′′(ǫ) are the second derivative, and the saddle-point approximation around its extremum u′ (ǫ∗ ) =
m m m m
v′ (ǫ∗ )=0 were used, respectively.
m m
(iii) A (ǫ) ǫa (a>0)
+
∼
Let us consider the gapless case with a power-law form. From Eq. (S19), the estimate will diverge since the gap is
actually zero, but it takes a finite value for the case of finite β:
W(a+3)
∆ˆ 1/β (β 1), (S30)
(n,β) ∼sW(a+1) ∼ ≫
where
ǫM ǫm ǫM
W(m) dǫ =: dǫw (ǫ) w (ǫ∗ ) β2(n+1)−m (β 1), (S31)
≡ n ǫ2+ω2 m ∼ m m ∼ ≫
Z0 j=0 j Z0
Q (cid:0) (cid:1)
because w′ (ǫ∗ )=0 ǫ∗ 1/β. Then ∆ˆ 1/β and ∆ˆ =0.
m m ⇒ m ∼ (∞,β) ∼ (n,∞)
(iv) A (ǫ) e−c/ǫa (a,c>0)
+
∼
At lastlet us considerthe case where the spectrum is gaplessandits function is exponential. The asymptotic form
gains non-trivial exponents:
X(a,c,2)
∆ˆ(∞,β) ∼sX(a,c,0) ∼β−a+11, (S32)
Y(a,c,2n 1)
∆ˆ(n,∞) ∼sY(a,c,2n−+1) ∼n−a1, (S33)
where
X(a,c,m) ǫM dǫ e−ǫcaǫm =: ǫM dǫe−xm,p(ǫ) 2π e−xm,p(ǫ∗m,p) e−βa+a1β−a2+(a2++21m) (β 1),
≡Z0 sinh β2ǫ Z0 ≈sx′m′,p(ǫ∗m,p) ∼ ≫
(cid:16) (cid:17) (S34)
Y(a,c,m) ǫM dǫe−ǫcaǫ−m ∞dǫe−ǫcaǫ−m = c−ma−1Γ m−1 m ma (m 1). (S35)
≡ ≈ a a ∼ ace ≫
Z0 Z0 (cid:18) (cid:19) (cid:16) (cid:17)
In Eq. (S34) x′′ (ǫ) is the second derivative and the saddle-point approximation around its extremum x′ (ǫ∗ ) = 0
m m m,p
was used. In Eq. (S35), Γ(t) is the gamma function.
10
Remarkably, the asymptotic behavior of the gap estimate is perfectly consistent also with the continuum spectrum
in all of the cases that we investigate here; that is
lim lim = lim lim ∆ˆ =∆ , (S36)
(n,β) 1
β→∞n→∞ n→∞β→∞
including the gapless (∆ =0) case.
1
The continuumspectrumthatwe haveinvestigatedappearsinthe thermodynamiclimit ofmanyrealisticquantum
systems. The asymptoticformshouldbe consideredtocheckthe convergenceofthe gapestimate andunderstandthe
spectrum structure. We note, nevertheless, that a gapless mode in critical phases obeys the asymptotic expression
obtained in the discrete formalism. It is because the ratios between the first gap and the higher gaps are kept even
though the system size increases. In other words, the low-energy-excitation spectrum is unchanged with the scale
transformation, which is certainly characteristic in critical phases. We hence securely applied the discrete formalism
to confirm the convergence in the present analysis of the spin-Peierls model.
III. COMPARISON OF ERROR-BAR ESTIMATIONS
We will show the comparison of the error-bar estimations for the relevant spin-Peierls model [Eq. (9) for L = 4,
ω =4,λ=1/2,D(cutoff)=4,inthemaintext]betweenthegapestimatorsandtheoptimalfit. Inthefittingmethod,
an optimal τ∗ is selected so that the gap error is minimized by the linear regressionof logC(τ) for τ τ β/4. To
1 1 ≤ ≤
checkthevalidityoftheerrorbar,weshowthemeanandthestandarderrorofthenormalizedquantityǫ¯=(∆ˆ ∆ )/σˆ
1
−
inTableII,where∆ 1.111388istheexactgapvalue,∆ˆ andσˆ arethegapanditserrorbar,respectively,estimated
1
from each simulation≃of 220( 106) Monte Carlo steps. Then the mean and the standard error of the normalized
∼
quantity were calculated from 2048 independent simulations. The mean should approach zero if the estimation is
unbiased, and the standard error should become one if the error bar is appropriately estimated. While the bias of
the optimal fit becomes smaller as temperature decreases, the standard error is significantly large as shown in the
table. It is because the correlation between data at different imaginary times is ignored and the estimated error bar
is improperlytoo small. This inappropriateerror-barestimationcausesthe deviationfromthe exactvalue to become
typically 20σ ( 2.22+18.93, which are in the table) even at T = 1/8. If the error bar is naively used for another
∼
analyses,e.g., the extrapolationto the thermodynamic limit as we showedin the main text, it might end up a wrong
conclusion. Ontheotherhand,ourhigher-ordergapestimatorisasymptoticallyunbiased,andtheerrorbarisreliably
estimated (see Table II). This is a clear advantage of the present approach, and it makes the precise analysis of the
criticality possible as demonstrated in the present paper.
We note that more careful statistical analyses like bootstrapping[S1] would improve the estimate of the error
bar even in the fitting method (also in the optimal fit). It is, however, necessary to run an additional (Monte
Carlo) simulation that is somewhat costly for users. Moreover, one has to make sure that the number of (almost)
independent bins is largeenough. Otherwise the varianceof the errorbar becomes improperly largeand the resultof
the bootstrapping is unreliable. The needed number of floating-point operations scales as O((nm+αm)M) once a
fitting region (τ) is fixed, where n ( 102–103 typically) is the number of bins, m ( 102–103) is the number of data
points(atdifferentτ), αmisthecom∼putationalcostforaregression,α( 10–103)d∼ependsonregressionschemesand
thenumberofparameters,andM ( 103–104)isthenumberofbootstra∼psamples. Thenthewholeprocessincluding
the fitting-region (τ ) optimization∼costs O((nm+αm2)M). It would take a few hours for large system sizes with
1
TABLE II. Mean µ¯(T) and standard error σ¯(T) of the normalized quantity ǫ¯=(∆ˆ −∆ )/σˆ estimated from 220(∼106)
1
Monte Carlo steps by the optimal fit and the gap estimators (1 ≤ n ≤ 5) at temperature T = 1/6 and 1/8. For the
calculationofeachvalue,2048independentsimulationswererun. Thestatisticalerrorsindicatedintheparenthesiswere
estimated by bootstrapping.
Estimator µ¯(1/6) σ¯(1/6) µ¯(1/8) σ¯(1/8)
optimal fit −33.75(48) 21.76(33) −2.22(59) 18.93(55
present approach
n=1 4.070(22) 1.008(15) 2.470(22) 1.012(16)
n=2 0.490(22) 0.992(15) 0.168(22) 1.005(16)
n=3 0.131(22) 0.993(15) 0.022(22) 1.003(16)
n=4 0.057(22) 0.993(15) −0.003(22) 1.002(16)
n=5 0.033(22) 0.993(15) −0.011(22) 1.002(17)
