Quantum Monte Carlo simulation in the canonical ensemble at finite temperature

0.2 MB·English

by  Kris Van Houcke

Quantum Monte Carlosimulationinthecanonical ensemble atfinite temperature
K. Van Houcke,1 S.M.A. Rombouts,1 and L. Pollet1
1Universiteit Gent - UGent, Vakgroep Subatomaire en Stralingsfysica
Proeftuinstraat 86, B-9000 Gent, Belgium
6 (Dated:February2,2008)
0
AquantumMonteCarlomethodwithanon-localupdateschemeispresented.Themethodisbasedonapath-
0
integraldecompositionandawormoperatorwhichislocalinimaginarytime. Itgeneratesstateswithafixed
2
numberofparticlesandrespectsotherexactsymmetries. Observablesliketheequal-timeGreen’sfunctioncan

n beevaluatedinanefficientway. Todemonstratetheversatilityofthemethod,resultsfortheone-dimensional
a Bose-Hubbardmodelandanuclearpairingmodelarepresented.WithinthecontextoftheBose-Hubbardmodel
J
theefficiencyofthealgorithmisdiscussed.

1
3 PACSnumbers: 02.70.Ss,05.10.Ln,21.60.Ka,71.10.Li

]
h
I. INTRODUCTION the(d+1)-dimensionalXY universalityclass[20]. Whenen-
c
tering the superfluid phase, it becomes difficult to keep the
e
m numberofbosonsconstantandtuningthechemicalpotential
QuantumMonteCarlo(QMC)simulationisapowerfuland
can become a very time-consuming task. When simulating
- versatilemethodfortheinvestigationofthermodynamicprop-
t mesoscopic systems like superconducting grains [21, 22] or
a erties of many-body systems. When generating a Markov-
t atomic nuclei [23, 24], it is primordial to keep the number
s chain of configurations using a Metropolis scheme [1], it is
ofparticlesconstant. Aworld-linealgorithmsatisfyingthese
t. knownthatupdatesbasedonlocalchangesareinefficient,par- conditionsispresented[25]anddiscussedindetailinthispa-
a
ticularly near critical points. At transitional points this type
m per. Besides particle-numberconservation, the algorithm al-
of algorithm gives very large autocorrelation times, a phe-
lowstoincludeothersymmetry-projections. Itisconstructed
- nomenon one refers to as ’critical slowing down’ [2]. By
d insuchawaythatallmovesareaccepted,whichmakesiteffi-
n developingnon-localupdateschemes, thisproblemhasbeen cienttorunandeasiertocode.Thoughworkinginthecanoni-
o overcomeforsecondorderphasetransitions.TheWolffsingle calensemble,ouralgorithmisstillabletogenerateconfigura-
c cluster algorithm[3] and the Swendsen-Wang multiple clus-
tionswithdifferentwindingnumbers,incontrasttothelocal
[ ter algorithm [4], both used to solve classical physics prob-
world-linemethodbyHirschetal. [26]. Wewilldemonstrate
2 lems, were successful applications of this idea. In the same the versatility of the method by applyingit to bosonsand to
v spirit, loopalgorithms[5, 6, 7]weredevelopedforthestudy pairedfermions.
1 ofdiscretequantumsystems.Newconfigurationsareobtained
3 by flipping clusters in the form of loops. The systematic er-
4
rorcausedbythetheSuzuki-Trotterapproximation[8,9]was
II. THEALGORITHM
8
eliminated by formulatingthe world-line algorithmsdirectly
0
in continuous imaginary time [10, 11]. In the worm algo-
5 PracticallyallQMCmethodsarebasedonadecomposition
0 rithm[12],thepartitionfunctionisembeddedinanextended of the evolutionoperatore−b H. The trick is to break up this
/ configuration space, allowing a direct and exact evaluation
t operatorintopieceswhichcanbehandledexactly[27]. Gen-
a of the one-bodyGreen’s function. The conceptof non-local
erallyonecanwritetheHamiltonianasconsistingofaneasy
m loopupdateshasalsobeenimplementedintheStochasticSe-
partH andaresidualinteractionV,
- ries Expansion (SSE) representation [13], leading to ’opera- 0
d tor loop’ update algorithms [14, 15] and ’directed loop’ al- H=H −V. (1)
n 0
gorithms [16, 17], which are a further optimization of the
o
loopconstruction.Thegeneralideabehindthisistoconstruct Theminussignhasbeenincludedtoeasenotationsfurtheron.
c
: movesinalocallyoptimalway[18]. For such a Hamiltonian, one can make an exactperturbative
v
expansioninV oftheevolutionoperator,usingthefollowing
i All the non-local update world-line algorithms which are
X integralexpression[11]:
mentioned above sample the grand-canonical ensemble. In
ar tmhiasgnweatiyzaotnioencoarnogcecnuepraattieonconnufimgbuerar.tioRnessuwltisthfoer.gth.evcaarnyoinng- e−b H = (cid:229)¥ b dt tmdt ··· t2dt V(t )
Z mZ m−1 Z 1 1
ical ensemble are obtained by using only the configurations m=0 0 0 0
with the right particle number [19] or by rejecting loop up- V(t2)···V(tm)e−b H0, (2)
dateswhichchangethisnumber[6]. Itisclearthatthisisan
inefficientwayofworking. Samplingthecanonicalensemble withV(t)=exp(−tH )Vexp(tH ) and b the inverse temper-
0 0
directlywouldbeadvantageouswhenstudyingsystemswhere ature (also called imaginary time). The basic idea of the
particle-conservationisimportant. Oneexampleisthetransi- continuous-time loop algorithm [6, 10] and the worm algo-
tion between the superfluid and the Mott phase in the Bose- rithm [12] is to insert two adjointworld-line discontinuities.
Hubbardmodelatconstantfilling. Thistransitionbelongsto Bypropagatingoneofthesediscontinuities(whicharecalled
2
the mobile ’worm head’ and stationary ’worm tail’) through for all basis states |ii), the sum of the weights of all diago-
latticespaceandimaginarytime,theconfigurationchangesin nalconfigurationsisproportionaltotheparticle-numberpro-
suchawaythatdetailedbalanceisfulfilled. Atthatpointone jected trace of the evolution operatorU(b ). This is nothing
isnotsamplingaccordingtothepartitionfunctionTr(e−b H), elsethanthecanonicalpartitionfunctionTr (e−b H),withTr
N N
butaccordingtoanextensionhereof, theparticle-numberprojectedtrace.Hence,samplingthecon-
figurationsproportionalto their weightsW(m,i,t,t ) leads to
Tr(W†e−t HWe−(b −t )H), (3) asamplingofthecanonicalensemble.TheMarkovprocessis
with t the imaginarytime intervalbetween the worm ’head’ setupusingtheMetropolis-Hastingsalgorithm[1,28],hereby
samplinginanextendedspaceaccordingtoTr [U′(b ,t )]. At
operatorW† and ’tail’ operatorW. The worm head can be N
eachMarkovsteponlyafewofthefactorsofEq. (5)areal-
creating or annihilating, dependingon the choice ofW. Af-
tered by the worm operator which moves to a new point in
ter some Markov steps, the worm head bites its tail and the
imaginary time. These worm operator moves will be con-
discontinuities are removed. Only configurations with con-
structed in a such a way that detailed balance is fulfilled
tinuous world-linescan contribute to the statistics according
to Tr(e−b H). In contrast to these algorithms, we will work locally at each Markov step. Therefore it is also fulfilled
when going from one diagonal configuration to another. It
withawormwhichislocalinimaginarytime. Theevolution
takesa numberofMarkovstepsbeforediagonalobservables
operatorextendedbysucha localworm(animaginarytime-
(i.e. observableswhich commute with H ) can be measured
independentoperatorAtobedefinedbelow)reads 0
again. While each Markovstep containsonlylocalchanges,
U′(b ,t )=e−t HAe−(b −t )H, (4) the chain of steps between two diagonal configurations cor-
respondsto a global update. Non-diagonaloperatorscan be
where t can be regardedas the imaginary time at which the measuredusingthefactthatonesamplesinanextendedspace.
worm operator is inserted. We will show that by working By keeping track of the worm moves between two diagonal
withalocalwormoperatoronecanconstructaveryefficient configurations, statistics for the expectation values of non-
samplingmethod,providedthatthewormoperatorcommutes diagonal operators can be collected, similar to the way one
with the residual interaction (AV =VA). If A furthermore evaluates the one-body Green’s function in the worm algo-
commutes with the generators of a symmetry of H and V, rithm [12]. Our method however will lead to much better
0
onecanrestrictthesamplingtoconfigurationswiththosespe- statistics for equal-time non-diagonaloperators, because the
cific symmetries, leading to symmetry-projected results. In wormoperatorisalwayslocalinimaginarytime.
particularonecansamplethecanonicalensemblewithaworm Before shifting the worm operator over some imaginary
operatorthatconservesparticlenumber.Wewouldliketoem- time interval Dt , a direction D has to be chosen. One can
phasizeatthispointthatthealgorithmstatedbelowdoesnot choosebetweenthedirectionsD=R(highervaluesoft )and
dependon the specific structureof A. The operatorA has to D=L(lowervaluesoft )with someprobabilityP(D), tobe
bechoseninsuchawaythatanergodicMarkovchaincanbe specified later. The presence of the exponentials in Eq. (5)
constructed,andthereforeitwilldependonthespecificform inspiresusto choosethe timeshiftDt proportionaltoanex-
oftheinteractionV. ponentialdistribution,
Ifonedecomposesthetrace(restrictedtothewantedparti-
cle numberand symmetry)ofU′(b ,t ) using Eq. (2) and in- P(Dt )dDt =e De−e DDt dDt , (6)
sertscompletesetsofeigenstatesofH atallimaginarytimes,
0
one obtains a set of integrals which can be evaluated using withe D anoptimizationparameter. Inshiftingthewormop-
MonteCarlosampling. TheMarkovprocesswillsamplethe eratorfromt toanewimaginarytimet ′ =t +Dt ,theworm
configurationsproportionaltotheweights operatorcanencounteraninteractionoperatorV atsometime
in between. Assume thedirectionR is chosenandtheworm
W(m,i,t,t )=hi0|V|i1ie−(t2−t1)Ei1hi1|V|i2ie−(t3−t2)Ei2··· operatormeetsaninteractionattimetR. Letusfirstconsider
hiL−1|V|iLie−(t −tL)EiLhiL|A|iRie−(tR−t )EiRhiR|V|iR+1i··· tthereascittiuoant,iownitwhohuetreanthneihwilaotrimngoipt,eratormovesthroughthisin-
him−1|V|imie−(tm−tm−1)Eimhim|V|i0ie−(b +t1−tm)Ei0. (5)
hi |A|i ihi |V|i i−→hi |V|i′ihi′ |A|i i. (7)
L R R R+1 L R R R+1
Each configuration is specified by an order m (the num-
ber of interactions), a set i of eigenstates of H0 (with When passing the interaction, the intermediate state can
i a shorthand notation for all the intermediate states change. A convenient way to pick one of these possible
|i0i,...,|iLi,|iRi,...,|imi),interactiontimest1,...,tm,andthe changesistochoosethenewconfigurationproportionaltoits
worm insertiontime t . We use the notationEij =hij|H0|iji. weight
The configuration|i i to the left of the worm is changedby
L
thewormoperatorintotheconfiguration|iRi. Weusethesub- P (i′ )= hiL|V|i′Rihi′R|A|iR+1i. (8)
scriptL (R=L+1)toindicatetheeigenstate|iLi(|iRi)and R+1,L R hiL|VA|iR+1i
interaction time t (t ) just before (after) the worm operator
L R
in imaginarytime. We willcallthe configurationsforwhich Part (a) of Figure 1 shows a diagrammatic representation of
i =i diagonalconfigurations.Bychoosingthewormopera- thedifferentwaysinwhichageneralone-bodywormoperator
L R
torsuchthatitsdiagonalelementsareconstant(i.e.hi|A|ii=c A=(cid:229) c a†a (for some constants c ) can pass a similar
i,j ij i j ij
3
interactionV. The worm operator is represented by a curly ofinteractions,
lineandtheinteractionbyaverticalstraightline.Forthistype
W(m,i′,t′,t ′)P(i′,t′,t ′→i,t,t )
ofwormandinteractionoperatortherearealwaysatmostfour q =
waysinwhichtheintermediatestatecanchange.Itshouldbe W(m,i,t,t )P(i,t,t →i′,t′,t ′)
notedthatourchoiceEq. (8)isnotuniqueandpossiblymore = (e D′)initial, (13)
optimal choices can be found [18]. Because of this choice (e )
D final
however,thereappearsafactor
where(e D)final ((e D′)initial) is the valueofe D (e D′) atthe end
hi |V|i′ihi′ |A|i iP (i ) (beginning)ofthewormoperatormoveintodirectionD,and
n = L R R R+1 L,R+1 R
R+1,L hi |A|i ihi |V|i iP (i′ ) D′denotestheoppositedirectionofD. Theactualacceptance
L R R R+1 R+1,L R
probabilityisgivenbymin(1,q),accordingtotheMetropolis-
hi |VA|i i
L R+1
= , (9) Hastingsalgorithm.
hi |AV|i i
L R+1 Let us now introduce a number of steps, which allow to
change the number of interactions in a reversible way. We
intheacceptancefactoroftheMetropolis-Hastingsalgorithm.
wanttheacceptancefactorofsuchupdatestobelocal,i.e. the
Everytime the worm operatorpasses an interaction an anal-
probabilitytopass, createorannihilateaninteractionshould
ogousfactorappears,dependingonthedirectionDofpropa-
only depend on the properties of the state at that point in
gation.Thereforeit’sadvantageoustoimposeonAthecondi-
imaginary space-time. We consider three extensions of the
tion,
procedureoutlinedabovewherenointeractionsarecreatedor
deleted:
AV−VA=0, (10)
• AtthebeginningoftheMarkovstep, we introducethe
because then nDD′ =1 in all cases, and we do not have to possibility that the worm operatorcreates a new inter-
worryaboutthisnormalizationfactoranymore. Furthermore,
actionwithprobabilityc ,whichdependsonthedirec-
D
Eq. (10)ensuresthatthewormoperatorcanalwayspassthe
tion D of propagation. This creation will also change
interactionitencounters. Ifonewouldchooseawormopera-
theintermediatestate:
torAthatdoesnotsatisfythiscondition,asingrand-canonical
algorithms, it is possible the worm operator cannot pass the hi |A|i i−→hi |V|i′ihi′|A|i i, (14)
L R L R
interactionandthedirectionofpropagationhastochange,in
that way undoing changes previously made. It is intuitively assuming again the worm operator is moving in the
clearthatthese’bounces’giverisetoaslowdecorrelationand D=R direction. The new intermediate state |i′i will
shouldbeavoided[16,17,18]. Inthedirectedloopalgorithm bechosenwithprobabilities
one increases the efficiency of the loop update by minimiz-
hi |V|i′ihi′|A|i i
ingtheappearanceofbounces,buttheycannotbeeliminated P (i′)= L R . (15)
RL
hi |VA|i i
completelybecauseofthereversibilitycondition.Inwhatfol- L R
lowswewillassumetheconditionEq. (10)isfulfilled,mak-
ForawormoperatormoveintheD=Ldirection,prob-
ingthealgorithmbounce-free. We willdropthefactorsnDD′ abilities PLR(i′) can be defined in an analogous way.
to ease the equations. After passing through the interaction
Figure1(parts(b)and(c))showsadiagrammaticrep-
attimet , one hasto choosea newimaginary-timeshiftDt .
D resentationoftheinsertionofaone-bodyinteractionat
However,onecanavoidgeneratinganewrandomnumberby
the beginningof the worm move. For a diagonalcon-
takingthenewshiftasfollows:
figurationonlythediagonalpartofAcontributestothe
matrixelementhi |A|i i. Inthissituationthewormop-
(e ) L R
Dt =(t ′−t ) D old , (11) eratorisrepresentedbylittlecirclesandallworld-lines
D (e )
D new arecontinuous.Wewillcallthisthe’diagonalworm’.
where the parameter e D has been updated after passing the • When the worm operatorarrivesat an interaction, one
interaction. also hasto considerthe possibility of annihilatingthat
Thechoiceoftheparameterse D followsfromdetailedbal- interaction.Supposetheinteractioncanbedeleted.Let
ance. Because the time shifts Dt of the worm operator are a betheprobabilitytoannihilatetheinteractionwhile
D
chosenbyEqs. (6)and(11),addingtheconstraint continuingthe worm update, and s the probabilityto
D
annihilate the interaction and stop the worm update.
ER−EL=e L−e R, (12) Then1−aD−sDistheprobabilitytopassthroughthat
interactionandcontinuethewormoperator.
ensures that all the exponentials which appear in the ac-
ceptancefactorofthe Metropolis-Hastingsalgorithmcancel, • To maintain reversibility, one also has to include the
leading to an efficient sampling method. So in practice one possibilitythattheconstructionoftheMarkovstepdoes
can choose any positive value for e and e , as long as Eq. nothaltatthe momentthe wormoperatorhasfinished
L R
(12)isfulfilledateachstepofthewormmovement. Tocon- a shiftDt withoutencounteringaninteraction. Atthat
clude,wewritedowntheacceptancefactorfortheabovepro- pointonehastochoosebetweenstoppingthewormop-
cedurewhenthewormoperatordoesnotchangethenumber erator,ortocontinue,withthepossibilityofinsertinga
4
newinteractionatthatpoint. Let f betheprobability where
D
tocontinuethewormoperatorwithoutinsertinganin-
teraction,gDtheprobabilitytoinsertaninteractionand qD=e D′+NDD′(sD′−aD). (22)
to continue the worm operator, then 1− f −g will
D D
Thefactor(q ) hastobeevaluatedatthebeginning
be the probability to stop the worm operator, without D initial/final
and the end of the Markov step in direction D (with D′ the
insertinganinteraction.
oppositeofD). ThecreationprobabilityisgivenbyEq. (17),
After creating, annihilating or passing an interaction, a new
time shift Dt should again be chosen according to Eqs. (6) c = NDD′sD′. (23)
or (11). Note that the parameter fD is redundant: jumping D qD
with aparametere andcontinuingthewormoperatorunal-
D
teredwithprobability f isstatisticallyequivalenttomaking WestillhavetodeterminehowtochoosethedirectionD. The
D
ajumpwithparametere (1−f ),andthenchoosingbetween acceptanceratioofEq.(21)inspiresustochoosethedirection
D D
eitherstoppingthe wormoperatoror insertingan interaction ofthemovewithprobabilities
and move on. Therefore, one can set f = 0 without loss
D q
R
of generality. We now determine how the other parameters P(D=R) = , (24)
R
shouldbechoseninordertosatisfydetailedbalance.Assume LR
q
L
adirectionDischosen.Whennointeractionisinsertedatthe P(D=L) = . (25)
R
beginningofthewormmove,afactor LR
q0 = e D′(1−gD′), (16) dwiistthriRbuLtRio=nqgRiv+enqbLy. Byacceptingallwormoperatormovesa
D 1−c
D
appearsintheacceptancefactor.Ifontheotherhandaninter- W′(m,i,t,t )=RLRW(m,i,t,t ), (26)
actioniscreatedattheinitialtimet ofthewormoperator,this
will be sampled, instead of the distributionW(m,i,t,t ). Be-
willleadtoafactor
causethefactorsR fluctuateonlymildlyinpractice,accept-
LR
qc = NDD′sD′, (17) ing allmovesstill leads to a a veryusefulsamplingmethod.
D c Itspeedsupthealgorithmandreducesthecomplexityofthe
D
code.Finitetemperatureobservablescanstillbeevaluatedby
with
takingtheextraweightingfactorintoaccount:
NDD′ = hhiDiD|V|AA|i|DiD′i′i. (18) hQib = TTrr[e[e−−b b(H(H0−0−VV)Q)]]
AnewintermediatestateischosenwithprobabilitiesEq.(15).
AttheendoftheMarkovstep,thewormoperatorwillannihi- = (cid:229) s∈S(Qs/(RLR)s). (27)
lateaninteractionornot,leadingtoextrafactorsintheglobal (cid:229) s∈S(1/(RLR)s)
acceptance factor which have the inverse form of Eqs. (17)
Eachtimethewormoperatorcreates,annihilatesorpasses
and(16),becauseoftheinversesymmetrybetweenbeginning
an interaction, the parameters e , c , a , s , g are deter-
andendofthemove.Atintermediatepoints,wecanencounter D D D D D
minedbyEqs. (12),(20),(22)and(23). Thisstillleavesalot
the followingsituations. The wormoperatorcan stopaftera
shift Dt between two interactions, insert an interaction and of freedom. We will focus here on two limiting cases. First
we will consider the case where one of the two parameters
move on. The inverse situation of this can also occur, when
e and e is always zero. In that way it can occur the time
aninteractionisannihilatedandthewormoperatormoveson. R L
shift Dt becomesinfinite. This amountsto the worm opera-
Orthewormoperatorcansimplypassaninteractionwithout
tordirectlyjumpingto the nextinteraction,whichspeedsup
annihilatingit. In orderto havea goodtotalacceptancefac-
decorrelationinimaginarytimedirection.Inordertoobtaina
tor, we will requirethatthese intermediatesteps do notcon-
wormmovethatchangestheconfigurationsasmuchaspossi-
tribute to the acceptance factor. This condition leads to the
ble,theparametersg ,g ,a anda aremaximized. Theset
constraints L R L R
ofparametersobtainedinthiswayisshowninTableIforthe
NDD′aD′ = e DgD, (19) case ER >EL. We will callthis solutionA. The solutionfor
thecaseE >E canbefoundbyinterchangingthesubscripts
aD+sD = aD′+sD′. (20) L R
LandR. Notethatinthissolutionthewormoperatoralways
Apart from that, we want the sampling to be as uniform as starts to move into the direction of the highest diagonal en-
possible, which suggests the condition q0 =qc. Putting all ergy. It is clear that wheneverthe worm operator is moving
D D
thistogether,thetotalacceptancefactorisgivenby in the direction of the highest diagonal energy or whenever
E =E ,thetimeshiftDt becomesinfinite. Thereareanum-
W(m′,i′,t′,t ′)P(m′,i′,t′,t ′→m,i,t,t ) R L
q = berofextraconditionsoneshouldkeepinmind. Assumethe
W(m,i,t,t )P(m,i,t,t →m′,i′,t′,t ′) wormoperatorstartstomoveinthedirectionD=R(because
= (qD)initial, (21) ER>EL). When the worm operatorarrivesat an interaction
(qD′)final thatcanbeannihilated,onehastodetermineEL,ER andNLR
5
parameters diagonalconfigurations diagonalconfigurations
(iL=iR) (iL6=iR)
(EL=ER) (EL<ER)
e R 0 0
e L 0 ER−EL
qR f ER−EL
qL f 0
N
cR 1 min(1,ER−LREL)
c 1 0
L
f
sR N 0
sL NfLR min(1,ERN−EL)
LR LR
aR 0 min(1,ERN−EL)
LR
a 0 0
L
N
RgLLR 20f minE(R1,−ERE−LLREL)
TABLEI: AsetofalgorithmparameterssatisfyingEqs. (12),(20),
(22)and(23)forthecasesEL=ER andEL<ER (otherwiseinter-
change thesubscriptsLandR). WecallthissolutionA,for which
oneoftheparameterse isalwayszero.
FIG.1: Adiagrammaticrepresentationforone-bodywormoperator D
moves. The worm operator is represented by a curly line and the
interactionV by a solid vertical line. We distinguish between the
parameters diagonalconfigurations diagonalconfigurations
followingupdates. (a)WhenthewormoperatormovesintheD=R
(iL=iR) (iL6=iR)
direction,itcanencountersomeinteraction. Thewormoperatorcan
(EL=ER) (EL<ER)
passtheinteraction,inthatwaychangingtheintermediatestate.For
e N N
generalone-bodywormandinteractionoperators,thereareatmost R LR LR
four possible ways in doing this. (b)-(c) At the beginning of the e L NLR NLR+ER−EL
wormmove,weintroducethepossibilityofinsertinganinteraction. qR NLR ER−EL
q N 0
Whentheinitialwormisdiagonal(representedbycircles),anumber L LR
of interaction insertions of the type shown in (b) are possible. In cR f 0
c f 0
part(c)theinitialwormoperatorisnotdiagonal,andaninteraction L
insertioncanmakethewormdiagonalornot. sR f 0
s f 0
L
aR f 1
aL f 1
g f 1
R
aaftefrrothmeaTnanbilheilIatairoen.thIefEcoRr>recEtLpirsosbtailblilsiatiteissfiteodstthoepn,osrRcaonnd- gL f NLR+NELRR−EL
R RLR 2NLR ER−EL
tinuethe wormoperator. Ifnow E >E onthe otherhand,
L R
onehastousethesolutions =min(1,EL−ER)anda =0,but
R N R TABLE II: An alternative set of parameters for which one of the
LR
the worm operatorkeepsmovingin the same direction. The parameterse DisalwaysNLR.WecallthissolutionB.Theparameter
timeshiftofthewormoperatorisonlyfinitewhenitmovesin f canbechosentooptimizethealgorithm.
thedirectionDandE <E′ . Notealsothatonlyg ismen-
D D L
tionedinTableI,becausetheparameterg hasonlymeaning
D
whenthetimeshiftisfinite. Inthepresentsolutionhowever, Thereforeweconsiderthecasewhereoneofthetwoparam-
aproblemariseswheneverEL=ER. Inthiscasee R=e L=0 eters e R and e L is NLR. As a consequence the time shift is
andthetimeshiftDt isalwaysinfinite. BecausesR =sL=0 alwaysfinite. ForER>ELsuchasetofparametersisgivenin
inaddition,thewormoperatorneverhalts. Asaconsequence Table II. Again, the case E =E needsan alternativesolu-
L R
configurationswith a diagonalworm will never be sampled. tion,sinceotherwiseR wouldbezerofordiagonalconfigu-
LR
This can be solved by proposing a small but finite stopping rations.WewillrefertothissolutionassolutionB.
probability. Thisalternativesolutionforthe case EL=ER is In short, the algorithm is based on a time-dependent per-
alsogiveninTableI. Theglobalparameterf shouldbetaken
turbation in the interaction V (see Eq. (2)), so there is no
small(suchthat0<f <NDD′ foralldiagonalconfigurations) systematicerrorarisingfromtimediscretization. Becausewe
butnotzero,andcanbechoseninordertooptimizethedecor-
choose time shifts of a worm operator according to Eq. (6)
relationbetweensuccessiveevaluationsofobservables. Note
thediagonalpartoftheHamiltonianishandledexactly.There
R ofEq. (26)takestheconstantvalue2f .
LR are a numberof algorithmparameters, which have to satisfy
Another possibility to find algorithm parameters follows Eqs. (12), (20), (22) and (23). We have derived two sets of
fromtheideathatwewantthestepsizeDt tobeoftheorder algorithmparameters, satisfying these equations. In the first
ofthemeanimaginarytimeintervalbetweentwointeractions. set (solutionA) one ofthe parameterse isalways zeroand
D
6
in the second (solution B) it is equal to N . Therefore the The one-body Green’s function G(r)=hb†b i is calcu-
LR i i+r
main differencebetween these two solutionswill be the size lated by updating the entry r in a histogram for G(r) at ev-
of the imaginary time shift Dt . Other algorithms where e ery Markov step. The function G(r) can be normalized di-
R
ore takevaluesbetween0andN canbeconstructedina rectly from the diagonal / non-diagonalworm fraction. The
L LR
similar way, takingfornowonlythese limitingcases. Inthe non-diagonal worm components b†b of Eq. (30) can be
i i+r
nextsectionourQMCalgorithmwillbeappliedtotheBose- givena differentweight, leadingto a wormmatrixrepresen-
Hubbardmodel. Theefficiencyofthetwodifferentsolutions tationofthesymmetricToeplitztype(i.e. asymmetricmatrix
leadingto differentalgorithmswillbe comparedin thiscon- withconstantvaluesalongnegative-slopingdiagonals). Such
text. awormoperatorstillcommuteswiththeinteractionpartV of
theHamiltonian. Bygivingsomewormcomponentsabigger
weight, the correspondingcomponentsG(r) will be updated
III. APPLICATIONTOTHEBOSE-HUBBARDMODEL more often, leading to a higher accuracy and mimicking flat
histograms[34].Thecondensedfractionr canbedetermined
c
Ultracoldbosonicatomsin anopticallattice aredescribed via
bytheBose-Hubbardmodel[29,30,31], r = 1 (cid:229)Nshb†b i. (31)
c NN i j
H=−t (cid:229)Ns b†b +U (cid:229)Ns n(n −1), (28) s i,j
i j 2 i i Thesuperfluidfractioncan bedeterminedusingthe winding
hi,ji i
number[35],
iw,inthi tbh†ie(nbui)mthbeerboopseornatcorreaotniosnit(eaninainhdilahit,iojni)doepneortaintogrnoenarseitset r s= hW2t2NiNb s2, (32)
neighbors. The lattice has N sites, occupied by N bosons.
s
The parametert is the tunneling amplitude andU is the on- wherehW2iisthemeansquareofthewindingnumberoper-
site repulsion strength. We will restrict the discussion here atorinonedimension. Figure2showsthecondensedandsu-
to the one-dimensional homogeneous case without trap. At perfluidfractionforauniformone-dimensionalsystemof128
low values of U/t the system forms a compressible super- sitesatadensityofexactlyoneparticlepersite. Calculations
fluid.Thisphaseischaracterizedbyagaplessexcitationspec- wereperformedatan inversetemperatureb =128t−1, using
trum and long-rangephase coherence. By increasingU/t, a thealgorithmbasedonsolutionA.Wehaveusedthealgorithm
quantum phase transition from a superfluid state to a Mott tostudythequantumcriticalbehavioroftheone-dimensional
insulating state is achieved at integer densities. In the pure Bose-Hubbard model with constant filling, using Renormal-
Mottphasethebosonsarelocalizedatindividuallatticesites izationGroupflowequations. StudyingtheBKTtransitionis
and all phase coherence is lost due to quantum fluctuations. notoriouslydifficultbecauseofthelogarithmicfinite-sizecor-
In addition, densityfluctuationsdisappearwhenenteringthe rections.Thepresentalgorithmhasthebigadvantageofkeep-
Mottphaseandagapappearsintheexcitationspectrum.This ingthedensityconstant,incontrasttothegrand-canonicalap-
phasedriventransitionbelongstotheBerezinskii-Kosterlitz- proaches.Theresultsofthisstudywillbepresentedelsewhere
Thouless(BKT)[32,33]universalityclassinonedimension. [36].
TheBose-HubbardHamiltoniancanberewrittenintheform We comparethe efficienciesof solutionsA and B. To this
Eq. (1), purpose, we simulate a one-dimensionallattice with 32 sites
ataninversetemperatureb =32t−1andafillingfactorofone
U (cid:229) bosonpersite.Thesimulationsconsistedof40Markovchains
H = n(n −1),
0 i i
2 thateachran600secondsafterthermalizationonaPentiumIII
i
V = t (cid:229) b†b . (29) processor.Thesamecodewasused,withonlyminorchanges
i j to go from algorithm A to B. We discuss the efficiency by
hi,ji
looking at the standard deviations on the expectation value
of V of Eq. (29) and on the average squared density. We
Asarguedabove,itisadvantageoustotakethewormoperator
calculatedthesquareddensityn2byaveragingoverallsites,
AsuchthatitcommuteswithV. Anoperatorthatsatisfiesthis
conditionisgivenby hn2i= 1 (cid:229)Nshn2i. (33)
A= 1 (cid:229) n +(cid:229) b†b , (30) Ns i i
N¯ i i i6=j i j The expectationvalue of the interactiontermV was notcal-
culatedviathecorrelationfunctionG(1),butbycountingthe
with N¯ a c-numberto be optimized. In our calculationsthis numberof interactionverticesin the configurationwhenever
parameterisalwayssetequaltothetotalnumberofparticles. the worm operator was diagonal [37]. When looking at the
Wehavecheckedourcodebycomparingwithexactdiagonal- standarddeviationonthesquareddensity(Figure3),onecan
izationresultsforsmalllattices. Ergodicitywastestednumer- concludesolutionAisthemostefficientone.Wefoundasim-
ically. Powerlawbehaviorofcorrelationfunctionscoincides ilarpicturewhenlookingatthestandarddeviationsontheex-
withpredictionsfrombosonizationtheoryforlargelattices. pectationvalueofH ofEq. (29). Theerrorsonthestandard
0
7
updates, and increases only slowly with increasingU/t. Of
1
r
0.9 r s course it should be noted the algorithm based on solution A
c hadtorunmuchlongerinordertogetthesamenumberofdi-
0.8
agonalmeasurements.Forallmeasuredobservableswefound
0.7 similarautocorrelationtimes.
0.6 We conclude that solution A, derived in the previoussec-
0.5 tion, is more efficient than solution B, except in the super-
r
fluid phase when looking at the interaction energyV. This
0.4
can be understood by remarking that the time shifts of the
0.3
worm operator are much larger for solution A. In the algo-
0.2
rithmbasedonsolutionB, the timeshiftsare ofthe orderof
0.1 themeanimaginarytimebetweentwo interactionverticesin
0 theconfigurations. Thisexplainswhythestandarddeviations
0 2 4 6 8 10 onthe interactionenergyaresmallest forthis solutionin the
U/t superfluidphase.Wealsoconcludethatouralgorithmismore
FIG. 2: Superfluid (r s) and condensed fraction (r c) for the one- efficientthanthedirectedloopSSE algorithmwhensimulat-
ing the one-dimensional Bose-Hubbard model. In the next
dimensional homogeneous Bose-Hubbard model as a function of
U/t. Thefractions havebeen calculated for auniform latticewith section, we will apply the algorithm to a pairing model. In
Ns=128ataninversetemperatureb =128t−1,usingEqs. (31)and what follows, all calculations are performed using the algo-
(32). The condensed fraction was calculated from the correlation rithmbasedonsolutionB.
functionhb†ibi+ri,forwhichwehavehighstatistics.
0.01
solution A
solution B
SSE
deviationsliewithintenpercent.Fromthestandarddeviation
ontheexpectationvalueoftheinteractiontermV (Figure4), 0.001
itfollowsthatsolutionBisbetterinthesuperfluidphase.The 2)
same conclusionfollowsfrom the total energy. For the con- (n
densedfractionthedeviationsaresmallestforsolutionA for s
1e-04
all values of U/t. Those on the superfluid fraction lie very
closeforsolutionAandB(seeFigure5).Wefoundthatvary-
ing the parameter f of Tables I and II does not change the
efficiency in a significant way, as long as f is not too small.
1e-05
Forfurthersimulationswe willalwayschoosetheparameter
1 2 3 4 5 6 7 8 9 10
f as big as possible, underthe constraintf ≤NDD′. Figures U/t
3, 4 and 5 show also standard deviations resulting from the
directedloopSSEalgorithm[16,17,18]. Onehastobevery FIG. 3: A comparison between the standard deviations on the
careful when comparing efficiencies of different algorithms. squareddensity(seeEq. (33))resultingfromthedirectedloopSSE
algorithmandthealgorithmsbasedonsolutionsAandB.Thehomo-
First, the SSE code works in the grand-canonicalensemble.
geneousBose-Hubbardmodelissimulatedforalatticewith32sites
In the SSE simulations, the chemical potential was changed
and32bosonsataninversetemperatureb =32t−1. Thedeviations
in such a way that N remained constant. Second, the effi-
resultedafteraQMCcalculationwith40independentMarkovchains
ciencydoesnotonlydependonthealgorithm,butalsoonthe
thateachran600secondsonaPentiumIIIprocessor.
used data structures. In a SSE approach, the decomposition
oftheevolutionoperatorcorrespondstoaperturbationexpan-
sioninalltermsoftheHamiltonian,whilethedecomposition
Eq. (2) perturbs only in the off-diagonal terms V. For the
Bose-Hubbardmodel,wherethecontributionofthediagonal
terms is large, the last approachis preferable. For all calcu- IV. APPLICATIONTOAPAIRINGMODEL
lated observables the standard deviations resulting from the
SSEcodewerethelargest. Figures3and4showtheSSEde- Inthenuclearshellmodel,quantumMonteCarlomethods
viationsincreaseratherrapidlywithincreasingU/t,whereas arevaluablebecausetheyofferthepossibilityofdoingcalcu-
thedeviationsresultingfromourmethodremainofthesame lationsinmuchlargermodelspacesthanconventionaldiago-
order. We also calculated autocorrelationtimes for different nalizationtechniques. Finitetemperatureshellmodelstudies
observables. Hereeachbinendedafteraconstantnumberof havebeendonewiththeaidofauxiliary-fieldQMCmethods
measurements. We noticed that for solution B the autocor- [23],andground-statepropertiesoflightnucleihavebeencal-
relation times became very big for high values ofU/t. For culatedusingvariationalanddiffusionQMCtechniques[38].
small U/t, the autocorrelation time for solution A is of the Furthermore,beingabletocalculatethermalpropertiesofnu-
orderofthenumberofMarkovstepsneededfor10diagonal clei makes it in principle possible to calculate nuclear level
8
k. Theoperatorn isthecorrespondingnumber-operatorand
0.0025 kt
solution A e =e thesingleparticleenergy-eigenvalue. G isthepair-
solution B kt kt t
ingstrengthforprotonsorneutrons.Proton-neutronpairingis
SSE
0.002
notincluded,butthis couplingcontributesonlyin an impor-
tantwayforN =Z nuclei[42]. Asaconsequence,theprob-
0.0015 lemdecouplesforprotonsandneutrons.Inthesequelonlythe
V) neutronpartofthemodelisconsidered,andtheisospinindex
(
s 0.001 t isdroppedtoeasethenotations.
Based on algebraic techniques developed by Richardson
[43],thepairingmodelcanbesolvedexactlyforanarbitrary
0.0005 setofsingleparticlelevelsatzerotemperature[44]. Inprac-
ticeitremainsdifficulttousetheseexactresultstostudythe
0 thermodynamicsof the model, because the numberof states
1 2 3 4 5 6 7 8 9 10 neededintheensembleincreasesveryrapidlywithincreasing
U/t temperature. Thermodynamicalpropertieshavebeenstudied
usingauxiliary-fieldQMC,whichisfreeofsignproblemsfor
FIG. 4: The standard deviation on the mean value of V (see Eq.
thepairingHamiltonianofEq.(34)whendealingwithaneven
29)forthehomogeneousBose-HubbardmodelasafunctionofU/t.
HeresolutionAisthemost efficientoneintheMott phase. Inthe numberofparticles[24]. However,thepresentalgorithmcan
superfluidphase,solutionBbecomesmoreefficient. considernucleiwithevenandoddnucleonnumbers.Notethe
auxiliary-field method scales as O(N3) with N the number
s s
consideredsingleparticlestates,whileaworld-linealgorithm
scaleslinearwithN .
0.006 s
solution A When a nucleon occupies a single particle state k and its
solution B
time-reversedstate k isunoccupied,thenucleonissaidtobe
0.005 SSE
’unaccompanied’. These states do notparticipate in the pair
scattering by H . The mean-field plus pairing Hamiltonian
0.004 P
canberewrittenasEq. (1),
)s 0.003
r( H = H −V, (35)
0
s
0.002 H = (cid:229) e n −G(cid:229) b†b , (36)
0 k k 2 k k
k k
0.001 V = G4 (cid:229) b†kbk′. (37)
0 k6=k′
1 2 3 4 5 6 7 8 9 10
U/t Theoperatorsb†=a†a†createapairofnucleonsintwotime-
k k k
reversedstatesandsatisfyhard-corebosoncommutationrela-
FIG.5: Thestandarddeviationonthemeanvalueofthesuperfluid
tions. Inordertogetthecorrectfinitetemperatureproperties,
fractionr s. Thedeviations result fromsolutions A and B,and the
thepossibilityofchangingthenumberofunaccompaniednu-
directedloopSSEmethod.Eachsimulationconsistedof40indepen-
dentMarkovchainsthateachran600seconds. cleonsduringthe simulationshouldbeincorporated. A path
integralMonteCarlomethodforthepairingHamiltonianhas
been developed by Cerf and Martin [45, 46], but there the
number of pairs remained fixed [24, 47]. This problem can
densities [39]. These densities are extremely important for
nowbeovercomeelegantlybyaddinganextrapairbreaking
makinggoodtheoreticalestimatesofnuclearreactionrates.
term
Thebasicassumptionintheshellmodelisthepresenceofa
mreseiadnufialelidntienrawcthiiocnhbthetewneuecnletohnesnmucolveeo.nTsoisiminptrroodvuecoend.thPias,ira- Vpert= G2g(cid:229) (cid:229) b†kak′ak′′+H.c. , (38)
k k′6=k′′(cid:0) (cid:1)
ing between nucleonsis the main short-rangecorrelationin-
duced by the residual interaction. Adding a simple pairing to the interaction part V of Eq. (37). We define the worm
Hamiltoniantothemean-fieldHamiltonian operatoras
Hmf+HP= (cid:229) (cid:229) ektnkt− (cid:229) G4t (cid:229) a†k′tak†′taktakt, (34) A= N1¯ (cid:229) nk+41 (cid:229) b†kbk′+g2(cid:229) (cid:229) b†kak′ak′′+H.c. ,
t=p,n k t=p,n k,k′ k k6=k′ k k′6=k′′(cid:0) (cid:1)
(39)
can account for this [40, 41]. The operators a† create a with two extra parametersN¯ and g to be optimized. A term
kt
particle in the single-particle eigenstate k of the mean-field proportionaltoV isincludedinthewormoperator,inorder
pert
Hamiltonian in the valence shell. The indext indicates pro- tosatisfyconditionEq.(10). Thistermwillgenerateconfigu-
tonorneutronstatesandk denotesthetime-reversedstateof rationswithpairbreakinginteractions. However,itcanoccur
9
thattoomanyofthese interactionsare generated,thoughwe Single-particleenergies(MeV)
areonlyinterestedingeneratingconfigurationswithadiffer- Orbital Protons Neutrons
entnumberof unaccompaniedparticles, butwithoutinterac- 1f -4.1205 -10.4576
7/2
tionsof the type Eq. (38). Thiscan be preventedby impos- 2p -2.0360 -8.4804
3/2
ingtheconstraintthataconfigurationcancontainatmosttwo 2p1/2 -1.2334 -7.6512
pair breakinginteractionsof this type. Observablesare only 1f5/2 -1.2159 -7.7025
updatedif there are no pair-breakinginteractionsin the con- 3s1/2 4.7316 -0.3861
figuration. This means that a number of Markov steps are 2d5/2 5.6562 0.2225
neededin orderto reach a new allowed configurationwith a 2d3/2 6.1324 0.9907
differentnumberofunaccompaniedparticles. WhengofEq. 1g9/2 6.6572 0.5631
(39) is put equalto one, the percentageof diagonalconfigu-
rations which contain no V interactions (see Eq. (38)) is TABLEIII: SingleparticleeigenstatesofaWoods-Saxonpotential,
pert
about15%.Thisisstillefficientenoughtosamplethepairing takenfromRef. [24]. Thechosenvalencespacecontains42states.
Theprotonandneutronsingleparticleenergies(inMeV)areshown
Hamiltonian. Thereareanumberofwaystoincreasetheef-
ontheright.
ficiency. Firstof allonecan changetheparameterg, hereby
influencingtheappearanceofpairbreakinginteractions. One
canalsorestrictthenumberoftimesthewormtriestoinsert -1
aV interactionbyallowingthisonlyafteracertainMarkov
pert
-2
timeinwhich’good’(i.e. withoutpairbreakinginteractions)
configurationsaresampled.Oneshouldalsokeepinmindthat
-3
whileaconfigurationcontainspairbreakinginteractions,the
V]
wormitselfisnotnecessarilyofthepairbreakingtype. Soa e -4
M
lotof Markovtime is spendto changethe configurationin a >[
P -5
global way without removing the pair breaking interactions, H
<
leadingtostrongdecorrelation.
-6
The main physical properties of nuclei in the Iron region
are modeled by a schematic mean-field plus pairing Hamil- -7
10 neutrons
tonian. For the mean-fieldpotential, we use a Woods-Saxon 11 neutrons
potential.SingleparticleenergiesaretakenfromRef. [24].A -8
0 0.5 1 1.5 2 2.5 3 3.5
full fp+sdg valence space is chosen. These single-particle
T[MeV]
statesandenergiesareshowninTableIII. Apairingstrength
G=16/56MeVisused. Duetothesizeofthemodelspace FIG. 6: Expectation value of the neutron part of the pairing-
a strength smaller than the suggested value of 20 MeV per interaction operator as a function of temperature. The pairing
nucleon is used [24]. We have tested our code by compar- strengthGnisequalto16/56MeV.Weconsider10and11neutrons
ing finite temperatureresults in a fp valence space with the inthe fp+sdg valence space (seeTabel III). At temperatures be-
onesobtainedviaanexactdiagonalizationtechnique[48].We low0.5MeVthepairingenergyismuchlowerfortheevenneutron
number.
show results of calculations with the valence shell given in
Table III occupied by 10 and 11 valence neutrons. Figure 6
showstheexpectationvalueoftheneutronpairing-interaction
operatorhHPiasafunctionoftemperature. Atlowtempera- commuteswiththeangularmomentumprojectionoperatorJz
ture,thepairingenergiesaremuchlowerfortheevennumber (butnotwithJ2),ourcurrentcodeallowsrestrictingthesim-
ofneutrons. Thiscanbeunderstoodbyremarkingthatforan ulationto configurationswith a fixedJz. Work on extending
oddnumberofneutronsthereisalwaysatleastoneunpaired thistechniquetofullJ-projectionisinprogress.
nucleon.Attemperatureshigherthan1MeV,thepairingener- When the projectionon Jz was turned on, we included an
giesdifferonlyslightly,becausethereisanincreasingnumber extra global step in order to get a good convergence at the
of unpaired nucleonsdue to thermal excitation. This is also lowesttemperatures. Thisextraglobalchangeallowsforone
reflected in the specific heat (see Figure 7). A peak appears or two unaccompaniednucleons (which block the state they
around0.8MeVduetothedevelopmentofpaircorrelations. occupy) to move to other states, and can occur whenever
Becausethewormoperatorconservesangularmomentum, the worm is diagonal. First an unaccompaniednucleon at a
one can restrict the intermediate states to a specific value of blockedstatel ischosenatrandom. A’non-blocked’pairof
the quantumnumbersJ andJ . Thisis notpossiblewith the states(k,k)isthenchosenwithprobability
z
auxiliary-fieldQMC method. In our currentimplementation
b
of the algorithm however, the occupation of each couple of P(k)=eR0(nk(t)+nk(t)−1)ekdt/Nl, (41)
time-reversedsingle-particlestates (k,k) is exactlyknownat
withn (t)theoccupationnumberofstatekatimaginarytime
alltimes. Becausetheunaccompaniedparticlenumberopera- k
t, and N a normalization factor. The subscript l indicates
tor l
that the norm is determined for a configurationcontaining a
Nu=(cid:229) n −b†b , (40) blockedstatel. TheideabehindEq. (41)istogetaprobabil-
k k k itydistributionP(k)whichispeakedaroundtheFermilevel,
k (cid:0) (cid:1)
10
10 -98
J
z
8 -99 7
-100
6
eV] -101 65
Cn 4 >[M 3,4
H -102
< 1,2
2
-103
0
-104
10 neutrons
11 neutrons 0
-2 -105
0 0.5 1 1.5 2 2.5 3 3.5 0 0.1 0.2 0.3 0.4 0.5 0.6
T[MeV] T[MeV]
FIG.7: TheneutronspecificheatCnasafunctionoftemperaturefor FIG.8: Jz projectedinternalenergiesasafunctionoftemperature.
10and11neutronsinthefull fp+sdgvalencespace(seeTableIII). ThevaluesofJz from0to7areindicatedontheleft. Aclearcon-
Calculationswereperformedat aconstant neutronpairing strength vergencetoexactzerotemperatureresultscalculatedviatechniques
Gn=16/56MeV. explainedinRef.[44]canbeseen.
butotherdistributionscanbechosenaswell. Theinterchange
of the occupationsof the blockedpair of states (l,l) and the
non-blockedpair(k,k)overthewholeimaginarytimeinterval erator in the path integral approach. This method allows to
b ,isacceptedwithprobability sampleconfigurationswithspecificsymmetriesand,inpartic-
ular,tosamplethecanonicalensemble. Itleadstoaveryef-
N
l
p=min(1, ). (42) ficientsamplingschemewithallmovesacceptedandwithout
N
k ’bounces’ or critical slowing down near second order phase
The acceptance factor for the case when the occupations of transitions. We have proven detailed balance and tested er-
twopairsofnon-blockedandblockedstatesareinterchanged, godicity. Our method opens new perspectives for the study
canbeconstructedinasimilarway. Theextrastephasahigh of quantum many-body systems where particle number and
acceptance rate, but is only necessary to enhance decorrela- other symmetries play an important role. It can be applied
tionatverylowtemperaturewhenaJz-projectionisincluded. to bosons, to fermions in absence of a sign problem and to
Athighertemperaturestheunaccompaniednucleonsmoveef- non-frustratedspinsystems atfixedmagnetization. We have
ficientlyfromstatetostateviathelastwormpieceinEq.(39). demonstratedthisbysimulatingtheBose-Hubbardmodeland
Figure 8 shows total energiesafter Jz-projectionat low tem- a nuclear pairing model. The equal-time one-body Green’s
perature. Calculationswereperformedfortenneutronsmov- function can be evaluated with high efficiency. When non-
ing in the modelspace listed in Table III. The neutronpair- equaltimeobservablesarerequired,thecurrentmethodcanin
ingstrengthisagainGn=16/56MeV.Thefigurealsoshows principlestillbecombinedwithconventionalnon-localworm
exact energy eigenvalues for Jz =0 to Jz =7. These were steps. Thereisstillalotoffreedominchoosingthealgorithm
calculatedviaatechniqueexplainedinRef. [44]. Thelowest parameters,whichcanbeusedtooptimizethealgorithm.For
Jz=1,2andthelowestJz=3,4statesaredegenerate.Forlow the Bose-Hubbardmodel we comparedthe efficiency of our
enoughvaluesofT thefinitetemperatureresultsclearlycon- algorithm(withdifferentparametersets)withadirectedloop
vergetothegroundstateswithintheconsideredensembles. SSEcode. Thoughoneshouldalwaysbecarefulwhencom-
Notethatwecancomparewithexactsolutionsbecausethe paring different algorithms, we have strong indications that
pairingstrengthwas taken constantfor all levels. Our QMC our method is very efficient. We have simulated a pairing
methodallows to solve pairingmodelswith a single particle modelforevenandoddparticlenumbers. Ourfinitetemper-
state dependentpairingstrengthGkk′, forwhichnoalgebraic atureresultsclearlysupplementalgebraicmethodsandother
solutions are available. Taking in mind the method is appli- QMCmethods.Furthermore,aprojectiononangularmomen-
cable for even and odd nucleon systems and allows angular tumsymmetriescanbeincluded. We havedemonstratedthis
momentum symmetry projections, this could greatly extend byshowingJ -projectedresults. A workonfullJ-projection
z
theapplicabilityofthepairingmodel. isinprogress.
The authorswish to thank K.Heyde, J. Dukelsky, S. Wes-
V. CONCLUSIONSANDOUTLOOK sel, M. Troyer and S. Trebst for interesting discussions and
the Fund for Scientific Research - Flanders (Belgium), the
WehavesetupaquantumMonteCarlomethodwithanon- ResearchBoardoftheUniversityofGhentandN.A.T.O.for
local loop updating scheme starting from a local worm op- financialsupport.