arXiv:0704.1469v1  [cond-mat.str-el]  11 Apr 2007
Monte Carlo Simulations of Quantum Spin
Systems in the Valence Bond Basis
Anders W. Sandvik1 and K. S. D. Beach1,2
1 Department of Physics, Boston University, Boston, Massachusetts, USA
2 Theoretische Physik I, Universität Würzburg, Würzburg, Germany
Abstract. We discuss a projector Monte Carlo method for quantum spin models
formulated in the valence bond basis, using the S = 1/2 Heisenberg antiferromagnet
as an example. Its singlet ground state can be projected out of an arbitrary basis
state as the trial state, but a more rapid convergence can be obtained using a
good variational state. As an alternative to first carrying out a time consuming
variational Monte Carlo calculation, we show that a very good trial state can be
generated in an iterative fashion in the course of the simulation itself. We also show
how the properties of the valence bond basis enable calculations of quantities that
are difficult to obtain with the standard basis of Sz eigenstates. In particular, we
discuss quantities involving finite-momentum states in the triplet sector, such as
the dispersion relation and the spectral weight of the lowest triplet.
1 Introduction
Quantum Monte Carlo (QMC) simulations of spin systems have traditionally
been carried out in the basis of eigenstates of the spin-z operators Szi , i =
1, . . . , N , i.e., the basis of “up” and “down” spins in the case of S = 1/2
(which is the case we consider here). For the prototypical model of interacting
quantum spins, the antiferromagnetic (J > 0) Heisenberg hamiltonian,
H = J
〈i,j〉
Si · Sj = J
〈i,j〉
[Szi S
(S+i S
j + S
j )], (1)
this basis is clearly natural and convenient, as an off-diagonal operator acting
on a basis state just flips two spins or destroys the state. Starting with the
work of Suzuki [1], finite-temperature simulation methods employing the spin-
z basis were developed in which a quantum mechanical expectation value for
a system in D dimensions is mapped onto an anisotropic classical statistical-
mechanics problem in D+1 dimensions—the discretized [2–4] or continuous
[5,6] imaginary-time path integral. There are now very efficient methods uti-
lizing loop-cluster [7,8,5] or “worm” [6] updates of the world-line spin config-
urations. These methods have enabled studies of systems with ≈ 104 − 105
spins in the low-temperature (ground-state) limit and much more at elevated
temperatures. Loop updates have been developed [9,10] also for the alterna-
tive and now frequently used power-series expansion representation [11–14]
of the partition function (stochastic series expansion; SSE), where the spin-z
http://arxiv.org/abs/0704.1469v1
2 A. W. Sandvik and K. S. D. Beach
|Vk> |Vl> <Vl|Vk>
Fig. 1. Two valence bond states |Vk〉, |Vl〉 in two dimensions and their overlap graph
corresponding to 〈Vl|Vk〉 = 2
N◦−N/2, where N◦ is the number of loops formed (in
this caseN◦ = 3 and the number of sitesN = 16). Filled and open circles correspond
to sublattices A and B. The sign convention in Eq (3) for a singlet valence bond
(i, j) dictates that spins i and j belong to A and B, respectively
basis is also normally used. It is in principle possible to adapt these ap-
proaches to other local bases, e.g., that of singlet and triplet states of spin
pairs on a dimerized lattice. This basis is often used in diagrammatic and
series-expansion calculations [15], but its implementation in QMC simula-
tions is typically rather cumbersome.
Zero-temperature (T = 0) simulations, in which the ground state is pro-
jected out of a trial wave function, are also normally carried out in the
spin-z basis [16,17]. Here we will discuss an alternative approach to ground-
state calculations which turns out to have some unique features enabling
access to quantities that are normally difficult to obtain with standard finite-
temperature or projector methods. We make use of the valence bond basis,
i.e., states in which the spins are paired up into singlets;
|V 〉 = |(i1, j1)(i2, j2) · · · (iN/2, jN/2)〉. (2)
Here (i, j) denotes a singlet formed by the spins at sites i and j;
(i, j) = (| ↑i↓j〉 − | ↓i↑j〉)/
2, (3)
and the total number of sites N is assumed to be even. While in principle
one can include all possible pairings of the spins, it is in most cases better to
consider a smaller basis in which the sites are first divided into two groups,
A and B, of N/2 spins each, and to only consider singlets (i, j) in which
the first index i ∈ A and the second j ∈ B [18–20]. In the case of a bipartite
lattice, these groups are naturally the two sublattices, as shown in Fig. 1. This
restricted VB basis has (N/2)! states and is still massively overcomplete—the
singlet space has N !/[(N/2)!]2(N/2+1) dimensions [18]. The VB basis states
are all non-orthogonal, overlapping with each other according to the simple
loop rule illustrated in Fig. 1.
The VB basis was introduced already in the early 1930s [18,21,22] and
has played an important role in exactly solvable models [18,23–25]. Later,
it became a tool for describing spin liquids—the resonating valence bond
(RVB) mechanism introduced by Fezekas and Anderson [26–28], in which the
ground state is dominated by short valence bonds. In exact diagonalization
Monte Carlo Simulations of Quantum Spin Systems 3
studies, the VB basis is useful in cases where it is a good approximation to
only consider a restricted (and incomplete) space of short bonds (spin liquids
and other states with no magnetic long-range order) [29–32]. Variational cal-
culations in the VB basis have been carried out for the 2D Heisenberg model
[20,33]. Furthermore, Liang realized that a variational VB state could be con-
siderably improved by stochastically projecting it with an operator (−H)m
for large m [34]. Later, Santoro et al. devised a Green’s function method for
calculating energies in the VB basis [35]. Despite the promising results ob-
tained in these studies, there was, to our knowledge, no further developments
of QMC methods in the VB basis until one of us recently introduced two
related projector algorithms [36], improving on the schemes of Liang [34] and
Santoro et al. [35]. These algorithms have already been applied in studies of
quantum phase transitions [37,38] and entanglement entropy [39].
Some previously unnoticed advantages of the VB basis in QMC algorithms
were pointed out in Ref. [36]. Here we summarize our recent work on VB
projector methods and highlight some of their unique features. We discuss
in particular a scheme for “self-optimizing” the trial state out of which the
ground state is projected, and also show how to study properties of triplet
excitations at finite momentum.
2 Ground state projection
Consider a state |Ψ〉 and its expansion in terms of eigenstates |n〉, n = 0, 1, . . .,
of some hamiltonian H ;
|Ψ〉 =
cn|n〉. (4)
With C a constant chosen such that the lowest eigenvalue E0 − C is the
largest in magnitude, a large number m of repeated operations with C −H
projects out the ground state,
(C −H)m|Ψ〉 → c0(C − E0)m
|0〉+ c1
C − E1
C − E0
|1〉+ . . .
, (5)
provided that the overlap c0 6= 0. Here we will first consider singlet eigenstates
of the Heisenberg model (1), which can be expanded in VB states;
|Ψ〉 =
fi|Vi〉. (6)
Because of the overcompleteness of the VB basis, this expansion is not unique.
That, however, does not prohibit the ground state |0〉 to be projected out
according to Eq. (5). The Heisenberg hamiltonian can be written in terms
of singlet projection operators Hb ≡ Hi(b),j(b) on the interacting spin pairs,
{i(b), j(b)}, b = 1, . . . , Nb = DN (for a periodic cubic D-dimensional lattice);
H = −J
Hb = −J
Hi(b),j(b), Hij = −(14 − Si · Sj). (7)
4 A. W. Sandvik and K. S. D. Beach
When a singlet projector Hij acts on a VB basis state, one of two things
can happen; 1) if i, j belong to the same bond the state is unchanged with a
matrix element 1, or 2) if they belong to different bonds these two bonds are
reconfigured (“flipped”) with matrix element 1/2;
Hij | · · · (i, j) · · ·〉 = | · · · (i, j) · · ·〉, (8)
Hij | · · · (i, k) · · · (l, j) · · ·〉 = 12 | · · · (i, j) · · · (l, k) · · ·〉. (9)
Here positive-definitness of (9) is directly related to the two sites i and j being
in different sublattices. For a frustrated system, where there are operators
with both sites i, j in the same sublattice, positive-definitness does not hold
[19,20]. For a non-frustrated system the simple bond flip (9) makes for a
convenient stochastic implementation of the ground state projection (5). We
write the projection operator as (with J = 1 henceforth)
(C −H)m =
Pr, Pr = Hbr
· · ·Hbr
, (10)
where we have introduced a compact notation Pr, r = 1, . . . , N
b , for the
different strings of singlet projectors. When a string Pr acts on a VB basis
state |V 〉 the result is another basis state, which we denote |V (r)〉, with
a prefactor (weight) Wr which is simply given by the number moff of off-
diagonal operations in the course of evolving |V 〉 to |V (r)〉;
Pr|V 〉 = Wr|V (r)〉, Wr = 2−moff . (11)
We here first consider projecting the ground state out of a single VB
basis state; later we will consider the use of a more complicated trial state.
We consider two ways to calculate expectation values:
〈H〉 =
r〈Ψ |HPr|V 〉
r〈Ψ |Pr |V 〉
r Wr〈Ψ |H |V (r)〉
r Wr〈Ψ |V (r)〉
, (12)
〈A〉 =
rl〈V |P ∗l APr|V 〉
rl〈V |P ∗l Pr|V 〉
rl WrWl〈V (l)|A|V (r)〉
rl WrWl〈V (l)|V (r)〉
. (13)
We will discuss how to estimate these using importance sampling; terms
(configurations) of the numerators are illustrated in Fig. 2. We will refer to
(12) and (13) as the single and double projection, respectively.
In (12), which is an exact (when m → ∞) expectation value only of the
hamiltonian (or other operators for which the ground state is an eigenstate)
the state |Ψ〉 is in principle arbitrary. It is very convenient to use a state
which has equal overlaps with all the VB basis states, e.g., the Néel state
|ΨN 〉 (all spins up on sublattice A and down on B). It is easy to see that
〈ΨN |V 〉 = 2−N/2 for any basis state |V 〉. Since H |V (r)〉 = −
b Hb|V (r)〉 is
a sum of basis states multiplied by factors −1 or −1/2, the overlaps with 〈ΨN |
drop out altogether and do not have to be considered further. If the projector
Monte Carlo Simulations of Quantum Spin Systems 5
|V ><V | A
(a) (b)
|V >H
Fig. 2. Propagation of a VB state on a 6-site chain. The horizontal bars represent
nearest-neighbor Heisenberg interactions (singlet projectors). In the single projec-
tion (a) the state is propagated from right to left, and an estimator for the ground
state energy is obtained by acting once more with all terms of the hamiltonian. In
the double projection (b) the state is projected from the right and the left, and
any operator expectation value can be estimated by calculating the corresponding
matrix elements between the propagated states
strings Pr in (12) are importance-sampled according to their weights Wr, the
estimator for the ground state energy is thus
E0 = 〈H〉 = −〈md + 12mo〉, (14)
where md and mo are, respectively, the number of diagonal and off-diagonal
operations Hb|V (r)〉 (and md +mo = Nb). It should be noted that although
this estimator is exact in the limit m → ∞, it is not variational. The correct
energy may thus be approached with increasing m from above or below.
Eq. (13) is valid for any expectation value and in the case of A = H gives
a variational estimate of the energy. Using WrWl〈V (l)|V (r)〉 as the sampling
weight, the estimator for any 〈A〉 is of the form
〈A〉 =
〈V (l)|A|V (r)〉
〈V (l)|V (r)〉
. (15)
In the case of a spin correlation function 〈Si·Sj〉, the matrix element is related
to the loop structure of the overlap graph [19,20] (illustrated in Fig. 1):
〈Vl|Si · Sj |Vr〉
〈Vl|Vr〉
+3/4, if i, j ∈ same loop, same sublattice,
−3/4, if i, j ∈ same loop, different sublattices,
0, if i, j ∈ different loops.
Measuring the spin correlations is hence straight-forward once the overlap-
loops have been constructed. Higher-order functions, e.g., dimer-dimer cor-
relations 〈(Si · Sj)(Sk · Sl)〉, are also related to the loop structure [40].
Note again that no bond operator Hb can destroy a VB state and that
all the states have non-zero overlap with each other. Thus all terms in (13)
contribute to the expectation value. This turns out to be an advantage in
constructing a Monte Carlo algorithm, as any change made in the operator
strings can be accepted with some probability. With an orthogonal basis, such
6 A. W. Sandvik and K. S. D. Beach
as the Szi eigenstates, there would be considerable constraints, both in terms
of individual operators in the projection [the spin flip operator in (1) can act,
without destroying the state, only on anti-parallel spins], and in ensuring a
non-zero overlap between the propagated states (the two propagated states
have to be identical). Note also that the singlet projectors are non-hermitian
in the VB basis. As indicated in (13), and illustrated in Fig. 2, we here
propagate two states, |V (l)〉 ∝ Pl|V 〉 and |V (r)〉 ∝ Pr|V 〉, and subsequently
compute their overlap and various matrix elements. Propagating |V 〉 with
P ∗l Pr and then taking the overlap with |V 〉 is not equivalent term-by-term.
To carry out the projection stochastically, the operator strings are stored
in arrays [Pα] = [b
1 ][b
2 ] · · · [bαm], where bαi ∈ {1, . . . , Nb} with α = 1, 2 cor-
responding to Pr and Pl, respectively, in the double projection; in the single
projection α is redundant. A table holds the site pairs i(b), j(b). The state
|V 〉 is stored in a list [V ] = [v1][v2] · · · [vN ] where vi = j and vj = i if there is
a valence bond at (i, j). Propagation with the bond flips (9) is easily carried
out in this representation. The state list [V ] is then first copied into two lists,
[V1] and [V2], in which |V (r)〉 and |V (l)〉 are constructed.
The simulation can be started with a randomly generated operator string.
The strings can can be updated in an trivial way, by changing a number R
of operators at random. In either [P1] or [P2], R positions pj , j = 1, . . . , R,
pj ∈ {1, . . . ,m} (all different) are generated. Their contents bαpj are picked
randomly from the set {1, . . . , Nb} (excluding the old value for each pj). To
calculate the Metropolis acceptance probability, the state is propagated with
the updated operator string and the number of off-diagonal operationsmoff in
(11) is counted. In the single expansion, the acceptance probability is simply
Paccept = min
W newr
W oldr
= min[2m
−mnew
off , 1], (17)
whereas in the double projection an overlap ratio appears as well. In the
double projection, we change operators only in one of the operator strings at
a time, so that only one state has to be propagated. It is customary to define
a size-normalized Monte Carlo step (or “sweep”). For projector length m we
do m replacement attempts, and so our sweep is independent of the number
of operators R replaced in each update.
Normally, in Monte Carlo simulations one does not compute the full
weight, because it is possible to more speedily calculate just the change in
the weight [the weight ratio in (17)]. In the present formulation of the VB
projector algorithm the weight is, however, recomputed from scratch each
time, because a better scheme is not yet known. Each update hence requires
on the order of m operations. In the double expansion, construction of the
loops needed to compute the overlap scales as N , but typically m > N and
the propagation of the state dominates the simulation. In spite of the need
to recalculate the weight, the scheme is sufficiently efficient to compete with
other ground state QMC methods. More importantly, as we will discuss in
Sec. 4, the VB basis offers access to quantities out of reach for other methods.
Monte Carlo Simulations of Quantum Spin Systems 7
1 2 3 4 5
1 2 3 4 5
(a) (b)
Fig. 3. (a) Acceptance rate in a double projection and (b) the number of bonds
changed in the propagated state in accepted updates for a 16× 16 lattice (symbols
with lines; the bare lines are for L = 8) versus the projection length (normalized
by the number of sites N). The trial state |V 〉 was a columnar dimer state
The optimum number R of operators to replace depends on the accep-
tance rate. In Fig. 3(a) we show the acceptance rate for a double-projected
2D system versus the length m of the projection for R = 1− 6. As expected,
the acceptance rate decreases with increasing R, but it does not change ap-
preciably with m. It also depends only weakly on the lattice size. Multiplying
the acceptance rate by R gives the average number of operators changed; it
initially increases with R but has a maximum for R = 7 (for L = 16). The
optimum R is clearly model/lattice dependent. Another characteristic of the
update is the number of bonds changed in the projected state |V (r)〉 as a
consequence of the modifications of Pr. This number is shown in Fig. 3(b).
It is seen to increase with R, as expected. As a function of m the num-
ber of changed bonds decreases. This behavior reflects a loop structure of
the singlet-projection operators [41,8], which implies that some changes “up-
stream” in the operator string may be healed further downstream in the
propagation. For a finite lattice in the limit m → ∞ one would expect a
substitution of an operator far upstream in Pr to have no effect on the final
propagated state |V (r)〉. This does not imply that this update is inconsequen-
tial, as the sampling is over paths, not just the final states in the propagations.
In principle, it should be possible to take advantage of the underlying loop
structure of the singlet projectors [41,8] to devise a loop update in the VB
basis, analogous to such updates in world-line [7,8] and SSE [9] simulations.
However, we have not yet been able to construct a scheme that is in practice
faster than the trivial random substitution with full state propagation.
3 Self-optimized trial state
So far, we have projected the ground state out of a single VB basis state |V 〉.
This works well [36], but the rate of convergence with m of course depends
8 A. W. Sandvik and K. S. D. Beach
on the state chosen; ideally one would like to maximize the overlap 〈V |0〉.
One way to obtain a typically good single-configuration trial state is to first
start with an arbitrary one; a regular bond pattern or a randomly generated
configuration. After carrying out some projection steps with this state, the
current propagated state |V (r)〉 is chosen as a new trial state. Since this state
has been generated in the projection it should contribute substantially to the
ground state and hence typically will be better than a completely arbitrary
one. However, as we will discuss next, we can do much better than this.
Liang’s original motivation for introducing a projector technique in the
VB basis was to improve on a variational calculation [34]. Previously, Liang,
Doucot, and Anderson had studied a variational amplitude-product state for
the 2D Heisenberg model of the form [20]
|Ψ〉 =
fk|Vk〉, fk =
h(xbk, ybk), (18)
where xbk and ybk are the x- and y-lengths of bond b in VB state k, as il-
lustrated in Fig. 4(a). Liang et al. tried power-law and exponential forms
depending only on the total length l of the bonds [the ”Manhattan” length
l = x + y was used, but defining l = (x2 + y2)1/2 should not change things
qualitatively], in addition to keeping several short-bond amplitudes as pa-
rameters to optimize [20]. More recently, all the amplitudes were optimized
without any assumed form on lattices with up to 32×32 sites, with the result
that h(l) ∼ l−3 for long bonds [33]. One of us has also recently showed the
more general result h(l) ∼ l−(D+1) within a mean-field approach for a D-
dimensional cubic lattice [42]. In 2D the fully-optimized amplitude-product
state turns out to be extremely good, with an energy deviating by only
≈ 0.06% from the exact ground state energy and with the long-distance
spin-spin correlations reproduced to within 2% [33].
With the state (18) an expectation value is given by
〈A〉 =
kp fkfp〈Vp|Vk〉
〈Vp|A|Vk〉
〈Vp|Vk〉
rl fkfp〈Vp|Vk〉
, (19)
which can be evaluated using importance sampling of the VB configurations
with weight fkfp〈Vp|Vk〉. Liang et al. introduced a very simple scheme for
updating the dimer configurations [20], which we here illustrate in Fig. 4(b).
Choosing two next-nearest-neighbor sites, i.e., ones on a diagonal of a 4-site
plaquette (or, in principle, any two sites on the same sublattice), the two
bonds connected to them are reconfigured in the only possible way which
maintains only bonds between the A and B sublattices, as shown in the figure.
Labeling the two initially chosen sites 1 and 2, and the bonds connected to
them b = 1, 2, the Metropolis acceptance probability is, assuming that the
bond update was made in |Vk〉, resulting in |Vk′ 〉,
Paccept = min
h(x1k′ , y1k′)h(x2k′ , y2k′)
h(x1k, y1k)h(x2k, y2k)
〈Vp|Vk′ 〉
〈Vp|Vk〉
. (20)
Monte Carlo Simulations of Quantum Spin Systems 9
(a) (b)
Fig. 4. (a) Definition of the size of a bond b. (b) Reconfiguration of two bonds in
an update of the trial state. The two sites marked 1, 2 are chosen at random among
all pairs of next-nearest neighbor sites.
We can use an amplitude-product state as the trial state in the projector
QMC method, using some set of amplitudes not necessarily originating from
a variational calculation. In updating the bond configurations, we then also
must compute the new weight of the propagation; Pr|Vk′ 〉 = Wk′r|Vk′ (r)〉, and
of course 〈Vp|Vk〉 in (20) is replaced by 〈Vp(l)|Vk(r)〉 (in the case of the double
projection; for the single projection there is no overlap). The acceptance rate
of a state update is similar to that of an operator string update with a smallR,
and often we find it advantageous to combine state and operator updates. To
save some time, one can tentatively accept/reject a bond update based solely
on an amplitude ratio—Eq. (20) without the overlaps—and then calculate
the overlap and the propagation weight for a final accept/reject probability
only for tentatively accepted bond updates.
In the variational calculation 〈H〉 is minimized with respect to all h(x, y).
With a recently developed stochastic optimization method [33], all the ∝ N
amplitudes can be minimized for moderate-size lattices (up to 32 × 32 sites
were considered in Ref. [33]). In principle we could follow Liang [34] and use
the best possible variational state as our trial state in the projector QMC
method—indeed this can be expected to be the optimum starting point.
However, we will now describe a scheme which delivers a trial state nearly as
good as the best variational state, at a smaller computational cost.
Consider the probability distribution P (x, y) of valence bonds. In an
amplitude-product state (18) we would have P (x, y) ∼ h(x, y), were it not
for the “hard-core” constraint of only one bond per site. Even with this con-
straint, it is clear that the probabilities and amplitudes are related in a mono-
tonic way; increasing h(x, y) for some given (x, y), while keeping the other
amplitudes fixed, will lead to a larger P (x, y). This fact can be exploited
in constructing a good trial state. We define two different probability dis-
tributions, P0(x, y) and Pm(x, y), the former being the just discussed bond
probability in a trial state of the form (18) and the latter the probability
distribution in the projected state. For sufficiently large m, Pm is an exact
property of the ground state, whereas P0 is a property of the trial state and is
in general different from Pm. However, for given m, we can adjust the ampli-
tudes h(x, y) of the trial state such that Pm(x, y) = P0(x, y) for all x, y. If this
is done for m sufficiently large, then our trial state has a bond distribution
identical to that of the exact ground state. Such a state is often almost as
good as the best variational state. The reason that this is useful in practice
10 A. W. Sandvik and K. S. D. Beach
0.0 0.5 1.0 1.5 2.0 2.5
1.0 1.5 2.0 2.5
0.650
0.655
0.660
0.665
0.670
 self-optimized
 projected random
 columnar dimer
0.0 0.5 1.0 1.5 2.0
0.6694
0.6695
0.6696
0.6697
0.6698
0.6699
0.6700
 Projected state
 Self-optimized trial state
Fig. 5. Left: Energy versus projection length for a 16 × 16 system, using three
different trial states (labeled in the inset). Right: The energy of the self-optimized
trial state compared with the projected energy using that trial state. The dashed
line shows the energy evaluated independently using the SSE method [43]
is that it is very easy to adjust the amplitudes to achieve self-consistency.
Because of the monotonous relationship between h(x, y) and P0(x, y), we can
simply increase h(x, y) by some amount if P0(x, y) < Pm(x, y) and decrease it
if P0(x, y) > Pm(x, y), and repeat this until self-consistency is achieved. We
use the following scheme to update the amplitudes after the kth iteration;
ln[h(x, y)] → ln[h(x, y)] + RAN · β(k) · sign[Pm(x, y) − P0(x, y)], (21)
where RAN is a random number in the range [0, 1) and β(k) decreases with
the iteration step k = 1, 2, . . ., according to β(k) ∝ k−α. For the exponent,
we typically use α = 3/4. To evaluate the probabilities P0 and Pm (in two
independent simulations; with and without projection of the trial state), the
number of Monte Carlo sweeps does not have to be very large, because we
only need the sign of the difference of the two probabilities. We normally
use on the order of 100-1000 sweeps per iteration. Even if the stochastically
evaluated sign in (21) occasionally may be wrong, it is correct on average
and the amplitudes typically converge to a self-consistent solution after a
few hundred iterations. Due to the stochastic nature of the procedure, self-
consistency of course obtains only to within some statistical error, which can
be reduced by increasing the number iterations and/or sweeps per iteration.
Fig. 5 shows results for the energy of a 16× 16 system obtained in double
projections with three different trial states; a columnar dimer state, a ran-
domly picked state generated while projecting the columnar state, as well as
the self-optimized state. Already for the shortest projection, m = N/8 = 32,
the self-optimized state gives a projected energy which deviates by only 0.06%
from the exact ground state energy, and for largerm the energy is exact within
statistical errors. The other two trial states also lead to the correct energy
but only for much larger m. The energy of the self-optimized trial state itself
is shown in the right panel—its error is as small as that of the best variational
amplitude-product state [33].
Monte Carlo Simulations of Quantum Spin Systems 11
0.0 0.5 1.0 1.5 2.0
 self-optimized
 random projected
 columnar dimer
0.0 0.5 1.0 1.5 2.0
0.116
0.118
0.120
0.122
0.124
0.126
 Projected state
 Self-optimized trial state
Fig. 6. Long-distance spin correlations in the same runs as in Fig. 5. A different
“random projected” state is picked for each m, which leads to an un-smooth curve.
Fluctuations beyond the small statistical errors in the right panel reflect differences
in how closely the rather short self-optimization runs have approached P0 = Pm
Fig. 6 shows the long-distance spin correlation calculated in the same
runs. Again, the self-optimized state delivers superior results, although here
the convergence is not as fast as for the energy. The error of the spin correla-
tion in the trial state (right panel) is about 4% for large m, which is twice the
error in the best variational state [33]. Thus, the self-optimized state is not
identical to the best variational state, and an even faster convergence could
be achieved by using the fully optimized variational state. However, the vari-
ational calculation is much more time consuming than the self-optimization.
We can go beyond the amplitude-product state by taking into account
bond correlations. We are currently exploring this with both variational and
self-optimized states.
4 Triplet excitations
A unique advantage of the VB basis is that an mz = 0 triplet state can
be projected simultaneously with the singlet, with essentially no additional
overhead. Any triplet can be expanded in VB states where one of the bonds
corresponds to a triplet [18]; (i, j) → [i, j], where
[i, j] = (| ↑i↓j〉+ | ↓i↑j〉)/
2. (22)
Formally, such a triplet can be created by acting on a singlet with Szi − Szj ;
(Szi − Szj )(i, j) = [i, j]. (23)
To create a triplet |τ(q)〉 with some momentum q, we can apply
Szq =
eiq·rjSzj (24)
12 A. W. Sandvik and K. S. D. Beach
to a singlet state |σ(0)〉 with zero momentum;
Szq|σ(0)〉 ∝ |τ(q)〉. (25)
The amplitude-product state (18) for a periodic-boundary system has q = 0
if the number of bonds, N/2, is even, whereas for odd N/2 it has q = (π, π).
This simply follows from the fact that sublattice A → B and B → A when
translating by one lattice constant, whence each singlet (3) acquires a minus
sign. We typically work with systems with even N/2 (e.g., L×L lattices with
even L) and so we will here consider q = 0 singlets.
At the antiferromagnetic wave-vector, Q = (π, π), SzQ acting on an ar-
bitrary VB basis state |V 〉 can be written as a sum of N/2 terms of the
form (Szi − Szj )|V 〉, with i, j corresponding to the sites connected by bonds.
Operating on a q = 0 singlet thus gives
SzQ|σ(0)〉 = SzQ
fk|(i1, j1)k(i2, j2)k · · · (iN/2, jN/2)k〉
fk|(i1, j1)k · · · [ib, jb]k · · · (iN/2, jN/2)k〉, (26)
where the unspecified coefficients fk, e.g., the amplitude products in (18),
have translational invariance built in. Thus, for a triplet with q = Q the
wave function phases are buried in the definition of the singlets. Often, the
lowest excitation of a Heisenberg system is a q = Q triplet, which we thus
can sample without any difficulties with signs or phases. We consider this
case first, before turning to triplets with arbitrary momentum.
There are two possible actions of a singlet projectorHij on a triplet bond;
Hij | · · · [i, j] · · ·〉 = 0, (27)
Hij | · · · [i, k] · · · (l, j) · · ·〉 = 12 | · · · (i, j) · · · [l, k] · · ·〉, (28)
i.e., a diagonal operation on a triplet bond destroys the state whereas an
off-diagonal operation on one triplet and one singlet bond creates a singlet
at the sites on which the operator acts and moves the triplet to the other
two sites involved. Importantly, the matrix element remains the same as in
the off-diagonal operation on two singlet bonds. The weight of a triplet path
is therefore the same as the corresponding singlet path (11), except that the
triplet dies (giving zero weight) if an operator in Pr acts diagonally on the
triplet bond. We can thus measure triplet properties using paths generated in
a singlet simulation, by considering only those triplet paths that do survive
the propagation. In the trial state we have N/2 possible locations of the
triplet. All of them can be attempted collectively in a single propagation, by
keeping counters t(i) for the number of surviving states in which the triplet is
connected to site i (with i on sublattice A). Initially t(i) = 1 for all i. During
the propagation, for each diagonal operation (27) t(i) → t(i) − 1, and for
each off-diagonal operation (28) t(i) → t(i) − 1, t(l) → t(l) + 1. Eventually,
Monte Carlo Simulations of Quantum Spin Systems 13
as m → ∞, all triplets die; t(i) = 0 for all i, but typically there are enough
survivors left at large enough m to compute converged triplet properties.
An added advantage of calculating singlet and triplet properties in the
same run is that there are error cancellations which in some cases can increase
the statistical precision of differences, e.g., the singlet-triplet gap [36],
∆ = ET (π, π)− ES(0, 0), (29)
by up to orders of magnitude relative to two independent calculations. The
triplet energy ET (π, π) can be estimated using (15), taking into account that
a diagonal operation on a triplet bond gives zero, i.e., nd → nd − nt, where,
for surviving triplet configurations, nt = 0, 1 is the number of triplet bonds
of length 1. Other triplet properties have been discussed in Ref. [40].
We now discuss calculations with triplets of arbitrary momentum q. The
energy can be evaluated according to
E(q) =
r〈σ(0)|Sz−qHPrSzq|σ(0)〉
r〈σ(0)|Sz−qPrSzq|σ(0)〉
. (30)
We want to evaluate this expression using the sampled q = 0 singlet bond
configurations, and so we rewrite it as
E(q) =
r〈σ(0)|Sz−qHPrSzq|σ(0)〉
r〈σ(0)|Pr |σ(0)〉
r〈σ(0)|Sz−qPrSzq|σ(0)〉
r〈σ(0)|Pr |σ(0)〉
. (31)
The two factors can be evaluated based on sampling the propagations Pr and
the amplitude-product state |σ(0)〉 [which we have not explicitly written as
a sum of bond configurations in (30) and (31)]. In the singlet energy (12), we
could pick the Néel state for 〈Ψ | and then obtained the very simple expression
(15). Now we must consider the overlap with a momentum q triplet state.
We use 〈σ(0)|Sz−q, but in (31) have rewritten E(q) so that only the overlap
with 〈σ(0)| has to be considered for the sampling weight. Phases arising from
Szq only appear in the measurements, but in the end we have to evaluate
the ratio of the two quantities in (31), which can be challenging in practice.
However, close to q = (π, π) and (0, 0) [q 6= (0, 0)] we find that it can be
done; in some cases the method works even far away from these momenta.
We define the dispersion relative to the gap (29) at (π, π);
ω(q) = ET (q)−∆. (32)
Fig. 7 shows results for∆ and ω(q1), where q1 is the momentum closest to but
not equal to (π, π); q1 = (π−2π/L, π). We show the convergence as a function
of N/m for 4× 4 and 16× 16 lattices, comparing with exact diagonalization
results in the former case. From (5) one would expect the convergence to
be asymptotically exponential, which is seen clearly for L = 4. For L = 16,
the three largest-m points are equal within statistical errors, suggesting that
these results are also close to converged. Using ω(q1) = 0.62 for the L = 16
14 A. W. Sandvik and K. S. D. Beach
0.578
0.580
0.582
0.584
0.586
0 0.5 1 1.5 2
1.760
1.765
1.770
1.775
1.780
0 0.5 1 1.5 2
L=4 L=16
Fig. 7. Projection-length convergence of the singlet-triplet gap and the excitation
energy at momentum q1 = (π − 2π/L, π) for L = 4 (left) and 16 (right) lattices.
The horizontal lines show the exact results for L = 4
system gives the spin-wave velocity c = 1.58, in very close agreement with
the known value [43].
Another useful quantity accessible with the VB projector is the matrix
element 〈τ(q)|Szq|σ(0)〉, the square of which gives the single-magnon weight
in the dynamic structure factor (which is experimentally measurable using
neutron scattering). It can be calculated in a way similar to E(q), and we have
done so successfully for q close to (π, π). Results will be presented elsewhere.
Finally, we also note that the triplet bond-length distribution gives a di-
rect, albeit basis dependent, window into the “spinon” aspects [44] of the
excitations. Spinon deconfinement should be manifested as a delocalized dis-
tribution function, whereas two spinons bound into a magnon or “triplon”
should be reflected in a well-defined peak in the distribution function. We are
currently exploring this.
Acknowledgments
We would like to thank H. G. Evertz for useful discussions. This work was sup-
ported by the National Science Foundation under grant No. DMR-0513930.
References
1. M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976); M. Suzuki, S. Miyashita, and
A. Kuroda, Prog. Theor. Phys. 58, 1377 (1977)
Monte Carlo Simulations of Quantum Spin Systems 15
2. M. Barma and B. S. Shastry, Phys. Rev. B 18, 3351 (1977)
3. J. E. Hirsch, R. L. Sugar, D. J. Scalapino and R. Blankenbecler, Phys. Rev. B
26, 5033 (1982)
4. J. J. Cullen and D. P. Landau, Phys. Rev. B 27, 297 (1983)
5. B. B. Beard and U. -J. Wiese, Phys. Rev. Lett. 77, 5130 (1996)
6. N. V. Prokofév, B. V. Svistunov, and I. S. Tupitsyn, Pisḿa Zh. Eks. Teor. Fiz.
64, 853 (1996) [JETP Lett. 64, 911 (1996)]; Zh. Eks. Teor. Fiz. 114, 570 (1998)
[JETP 87, 311 (1998)]
7. H. Evertz, G. Lana and M. Marcu, Phys. Rev. Lett. 70, 875 (1993)
8. H. G. Evertz, Adv. Phys. 52, 1 (2003)
9. A. W. Sandvik, Phys. Rev. B 59, R14157 (1999)
10. O. F. Syljůasen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002)
11. D. C. Handscomb, Proc. Cambridge Philos. Soc. 58, 594 (1962); 60, 115 (1964)
12. D. H. Lee, J. D. Joannopoulos, and J. W. Negele, Phys. Rev. B 30, 1599 (1984)
13. A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991)
14. A. W. Sandvik, Phys. Rev. B 56, 11678 (1997)
15. V. N. Kotov, O. Sushkov, Z. Weihong, and J. Oitmaa, Phys. Rev. Lett. 80,
5790 (1998)
16. N. Trivedi and D. M. Ceperley, Phys. Rev. B 40, 2737 (1989)
17. S. Sorella and L. Capriotti, Phys. Rev. B 61, 2599 (2000)
18. L. Hulthén, Ark. Mat. Astron. Fys. 26A, No. 11 (1938)
19. B. Sutherland, Phys. Rev. B 37, 3786 (1988)
20. S. Liang, B. Doucot, and P. W. Anderson, Phys. Rev. Lett. 61, 365 (1988)
21. G. Rumer, Göttingen Nachr. Tech. 1932, 377 (1932)
22. L. Pauling, J. Chem. Phys. 1, 280 (1933)
23. C. K. Majumdar and D. K. Ghosh, J. Math. Phys. 10, 1388 (1969)
24. B. S. Shastry and B. Sutherland, Physica B&C 108, 1069 (1981)
25. I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799
(1987)
26. P. Fezekas and P. W. Anderson, Philos. Mag. 30, 23 (1774)
27. P. W. Anderson, Science 235, 1196 (1987)
28. E. Fradkin, Field Theories of Condensed Matter Systems (Adison-Wesley 1991).
29. P. L. Iske and W. J. Caspers, Physica 142A, 360 (1987)
30. M. Kohmoto, Phys. Rev. B 37, 3812 (1988)
31. S. Tang and H.-Q. Lin, Phys. Rev. B 38, 6863 (1988)
32. M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422
(2006)
33. J. Lou and A. W. Sandvik, cond-mat/0605034
34. S. Liang, Phys. Rev. B 42, 6555 (1990); Phys. Rev. Lett. 64, 1597 (1990).
35. G. Santoro, S. Sorella, L. Guidoni, A. Parola, and E. Tosatti, Phys. Rev. Lett.
83, 3065 (1999)
36. A. W. Sandvik, Phys. Rev. Lett. 95, 207203 (2005)
37. A. W. Sandvik, cond-mat/0611343
38. A. W. Sandvik and K. S. D. Beach, cond-mat/0612126
39. F. Alet, S. Capponi, N. Laflorencie, and M. Mambrini, cond-mat/0703027
40. K. S. D. Beach and A. W. Sandvik, Nucl. Phys. B 750, 142 (2006)
41. M. Aizenman and B. Nachtergaele, Commun. Math. Phys. 164, 17 (1994)
42. K. S. D. Beach (unpublished).
43. A. W. Sandvik. Phys. Rev. B 56, 11678 (1997)
44. S. Sachdev. Rev. Mod. Phys. 75, 913 (2003)

  We discuss a projector Monte Carlo method for quantum spin models formulated
in the valence bond basis, using the S=1/2 Heisenberg antiferromagnet as an
example. Its singlet ground state can be projected out of an arbitrary basis
state as the trial state, but a more rapid convergence can be obtained using a
good variational state. As an alternative to first carrying out a time
consuming variational Monte Carlo calculation, we show that a very good trial
state can be generated in an iterative fashion in the course of the simulation
itself. We also show how the properties of the valence bond basis enable
calculations of quantities that are difficult to obtain with the standard basis
of Sz eigenstates. In particular, we discuss quantities involving
finite-momentum states in the triplet sector, such as the dispersion relation
and the spectral weight of the lowest triplet.

Introduction
Quantum Monte Carlo (QMC) simulations of spin systems have traditionally
been carried out in the basis of eigenstates of the spin-z operators Szi , i =
1, . . . , N , i.e., the basis of “up” and “down” spins in the case of S = 1/2
(which is the case we consider here). For the prototypical model of interacting
quantum spins, the antiferromagnetic (J > 0) Heisenberg hamiltonian,
H = J
〈i,j〉
Si · Sj = J
〈i,j〉
[Szi S
(S+i S
j + S
j )], (1)
this basis is clearly natural and convenient, as an off-diagonal operator acting
on a basis state just flips two spins or destroys the state. Starting with the
work of Suzuki [1], finite-temperature simulation methods employing the spin-
z basis were developed in which a quantum mechanical expectation value for
a system in D dimensions is mapped onto an anisotropic classical statistical-
mechanics problem in D+1 dimensions—the discretized [2–4] or continuous
[5,6] imaginary-time path integral. There are now very efficient methods uti-
lizing loop-cluster [7,8,5] or “worm” [6] updates of the world-line spin config-
urations. These methods have enabled studies of systems with ≈ 104 − 105
spins in the low-temperature (ground-state) limit and much more at elevated
temperatures. Loop updates have been developed [9,10] also for the alterna-
tive and now frequently used power-series expansion representation [11–14]
of the partition function (stochastic series expansion; SSE), where the spin-z
http://arxiv.org/abs/0704.1469v1
2 A. W. Sandvik and K. S. D. Beach
|Vk> |Vl> <Vl|Vk>
Fig. 1. Two valence bond states |Vk〉, |Vl〉 in two dimensions and their overlap graph
corresponding to 〈Vl|Vk〉 = 2
N◦−N/2, where N◦ is the number of loops formed (in
this caseN◦ = 3 and the number of sitesN = 16). Filled and open circles correspond
to sublattices A and B. The sign convention in Eq (3) for a singlet valence bond
(i, j) dictates that spins i and j belong to A and B, respectively
basis is also normally used. It is in principle possible to adapt these ap-
proaches to other local bases, e.g., that of singlet and triplet states of spin
pairs on a dimerized lattice. This basis is often used in diagrammatic and
series-expansion calculations [15], but its implementation in QMC simula-
tions is typically rather cumbersome.
Zero-temperature (T = 0) simulations, in which the ground state is pro-
jected out of a trial wave function, are also normally carried out in the
spin-z basis [16,17]. Here we will discuss an alternative approach to ground-
state calculations which turns out to have some unique features enabling
access to quantities that are normally difficult to obtain with standard finite-
temperature or projector methods. We make use of the valence bond basis,
i.e., states in which the spins are paired up into singlets;
|V 〉 = |(i1, j1)(i2, j2) · · · (iN/2, jN/2)〉. (2)
Here (i, j) denotes a singlet formed by the spins at sites i and j;
(i, j) = (| ↑i↓j〉 − | ↓i↑j〉)/
2, (3)
and the total number of sites N is assumed to be even. While in principle
one can include all possible pairings of the spins, it is in most cases better to
consider a smaller basis in which the sites are first divided into two groups,
A and B, of N/2 spins each, and to only consider singlets (i, j) in which
the first index i ∈ A and the second j ∈ B [18–20]. In the case of a bipartite
lattice, these groups are naturally the two sublattices, as shown in Fig. 1. This
restricted VB basis has (N/2)! states and is still massively overcomplete—the
singlet space has N !/[(N/2)!]2(N/2+1) dimensions [18]. The VB basis states
are all non-orthogonal, overlapping with each other according to the simple
loop rule illustrated in Fig. 1.
The VB basis was introduced already in the early 1930s [18,21,22] and
has played an important role in exactly solvable models [18,23–25]. Later,
it became a tool for describing spin liquids—the resonating valence bond
(RVB) mechanism introduced by Fezekas and Anderson [26–28], in which the
ground state is dominated by short valence bonds. In exact diagonalization
Monte Carlo Simulations of Quantum Spin Systems 3
studies, the VB basis is useful in cases where it is a good approximation to
only consider a restricted (and incomplete) space of short bonds (spin liquids
and other states with no magnetic long-range order) [29–32]. Variational cal-
culations in the VB basis have been carried out for the 2D Heisenberg model
[20,33]. Furthermore, Liang realized that a variational VB state could be con-
siderably improved by stochastically projecting it with an operator (−H)m
for large m [34]. Later, Santoro et al. devised a Green’s function method for
calculating energies in the VB basis [35]. Despite the promising results ob-
tained in these studies, there was, to our knowledge, no further developments
of QMC methods in the VB basis until one of us recently introduced two
related projector algorithms [36], improving on the schemes of Liang [34] and
Santoro et al. [35]. These algorithms have already been applied in studies of
quantum phase transitions [37,38] and entanglement entropy [39].
Some previously unnoticed advantages of the VB basis in QMC algorithms
were pointed out in Ref. [36]. Here we summarize our recent work on VB
projector methods and highlight some of their unique features. We discuss
in particular a scheme for “self-optimizing” the trial state out of which the
ground state is projected, and also show how to study properties of triplet
excitations at finite momentum.
2 Ground state projection
Consider a state |Ψ〉 and its expansion in terms of eigenstates |n〉, n = 0, 1, . . .,
of some hamiltonian H ;
|Ψ〉 =
cn|n〉. (4)
With C a constant chosen such that the lowest eigenvalue E0 − C is the
largest in magnitude, a large number m of repeated operations with C −H
projects out the ground state,
(C −H)m|Ψ〉 → c0(C − E0)m
|0〉+ c1
C − E1
C − E0
|1〉+ . . .
, (5)
provided that the overlap c0 6= 0. Here we will first consider singlet eigenstates
of the Heisenberg model (1), which can be expanded in VB states;
|Ψ〉 =
fi|Vi〉. (6)
Because of the overcompleteness of the VB basis, this expansion is not unique.
That, however, does not prohibit the ground state |0〉 to be projected out
according to Eq. (5). The Heisenberg hamiltonian can be written in terms
of singlet projection operators Hb ≡ Hi(b),j(b) on the interacting spin pairs,
{i(b), j(b)}, b = 1, . . . , Nb = DN (for a periodic cubic D-dimensional lattice);
H = −J
Hb = −J
Hi(b),j(b), Hij = −(14 − Si · Sj). (7)
4 A. W. Sandvik and K. S. D. Beach
When a singlet projector Hij acts on a VB basis state, one of two things
can happen; 1) if i, j belong to the same bond the state is unchanged with a
matrix element 1, or 2) if they belong to different bonds these two bonds are
reconfigured (“flipped”) with matrix element 1/2;
Hij | · · · (i, j) · · ·〉 = | · · · (i, j) · · ·〉, (8)
Hij | · · · (i, k) · · · (l, j) · · ·〉 = 12 | · · · (i, j) · · · (l, k) · · ·〉. (9)
Here positive-definitness of (9) is directly related to the two sites i and j being
in different sublattices. For a frustrated system, where there are operators
with both sites i, j in the same sublattice, positive-definitness does not hold
[19,20]. For a non-frustrated system the simple bond flip (9) makes for a
convenient stochastic implementation of the ground state projection (5). We
write the projection operator as (with J = 1 henceforth)
(C −H)m =
Pr, Pr = Hbr
· · ·Hbr
, (10)
where we have introduced a compact notation Pr, r = 1, . . . , N
b , for the
different strings of singlet projectors. When a string Pr acts on a VB basis
state |V 〉 the result is another basis state, which we denote |V (r)〉, with
a prefactor (weight) Wr which is simply given by the number moff of off-
diagonal operations in the course of evolving |V 〉 to |V (r)〉;
Pr|V 〉 = Wr|V (r)〉, Wr = 2−moff . (11)
We here first consider projecting the ground state out of a single VB
basis state; later we will consider the use of a more complicated trial state.
We consider two ways to calculate expectation values:
〈H〉 =
r〈Ψ |HPr|V 〉
r〈Ψ |Pr |V 〉
r Wr〈Ψ |H |V (r)〉
r Wr〈Ψ |V (r)〉
, (12)
〈A〉 =
rl〈V |P ∗l APr|V 〉
rl〈V |P ∗l Pr|V 〉
rl WrWl〈V (l)|A|V (r)〉
rl WrWl〈V (l)|V (r)〉
. (13)
We will discuss how to estimate these using importance sampling; terms
(configurations) of the numerators are illustrated in Fig. 2. We will refer to
(12) and (13) as the single and double projection, respectively.
In (12), which is an exact (when m → ∞) expectation value only of the
hamiltonian (or other operators for which the ground state is an eigenstate)
the state |Ψ〉 is in principle arbitrary. It is very convenient to use a state
which has equal overlaps with all the VB basis states, e.g., the Néel state
|ΨN 〉 (all spins up on sublattice A and down on B). It is easy to see that
〈ΨN |V 〉 = 2−N/2 for any basis state |V 〉. Since H |V (r)〉 = −
b Hb|V (r)〉 is
a sum of basis states multiplied by factors −1 or −1/2, the overlaps with 〈ΨN |
drop out altogether and do not have to be considered further. If the projector
Monte Carlo Simulations of Quantum Spin Systems 5
|V ><V | A
(a) (b)
|V >H
Fig. 2. Propagation of a VB state on a 6-site chain. The horizontal bars represent
nearest-neighbor Heisenberg interactions (singlet projectors). In the single projec-
tion (a) the state is propagated from right to left, and an estimator for the ground
state energy is obtained by acting once more with all terms of the hamiltonian. In
the double projection (b) the state is projected from the right and the left, and
any operator expectation value can be estimated by calculating the corresponding
matrix elements between the propagated states
strings Pr in (12) are importance-sampled according to their weights Wr, the
estimator for the ground state energy is thus
E0 = 〈H〉 = −〈md + 12mo〉, (14)
where md and mo are, respectively, the number of diagonal and off-diagonal
operations Hb|V (r)〉 (and md +mo = Nb). It should be noted that although
this estimator is exact in the limit m → ∞, it is not variational. The correct
energy may thus be approached with increasing m from above or below.
Eq. (13) is valid for any expectation value and in the case of A = H gives
a variational estimate of the energy. Using WrWl〈V (l)|V (r)〉 as the sampling
weight, the estimator for any 〈A〉 is of the form
〈A〉 =
〈V (l)|A|V (r)〉
〈V (l)|V (r)〉
. (15)
In the case of a spin correlation function 〈Si·Sj〉, the matrix element is related
to the loop structure of the overlap graph [19,20] (illustrated in Fig. 1):
〈Vl|Si · Sj |Vr〉
〈Vl|Vr〉
+3/4, if i, j ∈ same loop, same sublattice,
−3/4, if i, j ∈ same loop, different sublattices,
0, if i, j ∈ different loops.
Measuring the spin correlations is hence straight-forward once the overlap-
loops have been constructed. Higher-order functions, e.g., dimer-dimer cor-
relations 〈(Si · Sj)(Sk · Sl)〉, are also related to the loop structure [40].
Note again that no bond operator Hb can destroy a VB state and that
all the states have non-zero overlap with each other. Thus all terms in (13)
contribute to the expectation value. This turns out to be an advantage in
constructing a Monte Carlo algorithm, as any change made in the operator
strings can be accepted with some probability. With an orthogonal basis, such
6 A. W. Sandvik and K. S. D. Beach
as the Szi eigenstates, there would be considerable constraints, both in terms
of individual operators in the projection [the spin flip operator in (1) can act,
without destroying the state, only on anti-parallel spins], and in ensuring a
non-zero overlap between the propagated states (the two propagated states
have to be identical). Note also that the singlet projectors are non-hermitian
in the VB basis. As indicated in (13), and illustrated in Fig. 2, we here
propagate two states, |V (l)〉 ∝ Pl|V 〉 and |V (r)〉 ∝ Pr|V 〉, and subsequently
compute their overlap and various matrix elements. Propagating |V 〉 with
P ∗l Pr and then taking the overlap with |V 〉 is not equivalent term-by-term.
To carry out the projection stochastically, the operator strings are stored
in arrays [Pα] = [b
1 ][b
2 ] · · · [bαm], where bαi ∈ {1, . . . , Nb} with α = 1, 2 cor-
responding to Pr and Pl, respectively, in the double projection; in the single
projection α is redundant. A table holds the site pairs i(b), j(b). The state
|V 〉 is stored in a list [V ] = [v1][v2] · · · [vN ] where vi = j and vj = i if there is
a valence bond at (i, j). Propagation with the bond flips (9) is easily carried
out in this representation. The state list [V ] is then first copied into two lists,
[V1] and [V2], in which |V (r)〉 and |V (l)〉 are constructed.
The simulation can be started with a randomly generated operator string.
The strings can can be updated in an trivial way, by changing a number R
of operators at random. In either [P1] or [P2], R positions pj , j = 1, . . . , R,
pj ∈ {1, . . . ,m} (all different) are generated. Their contents bαpj are picked
randomly from the set {1, . . . , Nb} (excluding the old value for each pj). To
calculate the Metropolis acceptance probability, the state is propagated with
the updated operator string and the number of off-diagonal operationsmoff in
(11) is counted. In the single expansion, the acceptance probability is simply
Paccept = min
W newr
W oldr
= min[2m
−mnew
off , 1], (17)
whereas in the double projection an overlap ratio appears as well. In the
double projection, we change operators only in one of the operator strings at
a time, so that only one state has to be propagated. It is customary to define
a size-normalized Monte Carlo step (or “sweep”). For projector length m we
do m replacement attempts, and so our sweep is independent of the number
of operators R replaced in each update.
Normally, in Monte Carlo simulations one does not compute the full
weight, because it is possible to more speedily calculate just the change in
the weight [the weight ratio in (17)]. In the present formulation of the VB
projector algorithm the weight is, however, recomputed from scratch each
time, because a better scheme is not yet known. Each update hence requires
on the order of m operations. In the double expansion, construction of the
loops needed to compute the overlap scales as N , but typically m > N and
the propagation of the state dominates the simulation. In spite of the need
to recalculate the weight, the scheme is sufficiently efficient to compete with
other ground state QMC methods. More importantly, as we will discuss in
Sec. 4, the VB basis offers access to quantities out of reach for other methods.
Monte Carlo Simulations of Quantum Spin Systems 7
1 2 3 4 5
1 2 3 4 5
(a) (b)
Fig. 3. (a) Acceptance rate in a double projection and (b) the number of bonds
changed in the propagated state in accepted updates for a 16× 16 lattice (symbols
with lines; the bare lines are for L = 8) versus the projection length (normalized
by the number of sites N). The trial state |V 〉 was a columnar dimer state
The optimum number R of operators to replace depends on the accep-
tance rate. In Fig. 3(a) we show the acceptance rate for a double-projected
2D system versus the length m of the projection for R = 1− 6. As expected,
the acceptance rate decreases with increasing R, but it does not change ap-
preciably with m. It also depends only weakly on the lattice size. Multiplying
the acceptance rate by R gives the average number of operators changed; it
initially increases with R but has a maximum for R = 7 (for L = 16). The
optimum R is clearly model/lattice dependent. Another characteristic of the
update is the number of bonds changed in the projected state |V (r)〉 as a
consequence of the modifications of Pr. This number is shown in Fig. 3(b).
It is seen to increase with R, as expected. As a function of m the num-
ber of changed bonds decreases. This behavior reflects a loop structure of
the singlet-projection operators [41,8], which implies that some changes “up-
stream” in the operator string may be healed further downstream in the
propagation. For a finite lattice in the limit m → ∞ one would expect a
substitution of an operator far upstream in Pr to have no effect on the final
propagated state |V (r)〉. This does not imply that this update is inconsequen-
tial, as the sampling is over paths, not just the final states in the propagations.
In principle, it should be possible to take advantage of the underlying loop
structure of the singlet projectors [41,8] to devise a loop update in the VB
basis, analogous to such updates in world-line [7,8] and SSE [9] simulations.
However, we have not yet been able to construct a scheme that is in practice
faster than the trivial random substitution with full state propagation.
3 Self-optimized trial state
So far, we have projected the ground state out of a single VB basis state |V 〉.
This works well [36], but the rate of convergence with m of course depends
8 A. W. Sandvik and K. S. D. Beach
on the state chosen; ideally one would like to maximize the overlap 〈V |0〉.
One way to obtain a typically good single-configuration trial state is to first
start with an arbitrary one; a regular bond pattern or a randomly generated
configuration. After carrying out some projection steps with this state, the
current propagated state |V (r)〉 is chosen as a new trial state. Since this state
has been generated in the projection it should contribute substantially to the
ground state and hence typically will be better than a completely arbitrary
one. However, as we will discuss next, we can do much better than this.
Liang’s original motivation for introducing a projector technique in the
VB basis was to improve on a variational calculation [34]. Previously, Liang,
Doucot, and Anderson had studied a variational amplitude-product state for
the 2D Heisenberg model of the form [20]
|Ψ〉 =
fk|Vk〉, fk =
h(xbk, ybk), (18)
where xbk and ybk are the x- and y-lengths of bond b in VB state k, as il-
lustrated in Fig. 4(a). Liang et al. tried power-law and exponential forms
depending only on the total length l of the bonds [the ”Manhattan” length
l = x + y was used, but defining l = (x2 + y2)1/2 should not change things
qualitatively], in addition to keeping several short-bond amplitudes as pa-
rameters to optimize [20]. More recently, all the amplitudes were optimized
without any assumed form on lattices with up to 32×32 sites, with the result
that h(l) ∼ l−3 for long bonds [33]. One of us has also recently showed the
more general result h(l) ∼ l−(D+1) within a mean-field approach for a D-
dimensional cubic lattice [42]. In 2D the fully-optimized amplitude-product
state turns out to be extremely good, with an energy deviating by only
≈ 0.06% from the exact ground state energy and with the long-distance
spin-spin correlations reproduced to within 2% [33].
With the state (18) an expectation value is given by
〈A〉 =
kp fkfp〈Vp|Vk〉
〈Vp|A|Vk〉
〈Vp|Vk〉
rl fkfp〈Vp|Vk〉
, (19)
which can be evaluated using importance sampling of the VB configurations
with weight fkfp〈Vp|Vk〉. Liang et al. introduced a very simple scheme for
updating the dimer configurations [20], which we here illustrate in Fig. 4(b).
Choosing two next-nearest-neighbor sites, i.e., ones on a diagonal of a 4-site
plaquette (or, in principle, any two sites on the same sublattice), the two
bonds connected to them are reconfigured in the only possible way which
maintains only bonds between the A and B sublattices, as shown in the figure.
Labeling the two initially chosen sites 1 and 2, and the bonds connected to
them b = 1, 2, the Metropolis acceptance probability is, assuming that the
bond update was made in |Vk〉, resulting in |Vk′ 〉,
Paccept = min
h(x1k′ , y1k′)h(x2k′ , y2k′)
h(x1k, y1k)h(x2k, y2k)
〈Vp|Vk′ 〉
〈Vp|Vk〉
. (20)
Monte Carlo Simulations of Quantum Spin Systems 9
(a) (b)
Fig. 4. (a) Definition of the size of a bond b. (b) Reconfiguration of two bonds in
an update of the trial state. The two sites marked 1, 2 are chosen at random among
all pairs of next-nearest neighbor sites.
We can use an amplitude-product state as the trial state in the projector
QMC method, using some set of amplitudes not necessarily originating from
a variational calculation. In updating the bond configurations, we then also
must compute the new weight of the propagation; Pr|Vk′ 〉 = Wk′r|Vk′ (r)〉, and
of course 〈Vp|Vk〉 in (20) is replaced by 〈Vp(l)|Vk(r)〉 (in the case of the double
projection; for the single projection there is no overlap). The acceptance rate
of a state update is similar to that of an operator string update with a smallR,
and often we find it advantageous to combine state and operator updates. To
save some time, one can tentatively accept/reject a bond update based solely
on an amplitude ratio—Eq. (20) without the overlaps—and then calculate
the overlap and the propagation weight for a final accept/reject probability
only for tentatively accepted bond updates.
In the variational calculation 〈H〉 is minimized with respect to all h(x, y).
With a recently developed stochastic optimization method [33], all the ∝ N
amplitudes can be minimized for moderate-size lattices (up to 32 × 32 sites
were considered in Ref. [33]). In principle we could follow Liang [34] and use
the best possible variational state as our trial state in the projector QMC
method—indeed this can be expected to be the optimum starting point.
However, we will now describe a scheme which delivers a trial state nearly as
good as the best variational state, at a smaller computational cost.
Consider the probability distribution P (x, y) of valence bonds. In an
amplitude-product state (18) we would have P (x, y) ∼ h(x, y), were it not
for the “hard-core” constraint of only one bond per site. Even with this con-
straint, it is clear that the probabilities and amplitudes are related in a mono-
tonic way; increasing h(x, y) for some given (x, y), while keeping the other
amplitudes fixed, will lead to a larger P (x, y). This fact can be exploited
in constructing a good trial state. We define two different probability dis-
tributions, P0(x, y) and Pm(x, y), the former being the just discussed bond
probability in a trial state of the form (18) and the latter the probability
distribution in the projected state. For sufficiently large m, Pm is an exact
property of the ground state, whereas P0 is a property of the trial state and is
in general different from Pm. However, for given m, we can adjust the ampli-
tudes h(x, y) of the trial state such that Pm(x, y) = P0(x, y) for all x, y. If this
is done for m sufficiently large, then our trial state has a bond distribution
identical to that of the exact ground state. Such a state is often almost as
good as the best variational state. The reason that this is useful in practice
10 A. W. Sandvik and K. S. D. Beach
0.0 0.5 1.0 1.5 2.0 2.5
1.0 1.5 2.0 2.5
0.650
0.655
0.660
0.665
0.670
 self-optimized
 projected random
 columnar dimer
0.0 0.5 1.0 1.5 2.0
0.6694
0.6695
0.6696
0.6697
0.6698
0.6699
0.6700
 Projected state
 Self-optimized trial state
Fig. 5. Left: Energy versus projection length for a 16 × 16 system, using three
different trial states (labeled in the inset). Right: The energy of the self-optimized
trial state compared with the projected energy using that trial state. The dashed
line shows the energy evaluated independently using the SSE method [43]
is that it is very easy to adjust the amplitudes to achieve self-consistency.
Because of the monotonous relationship between h(x, y) and P0(x, y), we can
simply increase h(x, y) by some amount if P0(x, y) < Pm(x, y) and decrease it
if P0(x, y) > Pm(x, y), and repeat this until self-consistency is achieved. We
use the following scheme to update the amplitudes after the kth iteration;
ln[h(x, y)] → ln[h(x, y)] + RAN · β(k) · sign[Pm(x, y) − P0(x, y)], (21)
where RAN is a random number in the range [0, 1) and β(k) decreases with
the iteration step k = 1, 2, . . ., according to β(k) ∝ k−α. For the exponent,
we typically use α = 3/4. To evaluate the probabilities P0 and Pm (in two
independent simulations; with and without projection of the trial state), the
number of Monte Carlo sweeps does not have to be very large, because we
only need the sign of the difference of the two probabilities. We normally
use on the order of 100-1000 sweeps per iteration. Even if the stochastically
evaluated sign in (21) occasionally may be wrong, it is correct on average
and the amplitudes typically converge to a self-consistent solution after a
few hundred iterations. Due to the stochastic nature of the procedure, self-
consistency of course obtains only to within some statistical error, which can
be reduced by increasing the number iterations and/or sweeps per iteration.
Fig. 5 shows results for the energy of a 16× 16 system obtained in double
projections with three different trial states; a columnar dimer state, a ran-
domly picked state generated while projecting the columnar state, as well as
the self-optimized state. Already for the shortest projection, m = N/8 = 32,
the self-optimized state gives a projected energy which deviates by only 0.06%
from the exact ground state energy, and for largerm the energy is exact within
statistical errors. The other two trial states also lead to the correct energy
but only for much larger m. The energy of the self-optimized trial state itself
is shown in the right panel—its error is as small as that of the best variational
amplitude-product state [33].
Monte Carlo Simulations of Quantum Spin Systems 11
0.0 0.5 1.0 1.5 2.0
 self-optimized
 random projected
 columnar dimer
0.0 0.5 1.0 1.5 2.0
0.116
0.118
0.120
0.122
0.124
0.126
 Projected state
 Self-optimized trial state
Fig. 6. Long-distance spin correlations in the same runs as in Fig. 5. A different
“random projected” state is picked for each m, which leads to an un-smooth curve.
Fluctuations beyond the small statistical errors in the right panel reflect differences
in how closely the rather short self-optimization runs have approached P0 = Pm
Fig. 6 shows the long-distance spin correlation calculated in the same
runs. Again, the self-optimized state delivers superior results, although here
the convergence is not as fast as for the energy. The error of the spin correla-
tion in the trial state (right panel) is about 4% for large m, which is twice the
error in the best variational state [33]. Thus, the self-optimized state is not
identical to the best variational state, and an even faster convergence could
be achieved by using the fully optimized variational state. However, the vari-
ational calculation is much more time consuming than the self-optimization.
We can go beyond the amplitude-product state by taking into account
bond correlations. We are currently exploring this with both variational and
self-optimized states.
4 Triplet excitations
A unique advantage of the VB basis is that an mz = 0 triplet state can
be projected simultaneously with the singlet, with essentially no additional
overhead. Any triplet can be expanded in VB states where one of the bonds
corresponds to a triplet [18]; (i, j) → [i, j], where
[i, j] = (| ↑i↓j〉+ | ↓i↑j〉)/
2. (22)
Formally, such a triplet can be created by acting on a singlet with Szi − Szj ;
(Szi − Szj )(i, j) = [i, j]. (23)
To create a triplet |τ(q)〉 with some momentum q, we can apply
Szq =
eiq·rjSzj (24)
12 A. W. Sandvik and K. S. D. Beach
to a singlet state |σ(0)〉 with zero momentum;
Szq|σ(0)〉 ∝ |τ(q)〉. (25)
The amplitude-product state (18) for a periodic-boundary system has q = 0
if the number of bonds, N/2, is even, whereas for odd N/2 it has q = (π, π).
This simply follows from the fact that sublattice A → B and B → A when
translating by one lattice constant, whence each singlet (3) acquires a minus
sign. We typically work with systems with even N/2 (e.g., L×L lattices with
even L) and so we will here consider q = 0 singlets.
At the antiferromagnetic wave-vector, Q = (π, π), SzQ acting on an ar-
bitrary VB basis state |V 〉 can be written as a sum of N/2 terms of the
form (Szi − Szj )|V 〉, with i, j corresponding to the sites connected by bonds.
Operating on a q = 0 singlet thus gives
SzQ|σ(0)〉 = SzQ
fk|(i1, j1)k(i2, j2)k · · · (iN/2, jN/2)k〉
fk|(i1, j1)k · · · [ib, jb]k · · · (iN/2, jN/2)k〉, (26)
where the unspecified coefficients fk, e.g., the amplitude products in (18),
have translational invariance built in. Thus, for a triplet with q = Q the
wave function phases are buried in the definition of the singlets. Often, the
lowest excitation of a Heisenberg system is a q = Q triplet, which we thus
can sample without any difficulties with signs or phases. We consider this
case first, before turning to triplets with arbitrary momentum.
There are two possible actions of a singlet projectorHij on a triplet bond;
Hij | · · · [i, j] · · ·〉 = 0, (27)
Hij | · · · [i, k] · · · (l, j) · · ·〉 = 12 | · · · (i, j) · · · [l, k] · · ·〉, (28)
i.e., a diagonal operation on a triplet bond destroys the state whereas an
off-diagonal operation on one triplet and one singlet bond creates a singlet
at the sites on which the operator acts and moves the triplet to the other
two sites involved. Importantly, the matrix element remains the same as in
the off-diagonal operation on two singlet bonds. The weight of a triplet path
is therefore the same as the corresponding singlet path (11), except that the
triplet dies (giving zero weight) if an operator in Pr acts diagonally on the
triplet bond. We can thus measure triplet properties using paths generated in
a singlet simulation, by considering only those triplet paths that do survive
the propagation. In the trial state we have N/2 possible locations of the
triplet. All of them can be attempted collectively in a single propagation, by
keeping counters t(i) for the number of surviving states in which the triplet is
connected to site i (with i on sublattice A). Initially t(i) = 1 for all i. During
the propagation, for each diagonal operation (27) t(i) → t(i) − 1, and for
each off-diagonal operation (28) t(i) → t(i) − 1, t(l) → t(l) + 1. Eventually,
Monte Carlo Simulations of Quantum Spin Systems 13
as m → ∞, all triplets die; t(i) = 0 for all i, but typically there are enough
survivors left at large enough m to compute converged triplet properties.
An added advantage of calculating singlet and triplet properties in the
same run is that there are error cancellations which in some cases can increase
the statistical precision of differences, e.g., the singlet-triplet gap [36],
∆ = ET (π, π)− ES(0, 0), (29)
by up to orders of magnitude relative to two independent calculations. The
triplet energy ET (π, π) can be estimated using (15), taking into account that
a diagonal operation on a triplet bond gives zero, i.e., nd → nd − nt, where,
for surviving triplet configurations, nt = 0, 1 is the number of triplet bonds
of length 1. Other triplet properties have been discussed in Ref. [40].
We now discuss calculations with triplets of arbitrary momentum q. The
energy can be evaluated according to
E(q) =
r〈σ(0)|Sz−qHPrSzq|σ(0)〉
r〈σ(0)|Sz−qPrSzq|σ(0)〉
. (30)
We want to evaluate this expression using the sampled q = 0 singlet bond
configurations, and so we rewrite it as
E(q) =
r〈σ(0)|Sz−qHPrSzq|σ(0)〉
r〈σ(0)|Pr |σ(0)〉
r〈σ(0)|Sz−qPrSzq|σ(0)〉
r〈σ(0)|Pr |σ(0)〉
. (31)
The two factors can be evaluated based on sampling the propagations Pr and
the amplitude-product state |σ(0)〉 [which we have not explicitly written as
a sum of bond configurations in (30) and (31)]. In the singlet energy (12), we
could pick the Néel state for 〈Ψ | and then obtained the very simple expression
(15). Now we must consider the overlap with a momentum q triplet state.
We use 〈σ(0)|Sz−q, but in (31) have rewritten E(q) so that only the overlap
with 〈σ(0)| has to be considered for the sampling weight. Phases arising from
Szq only appear in the measurements, but in the end we have to evaluate
the ratio of the two quantities in (31), which can be challenging in practice.
However, close to q = (π, π) and (0, 0) [q 6= (0, 0)] we find that it can be
done; in some cases the method works even far away from these momenta.
We define the dispersion relative to the gap (29) at (π, π);
ω(q) = ET (q)−∆. (32)
Fig. 7 shows results for∆ and ω(q1), where q1 is the momentum closest to but
not equal to (π, π); q1 = (π−2π/L, π). We show the convergence as a function
of N/m for 4× 4 and 16× 16 lattices, comparing with exact diagonalization
results in the former case. From (5) one would expect the convergence to
be asymptotically exponential, which is seen clearly for L = 4. For L = 16,
the three largest-m points are equal within statistical errors, suggesting that
these results are also close to converged. Using ω(q1) = 0.62 for the L = 16
14 A. W. Sandvik and K. S. D. Beach
0.578
0.580
0.582
0.584
0.586
0 0.5 1 1.5 2
1.760
1.765
1.770
1.775
1.780
0 0.5 1 1.5 2
L=4 L=16
Fig. 7. Projection-length convergence of the singlet-triplet gap and the excitation
energy at momentum q1 = (π − 2π/L, π) for L = 4 (left) and 16 (right) lattices.
The horizontal lines show the exact results for L = 4
system gives the spin-wave velocity c = 1.58, in very close agreement with
the known value [43].
Another useful quantity accessible with the VB projector is the matrix
element 〈τ(q)|Szq|σ(0)〉, the square of which gives the single-magnon weight
in the dynamic structure factor (which is experimentally measurable using
neutron scattering). It can be calculated in a way similar to E(q), and we have
done so successfully for q close to (π, π). Results will be presented elsewhere.
Finally, we also note that the triplet bond-length distribution gives a di-
rect, albeit basis dependent, window into the “spinon” aspects [44] of the
excitations. Spinon deconfinement should be manifested as a delocalized dis-
tribution function, whereas two spinons bound into a magnon or “triplon”
should be reflected in a well-defined peak in the distribution function. We are
currently exploring this.
Acknowledgments
We would like to thank H. G. Evertz for useful discussions. This work was sup-
ported by the National Science Foundation under grant No. DMR-0513930.
References
1. M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976); M. Suzuki, S. Miyashita, and
A. Kuroda, Prog. Theor. Phys. 58, 1377 (1977)
Monte Carlo Simulations of Quantum Spin Systems 15
2. M. Barma and B. S. Shastry, Phys. Rev. B 18, 3351 (1977)
3. J. E. Hirsch, R. L. Sugar, D. J. Scalapino and R. Blankenbecler, Phys. Rev. B
26, 5033 (1982)
4. J. J. Cullen and D. P. Landau, Phys. Rev. B 27, 297 (1983)
5. B. B. Beard and U. -J. Wiese, Phys. Rev. Lett. 77, 5130 (1996)
6. N. V. Prokofév, B. V. Svistunov, and I. S. Tupitsyn, Pisḿa Zh. Eks. Teor. Fiz.
64, 853 (1996) [JETP Lett. 64, 911 (1996)]; Zh. Eks. Teor. Fiz. 114, 570 (1998)
[JETP 87, 311 (1998)]
7. H. Evertz, G. Lana and M. Marcu, Phys. Rev. Lett. 70, 875 (1993)
8. H. G. Evertz, Adv. Phys. 52, 1 (2003)
9. A. W. Sandvik, Phys. Rev. B 59, R14157 (1999)
10. O. F. Syljůasen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002)
11. D. C. Handscomb, Proc. Cambridge Philos. Soc. 58, 594 (1962); 60, 115 (1964)
12. D. H. Lee, J. D. Joannopoulos, and J. W. Negele, Phys. Rev. B 30, 1599 (1984)
13. A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991)
14. A. W. Sandvik, Phys. Rev. B 56, 11678 (1997)
15. V. N. Kotov, O. Sushkov, Z. Weihong, and J. Oitmaa, Phys. Rev. Lett. 80,
5790 (1998)
16. N. Trivedi and D. M. Ceperley, Phys. Rev. B 40, 2737 (1989)
17. S. Sorella and L. Capriotti, Phys. Rev. B 61, 2599 (2000)
18. L. Hulthén, Ark. Mat. Astron. Fys. 26A, No. 11 (1938)
19. B. Sutherland, Phys. Rev. B 37, 3786 (1988)
20. S. Liang, B. Doucot, and P. W. Anderson, Phys. Rev. Lett. 61, 365 (1988)
21. G. Rumer, Göttingen Nachr. Tech. 1932, 377 (1932)
22. L. Pauling, J. Chem. Phys. 1, 280 (1933)
23. C. K. Majumdar and D. K. Ghosh, J. Math. Phys. 10, 1388 (1969)
24. B. S. Shastry and B. Sutherland, Physica B&C 108, 1069 (1981)
25. I. Affleck, T. Kennedy, E. H. Lieb, and H. Tasaki, Phys. Rev. Lett. 59, 799
(1987)
26. P. Fezekas and P. W. Anderson, Philos. Mag. 30, 23 (1774)
27. P. W. Anderson, Science 235, 1196 (1987)
28. E. Fradkin, Field Theories of Condensed Matter Systems (Adison-Wesley 1991).
29. P. L. Iske and W. J. Caspers, Physica 142A, 360 (1987)
30. M. Kohmoto, Phys. Rev. B 37, 3812 (1988)
31. S. Tang and H.-Q. Lin, Phys. Rev. B 38, 6863 (1988)
32. M. Mambrini, A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422
(2006)
33. J. Lou and A. W. Sandvik, cond-mat/0605034
34. S. Liang, Phys. Rev. B 42, 6555 (1990); Phys. Rev. Lett. 64, 1597 (1990).
35. G. Santoro, S. Sorella, L. Guidoni, A. Parola, and E. Tosatti, Phys. Rev. Lett.
83, 3065 (1999)
36. A. W. Sandvik, Phys. Rev. Lett. 95, 207203 (2005)
37. A. W. Sandvik, cond-mat/0611343
38. A. W. Sandvik and K. S. D. Beach, cond-mat/0612126
39. F. Alet, S. Capponi, N. Laflorencie, and M. Mambrini, cond-mat/0703027
40. K. S. D. Beach and A. W. Sandvik, Nucl. Phys. B 750, 142 (2006)
41. M. Aizenman and B. Nachtergaele, Commun. Math. Phys. 164, 17 (1994)
42. K. S. D. Beach (unpublished).
43. A. W. Sandvik. Phys. Rev. B 56, 11678 (1997)
44. S. Sachdev. Rev. Mod. Phys. 75, 913 (2003)
