### The stochastic series expansion (SSE) method

Historical backround; Handscomb's method

In the late 1970's and early 1980's, several types of Monte Carlo methods for quantum lattice models at finite temperature were developed based on the Trotter decomposition of the operator exp(-H/T), where T is the temperature and H the Hamiltonian. An earlier proposal by Handscomb, dating back to the early 1960's, to instead use the power-series expansion of exp(-H/T), was largely neglected because its applicability appeared to be very limited in practice. Handscomb's method was originally developed for the ferromagnetic S=1/2 Heisenberg model and was later extended also to the S=1/2 antiferromagnet. Further generalizations were not successful, however, because the scheme relied on the ability to calculate exactly the traces of the operator products appearing when the powers of H are written as sums of products of the individual bond operators in H (e.g., S_i*S_j in the case of the Heisenberg model). In practice the traces can only be calculated exactly in a few special cases.

Generalization of Handscomb's method; stochastic series expansion

The problem of generalizing Handscomb's method to Heisenberg models with higher spins was solved by Sandvik and Kurkijärvi in 1989 . The solution is in hindsight quite simple: Instead of attempting to calculate the traces analytically, they are written as a sum of diagonal matrix elements in a suitably chosen basis. The matrix elements can be very easily calculated and the sum can be sampled along with the operator sequences. Such an approach can in principle be used for any model written in a discrete basis. An implementation for the 1D Hubbard model was presented in , where a more efficient sampling scheme was also developed. In practice, sign problems appear for the same models as in methods based on direct sampling of the path integral (world-line methods), and hence one is in practice limited to non-frustrated spin and boson systems and one-dimensional fermion systems.

The term "stochastic series expansion", or SSE, refers to this generalization of Handscomb's scheme, where individual terms of the traces of the power series expanded exp(-H/T) are sampled. Handscomb's method can then be considered as a special case of SSE, applicable in cases where the full traces can be calculated exactly. However, even in the cases where Handscomb's method can be used, the isotropic S=1/2 Heisenberg antiferromagnet in particular, the algorithm based on the full SSE scheme is preferrable, due to the more powerful updating schemes available.

The SSE representation of quantum statistical mechanics

The SSE approach constitutes an alternative to the standard path integral for representing a quantum statistical mechanics problem in terms of classical variables that can be sampled using Monte Carlo techniques. There is a close relationship between the SSE representation and the path integral in continuous imaginary time, which was explored in detail in . In the SSE representation the continuum limit is achieved in practice although the space is discrete. SSE was the first widely applicable finite-temperature quantum Monte Carlo method that was directly formulated in the continuum limit, so that numerically exact results could be obtained without having to extrapolate to zero discretization. Although path integral (world-line) methods have also been formulated directly in the continuum limit [B. Beard and U. J. Wiese, Phys. Rev.Lett. 77, 5130 (1996), N. V. Prokof'ev, B. V. Svistunov, and I. S. Tupitsyn, Sov. Phys JETP 87, 310 (1998); cond-mat/9703200], the discreteness of the SSE representation is often an advantage in practice.

Loop-cluster algorithms for SSE and path integral methods

The SSE method was originally based on local updates for sampling the states and operator sequences. In  it was shown that cluster-type updates, that previously had been developed in the context of world-line quantum Monte Carlo (called loop or worm algorithms) could also be constructed for sampling in the SSE representation. In the SSE formulation, the clusters are constructed in a linked list of operators and were therefore called operator-loops. A systematic way to find efficient operator-loop algorithms was presented in , and it was also shown that the scheme could be directly translated into the path integral formalism. The generalized clusters are called "directed loops", as the detailed balance equations derived for them are associated with a directionality when the loops are constructed. The directed loop formalism allows for highly efficient simulations also in the presence of external fields.