The sign problem in Monte Carlo simulations of frustrated quantum spin systems

13 Pages·0.25 MB·English

by  Patrik Henelius

The sign problem in Monte Carlo simulations of frustrated quantum spin systems
Patrik Henelius1 and Anders W. Sandvik2
1 National High Magnetic Field Laboratory, 1800 East Paul Dirac Dr., Tallahassee, Florida 32310
2 Department of Physics, University of Illinois at Urbana-Champaign, 1110 West Green Street, Urbana, Illinois 61801
(February 1, 2008)
We discuss the sign problem arising in Monte Carlo simulations of frustrated quantum spin sys-
0 tems. Weshowthatforaclassof“semi-frustrated”systems(Heisenbergmodelswithferromagnetic
0
couplings Jz(r) < 0 along the z-axis and antiferromagnetic couplings Jxy(r) = Jz(r) in the xy-
0 −
plane, for arbitrary distances r) the sign problem present for algorithms operating in the z-basis
2
canbesolvedwithinarecent“operator-loop”formulationofthestochasticseriesexpansionmethod

n (a cluster algorithm for sampling the diagonal matrix elements of the power series expansion of
a exp( βH) to all orders). The solution relies on identification of operator-loops which change the
J −
configuration sign when updated(“merons”) and is similar tothemeron-cluster algorithm recently

4 proposed byChandrasekharanandWieseforsolvingthesign problemforaclassoffermion models
2 (Phys. Rev. Lett. 83, 3116 (1999)). Some important expectation values, e.g., the internal energy,

can be evaluated in the subspace with no merons, where the weight function is positive definite.
]
l Calculationsofotherexpectationvaluesrequiresamplingofconfigurationswithonlyasmallnumber
e
ofmerons(typicallyzeroortwo),withanaccompanyingsignproblemwhichisnotserious. Wealso
-
r discuss problems which arise in applying the meron concept to more general quantum spin models
t
s with frustrated interactions.
.
t
a PACS numbers: 75.40.Mg, 75.10.Jm, 02.70.Lq
m
- I. INTRODUCTION tum problem into a form resembling a classical statisti-
d calmechanics problem. There are two classes of systems
n
Recently, there have been several significant develop- forwhichthisissueisparticularlypressing—interacting
o
c mentsofmoreefficientMonteCarlomethodsforinteract- fermions in more than one dimension and quantum spin
[ ing quantum many-body systems.1 The Trotter decom- systems with frustrated interactions (in any number of

positionformula2,3 hastraditionallybeenusedasastart- dimensions). For fermions in one dimension, and hop-
1
v ing point for finite-temperature simulation algorithms, ping between nearest-neighbor sites only, the sign prob-
1 such as the worldline3 and fermion determinant4 meth- lemcanbeavoidedbecausethefermionanticommutation
5 ods. Itintroducesasystematicerrorthatcanberemoved relations do not come into play (other than introducing
3 only by carrying out simulations for several different ahard-coreconstraint)intheone-dimensionalreal-space
1 imaginary time discretizations ∆τ and subsequently ex- path integral. In two or more dimensions (or even in
0 trapolating to ∆τ =0. Such extrapolations are not nec- one dimension if hopping further than between nearest
0 essary with the stochastic series expansion method,5–7 neighbors is included), permutations of fermions during
0
/ whichisbasedonsamplingthepowerseriesexpansionof the propagation in imaginary time lead to a mixed-sign
at exp( βH) to all orders and is related to a less general path integral which typically cannot be efficiently evalu-
m meth−od proposed much earlier by Handscomb.8–10 Re- ated using Monte Carlo methods. The sign problem can
sultsthatareexacttowithinstatisticalerrorscanalsobe be avoided with the fermion determinant algorithm in
-
d directly obtained with recent worldline11–13 and fermion special cases, such as the half-filled Hubbard model (be-
n determinant14algorithmsformulatedincontinuousimag- causeofparticle-holesymmetry),4butinothercasessim-
o inary time. Even more significant than the elimination ulationsarerestrictedtohightemperaturesand/orsmall
c of the Trotter decompositionerror are generalizationsto system sizes.18 For frustrated spin systems the source of
:
v thequantumcase15,16,12,7ofclusteralgorithmsdeveloped the sign problem is different. A minus sign appears for
Xi for classical Monte Carlo.17 These “loop algorithms” (so every event in the path integral in which two antiferro-
calledbecause the clusters areloops on a space-time lat- magnetically interacting spins are flipped.19 This causes
r
a tice) can reduce the autocorrelation times by orders of an over-allminus sign if the total number of spin flips is
magnitudeandenablehighlyaccuratestudiesofsystems odd, which can be the case, e.g., for a triangular lattice
inparameterregimeswherepreviousalgorithmsencoun- or a square lattice with both nearest- and next-nearest-
tered difficulties due to long autocorrelation and equili- neighbor interactions. Simulations of quantum spin sys-
bration times. tems are therefore restricted to models with no frustra-
In spite of these developments, the class of mod- tion(intheoff-diagonalpartoftheHamiltonian),suchas
els which can be studied using quantum Monte Carlo ferromagnets, or antiferromagnets on bipartite lattices.
methods is still severely restricted due to the “sign A promising approach to solving the sign problem
problem”,18,19 i.e., the non-positive-definiteness of the wasrecentlysuggestedbyChandrasekharanandWiese.20
weight function that can arise in transforming a quan- They considered a system of spinless fermions on a two-
1
dimensional square lattice within the context of the forthismodelhaveaseveresignproblemwhenusingthe
worldline loop algorithm.15 They showed that, for this z-axisasthequantizationaxis,however,itcanbeavoided
particular model and for a certain range of nearest- by a simple rotation to the x-representation. Using the
neighbor repulsion strengths, the properties of the loops SSE algorithm and the meron concept, the sign prob-
can be used to eliminate the sign problem. Flipping a lemcanbe eliminatedalsointhe z representation. With
loop can change the number of fermion permutations both representations accessible in simulations, correla-
from odd to even, or vice versa, thereby also changing tion functions both parallel and perpendicular to the z
theover-allsignoftheconfiguration. Suchsign-changing direction can be easily evaluated.
loops are called “merons”. The magnitude of the weight The model, Eq. (1), can be mapped onto a hard-core
is not affected by flipping a meron and therefore all con- boson model with attractive interactions and frustrated
figurations with one or more merons cancel in the par- hopping. Frustrationinthepotentialenergyhasbeenin-
tition function. The subspace of zero merons is posi- vestigatedinthiscontextasapossiblemechanismtoren-
tivedefiniteandcanbesampledwithoutasignproblem. der a disordered bosonic ground state.21 Frustration in
Typical operator expectation values of interest also con- the hopping [the xy-term in Eq.(1)] should decrease the
tain contributions from configurations with two merons tendency to forming off-diagonal long-range order and
which therefore also have to be included in the simu- couldthenleadto a normalfluid (non-superfluid). How-
lation and introduce a “mild” sign problem. The rela- ever, the highly symmetric case considered here has a
tive weights of the zero- and two-meron subspaces to be trivial, ordered ground state; the fully polarized ferro-
sampled canfurther be chosenin an optimum way using magnetic state (corresponding to a completely filled lat-
a reweighting technique, which further reduces the sign tice of hard-core bosons; a trivial case of normal solid).
problem. Effectsoffrustrationonlycomeintoplayatfinitetemper-
Inthispaperweexploreananalogousmethodforsolv- ature, where the model is different from the correspond-
ingthesignproblemforfrustratedquantumspinmodels. ing isotropic ferromagnet (on non frustrated, bipartite
Weconsidertheoperator-loopformulationofthestochas- lattices the two models are equivalent, since the sign of
tic series expansionmethod,7 in which sequences of two- the xy-termcanbeswitchedbyaspinrotationononeof
spinoperatorsaresampledbyformingclusters(loops)of the sublattices).
operators that can be simultaneously updated without Although we have not been able to solve the sign
changes in the weight function. The updated clusters problem for other cases, such as the Heisenberg model
contain operators acting on the same spins, but diago- with completely antiferromagnetic interactions [J (1) =
z
nal operators may be changed to off-diagonal ones and J (1)>0 and J (√2)=J (√2)>0], our work never-
xy z xy
vice versa. For a model with frustrated interactions an thelessgivessomefurtherinsightsintothemeronconcept
operator-loopupdate canlead to a signchange. In anal- and what is required in order to solve the sign problem
ogy with Chandrasekharan and Wiese20 we will refer to for arbitrary couplings.
suchsign-changingloopsas“merons”. Thesignproblem Theoutlineoftherestofthepaperisthefollowing: In
can be solvedif the operator-loopsfor a givenconfigura- Sec.IIwereviewthebasicsofthestochasticseriesexpan-
tion can be uniquely defined and the weight function is sionmethodanddiscussoperator-loopupdatingschemes
positivedefiniteintheconfigurationsubspacecontaining for both ferromagnetic and antiferromagnetic couplings.
nomerons. Unfortunately,wefindthatthesecriteriaare InSec.IIIwepresentthesolutionofthesignproblemfor
in general difficult to satisfy. Operator-loop algorithms the J (r) = J (r) model. The reweighting technique
z xy
−
withuniquelydeterminedloopsaretypicallynon-ergodic isanalyzedinsomedetailinSec.IV.InSec.Vwediscuss
for frustrated systems, and with supplemental local up- somesimulationresultsforthesemifrustratedmodeland
dates for ergodicity there are mixed signs in the zero- make comparisons with the isotropic Heisenberg ferro-
meron subspace. In fact, in such cases merons typically magnet. We summarize our work in Sec. VI.
do not even exist, i.e., none of the operator-loops can
change the sign when flipped. We have found only one
spinsystemforwhichthesignproblemcanbeeliminated II. OPERATOR-LOOP ALGORITHM
using merons — the Heisenberg model with ferromag-
netic couplingsJ (r)<0 alongthe z-axisandfrustrated
z Inthissectionwefirstbrieflyreviewtheexpansionun-
antiferromagneticcouplingsJ (r)= J (r)intheplane
xy z derlying the SSE method and then discuss the operator-
−
perpendicular to this axis, i.e., the Hamiltonian
loop updates used to efficiently sample the expansion.
We here assume a non-frustrated case and postpone the
H =− Jij[SizSjz − 21(Si+Sj−+Si−Sj+)], (1) discussion of the sign problem for frustrated models to
i,j Sec.III.FordefinitenessweconsidertheS =1/2Heisen-
X
berg model
where J > 0 and the range of the couplings is arbi-
ij
trary. We have implemented a meron algorithm for this
H = J S S , (2)
modelonasquarelatticewithnearest-andnext-nearest- ± i· j
neighborcouplingsJ(1)andJ(√2). Standardalgorithms hXi,ji
2
p
where ij denotes a pair of nearest-neighbor spins on
h i α(p) = H α , (9)
a cubic lattice (in an arbitrary number of dimensions), | i ai,bi| i
and J > 0. Depending on the sign, the model is an iY=1
antiferromagnet( ) or a ferromagnet(+). To construct
wherefor anallowedconfiguration α(0) = α(L) = α
−
theSSEconfigurationspacetheHamiltonianisrewritten | i | i | i
and the weight function corresponding to (7) is given by
as a sum of diagonal and off-diagonal operators
(Jβ)n(L n)!
J M W(α,SL)= 2nL!− . (10)
H = (H H )+C, (3)
1,b 2,b
−2 ∓
b=1 Having established the framework we will proceed
X
to describe the procedures for importance sampling
wheretheindexbdenotesaninteractingspinpair(bond)
of the terms (α,S ) according to the weight (10).
i(b),j(b) andC isanirrelevantconstantequaltoMJ/4, L
h i The initial state can be a sequence of the form
whereM isthetotalnumberofpairsofinteractingspins.
(0,0) ,(0,0) ,...,(0,0) (subscripts on the index pairs
The bond-indexed operators are given by 1 2 L
will sometimes be used to denote the position in S )
L
and a random state α . An ergodic procedure for sam-
H =2 1 Sz Sz (4) | i
1,b 4 ∓ i(b) j(b) pling the terms is achieved using two types of basic up-
H =S+(cid:16) S− +S− (cid:17)S+ . (5) dates; a simple substitution of single diagonal operators
2,b i(b) j(b) i(b) j(b) and the operator-loop update which involves simultane-
ous updates of a number (in principle varying between 1
Notethattheeigenvaluesofboththediagonal(H )and
1,b
and n) of diagonal and off-diagonal operators.
the off-diagonal (H ) operators are 0 and 1, both for
2,b
The diagonal update is carried out by traversing the
the antiferromagnet and the ferromagnet. The partition
operator-index sequence S from beginning (p = 1)
function Z =Tr exp( βH) is expanded as5 L
{ − } to end (p = L). Operator substitutions of the form
Z = ∞ (−β)n αHn α , (6) s(0to,r0e)dp s↔tat(e1,bα)piasruepadtatteemdpetveedrywthiemreeapnososffib-ldei,awgohnilaeltohpe-
n! h | | i | i
α n=0 eratoris encounteredso that the state α(p) is available
XX | i
whenneeded. Withtheweightfunction(10)detailedbal-
in the basis {|αi} = {|S1z,S2z,...,SzNi}, where N is the ance can be seen to be satisfied if the acceptance proba-
number of spins. Terms of order greater than n Nβ bilities are taken to be
∼
give an exponentially vanishing contribution and for the
purposeofastochasticsamplingtheexpansioncanthere- Mβ αb(p)H1,b αb(p)
P [(0,0) (1,b) ] = h | | i, (11)
p p
fore be truncated at some n = L of this order without → L n
−
loss of accuracy (see, e.g., Ref. 6 for details on how to L n+1
P [(1,b) (0,0) ] = − . (12)
choose a sufficiently high truncation power). Additional p → p Mβ α (p)H α (p)
b 1,b b
unitoperatorsH 1areintroducedtorewriteEq.(6) h | | i
0,0
≡
as Note that the diagonal update changes the expansion
power n by 1. Off-diagonal operators cannot be in-
( 1)n2(Jβ)n(L n)! L troducedone-±by-onebecauseoftheperiodicitycondition
Z = ∓ − α H α ,
Xα XSL 2nL! * (cid:12)(cid:12)iY=1 ai,bi(cid:12)(cid:12) + |cαa(nLb)ie=use|αd(f0o)ri.thLisocpaulrpuposdea,6tebsuitnvaorlevminogrtewcoomoppelircaattoerds
(cid:12) (cid:12)
(cid:12) (cid:12) (7) and far less efficient than the operator-loop update,7
(cid:12) (cid:12)
which is discussed next.
where SL denotes a sequence of operator-indices: We use a largely pictorial description of the operator
loops. First we consider the antiferromagnet. Note that
SL =(a1,b1)1,(a2,b2)2,...,(aL,bL)L, (8) inthiscasetheonlynon-zeromatrixelementsofthebond
operators are
with a 1,2 and b 1,...,M , or (a ,b )= (0,0).
i i i i
∈{ } ∈{ }
The number of non-(0,0) elements in SL is denoted n, H = H =1,
1 1
while n denotes the number of off-diagonal operators h↓↑| |↓↑i h↑↓| |↑↓i
2 H = H =1, (13)
in the sequence. Note that since the expectation value h↓↑| 2|↑↓i h↑↓| 2|↓↑i
in Eq. (7) is always equal to zero or one, the sign of a
i.e., they can act only on antiparallel spins. An example
term is negative only if n is odd. This sign problem
2 ofatermintheexpansionforafour-siteantiferromagnet
occurs(only)whenfrustrationispresentandisthemain
isdepictedinFig.(1). Thisrepresentationmakesevident
topic of this paper. However, for the discussion of the
thecloserelationshipbetweentheSSEexpansionandthe
samplingprocedures,inthissectionweassumeapositive
Euclideanpathintegral. Animaginarytimeseparationτ
definite expansion. We introduce the notation α(p) for
| i correspondstoadistributionofpropagations∆pbetween
a propagated state states, centered around ∆p = (τ/β)n.5,13 We will for
3
|α(0)>
H
1,1
|α(1)>
H
2,3
|α(2)>
H
1,1
|α(3)>
H
2,3
|α(4)>
H
1,2
|α(0)>
FIG.1. Representation of a term in the SSE expansion of
a four-site antiferromagnet. Up and down spins are repre-
FIG. 3. The SSE space-time configuration of Fig. 1
sentedassolid andopentriangles, respectively. Thehorizon-
uniquelydividedupintoloops(left). Therighthandconfigu-
tal bars indicate the presence of diagonal (narrow lines) and
ration results from flipping the loop indicated by the dashed
off-diagonal (thick lines) operators.
line. Note the periodic boundary conditions in the vertical
(“time”) direction.
convenience here refer to the propagation as the time
dimension.
only one cluster; see Fig. (3), where our example con-
The general idea15 behind the loop update is to flip a
figuration has been divided up into three clusters. All
clusterofspinsintheconfigurationinsuchawaythatthe
loops can be flipped independently with probability 1/2
weight is not changed. With the SSE method there will
— in Fig. (3) we depict the result of flipping the largest
also have to be changes made to the operatorsacting on
loop of the example configuration. A full operator-loop
thespins,sinceotherwiseoperatorsH orH mayacton
1 2 update amounts to dividing up the configuration into
parallel spins, resulting in zero-valued matrix elements.
all of its loops which are flipped with probability 1/2.
Since one of the states α(p) and the operator sequence
| i The (random) decision of whether or not to flip a loop
S uniquelydefinethewholespinconfiguration,theSSE
L can be made before the loop is constructed, so that each
loops are in practice treated as loops of operators, the
loop has to be traversed only once. Spin “lines” Sz(p),
exact meaning of which will be made clear below. i
p=0,...,n, which are not acted upon by any of the op-
eratorsinS willnotbeincludedinanyoftheoperator-
L
loops. They correspond to “free” spins which can be
flippedwithprobability1/2. Suchalinecanalsobecon-
sideredaloop,andthenit willalwaysbe truethatevery
FIG.2. Spinflipandaccompanyingoperatorexchangedur- spin Sz(p) belongs to one loop. Free spins appear fre-
i
ing a loop update for an antiferromagnet. The dashed line quently at high temperatures, when the total number of
indicates part of the loop. operators n<N, but are rare at low temperatures.
Note that∼the spin states at the four legs of the
ConsideroneoftheoperatorsH inFig.(1). Itcanbe operator-verticescompletely determine the full spincon-
1,1
depictedasa“vertex”withfourlegsassociatedwithspin figuration, except for free spins that happen not to be
states or . If we flip the upper left spin, a vanishing acted upon by any of the operators in S . Hence, the
L
↓ ↑
matrix element results. But if we flip both upper (or operator-loop update can be carried out using only a
lower) spins and also change the operator type to off- linked list of of operators, i.e., an array of vertices with
diagonal H , an allowed matrix element is generated; four spin states and associated pointers to the “previ-
2,1
seeFig.(2). Usingthisideawecanformaclusterofspins ous” and “next” vertex associated with the same spin.
bychoosingarandomspinSz(p)intheconfigurationand The storage requirements and the number of operations
i
traversing up or down until we encounter an operator needed for carrying out a full operator-loopupdate then
(bond) acting on that spin, then switch to the second scale as Nβ instead of N2β if the full spin configu-
∼ ∼
spin of the bond and change the direction of traversing ration were to be used.
the list. Eventually we will necessarily arrive back at In a simulation we first make a full cycle of diagonal
the initial starting point, whereupon a closed loop has updates in the sequence S and then create the linked
L
formed. All the spins on this loop can be flipped if the list of vertices in which the operator-loop updates are
operators encountered are also switched [(1,b) (2,b)]. carried out. The vertex list is then mapped back onto
↔
Note that the same operator can be encountered twice, the sequence S and the state α(0) . Alternatively, the
L
| i
which results in no net change of operator type (but the linked list can be updated simultaneously with each di-
spins at all four vertex legs are flipped). agonaloperator substitution, so that it does not have to
The whole configuration can be uniquely divided up bere-createdforeachMonteCarlostep—dependingon
into a set of loops, so that each spin belongs to one and the model studied there may be significantdifferences in
4
execution time between the two approaches. linked list of vertices must to be bi-directional, but for
theferromagnetitissufficienttokeepasingly-directional
list.
The diagonal and operator-loop updates satisfy de-
tailed balance and the combination of them leads to er-
godic sampling for a ferromagnet on any lattice, and for
FIG.4. Spinflipandaccompanyingoperatorexchangedur- antiferromagnets on bipartite lattices — for frustrated
ingaloopupdateforaferromagnet. Thedashedlineindicates antiferromagnets there are complications, in addition to
part of theloop. the sign problem, which will be discussed further in the
next section. The operator-loop sampling is highly ef-
ficient, with integrated auto-correlation times typically
less than one Monte Carlo step.
Theloopconstructiondescribedherereliesontherota-
tionalinvarianceofthemodels,i.e.,thefactthatboththe
diagonal and off-diagonal matrix elements in Eqs. (13)
and (14)are equal to one. For a generalanisotropiccase
the loops will lead to weight changes when flipped and
must then be assigned “a-posteriori” acceptance prob-
abilities which typically become small for large lattices
at low temperatures. Other types of loops avoiding this
problem have been proposed,7 but will not be discussed
here.
FIG. 5. SSE space-time configurations uniquely divided
III. THE SIGN PROBLEM
into loops for the case of a ferromagnet. The left and right
handconfigurations differbyflippingtheloop indicated bya
dashed line. The notorious sign problem arises in stochastic sam-
plingwhenthe functionusedtoweightthe differentcon-
For the ferromagnet we can construct the loops in a figurations in not positive definite. A typical quantity
similar manner, but the non-zero matrix elements are thatcanbecalculatedbyMonteCarlo(importancesam-
now pling) is of the form
h↓↓|H1|↓↓i=h↑↑|H1|↑↑i=1, A = iA(xi)W(xi) = A(x) W, (15)
H2 = H2 =1, (14) h i P iW(xi) h i
h↓↑| |↑↓i h↑↓| |↓↑i
where W is the wPeight function and A the measured
i.e., the off-diagonal operators act only on antiparal-
quantity, which both depend on the general coordinate
lel spins, as before, whereas the diagonal ones can act
x of the configuration space sampled. When the coordi-
only on parallel configurations. This implies qualita-
nates are sampled according to relative weight, the de-
tive changes in the structure of the loops, as depicted
sired quantity is simply the arithmetic average of the
in Fig. (4). If we again consider an operator H and flip
1 measured function A(x), as indicated by the notation
theupperleftspin,wenotethatweneedtoflipthelower
A(x) above. If the weight function is not positive
rightspinandchangetheoperatortoH2. Hence,instead dhefinitieW, the sampling can be done using the absolute
of changing the direction of traversing the configuration
value of the weight, and the expectation value can be
every time an operator is encountered we now continue
calculated according to
in the same direction after switching to the second spin
of the bond. Any configuration can still be uniquely di- A(x)s(x)
|W|
A = h i , (16)
vided up into loops that can be flipped with probability
h i s(x)
|W|
1/2. An example of a ferromagnetic configuration with h i
its loops is shown in Fig. (5). where s(x) equals 1, depending on whether the sign
±
Note that since the loops for the ferromagnet never of the weight function is positive or negative. For most
change direction as they go through the lattice, every models,whereasignproblemispresent,theaveragesign
singleloophastotraverseeachstate α(p) atleastonce. s(x) decreases exponentially to zero as a function
|W|
| i h i
It follows that the number of sites N is an upper bound of inverse temperature and system size, and the relative
of the number of loops. The antiferromagnetic loops, on statistical errors of calculated quantities increase expo-
theotherhand,traversethelatticeinbothdirectionsand nentially.
the number of loops is therefore instead limited by the A meron-cluster solution to the sign problem using
total number of operators n Nβ. As a consequence loopupdatesoffermionworld-lineconfigurationswasre-
of the change of direction, fo∼r the antiferromagnet the cently proposed by Chandrasekharanand Wiese.20 This
5
approach is based on the idea that if it is possible to where s(δ)= s(δ¯), where δ¯denotes a flipped loop, and
±
map every configuration with negative weight uniquely the sign is negative for merons and positive otherwise.
to a correspondingconfigurationwith equalweightmag- Since flipping any meron leads to two terms that cancel,
nitude butopposite sign,thenthe partitionfunction can it follows that
be sampled without a sign problem simply by not in-
cluding any configuration which is a member of such a 1 2NL
s(x )= δ , (20)
cancelingpair. Inthemeron-clusteralgorithm,flippinga 2NL l ± nM,0
loopofspinscanleadtoasignchangewithoutchangein Xl=1
the magnitude of the weight, and such a “meron” hence
wheren denotes the totalnumber ofmerons. The sign
M
identifies a canceling pair ofconfigurations. Here we will
infrontofthedeltafunctionisthe“inherent”signofthe
presenta similar approachwithin the SSE operator-loop
configuration, independent of the loop orientation when
method for frustrated quantum spins.
therearenomeronspresent. Thissignhastobe positive
Let us consider the Heisenberg antiferromagnet dis-
for the meron solution to be applicable in practice, and
cussed in the previous section. From Eq. (7) we see that
thenthepartitionfunctioncanbesampledinthepositive
aconfigurationhasnegativeweightifthetotalnumbern
2 definite subspace of configurations with no merons.
ofoff-diagonaloperatorsisodd. Thiscanonlybethecase
Having found an expression for the denominator in
on frustrated lattices. As described in the previous sec-
Eq. (17) we need to consider the numerator for cases of
tion, any SSE configuration can be divided up uniquely
interest. SSE estimators for a number of important op-
into a number of loops. Flipping a loop interchangesthe
erators have been discussed, e.g., in Ref. 6. The internal
diagonal and off-diagonal operators, but leaves the total
energy is given by
number of operators unchanged. It follows that the sign
will change if and only if a loop passing through an odd 1
E = n , (21)
number of operators is flipped (two passes through the −βh iW
same operator is counted as two operators). Since the
total weight remains unchanged we have thus found the where n denotes the total number of operators in the
desired mapping between positive and negative config- sequence SL. This number is not affected by the loop
urations (assuming that there exist loops which change updates, and hence it follows that
the sign when flipped, which in fact is not always the
case). Inanalogywithpreviousworkwecallsuchasign- E = 1hhn(x)s(x)ii|W| = 1hn(x)δnM,0i|W|. (22)
flipping loop a “meron”. Let us now see how we can use −β hhs(x)ii|W| −β hδnM,0i|W|
this concept to measure observables.
As in Ref. 20 we consider improved estimators that Assumingthatthissectorhaspositivedefiniteweightwe
average over all loop configurations. Denote the number havethereforecompletelyeliminatedthesignproblemby
ofloopsinthe systemN . Since eachloopcanbe inone restricting the simulation to the zero-meron sector. The
L
of two states there is a total of 2NL configurations that energy is then simply given by
can be reached by flipping all combinations of the loops
n(x0)
present. The improvedestimate thereforetakesthe form E = h iW, (23)
− β
A(x)s(x)
|W|
A = hh ii , (17)
where the superscript 0 indicates the restriction of the
h i s(x)
|W|
hh ii simulation to the zero-meron sector.
where the double bracketsdenote anaverageoverallthe Next we will consider the magnetic susceptibility,
loop states for each generated SSE configuration, e.g.,
2
s(x) = 1 2NLs(x ) . (18) χ=β* Siz! +/N =βhM2i/N. (24)
hh ii *2NL l + Xi
l=1
X
Since M is conserved by the Hamiltonian its value is
The generalcoordinate x here refers to the SSE configu-
the same in all propagated states; M(p) = M(0) M
ration space (α,S ) and x refers to one out of the 2NL ≡
L l (p = 1,...,n). In a configuration uniquely divided up
possible outcomes of “flipping” a number of loops. into loops, every spin Sz(p) belongs to one and only one
i
Let us consider this average. Denote the state of a
loop, if we count as loops also all “lines” of free spins,
loop with δ, with two possible “orientations” δ , . i.e., the spins Sz(p), p=1,...,n for all sites i which are
∈ {↑ ↓} i
Sinceflippingoneloopdoesnotaffectanyotherloops(in
not associated with any operator in the sequence (and
terms of their paths taken), the sign of a configuration
therefore can be flipped). It follows that the change in
factors according to
M(p) when flipping a loop must be the same for all p,
andhenceonlyloopsthatgothroughallstates α(p) (at
NL | i
s(δ ,δ ,...,δ )= s(δ ), (19) least once) can change M when flipped. We can there-
1 2 NL i
fore introduce a loopmagnetizationm , whichis simply
i=1 L
Y
6
equal to the sum of the spins traversed by the loop for evento odd, or vice versa,in loop updates). These simi-
an arbitrary α(p) . In the estimator (17) corresponding larities are not surprising,considering the close relation-
| i
to the susceptibility the numerator can hence be written ship between the SSE expansionandthe Euclideanpath
as integral.13
2NL
1
(m (x )+...+m (x ))2s(x ) , (25)
*2NL 1 l NL l l +
l=1
X
The magnetization on a loop always changes sign when
aloopisflipped; the over-allsigns(x )onlychangessign
l
when a meron is flipped. Therefore, in summing over all
loops in (25), a non-zero value results only if the config-
uration has zero or two merons. The full susceptibility
estimator therefore takes the form
FIG. 6. Three off-diagonal operators acting on a triangle
2NL m 2δ +2m m δ of spins. The left and right hand configurations differ by
l=1 | l| nM,0 | M1|| M2| nM,2 flippingtheloopindicatedbyadashedline. Notetheperiodic
χ= , (26)
DP hδnM,0i E boundary conditions.
whereM1andM2aretheindicesoftheloopscorrespond- Nowconsiderthe applicationofthe aboveideas to the
ing to merons in a two-meron configuration. Hence, un- Heisenberg model on a square lattice with nearest- and
like in the case of the total energy, the sign problem has next-nearest-neighborcouplingsJ(1)>0andJ(√2)>0
here not been completely eliminated, since the zero- and (antiferromagnetic). Thismodelhasasignproblemsince
two-meronconfigurationcontribute1and0,respectively, the total number n of spin-flipping operators in a con-
2
totheaveragesign. Whenthe SSEconfigurationvolume figuration can be odd, e.g., three operators on a triangle
growstherelativeweightofthezero-meronsectorshould ofspins,as showninFig.(6). Already this simple exam-
diminish, leading to a decreasing average sign. Chan- pleillustratesthattheoperator-loopalgorithmdiscussed
drasekharanandWiese statedthatthe statisticalfluctu- above does not sample the full configuration space and
ation in the improved susceptibility estimator increases that the meron concept therefore cannot be applied to
quadraticallywithNβ,i.e.,muchslowerthantheconven- solving the sign problem. Since all three bonds are anti-
tional exponential increase. They also pointed out that ferromagnetic,aloopwillchangedirectioneverytimean
this remaining sign problem can be solved by reweight- operator is encountered. In order for the loop to close,
ingthezero-andtwo-meronsectorswithexternalweight it therefore has to pass through an even number of op-
factors w(0) and w(2). This changes the above formula erators and hence flipping the loop cannot change the
to numbern ofoff-diagonaloperatorsfromoddtoeven,or
2
2NL m 2δ w(2)+2m m δ w(0) vice versa. This is illustrated in Fig. (6), where the only
l=1 | l| nM,0 | M1|| M2| nm,2 effect of flipping the single loop in the system is to flip
χ=
DP hδnM,0w(2)i E all the spins; the operators remain unchanged. Hence,
merons do not even exist within the operator-loop algo-
(27)
rithmforthismodel,andthesamplingisnon-ergodic. A
localupdatecaninprinciplebeusedincombinationwith
In the next section we will say more about reweighting.
theoperator-loopsinordertomakethesamplingergodic.
As we have shown above, the meron concept within
However, a configuration can then have a negative sign
theSSEmethodformallyleadstoexactlythesameequa-
(ofwhichFig.(6)isanexample)eventhoughthereareno
tions as in the world-line simulations of fermion systems
merons present. The principal requirement of positive-
consideredinRef.20. Thedifferenceisonlyinthestruc-
definiteness in the zero-meron subspace (which in this
tureofthemeronitself; thefermionicmeronchangesthe
case is the full space) is hence not fulfilled. A similar
number of particle permutations from even to odd, or
problemseems toaffectallmodels withfrustrationinall
vice versa, whereas the SSE meron in the case of a frus-
of the spin components. It is possible that some other
trated spin system instead changes the number of anti-
wayof constructing the loops couldremedy this, e.g.,by
ferromagnetic spin flips from even to odd or vice versa.
switching to some other, non-trivial basis in which the
Applying the SSE operator-loop algorithm to a fermion
SSE expansion could be carried out. Other ways pro-
system would lead to merons of exactly the same kind
posed for constructing loops in the standard basis con-
as those existing within the world-line framework, and,
sidered here do not uniquely define the set of loops7 and
conversely,applyingaworld-lineloopalgorithmtoafrus-
cantherefore noteasilybe usedwith the meronconcept.
tratedspinsystemshouldleadtomeronssimilartothose
We have found one class of spin models for which the
discussed here (there are no diagonal operators in the
meronideas can be successfully applied to solve the sign
world-lineconfigurations,butspinflipeventscorrespond
problem: Heisenbergmodelswhichareantiferromagnetic
to the SSE off-diagonal operators and can change from
7
in the xy-plane but ferromagnetic along the z-axis., i.e., tating the ferromagnetic component to the y-direction.
the“semi-frustrated”model(1),whichcanbewrittenas The Hamiltonian then takes the form
M M
H =− J2b (H1,b−H2,b)+C, (28) H =−J2 H1′,b+H2′,b +C, (31)
Xb=1 Xb=1(cid:0) (cid:1)
where the bond-indexed operators are given by
where the bond-indexed operators are given by
H =2 1 +Sz Sz , (29)
1,b 4 i(b) j(b) H′ =2 1 Sz Sz , (32)
H =S(cid:16)+ S− +S− (cid:17)S+ . (30) 1,b 4 − i(b) j(b)
2,b i(b) j(b) i(b) j(b) H′ =S+(cid:16) S+ +S− (cid:17)S− , (33)
2,b i(b) j(b) i(b) j(b)
Onanon-frustratedlatticethismodelisequivalenttoan
isotropicHeisenbergferromagnet,sincen isalwayseven and the fundamental spin-flips and operator exchange
2
and the sign in front of the operators H in Eq. (28) is during a loop update is shown in Fig. (7). Being able to
2,b
irrelevant as ( 1)n2 = 1 in Eq. (7) — the sign can also workinbothbaseswecaneasilymeasureallcomponents
be transforme−d away by a spin rotation on one of the of the susceptibility. Our main motivation for studying
two sublattices. On a frustrated lattice, on the other this model is to illustrate how the sign problem can be
hand, n can be odd (the lattice is no longer bipartite removed in the z-basis. Nevertheless, we will also show
2
so thatthe transformationmentionedabovedoes notre- some results calculated in the x-basis.
move all signs), and the system is no longer equivalent
to the isotropic ferromagnet. The model has a classi-
IV. REWEIGHTING
cal two-fold degenerate ferromagnetic ground state, but
at finite temperatures the transverse spin components
are frustrated, and the behavior will be different from Animportanttechnicalaspectofthemeron-solutionis
the isotropic ferromagnet. When simulated in the z- thereweightingofthezero-andtwo-meronsectors,which
basisusingstandardalgorithmsthesemifrustratedmodel was briefly mentioned in the previous section. Eq. (27)
has a severe sign problem, but the zero-meron sector is gives the correct estimator for the susceptibility after
positive-definite and the SSE meron solution can be ap- reweighting, but it gives no information on how to do
plied. The structure ofthe operator-loopsis the same as the reweighting in practice. How to determine the op-
for the ferromagnet and the loop algorithm is therefore timal reweighting and whether reweighting changes the
ergodic. scalingoftherelativeerrorareimportantquestionstobe
Themeronsolutioncanbe implementedinseveraldif- considered in this section.
ferent ways and we briefly describe how it was done in
this work: During the sequential diagonal updates the
1
linkedlist,representingtheloopstructure,isupdatedsi-
multaneously with each accepted diagonal update. The
loopsarenumberedandinformationisstoredonwhether δn,0
eachloopisameronornot. Duringanattempteddiago-
nal update only the loops directly affected by the opera-
0
tor substitutionareupdated. This permits easyandfast
1
checkingof whether the number of meronsin the system
has changed or not. If the new number of merons is dif-
ferentformzeroortwothe update is rejected,whereasif
the number of merons changes from i to j it is accepted δn,0
with probability w(j)/w(i), where i,j 0,2 , and w(i)
∈{ }
is the re-weighting factor assigned to meron sector i. 0
0 50 100 150 200
MC steps
FIG. 8. Fluctuations of n during a random process with
two outcomes (n = 0,2). The upper graph shows results for
W0 =0.5 and W2 =0.5, while in the lower graph W0 =0.01
and W2 = 0.99 and a reweighting factor W = 99 as used
duFriInGg. 7a. lSopopin uflpipdaatnedfoarccommopdaenlyiwngithopaerafteorrromexacghnaentgice (resulting in W0′ =W2′ =0.5).
x-component, but antiferromagnetic y and z-components.
The dashed line indicates part of theloop. As a first, simple example of how reweighting affects
the statistics of a simulation we will discuss a simple
Note that the sign problem for the semifrustrated random process. Consider a random variable n, which
model can also be very easily transformed away by ro- can take two different values, 0 or 2. Let W0 and W2
8
designate the probability for these two outcomes. The 1 N
x¯= x (41)
expected fractions and standard deviations of these out- i
N
comes from N random selections are given by Xi=1
and
√W W
0 2
δ =W (34)
n,0 0
h i ± √N x¯2 x¯2
δ =W √W0W2. (35) σx¯ =s N− . (42)
n,2 2
h i ± √N
When studying the behavior of the standard deviation
Weconsideranexpectationvalueofaformsimilartothe itself, we also want to obtain an estimate of its accu-
QMC susceptibility, Eq. (26); racy. This can be done by dividing the N bins into M
sets containing N/M bins each. For each set a standard
δ +δ 1 1 W 1
n,0 n,2 2 deviation σ can be calculated according to
f = h i = , (36) x
h i hδn,0i W0 ± W0rW0√N
σ = x¯2 x¯2, (43)
x
with a relative standard deviation −
p
wherethebardenotesanaverageoftheN/M binswithin
σ W 1
f = 2 . (37) the set. The final standard deviation and its statistical
f rW0√N fluctuation are then given by
This formula becomes valid for large N, when the stan-
M
1
dard deviation is small. As W decreases the standard
0 σ¯ = σ (44)
deviation increases, but one can reweight the two out- x M xi
i=1
comes by assigning an additional weight W to the n=0 X
outcome such that a transition from n = 0 to n = 2 is and
accepted with probability 1/W, while a transition from
n=2ton=0isalwaysaccepted. Aftersuchareweight- σ¯2 σ¯2
σ = x− x . (45)
ing the probabilities of obtaining n = 0 and n = 2 are σ¯x s M
given by
Eq. (44) represents the standard deviation of the dis-
W W
W′ = 0 (38) tribution of the binned values x, and not the standard
0 W0W +W2 deviation of an average of these. It does not decrease
W′ = W2 , (39) as the number of measurements N is increased, but it is
2 W0W +W2 dependent onthe number ofMC steps ineachbin, Nbin,
and will decrease as 1/√N . Hence it is important to
bin
and f is given by
state the number of MC steps in the bins for which the
deviationis calculated. The statisticalerrorofthis stan-
δ +Wδ
hfi= h n,0 δ n,2i. (40) dard deviation, σσ¯x, will, on the other hand, decrease as
h n,0i 1/√N.
When calculating the standard deviation for this case In this manner we can calculate the standard devia-
we have to be careful since the reweighting introduces tion for f. In order to show the standard deviation as a
correlationsinto the system. This is clearly visualizedin function of the reweighted probability W0′, Eq. (39), can
Fig. 8, where in the upper graph a series of independent be inverted to express the necessary weight factor that
outcomes with equal probability (W0 = W2 = 0.5) are causes the average to change from W0 to W0′;
shown , while in the lower graph a case with W = 0.01
0 W′(1 W )
isshownwithareweightingfactorofW =99(leadingto W = 0 − 0 . (46)
W0′ =W2′ =0.5). W0(1−W0′)
Let us now calculate the standard deviations for this
Simulation results for the standard deviation of f as a
case. We are interested in minimizing the standard de-
function of W′ is shown in Fig. 9. We see that the
viation of a variable evaluated in a simulation and also 0
reweighting actually increases the standard deviation.
need to calculate its statistical error. In a standard MC
This is due to the rapidly increasing auto correlation
simulationoneusuallywantstocalculatetheaverageand
times. The autocorrelationfunction
the standard deviation of the average for some quantity
x. This is typically achieved by dividing the run into a δ (i+t)δ (i) δ (i) 2
n,0 n,0 n,0
C (t)= h i−h i (47)
number of bins, N, and saving the average of x for each δ δ (i)2 δ (i) 2
n,0 n,0
bin. Ifthebins arestatisticallyindependentthefinalav- h i−h i
erageandstandarddeviationcanbecalculatedaccording isshowninFig.10,andonecanseethattheautocorrela-
to tiontimes(inverselyproportionaltotheslopesinFig.10)
9
0.5 100
10−1
0.4
σ/ff <s>10−2
0.3
10−3
0.2 10−4
0.0 0.1 0.2 0.3 0.4 0.5 1 10 100 1000
w’ V
0
FIG. 9. The relative standard deviation deviation of f FIG. 11. The average sign as function of system volume
′
as a function of W0, Calculations are performed over bins at temperature T/J = 1. Shown are results of an unre-
containing N =2000 MC steps. strictedsimulation(circles),andasimulationrestrictedtothe
bin
0- and 2-meron sectors without reweighting (squares). The
line shows a slope of -2.
1.00
averagesign in a simulation of the semifrustrated model
with J(1)=J(√2)=J is shown as a function of lattice
volume V = L L at a temperature T/J = 1.0. For
×
comparison we first performed a standard simulation by
Cδ 0.10 sampling all meron sectors, which leads to a severe sign
problem with a sign that decreases exponentially in sys-
tem volume. Next we sampled only the zero- and two-
meron sectors without reweighting, which dramatically
increases the averagesign. The scaling changes from ex-
0.01 ponential to quadratic in the volume, as can clearly be
0 50 100 150 200 seen from the graph.
MC steps
FIG. 10. Autocorrelation function, Eq.(47), for different
′ ′
valuesofW0. Shownareresultsfor W0=0.02, 0.05, 0.1, 0.2, 0.3
0.3, 0.4, and 0.5 (curvesleft to right).
are proportional to W′. Notice that the longest auto- 0.2
0
correlation times are significantly shorter than the indi- χ
vidual bins (consisting of 2000 MC steps) used above, a σ/χ
criteria for the analysis to be valid.
0.1
This simple example seems to indicate that reweight-
ing does not decrease the statistical errors. However, in
a standard Monte Carlo simulation the measured quan-
tities are notindependent evenwith no reweighting,and 0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
formula Eq. (27) contains measured quantities different
w’
from Eq. (36) consideredabove. Therefore the reweight- 0
ing will affect autocorrelation times differently than in FIG.12. Relativestandarddeviationofthez-susceptibility
theaboveexample,andreweightingcanactuallydecrease of the semifrustrated model as a function of reweighting.
the standard error.20 Shown are results for systems of linear size 8 (circles), 12
Using the above technique we can study how the rel- (squares), 16 (diamonds), and 20 (triangles). The standard
ative error in the susceptibility of the semi-frustrated deviation is calculated for bins containing Nbin = 103 MC
model is affected by reweighting. An initial run without steps at temperatureT/J =1.0.
reweighting has to be done first to determine the aver-
age sign δ = W . Thereafter Eq. (46) can be used Having determinedthe averagesign without reweight-
n,0 0
to determhine ithe desired weight factors. In Fig. 11 the ing we now use Eq. (46) to determine the desired weight
factors. In Fig. 12 the standard deviation (44) of the
10