Efficient Continuous-time Quantum Monte Carlo Method for the Ground State of Correlated Fermions

0.81 MB·English

by  Lei Wang

Efficient Continuous-time Quantum Monte Carlo Method for
the Ground State of Correlated Fermions
Lei Wang1, Mauro Iazzi1, Philippe Corboz2 and Matthias Troyer1
1Theoretische Physik, ETH Zurich, 8093 Zurich, Switzerland and
2Institute for Theoretical Physics, University of Amsterdam,
Science Park 904 Postbus 94485, 1090 GL Amsterdam, The Netherlands
We present the ground state extension of the efficient quantum Monte Carlo algorithm for lat-
ticefermionsofarXiv:1411.0683. Basedoncontinuous-timeexpansionofimaginary-timeprojection
operator,thealgorithmisfreeofsystematicerrorandscaleslinearly withprojectiontimeandinter-
actionstrength. ComparedtotheconventionalquantumMonteCarlomethodsforlatticefermions,
this approach has greater flexibility and is easier to combine with powerful machinery such as his-
5 togram reweighting and extended ensemble simulation techniques. We discuss the implementation
1 ofthecontinuous-timeprojectionindetailusingthespinlesst−V modelasanexampleandcompare
0 thenumericalresultswithexactdiagonalization,density-matrix-renormalization-groupandinfinite
2 projected entangled-pair states calculations. Finally we use the method to study the fermionic
quantum critical point of spinless fermions on a honeycomb lattice and confirm previous results
n
concerning its critical exponents.
a
J
PACSnumbers: 02.70.Ss,71.10.Fd,71.27.+a
5

]
l I. INTRODUCTION instabilities of the original approach have been remedied
e
in Refs. 23 and 24. The BSS algorithm has become the
-
r Quantum Monte Carlo (QMC) methods are power- methodofchoiceofmanylatticefermionsimulationsdue
t
s ful and versatile tools for studying quantum phases to its linear scaling in the inverse temperature β. We re-
.
t and phase transitions. Algorithmic development in past fer to Refs. 25 and 26 for pedagogical reviews.
a
m two decades including the nonlocal updates1–5 and the Closely related is the Hirsch-Fye algorithm,27 which
continuous-timeformulations6,7 havegreatlyboostedthe is numerically more stable and is more broadly applica-
-
d power of QMC methods, even surpassing the hardware ble because it is formulated using a (potentially time-
n improvements following Moore’s law. Using the mod- dependent) action rather than a Hamiltonian. How-
o ern QMC methods, the simulation of bosons and unfrus- ever, its computational effort scales cubically with the
c trated spin models is considered as a solved problem. inverse temperature and the interaction strength there-
[
QMCsimulationsthereforecanbeusedtotestnovelthe- fore is much less efficient than the BSS method for the

1 oretical scenarios8–12 and to verify experimental realiza- caseswherebothmethodsareapplicable. TheHirsch-Fye
v tions.13 methodthushastypicallybeenusedinthestudyofquan-
6 While efficient algorithms exist for the simulation of tum impurity problems and as impurity solvers in the
8 bosonsandunfrustratedspinmodels,1–5,14,15 simulations framework of dynamical mean field theory (DMFT),28
9
of fermions are more challenging than bosons and quan- where time-dependent actions need to be simulated.
0
0 tum spins because of the infamous fermionic sign prob- Both the BSS and the Hirsch-Fye algorithm are based
. lem.16,17 It causes exponential growth of computational on a discretization of imaginary time, thus introducing
1
effort as system size or inverse temperature increases. a systematic time step error, called the Trotter error.
0
5 Even for systems without a sign problem, the phase di- Nearly twenty years ago it was realized that the time-
1 agram of correlated fermions can be nontrivial to estab- discretization is not necessary for the simulation of lat-
: lish,18,19 nottomentiontoaccuratelydeterminetheuni- tice models.6,7 Besides increased accuracy due to the
v
versalityclassandassociatedcriticalexponents.20,21 The absence of a Trotter error, continuous imaginary time
i
X mainreasonforthisdifficultyistheunfavorablesuperlin- formulations often results in a more efficient and flexi-
r ear scaling with system size and/or inverse temperature ble algorithm.3 In Ref. 29 a first continuous-time QMC
a
of determinantal quantum Monte Carlo methods, which method for lattice fermions has been proposed. How-
are the workhorse of correlated lattice fermion simula- ever the scaling of this algorithm and numerical stabi-
tions. lization have not been discussed in this paper and we
Determinantal QMC method sums a factorially large are not aware of any application of the algorithm. Fur-
number of fermion exchange processes into a matrix de- ther development on fermionic continuous-time QMC
terminant,therebyavoidingthefermionsignproblemsin algorithms38 have focused on quantum impurity prob-
certaincases. Thefirstalgorithmbasedonthisideaisthe lems: the continuous-time interaction expansion (CT-
Blankenbecler-Scalapino-Sugar(BSS)method.22 Itmaps INT) algorithm31, continuous-time hybridization expan-
aninteractingfermionicsystemtofreefermionsinaspa- sion (CT-HYB) algorithm33 and the continuous-time
tially and temporally fluctuating external field and then auxiliaryfield(CT-AUX)32 algorithm. CT-INTandCT-
performs Monte Carlo sampling of this field. Numerical AUXarebasedonweak-couplingexpansionoftheaction
2
TABLE I. Comparison between various determinantal QMC methods for fermions. The ground state methods are extensions
of the corresponding finite temperature methods. They have similar scalings when replace the inverse temperature β by the
projection time Θ. N denotes the number of correlated sites and V denotes the interaction strength.
Lattice Models Impurity Models
Method name BSS — LCT-INT LCT-AUX Hirsch-Fye CT-INT CT-AUX CT-HYB
Finite temperature Ref. 22 Ref. 30 Ref. 27 Ref. 31 Ref. 32 Ref. 33
Ref. 29
Ground state Ref. 23, 34, and 35 This paper — Ref. 36 Ref. 37 — —
Trotter error Yes No No No Yes No No No
Auxiliary field Yes Yes No Yes Yes No Yes No
Scaling βVN3 a b βVN3 βVN3 (βVN)3 (βVN)3 (βVN)3 eN
a Althoughthenumberofoperationsdoesnotexplicitly dependontheinteractionstrengthV,oneneedstoincreasethenumberoftime
slicesproportionaltoV tokeepaconstantcouplingstrengthwiththeauxiliaryfield,i.e. toretainthesameleveloffluctuations.
b ThescalingofthiscodeisunclearsinceitisnotdiscussedinRef.29andimportantimplementationdetailsaremissing.
and share the same scaling as the Hirsch-Fye method.39 inating the time discretization error. Moreover, the
These methods have revolutionized the simulation of continuous-time formulation has greater flexibility for
quantum impurity problems and DMFT calculations.38 measuring observables and can easily be combined with
However, for lattice models they remained suboptimal histogram reweighting44,45 and extensive ensemble simu-
compared to the BSS method due to their cubic scal- lation14,46,47 techniques.
ingintheinversetemperature. Veryrecentlyanefficient The organization of this paper is as follows, in Sec-
continuous-time algorithm has been developed by two of tion II we introduce a model system of spinless fermions
the authors that scales identically to the time-honored that we will use to explain the algorithm in Section III.
BSS method30 and can be used both with an auxiliary Section IV contains comparisons of the method with
field (LCT-AUX) and without (LCT-INT). The prefix other numerical approaches and results on the quantum
“L” indicating both their linear scaling and their appli- criticalpointofspinlessfermionsonahoneycomblattice.
cabilitytolattice models. InTableIwesummarizesome We end with discussions of future prospects in Sec. V.
properties of these determinantal QMC methods.
Finite-temperature determinantal QMC methods can
II. MODEL
be extended to projector formulations,23,34–37 where the
ground state is obtained from imaginary time projection
To make the presentation of our algorithm more con-
ofatrialwavefunction. Inadditiontobeingmoredirect
crete, we will consider the following spinless fermion
toaddressquantumphasesatzerotemperature,thepro-
model at half-filling:
jectorformalismoftenallowsforadditionaloptimizations
such as symmetry and quantum number projections40,41
and combinations with fixed-node ideas in the presence Hˆ =Hˆ +Hˆ , (1)
omfearicsiaglnstparboiblilzeamti.o42nIanlstohebeccaosemoefstehaesiBerSSinmtehtehogdro,unnud- Hˆ0 =−t0(cid:88)(cid:16)1cˆ†icˆj +cˆ†jcˆi(cid:17)≡(cid:88)cˆ†iKijcˆj, (2)
state formulation.24,26 On the other hand, for projection (cid:104)i,j(cid:105) i,j
mjecettihoondtsimiteisscirnucceiatlhteoraecshuiletvseaarelienxeaacrtscoanlliynginintthheelpimroit- Hˆ1 =V (cid:88)(cid:18)nˆi− 21(cid:19)(cid:18)nˆj − 12(cid:19), (3)
of infinite projection time.43 The ground state variants (cid:104)i,j(cid:105)
of the Hirsch-Fye and the CT-INT methods36,37 exhibit
where cˆ is the fermion annihilation operator. t denotes
cubic scalingandthusarenotidealforlatticemodelsim- i
the nearest-neighbor tunneling, V > 0 denotes the ex-
ulations.
tendedHubbardrepulsiveinteraction,andwehaveintro-
In this paper we present details of the projection duced the matrix K to denote the single particle matrix
version of the LCT-INT method whose feasibility has elements.
already been mentioned in Ref. 30. This algorithm Quantum Monte Carlo studies of this model on a
provides an efficient continuous-time projection QMC square lattice date back to the early days of the BSS
approach for ground state simulations of correlated method.48,49 However, these simulations suffer from the
fermions. It retains the linear scaling with projection fermion sign problem because the Monte Carlo weight is
time and matches the one of the widely applied pro- a single determinant which is not guaranteed to be posi-
jector BSS method18–20,23,34,35 while completely elim- tiveingeneral. Recentlyitwasdiscoveredthatthemodel
3
(1) is naturally free from the sign problem on bipartite III. ALGORITHM
latticesathalf-fillingintheCT-INTformulation,50,51 be-
cause the Monte Carlo weight can be expressed as the A. General Description
determinant of a real skew-symmetric matrix. This de-
terminantequalsthesquareofthematrixPfaffianandis
InaprojectorQMCcalculationoneobtainstheground
thus nonnegative. A conventional auxiliary field decom-
state wave function using imaginary time projection of
position,ontheotherhand,breaksthissymmetry. Itwas
a trial wave function |Ψ (cid:105) and calculate ground state
T
shown that this model also allows sign problem free sim-
observables as34,35
ulation in the BSS formalism if one works in a Majorana
fermion representation,52,53 i.e. performs the auxiliary
fielddecompositionnotinthedensitychannelbutinthe (cid:104)Ψ |e−ΘHˆ/2Oˆe−ΘHˆ/2|Ψ (cid:105)
pairing channel. The idea applies not only to the BSS (cid:104)Oˆ(cid:105)= T T . (4)
(cid:104)Ψ |e−ΘHˆ|Ψ (cid:105)
algorithm but can be generalized to the continuous-time T T
QMC algorithm.30
For any |Ψ (cid:105) with non-vanishing overlap with the true
T
groundstate,Eq.(4)approachesthegroundstateexpec-
tation in large Θ limit. In this paper, we choose |Ψ (cid:105) as
T
a single Slater determinant,
(cid:89)NP (cid:32)(cid:88)N (cid:33)
Onthehoneycomblattice,thismodelexhibitsaquan- |Ψ (cid:105)= P cˆ† |0(cid:105), (5)
T ij i
tumphasetransitionfromaDiracsemimetaltoacharge-
j=1 i=1
density-wave (CDW) phase. The quantum critical point
is unconventional because of the coupling of CDW or- whereN isnumberofsites,NP isthenumberofparticles,
der parameter to the low-energy Dirac fermions.54,55 and P is a N ×NP rectangular matrix.56
Simulations using CT-INT found a critical point at Instead of breaking the projection operator into
V /t = 1.356(1) with critical exponents η = 0.302(7) small time steps as done in discrete time algo-
c
and ν = 0.80(3).51 Although CT-INT is free from the rithms,22–24,26,34,35 the continuous-time QMC formalism
time-discretization error, its cubic scaling with inverse writestheprojectionoperatorinaninteractionrepresen-
temperature β limited these simulations to inverse tem- tation
peratures βt≤20. To access the quantum critical point
fromafinitetemperaturesimulationβwasscaledlinearly (cid:34) (cid:35)
(cid:90) Θ
withthelinearextentofthesystem,assumingadynami- e−ΘHˆ =e−ΘHˆ0T exp − eτHˆ0Hˆ e−τHˆ0dτ . (6)
1
calcriticalexponentz =1. InSec. IVBwewill,asafirst
0
application of the projector LCT-INT algorithm, use it
to directly address the quantum critical point of model After a Taylor expansion of the exponential and time
(1) at zero-temperature and check our previous findings. ordering the terms57 the denominator of Eq. (4) reads
∞ (cid:90) Θ (cid:90) Θ (cid:90) Θ
(cid:104)ΨT|e−ΘHˆ|ΨT(cid:105)=(cid:88)(−1)k dτ1 dτ2... dτk(cid:104)ΨT|e−(Θ−τk)Hˆ0Hˆ1...Hˆ1e−(τ2−τ1)Hˆ0Hˆ1e−τ1Hˆ0|ΨT(cid:105). (7)
k=0 0 τ1 τk−1
In the CT-INT and the CT-AUX methods,31,32,37 one respect to the average expansion order. The algorithm
applies Wick’s theorem to the integrand of Eq. (7) and scales as ΘVN3, similar to the BSS algorithm.22,24,26
expresses it as a determinant of a matrix whose size is To proceed, we first express the interaction term
proportional to the expansion order k. The subsequent Hˆ through exponentials of bilinear fermion operators.
1
simulationmodifiesthematrixwithO(k2)operationsper Traditionally this is done via Hubbard-Stratonovich
Monte Carlo step. Since it takes k Monte Carlo steps to transformation, at the cost of introduction of auxiliary
generate an uncorrelated sample, these methods31,32,37 fields.29,32 Here we adopt a simpler approach based on
scale cubically with the average expansion order (cid:104)k(cid:105). As the operator identity nˆi = 21(1−eiπnˆi). The interaction
(cid:104)k(cid:105) increases linearly with inverse temperature β (or Θ term Eq. (3) can then be expressed as
in the ground state projection scheme) and interaction
strength,31,32,37 this unfavorable cubic scaling limits the
applicability of these methods for lattice models at low Hˆ1 = V4 (cid:88)eiπ(nˆi+nˆj). (8)
temperatureandstronginteractions. Hereweinsteaduse (cid:104)i,j(cid:105)
the LCT-INT algorithm30 to achieve linear scaling with
This reformulation works for any density-density
4
interaction58 and is crucial to respect the symmetry of lation. Substituting Eq. (8) into Eq. (7), the integrand
the model (1), ensuring a sign problem free QMC simu- is recognized as a sum of determinants
(cid:104)Ψ |e−ΘHˆ|Ψ (cid:105)=(cid:88)∞ (cid:18)−V (cid:19)k (cid:88) (cid:88) ... (cid:88) (cid:90) Θdτ (cid:90) Θdτ ...(cid:90) Θ dτ ×
T T 4 1 2 k
k=0 (cid:104)i1,j1(cid:105)(cid:104)i2,j2(cid:105) (cid:104)ik,jk(cid:105) 0 τ1 τk−1
(cid:104) (cid:105)
det P†e−(Θ−τk)KX(ik,jk)...X(i2,j2)e−(τ2−τ1)KX(i1,j1)e−τ1KP , (9)
where K is defined in Eq. (2). The vertex matrix X(i,j) consists of ordered times 0 ≤ τ < τ < ... < τ < Θ
1 2 k
isanN×N diagonalmatrixdefinedfornearestneighbors andcorrespondingpairsofnearestneighborsites(i ,j ),
1 1
i and j whose non-zero elements are (i ,j ) ...(i ,j ). An example of a configuration with
2 2 k k
three vertices is shown in Fig. 1.
(cid:40)
−1, l=i or l=j, Using Eq. (11) the weight of a configuration is ex-
X(i,j) = (10)
ll 1, otherwise. pressed as
Its form follows immediately from Eq. (8). In the follow- w(C)=(cid:18)−V (cid:19)kdet[P†B(Θ,0)P]dτ ...dτ . (14)
ingwedenotethematrixproductfromimaginarytimeτ 4 1 k
to τ(cid:48) >τ as a propagator
Since this weight is, up to a constant factor, identical
B(τ(cid:48),τ)=e−(τ(cid:48)−τm)KX(im,jm)...X(il,jl)e−(τl−τ)K, to the weight of such a configuration in a CT-INT cal-
(11) culation50,51 at zero temperature,37 these methods have
where τ and τ are the imaginary time locations of the identical sign problems.
first andl the lamst vertices in the time interval. The cor- Thequantummechanicalaverage(cid:104)Oˆ(cid:105)C,τ ofanoperator
responding vertex matrices are X(i ,j ) and X(i ,j ) Oˆ inserted into a configuration C at imaginary time τ:
l l m m
respectively. If there is no vertex in the time interval,
the propagator simply reads B(τ(cid:48),τ)=e−(τ(cid:48)−τ)K. (cid:104)Oˆ(cid:105)C,τ =
wrWiteetthheenexepxpecatnadtitohnevnaolumeininataorfoorfmEqs.u(it4a)bsliemfiolarrMlyoanntde (cid:104)ΨT|e−(Θ−τk)Hˆ0eiπ(nˆik+nˆjk)...Oˆ...eiπ(nˆi1+nˆj1)e−τ1Hˆ0|ΨT(cid:105)
Carlo sampling, (cid:104)ΨT|e−(Θ−τk)Hˆ0eiπ(nˆik+nˆjk)...eiπ(nˆi1+nˆj1)e−τ1Hˆ0|ΨT(cid:105)
(cid:80) w(C)(cid:104)Oˆ(cid:105) (cid:68) (cid:69) can be evaluated using Wick’s theorem since Oˆ is sand-
(cid:104)Oˆ(cid:105)= C (cid:80)Cw(C)C,Θ/2 = (cid:104)Oˆ(cid:105)C,Θ/2 MC, (12) wiched between two Slater determinants.
where the configuration C denotes a point in the sum-
mation and integration domain of Eq. (9) and (cid:104)...(cid:105) B. Monte Carlo Sampling
MC
denotes Monte Carlo averaging according to the configu-
ration weights w(C). Inthissectionwefirstexplainhowtosampleusingthe
For k vertices a configuration weights Eq. (14) and then discuss efficient ways to per-
form update and measurement using equal-time Green’s
C ={(τ1;i1,j1),(τ2;i2,j2),...,(τk;ik,jk)}. (13) function.
i3 i2 i1 ∆ 1. General Procedure
To sample configurations C according to the weight
Θ τ3 τ2 τ1 0 w(C) we use the Metropolis-Hastings algorithm.59,60
Starting from a configuration C we propose to move
j j j
3 2 1 to a new configuration C(cid:48) with an a-priori probability
A(C →C(cid:48)). The new configuration is then accepted with
FIG. 1. A configuration with k = 3 vertices. We divide the probability p(C → C(cid:48)) = min{1,r(C →C(cid:48))}, where the
imaginarytimeaxisintointervalsofsize∆andsweepthrough acceptance ratio r is
them sequentially. In each interval we propose updates that
eitherinsertortoremoveavertex. Measurementisperformed w(C(cid:48))A(C(cid:48) →C)
at the center of the imaginary time τ =Θ/2. r(C →C(cid:48))= . (15)
w(C)A(C →C(cid:48))
5
To facilitate fast computation of the acceptance rate, where N is the number of interacting bonds of the sys-
b
we divide the imaginary time axis into intervals of size tem, ∆ is the length of the time interval on which we
∆, as shown in Fig. 1. We focus our updates on one in- propose updates and n is the number of vertices in this
terval at a time, proposing several times to either insert interval. While insertion and removal updates are suffi-
or to remove an existing vertex at time τ for site indices cient to ensure ergodicity of the sampling, one can nev-
(cid:104)i,j(cid:105). Sweeping through the intervals we achieve ergod- ertheless implement additional updates to improve the
icity. While such sequential updates violate the detailed sampling efficiency, such as change the site index of a
balance condition, a global balance condition is still re- vertex (see Appendix A).
stored as long as the updates within each interval satisfy Afterafullsweepthroughalltheintervals,wemeasure
local detailed balance.61 the expection values of observables close to the center
Using short hand notations,24,26 τ =Θ/2.
L(τ)=P†B(Θ,τ) and R(τ)=B(τ,0)P, (16)
2. Fast Update Using Equal-time Green’s Function
the insertion or removal changes R to R± =X(i,j)±1R.
Note that for the model Eq. (1) studied here R+ = R− Crucial for the performance of the algorithm is a fast
becauseX(i,j)−1 =X(i,j),butthegeneralMonteCarlo calculation of the acceptance ratios Eqs. (17-18). They
scheme does not rely on this property. The acceptance can be efficiently computed from the equal time Green’s
ratios are functions G (τ) = (cid:104)cˆcˆ† (cid:105) , which in matrix form
lm l m C,τ
reads24,26
r =−det(LR+) × VNb∆ , (17) G(τ)=I−R(τ)[L(τ)R(τ)]−1L(τ). (19)
add det(LR) 4(n+1)
The determinant ratio in Eqs. (17-18) can be expressed
det(LR−) 4n
r =− × , (18) using the Green’s function as24,26
remove det(LR) N V∆
b
− det(LR±) =−det(cid:8)I+[X(i,j)±1−I](I−G)(cid:9)=−det(cid:18)1−2(1−Gii) 2Gij (cid:19)
det(LR) 2Gji 1−2(1−Gjj)
=4G G (20)
ij ji
The second equality follows from [X(i,j)±1 −I] = imaginary time τ(cid:48), which can be done by the following
lm
−2δ δ −2δ δ . Withappropriatelychosentrialwave similarity transformation (assuming τ(cid:48) >τ)
li im lj jm
function, the equal-time Green’s function of our model
(1)hasanimportantsymmetrypropertywhichweprove G(τ(cid:48))=B(τ(cid:48),τ)G(τ)[B(τ(cid:48),τ)]−1. (23)
in Appendix B:
G (τ)=δ −η η G (τ), (21) Propagating the Green’s function using the tricks dis-
ji ij i j ij
cussedinSec.IIIC,Eq.(23)ismoreefficientthancalcu-
where η = ±1 for site i belongs to the A(B) sublattice. lating G(τ(cid:48)) from scratch using Eq. (19).
i
SimilartothecaseofCT-INT,50,51Eq.(21)impliesG =
ii
1 and G = G for nearest neighbors. With this we
2 ij ji
canfurthersimplifythedeterminantratiosto4G G =
ij ji
4G2 . Since the remaining factors in Eqs. (17-18) are all 3. Observables
ij
positive,thereisnosignprobleminthesimulationofthe
model (1).50,51 All expectation values (cid:104)Oˆ(cid:105) can be related to the
C,τ
IfaproposedMonteCarlomoveisaccepted,weupdate Green’s function using Wick’s theorem. For example,
the Green’s function to G± =I−R±(LR±)−1L using thedensity-densitycorrelationfunctioncanbeexpressed
as (cid:104)nˆ nˆ (cid:105) = (1−G )(1−G )+(δ −G )G .
G (G −δ ) G (G −δ ) l m C,τ ll mm lm ml lm
G± =G − lj im im − li jm jm . (22) Using Eq. (21), the density-density correlation functions
lm lm Gij Gji is measured as
For a proof see Appendix C. (cid:28)(cid:18) (cid:19)(cid:18) (cid:19)(cid:29)
1 1
In the Monte Carlo updates we always keep track of C(l,m)= nˆ − nˆ −
l 2 m 2
the Green’s function at the imaginary time τ for which
(cid:28) (cid:18) (cid:19)(cid:29)
we propose an update. For the next Monte Carlo move Θ
=η η G2 , (24)
we need to propagate the Green’s function to a different l m lm 2
MC
6
the staggered CDW structure factor as diag(E ,E ,...,E ). The basis change modifies the
1 2 N
(cid:28)(cid:18) (cid:19)(cid:18) (cid:19)(cid:29) propagators to
1 (cid:88) 1 1
M = η η nˆ − nˆ −
2 N2 l m l 2 m 2
1 (cid:88)l,m(cid:28) (cid:18)Θ(cid:19)(cid:29) (cid:0)U†e−τKU(cid:1)lm =e−Elτδlm, (29)
= N2 G2lm 2 , (25) (cid:0)U†X(i,j)U(cid:1)lm =δlm−2Ul†iUim−2Ul†jUjm. (30)
l,m MC
and the kinetic energy and interaction energy as
In this basis the multiplication of G˜ by either Eq. (29)
(cid:18) (cid:28) (cid:18)Θ(cid:19)(cid:29) (cid:19) or Eq. (30) requires only O(N2) operations instead of
(cid:104)Hˆ (cid:105)=−Tr K G , (26)
0 2 O(N3).30 The disadvantage is that now the calcula-
MC
(cid:28) (cid:18) (cid:19)(cid:29) tion the determinant ratio Eq. (20) is slightly more ex-
(cid:104)Hˆ (cid:105)=−V (cid:88) G2 Θ . (27) pensive. However since we only need one matrix ele-
1 (cid:104)l,m(cid:105) lm 2 MC ment Gij = (UG˜U†)ij which can be calculated using
matrix-vector multiplication and vector inner-products,
Anotherusefulobservableistheaverageexpansionor-
this O(N2) overhead will not affect the overall scaling of
der
the algorithm. Similarly, updating of the Green’s func-
(cid:42)(cid:90) Θ (cid:43) tion in the eigenbasis also keeps an O(N2) scaling (see
(cid:104)k(cid:105)=− (cid:104)Hˆ1(cid:105)C,τdτ . (28) Appendix D).
0
MC WorkingintheeigenbasisofthenoninteractingHamil-
Sincethereisnotranslationalinvariancealongtheimag- tonian does not increase the complexity of the measure-
inary time axis, (cid:104)k(cid:105) is not directly related to the interac- ments either. One can choose to measure single par-
tionenergy(cid:104)Hˆ (cid:105)asitisinthefinitetemperaturecase.31 ticle observables in the eigenbasis and perform a ba-
1
Nevertheless, Eq. (28) still suggests (cid:104)k(cid:105)∼ΘVN, i.e. the sis rotation afterwords. Alternatively one can rotate G˜
average number of vertices scales linearly with the pro- backtoGwithO(N3)operationsforeachmeasurement.
jectiontime,theinteractionstrengthandthesystemsize. Sincemeasurementsareperformedonlyafterafullsweep
The fact that we are dealing with k of N ×N matrices throughallintervals,thisdoesnotaffecttheoverallscal-
compared to the single 2k ×2k matrix of the CT-INT ing of the algorithm. For many physical observables of
case51 allows LCT-INT to achieve an O(ΘVN3) scaling, interest we only need Gij for neighboring sites (cid:104)i,j(cid:105) or
as we will discuss in the next section. for fixed site i because of translational invariance, which
would reduce the required basis transformation to just
O(N2) operations.
C. Algorithm Optimization and Complexity
AchievingthesamescalingofO(ΘVN3)asintheBSS
algorithm requires a careful implementation, for which 2. Optimal Interval Size
an optimal choices of single particle basis and splitting
imaginary time into intervals is crucial.
Finally we show that by choosing the number of in-
tervals M = Θ/∆ proportional to the average number
1. Optimal Single-Particle Basis of vertices one can achieve an overall O(ΘVN3) scal-
ing in the algorithm. For each MC update, we need
to propagate the Green’s function from some time τ
The main computational effort in performing the
to another time τ(cid:48) in the same interval. This will on
Monte Carlo updates is the propagation of the Green’s
average pass through |τ−τ(cid:48)|(cid:104)k(cid:105) existing vertices, which
function to a new imaginary-time using Eq. (23). Imple- Θ
is of order O((cid:104)k(cid:105)/M). As we need O(N2) operations
mented naïvely this involves dense matrix-matrix multi-
plication and requires O(N3) operations, while the cost to pass through each vertex and O(N2) for calculat-
ing the acceptance rate and for the actual update, one
of the calculation of the determinant ratio in Eq. (20)
Monte Carlo step requires O(max{1,(cid:104)k(cid:105)/M}N2) opera-
is O(1) and update of the Green’s function Eq. (22) is
O(N2). tions,wherethemaxfunctionaccountsforthecaseofan
empty interval.
This unfavorable scaling can be circumvented by
working in the eigenbasis of the noninteracting Hamil- Asweepthroughallintervalsandupdating(cid:104)k(cid:105)vertices
tonian.30 In this way all the computations for MC results in an overall number of O(max{1,(cid:104)k(cid:105)/M}N2(cid:104)k(cid:105))
steps Eqs. (23,20,22) can be performed with complexity operations. By choosing the number of intervals M ∼
O(N2). (cid:104)k(cid:105) we can achieve an optimal scaling O((cid:104)k(cid:105)N2) =
For this we have to use basis-transformed Green’s O(ΘVN3). This should be compared to the scaling of
functions G˜ = U†GU, where U are the eigen- other continuous-time methods which scale as (cid:104)k(cid:105)3 ∼
vectors of the single-particle Hamiltonian U†KU = Θ3V3N3.31,32,37,38
7
D. Numerical Stabilization V = S(cid:96)+1B˜(τ(cid:96)+1,τ). The Green’s function is then re-
L
computed by G˜ =I−U (V U )−1V .
R L R L
As in the BSS algorithm the multiplication of the In the simulation we monitor the difference of the sta-
Green’s function with the propagator B(τ(cid:48),τ) for large bilized G˜ and the old one to dynamically adjust the fre-
imaginary-time suffers from numerical instabilities be- quencies of the SVD and the recomputation of G˜. It
cause the matrix multiplication mixes large and small turns out both frequencies are mainly set by the inverse
scalesinthepropagator. Westabilizethecalculationfol- bandwidthandareindependentofthesystemsizeorthe
lowingasimilarapproachasusedfortheBSSmethod.23 total projection time. Typically we need to perform one
The following discussion largely follows Refs. 24 and 26 of such stabilizations for a propagation time τ ∼ 1/t.
with the difference that our stabilization is done in con- Since each of these stabilization steps costs O(N3) due
tinuous time and in the eigenbasis of the single-particle to the SVD or matrix inverse and we need to perform
Hamiltonian. O(Θ) of them per sweep, it ends up with a scaling of
To avoid accumulation of numerical errors we need to O(ΘN3), conforming with the overall scaling of the al-
regularly recompute the Green’s function using G˜ =I− gorithm.
U†R(LR)−1LU, which requires us to have fast access to
thematricesU†RandLU inanumericalstableway. We
thusdividetheimaginarytimeaxisintoI intervalswhere E. Calculation of the Renyi Entanglement Entropy
within each interval the propagation is well-conditioned.
The interval length is set by the inverse bandwidth and Quantuminformationbasedmeasuresplayanincreas-
is independent of the total projection time Θ. These ingroleintheidentificationofquantumphasesandphase
intervals are different from the (shorter) intervals used transitions.62–65 In particular, Refs. 66–68 devised mea-
for MC updates discussed above. surements of the Renyi entanglement entropy in deter-
At the interval boundaries (corresponding to imagi- minantal QMC simulations.
nary times τ(cid:96)=0 = 0,...,τ(cid:96)=I = Θ) we store I +1 ma- Since in the present algorithm the many-body ground
trices S(cid:96). Depending on the current imaginary time τ statewavefunctionisexpressedasasumofSlaterdeter-
of the Monte Carlo sweep, they hold the matrix product minants, the derivations of Ref. 66 concerning the re-
either to the right or to the left, duced density matrix hold. In particular, the rank-2
U†R(τ(cid:96)), if τ >τ(cid:96), gRieonnyAi ecnatnanbgelecmalecnutlaetnetdroupsyingS,2 = −lnTr(ρˆ2A) of a re-
S(cid:96) = (31)
L(τ(cid:96))U, otherwise. (cid:80) w(C)w(C(cid:48))det[G G(cid:48) +(I−G )(I−G(cid:48) )]
e−S2 = C,C(cid:48) (cid:80) wA(CA)w(C(cid:48)) A A ,
On the rightmost and leftmost boundaries S(cid:96)=0 = U†P C,C(cid:48)
and S(cid:96)=I = P†U. The matrix S(cid:96) is updated whenever (34)
we cross the interval boundary τ(cid:96) in the sweep along the where C and C(cid:48) are configurations of two independent
replicasandG ,G(cid:48) arethecorrespondingGreen’sfunc-
imaginarytimeaxis. Duringasweepfromτ =0toΘ,we A A
multiply the propagator B˜(τ(cid:96),τ(cid:96)−1) ≡ U†B(τ(cid:96),τ(cid:96)−1)U tion restricted to the region A.
with S(cid:96)−1 to update S(cid:96) = B˜(τ(cid:96),τ(cid:96)−1)S(cid:96)−1. In the Although the estimator Eq. (34) is easy to implement,
we observe it shows large fluctuations at strong interac-
backward sweep from τ = Θ to 0, we update S(cid:96) to
S(cid:96) =S(cid:96)+1B˜(τ(cid:96)+1,τ(cid:96)). tion.67,69,70 We leave a discussion of extended ensemble
simulations67,68,71intheLCT-INTformalismforafuture
We still need to stabilize the calculation of S(cid:96)=0,...,I
study.
themselves. Performing a singular-value-decomposition
(SVD)
U†R=U D V , (32) F. Direct Sampling of Derivatives
R R R
LU =U D V , (33)
L L L
An advantage of continuous-time algorithms over
the different scales only appear in the diagonal matri-
discrete-time algorithms is that the Monte Carlo weight
ces of singular values D and D . Since G˜ = I −
L R Eq. (14) are homogeneous functions of the interaction
U (V U )−1V only depends on the well-conditioned
R L R L strengthV. Thisallowstodirectlysamplethederivatives
matricesU andV ,itissufficienttokeeptrackofthem
R L of any observable with respect to V using its covariance
instead of the full matrix products. Therefore before up-
with the expansion order k,
datingS(cid:96) wecanperformanSVDonthematrixproduct
BV˜(.τ(cid:96),τ(cid:96)−1)S(cid:96)−1 orS(cid:96)+1B˜(τ(cid:96)+1,τ(cid:96))andonlystoreUR or ∂(cid:104)Oˆ(cid:105) =(cid:42)∂Oˆ(cid:43)+ 1 (cid:16)(cid:104)Oˆk(cid:105)−(cid:104)Oˆ(cid:105)(cid:104)k(cid:105)(cid:17). (35)
L ∂V ∂V V
Using these stored matrices S(cid:96)=0,...,I we can easily re-
computetheGreen’sfunctionatanyimaginarytime. To Higherorderderivativescanbesampledinasimilarway.
compute G˜(τ) for τ(cid:96)+1 > τ > τ(cid:96) we can use the ma- Derivatives are useful for discovering quantum phase
trices S(cid:96)+1 and S(cid:96) and calculate U = B˜(τ,τ(cid:96))S(cid:96) and transitions and locating critical points.
R
8
100
10-1
)
r
(
C
r)
1
−
(
10-2
V/t=1
V/t=2
V/t=3
V/t=4
10-3
0 2 4 6 8 10 12 14 16
r
FIG.2. Density-densitycorrelationfunctionofa32-sitechain FIG. 3. Interaction energy (red dots), ground state energy
with periodic boundary condition. Solid lines are DMRG re- (yellow squares) and CDW structure factor M (blue trian-
2
sults. gles)ofanN =18sitehoneycomblatticeshownintheinset.
Solid lines are exact diagonalization results.
In discrete time approaches calculation of such ob-
servables either relies on the Hellmann-Feynman theo-
V/t=1.0
rem,72,73 which is limited to the first order derivative of 0.88
thetotalenergy(cid:104)Hˆ(cid:105),74,75 orrequiresnoisynumericaldif- iPEPS a
ferentiation of Monte Carlo data.18 0.89 QMC-PBC
QMC-APBC
0.90
IV. RESULTS
0.91
e
Wefinallypresentresultsobtainedwithouralgorithm, t
si 0.92
starting with benchmarks that demonstrate the correct- r
nessbeforepresentingnewresultsregardingthequantum pe V/t=1.4
critical point of the model Eq. (1). gy 0.95 iPEPS b
er QMC-PBC
n QMC-APBC
E 0.96
A. Benchmarks
0.97
For all of our results we use a projection time Θt=40
and use ground state of the noninteracting Hamiltonian
Hˆ0 asthetrialwavefunction. Incaseofdegeneratenon- 0.908.00 0.05 0.10 0.15 0.20 0.25
interacting ground states, we take as trial wave function 1/L or 1/D
thegroundstateofasystemwithanti-periodicboundary
conditioninx-directionandperiodicboundarycondition
FIG. 4. Ground state energy per site of a honeycomb lattice
along the y-direction.
versus inverse system length 1/L (QMC) or inverse bond di-
We start by showing results for a periodic chain. Fig- mension1/D(iPEPS).TheQMCresultsofperiodicboundary
ure 2 shows the density-density correlation function of a conditionsandanti-periodicboundaryconditionsapproachto
periodic chain compared with results from density ma- the thermodynamic limit value from different sides.
trix renormalization group (DMRG) calculations, where
C(r) is averaged over all pairs in Eq. (24) with the same
distancesr. Atmoderatecomputationalcostwecanper- ground state results for the total energy, interaction en-
fectly reproduce the exact ground state quantities using ergy as well as the staggered density structure factor.
projection LCT-INT (filled dots). Finally, we compare with infinite projected entangled-
Figure 3 compares the results obtained with our algo- pair states (iPEPS) results obtained for the honeycomb
rithmtoexactdiagonalization(solidlines)foranN =18 lattice.51,76,77iPEPSisavariationalmethodwhichworks
site honeycomb lattice. Our method correctly produces in thermodynamic limit, whose accuracy can be system-
9
V/t=1.0
0.060 0.60
iPEPS a L=6 a
0.045 QMC-PBC L=9
QMC-APBC L=12
0.030 +η 0.45 L=15
M2 zL L=18
2
0.015 M
0.30
0.000
1.30 1.32 1.34 1.36 1.38 1.40
V/t=1.4 V
0.060 0.60
iPEPS b L=6 b
0.045 QMC-PBC L=9
QMC-APBC L=12
0.030 +η 0.45 L=15
M2 zL L=18
2
0.015 M
0.30
0.000
0.00 0.05 0.10 0.15 0.20 0.25 2.5 2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
1/L or 1/D (V V )L1/ν
− c
FIG.5. TheQMCresultsfortheCDWstructurefactorcom- FIG. 6. (a) Scaled CDW structure factor of different system
pared with the square of CDW order parameter calculated sizescrossatthetransitionpoint(b)ScaledCDWstructures
using iPEPS. For V/t = 1 the CDW order parameter van- factorcollapseontoasinglecurvewhenplottedagainstscaled
ishes for bond dimensions D > 8 in iPEPS. For V/t = 1.4 interaction strength.
we used a linear fit in 1/D of the CDW order parameter to
obtainanestimateintheinfiniteDlimit(seeRef.51formore
details). up to L=18. Since a detailed finite size scaling study is
beyondthescopeofthispaper, weusethecriticalvalues
obtainedinRef.51andcheckforconsistency. TheCDW
atically improved by increasing the bond dimension D. structure factor should follow the scaling ansatz
Figure 4 shows the ground state energy per site versus (cid:16) (cid:17)
1/L together with iPEPS results versus 1/D. QMC re- M2Lz+η =F (V −Vc)L1/ν , (36)
sults for systems with periodic boundary conditions and
where we previously found z +η = 1.302, ν = 0.8 and
thoseanti-periodicboundaryconditionalongx-direction
V /t=1.356.51 Figure6(a)showsthescaledCDWstruc-
approach the L → ∞ limit from different sides, thus c
ture factor M Lz+η where all curves cross around V
bracketingthegroundstateenergyinthethermodynamic 2 c
whenusingthesecriticalexponents. Scalingofthex-axis
limit. Extrapolation of all data yields consistent results.
using (V −V )L1/ν yields good data collapse, shown in
Figure5showstheCDWstructurefactorM versus1/L, c
2
Figure6(b). Weconcludethatthenewzerotemperature
which extrapolates to the square of the CDW order pa-
resultsonlargersystemsizeareconsistentwithprevious
rameter. iPEPS on the other hand can directly mea-
findings concerning critical point and critical exponents
sure the order parameter since the symmetry is spon-
in Ref. 51.
taneously broken for an infinite system. Extrapolation
again yields consistent results and shows the system or-
ders for V/t=1.4 but not at V/t=1.0.
V. DISCUSSION
In this paper we presented details of the ground-state
B. Fermionic Quantum Critical Point
version of the LCT-INT algorithm of Ref. 30. As a
continuous-time QMC algorithm it eliminates the Trot-
We finally apply the projector LCT-INT to study the ter error due to time discretization of the BSS algorithm
quantum critical point of spinless t-V model on a hon- while still keeping the favorable linear scaling with pro-
eycomb lattice, which we previously studied by CT-INT jection time and interacting strength. It is therefore well
simulations.51 Our calculations go beyond the previous suitedforsimulationsofthegroundstateofstronglycor-
resultsintwoaspects. WecandirectlyaddresstheT =0 related lattice fermions.
quantum critical point using the projection version of Although the LCT-INT algorithms30 and the projec-
LCT-INT and we are able to reach larger system sizes tion version described here share operational similari-
10
ties with the BSS algorithm,22–24,26,34 there are impor- VI. ACKNOWLEDGMENTS
tant differences. In the BSS formalism, one breaks the
projection operator e−ΘHˆ into small discrete time steps The authors thank Fahker Assaad, Jakub Imriška and
andperformsTrotter-Suzukidecompositionforeachtime HiroshiShinokaforhelpfuldiscussions. Simulationswere
step, which leads to a systematic time-discretization er- performed on the Mönch cluster of Platform for Ad-
ror. The BSS algorithm then decouples the interaction vanced Scientific Computing (PASC), the Brutus clus-
terms using auxiliary fields. A typical update scheme is ter at ETH Zurich and the “Monte Rosa” Cray XE6
tosweepthroughthesetimeslices24,26 andfliptheauxil- at the Swiss National Supercomputing Centre (CSCS).
iaryfields,similartoourschemeofsweepingthroughthe We have used ALPS libraries84 for Monte Carlo simula-
intervals. However, the time slices of the BSS algorithm tions and data analysis. DMRG results in Fig. 2 have
are fixed in time and their number is proportional to the been obtained using mps_optim application85 of ALPS
projection time. Each time slice contains O(N) auxil- project. This work was supported by ERC Advanced
iary fields, therefore even with a brute force propagation Grant SIMCOFE, by the Swiss National Science Foun-
of the Green’s function on the site basis (Eq. (23)) one dation through the National Competence Centers in Re-
canachieveO(N3)scaling. Whileinourcasethenumber search QSIT and MARVEL, and the Delta-ITP consor-
and positions of vertices are allowed to fluctuate so we tium(aprogramoftheNetherlandsOrganisationforSci-
need to propagate in the eigenbasis (Eqs. (29-30)) and entificResearch(NWO)thatisfundedbytheDutchMin-
use M ∼ (cid:104)k(cid:105) intervals such that on average each inter- istry of Education, Culture and Science (OCW)).
val contains a single vertex to achieve a similar O(N3)
scaling.
Appendix A: Site-Shift Update
Formally the Monte Carlo weight in Eq. (7) is similar
to the local weight of the CT-HYB method.33,78 In par- An optional additional update shifts a vertex between
ticular the matrix version of CT-HYB78 also evaluates X(i,j) by moving the site j to another neighbor of the
the Monte Carlo weight in the eigenbasis of a propaga- site i denoted by j(cid:48). This will change the vertex matrix
tor. However, our case is simpler because e−τHˆ0 is a to X(i,j(cid:48)). It amounts to insert a vertex matrix X(j,j(cid:48))
at the same imaginary time without changing the per-
single particle propagator and the Monte Carlo weight
turbation order. The acceptance ratio is
simplifies to a determinant instead of than a trace in
the many-body Hilbert space in CT-HYB. The present
det(LX(j,j(cid:48))R)
method can still benefit from algorithmic developments r = =−4G G . (A1)
for the CT-HYB method. In particular a Krylov ap- shift det(LR) jj(cid:48) j(cid:48)j
proach for imaginary time propagation79 may bring the
cost of propagation of the Green’s function to O(N2). Since the site j and j(cid:48) belongs to the same sublattice,
Our “sweep through intervals” scheme is also similar to Gjj(cid:48) = −Gj(cid:48)j (see Eq. (21)) ensures the acceptance ra-
theslidingwindowapproachoftheCT-HYBalgorithm.80 tio is positive. The formula for updating the Green’s
One may alternatively consider using trees or skip list function after a shift move is identical to Eq. (22), with
data structure to store partial matrix products.38,81 indices i,j replaced by j,j(cid:48).
Besides being free of the discretization error, the Appendix B: Proof of Eq. (21)
continuous-time QMC approach provides a direct means
to compute quantities such as observable derivatives
Equation (21) is easiest to prove in the finite temper-
Eq. (35), that are harder to obtain in discrete time
ature formalism. Suppose the trial wave function |Ψ (cid:105)
simulations. These in turn may be used to locate in- T
is the ground state of a noninteracting trial Hamiltonian
teresting phase transitions with an accuracy that can-
K . The equal time Green’s function can be formally
not easily be reached by standard discrete-time algo- T
written as
rithms. Furthermore, the simple interaction dependency
of the Monte Carlo weight w(C) ∼ Vk allows straight- G(τ)= lim (cid:2)I+B(τ,0)e−βKTB(Θ,τ)(cid:3)−1. (B1)
forwardcombinationwiththehistogramreweighting44,45 β→∞
or the Wang-Landau sampling14,46 techniques. Both ap-
proachescanproduceresultsinacontinuous range ofin- We introduce a diagonal matrix D = η δ and the bi-
ij i ij
teraction strength by recording histograms over the per- partite conditions implies DKD = −K. Together with
turbation order k. Combined with Eq. (35), the method X(i,j)−1 = X(i,j) this shows that (B(τ ,τ )−1)T =
1 2
offers new exciting opportunities to bring the study of DB(τ ,τ )D. Similarly, assuming the trial Hamiltonian
1 2
quantumcriticalityofcorrelatedfermionstoanewlevel, alsosatisfiesDKTD =−KT,onehaseβKT =De−βKTD.
approaching to what was achieved in the simulations of Combing these facts it is then straightforward to show
the classical82 and quantum spin systems.83 that