Introduction to the
Stochastic Series Expansion method
Anders Sandvik, Boston University• Illustration of concept; classical Monte Carlo example• Detailed account of SSE for the S=1/2 Heisenberg modelThis presentation is based on material available at
http://physics.bu.edu/~sandvik/programs/

A simple SSE program (Fortran90) for the2D Heisenberg model can be downloaded from this siteWarm-up: SSE for a classical problem
Classical thermal expectation value

Classical (e.g., Ising) spins:
Classical Monte Carlo: Importance sampling of spin configurations
Probability of generating a configuration

Estimate of expectation value based on sampled configurations

Imagine that we are not able to evaluate the exponential function
How could we proceed then?

Use Taylor expansion of the exponential function

Expansion power n is a new “dimension” of the configuration space
To ensure positive-definitness we may have to shift E (must be < 0)
The sampling weight for the configurations (s,n) is

The function to be averaged (estimator) f(s) is the same as before;
it does not depend on n

However, if f(s) is a function of the energy it can be rewrittenas a function of n only!

Define:Shift summation index: m=n+1

Therefore the energy expectation value is

We can also easily obtain

And thus the specific heatis

What range of expansion orders n is sampled?
From the preceding results we obtain

Consider low T; C Æ 0Thus, for a system with N spins:
Average expansion order
Width of distribution
These results hold true for quantum systems as wellIn the quantum case H consists of non-commuting operators:
requires more complicated treatment

Quantum-mechanical SSE
Thermal expectation value

Choose a basis and Taylor expand the exponential operator

Write the hamiltonian as a sum of local operatorsa = operator type (e.g., 1=diagonal, 2=off-diagonal)b = lattice unit (e.g., bond connecting sites i,j)such that for every a, b:(no branching)Write the powers of H in terms of “strings” of these operators

Operator strings of varying length n• as in the classical case
Fixed-length operator strings: introduce unit operator:
Expansion cut-off M: add M-n unit operators to each string• there are M!/n!(M-n)! ways of doing this fin = numberof non-[0,0]operators
The truncation should not be considered an approximation• M can be chosen such that the truncation error is negligible

The terms
are sampled according to weight in this sum• requires positive-definiteness• to this end, a constant may have to be added to diagonal Hab• there can still be a “sign problem” arising from off-diagonal Hab

SSE algorithm for the S=1/2 Heisenberg model• The algorithm for this model is particularly simple and efficient

Consider bipartite lattice (sign problem for frustrated systems)

Standard z-component basis:
Bond operators: bond b connects sites i(b),j(b)Diagonal and off-diagonal bond operators

A minus sign in front of the off-diagonal H2b is neglected• this corresponds to a sublattice rotation; 180 degree rotation in
the xy-plane of the spin operators on sublattice B• The sign is irrelevant for a bipartite lattice (will be shown later)
SSE operator string
Represented in the computer program by
Spin state |a> represented by

SSE partition function

Both H1b and H2b give 0 when acting on parallel spins• non-zero matrix element = 1/2 in both cases
Define propagated states

For a contributing configuration:
The configuration weight is then(periodic)

Periodicity requires an even number of spin flips• This is why the sign of H2b is irrelevant for a bipartite lattice• For a frustrated lattice an odd number of flips is possible

Graphical representation• 1D example; 8 spins, M=12

1D: bond b connects sites b and b+1

Linked-list representation• vertex: operator and spins before and after the operator has acted• replace spins between vertices by links• linked vertex list used in some parts of the programA vertex has 4 “legs”, numbered l=0,1,2,3:

position p of operator in operator string opstring[p], vertex leg lfi position v in linked vertex list: v=1+l+4*(p-1)
vertexlist[v] contains the element # to which v is linked
l =[v] vertexlist[v]:01[ 1] 31[ 2] 32[ 5] 0[ 6] 0[ 9] 43[10] 44[13] 35[14] 47[17] 4[18] 11[21] 0[22] 0[25] 0[26] 0[29] 3[30] 18[33] 15[34] 16[37] 0[38] 0[41] 20[42] 12[45] 36[46] 48

23p[ 3] 29[ 7] 0[11] 18[15] 33[19] 30[23] 0[27] 0[31] 1[35] 13[39] 0[43] 9[47] 14[ 4] 17[ 8] 0[12] 42[16] 34[20] 41[24] 0[28] 0[32] 2[36] 45[40] 0[44] 10[48] 46

1
2
34
56
7
8
910
1112

Sampling the SSE configurations; updates
1) Diagonal update- replace unit operator by diagonal operator, and vice versa2) Off-diagonal update (local or loop)- change the operator type, diagonal
off-diagonal,
for two (local) or several (loop) operators
3) Flip spins in the state |a>- unconstrained “free” spins; weight unchanged after flip
- only possible at high temperatures; strictly not necessary
Updates satisfy detailed balance:Diagonal update• Carried out in opstring[p] for p=1,...,M• State |a(p-1)> stored in spin[]
insertion ofdiagonal operator

removal ofdiagonal operator

off-diagonal
no change,
propagate stateInsertion of a diagonal operator if opstring[p]=0Generate bond index b at random, attempt opstring[p]=2*b• can only be done if spin[i(b)]≠ spin[j(b)]• n increases by 1; weight ratio

Removal of a diagonal operator if opstring[p]≠0• n decreases by 1; weight ratio
B ways of selecting b but only one way of removing an operator;

Accept probabilities:

Local off-diagonal update (obsolete)
Change type of 2 operators on the same bond
• cannot always be done; check for constraining operators• no weight change; accept with fixed probability (e.g., P=1)

Note: periodic boundary conditions in the “propagation” direction• update spanning across the boundary affects the stored state |a>Local updates typically are not very efficient• critical slowing-down• no winding-number or particle-number fluctuations

Loop update• carried out in the linked-vertex-list representation• move “vertically” along links and “horizontally” on the same operator• spins flipped at all vertex-legs visited; operator type changes; weight unchanged• construct all loops, flip with probability 1/2 (as in Swendsen-Wang)

Monte Carlo step
• a cycle of diagonal updates (p=1,...,M in opstring[p])• construction of the linked vertex list
• construct all loops, flip each with probability 1/2• map updated vertex list back to opstring[], spin[]
Starting the simulation• “empty” perator string, opstring[p]=0, p=1,...,M• M is arbitrary, e.g., M=20
• random spin state; spin[p]=+1,-1
Determining the cut-off M• after each, MC step, compare expansion order n with M• if M-n<n/a, with, e.g., a=3, then M=n+n/a

Generalization of loop update; directed loops
In the case of the isotropic S=1/2 model• There are only 4 non-0 vertices• The operators uniquely define all loops• Loops are non-self-intersecting
Directed loops• In general, there are more than 4 allowed vertices• A vertex is entered at some entrance leg• The path can proceed (exit) through any of the 4 legs
• Exit probabilities are obtained from directed-loop equations• Loops can back-track (“bounce”) and self-intersect• Bounces can be avoided for some models (more efficient)