Quantum Phase Transition in the Four-Spin Exchange Antiferromagnet
Valeri N. Kotov, Dao-Xin Yao, A. H. Castro Neto, and D. K. Campbell
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215
We study the S=1/2 Heisenberg antiferromagnet on a square lattice with nearest-neighbor and plaquette four-
spin exchanges (introduced by A.W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).) This model undergoes a
quantum phase transition from a spontaneously dimerized phase to Néel order at a critical coupling. We show
that as the critical point is approached from the dimerized side, the system exhibits strong fluctuations in the
dimer background, reflected in the presence of a low-energy singlet mode, with a simultaneous rise in the triplet
quasiparticle density. We find that both singlet and triplet modes of high density condense at the transition,
signaling restoration of lattice symmetry. In our approach, which goes beyond mean-field theory in terms of the
triplet excitations, the transition appears sharp; however since our method breaks down near the critical point,
we argue that we cannot make a definite conclusion regarding the order of the transition.
PACS numbers: 75.10.Jm, 75.30.Kz, 75.50.Ee
I. INTRODUCTION
Problems related to quantum criticality in quantum spin
systems are of both fundamental and practical importance1.
Numerous materials, such as Mott insulators, exhibit either
antiferromagnetic (Néel) order or quantum disordered (spin
gapped) ground state depending on the distribution of Heisen-
berg exchange couplings and geometry. External perturba-
tions (such as doping or frustration) can also cause quantum
transitions between these phases. Systems with spin 1/2 are
indeed the most interesting as they are the most susceptible to
such transitions. It is well understood that the quantum transi-
tion between a quantum disordered and a Néel phase is in the
O(3) universality class1, where a triplet state condenses at the
quantum critical point (QCP).
A recent exciting development in our theoretical under-
standing of QCPs originated from the proposal that if the
quantum disordered (QD) phase spontaneously breaks lattice
symmetries (e.g. is characterized by spontaneous dimer or-
der), and the transition is of second order, then exactly at
the QCP spinon deconfinement occurs, i.e. the excitations
are fractionalized2. It is assumed that the Hamiltonian itself
does not break the lattice symmetries (i.e. does not have “triv-
ial” dimer order caused by some exchanges being stronger
than the others). We use the terms “dimer order” and “va-
lence bond solid (VBS) order” interchangeably. It is expected
that the dimer order vanishes exactly at the point where Néel
order appears, i.e. there is no coexistence between the two
phases. Deconfinement thus is intimately related to disap-
pearance of VBS order; indeed if the latter persisted in the
Néel phase it would be impossible to isolate a spinon, as
“pairing” would always take place. Spontaneous VBS or-
der driven by frustration has been a common theme in quan-
tum antiferromagnetism3, although its presence and the na-
ture of criticality in specific models, such as the 2D square-
lattice frustrated Heisenberg antiferromagnet, is still some-
what controversial4. It would be particularly useful to apply
unbiased numerical approaches, such as the Quantum Monte
Carlo (QMC) method, to study frustrated spin models; how-
ever due to the fatal “sign” problem5, frustrated Heisenberg
systems are beyond the QMC reach.
In a recent study, the QMC method was applied to a four-
spin exchange quantum spin model without frustration, which
was shown to exhibit columnar dimer VBS order and a mag-
netically ordered phase with a deconfined QCP separating
them6. These conclusions were later confirmed by further
QMC studies7. Extensions of the model, which include for
example additional (six-spin) interactions, provide additional
support for a continuous QCP8. A different VBS pattern (pla-
quette order) was also proposed for the four-spin exchange
model9. At the same time, the nature of the quantum phase
transition was challenged in Refs.[10,11], where arguments
were given that the transition is in fact of (weakly) first order.
It is the objective of the present work to study the Sand-
vik model6, by approaching the quantum transition from the
dimer VBS phase. Our approach uses as a starting point a
symmetry broken state (i.e. one out of four degenerate VBS
configurations), and we thus must search for signatures that
the system attempts to restore the lattice symmetry at the QCP.
Even though full restoration is impossible within the present
framework, we find a QCP characterized by condensation of
triplet modes of high density; this is in contrast to the conven-
tional situation when the condensing particles are in the dilute
Bose gas limit. The high density itself is due to the presence
of a singlet mode that condenses at the QCP, and reflects the
strong fluctuations of the background dimer order. The above
effects lead to the vanishing of the VBS order parameter; at
the same time our method, which accounts for the strong fluc-
tuations, leads to a rather sharp phase transition. It appears
that we cannot draw a definite conclusion about the order of
the transition because in the vicinity of the QCP the triplon
density increases uncontrollably, suggesting that other states
(such a plaquette states and larger clusters) are strongly ad-
mixed into the ground state. This is generally expected in a
situation where the lattice symmetry is restored at the quan-
tum critical point.
The model under consideration is
H = J
〈a,b〉
Sa.Sb −K
a,b,c,d
(Sa · Sb)(Sc · Sd), (1)
where J > 0,K > 0, and all spins are S = 1/2. Consider
the numbers 1,2,3,4 in Fig. 1. The summation in the four-
spin term is over indexes (a, b) = (1, 2), (c, d) = (3, 4) and
http://arxiv.org/abs/0704.0114v3
FIG. 1: (Color online) Dimer pattern in the quantum disordered
(VBS) phase, K/J > (K/J)c.
(a, b) = (1, 4), (c, d) = (2, 3) on a given plaquette, and then
summation is made over all plaquettes12. The range of pa-
rameters explored in Ref.[6] is K/J ≤ 2, and the QCP is at
(K/J)c ≈ 1.85. Our coupling notation is slightly different
from the one used in Refs.[6,7]; the coupling K is related to
the parameter Q6,7 via K = Q/(1 + Q/(2J)), and the criti-
cal point in that notation is (Q/J)c ≈ 25. The dimerization
pattern is proposed to be of the “columnar” type, as shown in
Fig. 1. Four such configurations exist. We will assume a con-
figuration of this type, will show that it is stable at K/J ≫ 1,
and will then search for an instability towards the Néel state
as K/J decreases.
The rest of the paper is organized as followed. In Sec-
tion II we present results based on the mean-field approach in
terms of the dimer (triplon) operators. In Section III we extend
our treatment beyond mean-field, and even further in Section
IV, where we also take into account low-energy singlet two-
triplon excitations. Section V contains our conclusions.
II. MEAN-FIELD TREATMENT
We start by rewriting Eq. (1) in the the bond-operator
representation13, where on a dimer i, the two spins forming it
are expressed as: Sα1,2 =
tiα±t†iαsi−iǫαβγt
tiγ), and
, α = x, y, z create a singlet and triplet of states. We re-
fer to the triplet (S=1) quasiparticle, t†
, as “triplon”. The
bold indexes i, j,m, l label the dimers (see Fig. 1). Summa-
tion over repeated Greek indexes is assumed, unless indicated
otherwise.
The hard-core constraint, s†
si + t
tiα = 1, must be en-
forced on every site, which at the mean-field (MF) level can be
done by introducing a term in the Hamiltonian, −µ
(s2 +
tiα − 1). Then µ and the (condensed) singlet amplitude
s ≡ 〈si〉, are determined by the MF equations13. We obtain at
the quadratic level, in momentum representation:
tkα +
−kα + h.c.
where
Ak = J/4− µ+ s2(ξ−k +K/2) + s
4Σ(k) ,
Bk = s
+ s4Σ(k) ,
= −(J/2) coskx + (J ±K/4) cosky . (3)
The four-spin interaction from (1) acting between two
dimers (e.g. i, j in Fig. 1) contributes to the “on-site” gap and
hopping (ξ−
) via Ak, as well as to the quantum fluctuations
term Bk. The part involving four dimers has been split in a
mean-field fashion, leading to the Hartree-Fock self-energy
−Σ(k)/K = 2Σx cos kx + 2Σy cos ky +Σxy cos kx cos ky,
tmα + t
t†mα〉 , (5)
where i,m are neighboring dimers in the x (horizontal) direc-
tion (Fig. 1), and similarly for the y and the diagonal contribu-
tions. The triplon dispersion is ω(k) =
, and has
a minimum at the Néel ordering wave-vector kAF = (0, π)
(since we work on a dimerized lattice). The ground state en-
ergy is then easily computed,
EGS = E0 + 〈H2〉 , (6)
where
E0/N = −
(Js2 +Ks4) + µ(−s2 + 1) + (7)
Σ2x +Σ
, (8)
〈H2〉 =
(ω(k)−Ak) . (9)
The mean-field equations require a numerical minimization
of EGS with respect to the parameters {µ, s,Σx,Σy,Σxy}.
This amounts to the self-consistent Hartree-Fock approxima-
tion for Σ(k). The result for the triplon gap ∆ = ω(kAF ) is
presented in Fig. 2 (black curve).
The MF result (K/J)c ≈ 0.6 substantially underestimates
the location of the critical point, compared to the the QMC
calculations, where (K/J)c ≈ 1.856,7. Interestingly, if one
solves the MF equations ignoring both the hard core and the
Σ(k), one finds (K/J)c = 1. Physically, in the full MF, the
hard core contribution increases the gap (and hence the stabil-
ity of the dimer phase) while at the same time suppressing the
antiferromagnetic fluctuations (which favor the Néel state).
We also note that a recent (hierarchical) MF treatment
based on the plaquette ground state also underestimates very
strongly the QCP location ((K/J)c ≈ 19), similarly to our
result. In our view this means that both mean field approaches
are not sufficient to attack the present problem, where fluctu-
ations are apparently very strong. We choose to accept that
the numerical QMC result gives the most accurate determina-
tion of the QCP, and therefore in what follows we extend our
treatment in several directions beyond mean-field theory.
0.5 1 1.5 2 2.5
Brueckner field theory (II)
Brueckner +Singlet fluctuations (III)
Mean-Field Theory (I)
(III)
FIG. 2: (Color online) Triplon excitation gap ∆ = ω(kAF ) in vari-
ous approximations. The point ∆ → 0 corresponds to transition to
the Néel phase.
III. BEYOND MEAN-FIELD: THE DILUTE TRIPLON GAS
APPROXIMATION
A more accurate treatment of fluctuations is possible by
taking into account the hard-core constraint beyond mean-
field. One can set the singlet amplitude s = 1 in the pre-
vious formulas, but introduce an infinite on-site repulsion be-
tween the triplons, U
tβitαi, U → ∞. As long as
the triplon density (determined by the quantum fluctuations) is
low, an infinite repulsion corresponds to a finite scattering am-
plitude between excitations and can be calculated by resum-
ming ladder diagrams for the scattering vertex14. This leads
to the effective triplon-triplon vertex Γ(k, ω) which was pre-
viously calculated15:
Γ−1(k, ω)=
ω(q) + ω(k− q)− ω
u → v
ω → −ω
This vertex in turn affects the triplon dispersion via (what we
call) the Brueckner self-energy15:
ΣB(k, ω) = 4
v2qΓ(k+ q, ω − ω(q)). (11)
The corresponding parameters in the quadratic Hamiltonian
(2) in this case are
Ak = J + 2K(1− 4nt/3) + ξ−k +Σ(k) + ΣB(k, 0),
Bk = ξ
+Σ(k). (12)
The Bogolubov coefficients are defined in the usual way
= 1/2 + Ak/(2ω(k)) = 1 + v
. The various terms
in Σ(k) can be expressed through them: for example Σx =
+ vkuk) cos kx, and so on. The density of triplons
is nt = 〈t†iαtiα〉 = 3
v2k. In addition, the renormaliza-
tion of the quasiparticle residue, Z−1
= 1 − ∂ΣB(k, 0)/∂ω,
−k y y
FIG. 3: Renormalization of quantum fluctuations by resummation of
a ladder series, with (13) at the vertices.
implies the replacement uk →
Zkuk, vk →
Zkvk in
all the formulas15, and the renormalized spectrum ω(k) =
An iterative numerical evaluation of the spectrum using
the above equations, which amounts to solution of the Dyson
equation, leads to the result shown in Fig. 2 (blue curve). The
above approach appears to be well justified since the quasi-
particle density nt < 0.1. The resulting critical point is still in
the “weak-coupling” regime K/J < 1, with about 100% de-
viation from the QMC result ((K/J)c ≈ 1.85). This suggests
that the on-site triplon fluctuations are not the dominant cause
for the disagreement with the QMC results; thus we proceed
to include two-particle fluctuations (in the triplon language),
which amounts to including dimer-dimer correlations.
IV. STRONG FLUCTUATIONS IN THE SINGLET
BACKGROUND: QCP BEYOND THE DILUTE TRIPLON
GAS APPROXIMATION
It is clear that “non-perturbative” effects are responsible for
driving the QCP towards the “strong-coupling” regionK/J ∼
2. To proceed we make two improvements to the previous
low-density, weak-coupling theory.
First, we take into account fluctuations in the singlet back-
ground, i.e. the manifold on which the triplons are built and
interact. The main effect originates from the action of the
four-spin K-term from (1) on two dimers, e.g. i, j in Fig. 1.
Part of this action has led to the on-site gap 2K in (12), fa-
voring dimerization. However, a strong attraction between the
two dimers is also present, since theK-term is symmetric with
respect to the index pair exchange (1, 2)(3, 4) ↔ (1, 4)(2, 3),
leading to a “plaquettization” tendency as well. In the triplon
language this is manifested by formation of bound states of
two triplons, due to their nearest-neighbor interactions
H4,y =
〈i,j〉y,αβ
tβitαj + γ2t
tβitβj
+ γ3t
tαitβj
, (13)
γ1 = −
, γ2 = −
, γ3 = −
We also checked that on the perturbative (Hartree-Fock) level,
the effect of this term on equations (3) and (12) was negligible
(and we did not write it explicitly).
An intuitive way of taking into account the effect of two-
triplon bound states (with total spin S=0) on the one-triplon
2.1 2.2 2.3 2.4 2.5
FIG. 4: (Color online) (a.) Singlet bound state energy Es (black),
binding energy ǫ = 2∆ − Es (blue), and the triplon gap ∆ (red).
(b.) Triplon density nt. (c.) Dimer order parameters. Dashed parts
of the lines represent points corresponding to rapid growth of the
quasiparticle density.
spectrum, is to work in the “local” approximation. This means
effectively neglecting the triplon dispersion and directly eval-
uating the ladder series that renormalizes the quantum fluctu-
ation term Bk in (2), corresponding to emission of a pair of
triplons with zero total momentum. This is illustrated graphi-
cally in Fig. 3, with the result
Bk = −
cos kx +
J +K/4
1− |γ|
 cos ky +Σ(k) ,
γ ≡ γ1 + 3γ2 + γ3 = −J −
K, (14)
where γ is the effective attraction of two triplons with total
S = 0, and
∆E = 2J +
K (15)
is the energy of two (non-interacting) triplons on adjacent
sites. This calculation is justified for K/J ≫ 1 and leads
to an increase of the quantum fluctuations, and from there to
almost doubling of the triplon density nt (see Fig. 4 below).
It contributes significantly to the shift of the QCP.
We can go beyond the “local” approximation by solving
the Bethe-Salpeter equation for the bound state, formed due
to the attraction (13), and taking into account the full triplon
dispersion. The equation for the singlet bound state energy
Es(Q), corresponding to total pair momentum Q is
1 = 2γ
u4q cos
Es(Q)− ω(Q/2 + q)− ω(Q/2− q)
. (16)
Here we have, for simplicity, written only the main contribu-
tion to pairing (Eq. (13)) in the limit K/J ≫ 1, and have
neglected the on-site repulsion (which leads to slightly di-
minished pairing), as well as small pairing due to the ex-
change J from dimers in the x-direction on Fig. 1. It is
easily seen that the lowest energy corresponds to Q = 0;
we define from now on Es ≡ Es(Q = 0). The bind-
ing energy is ǫ = 2∆ − Es, where ∆ is the one-particle
gap. The bound state wave-function corresponding to Es is
|Ψ〉 =
α,i,j,qy
iqy(i−j)t
|0〉. In the “local” limit
(nearest-neighbor pairing), Ψqy =
2 cos qy .
Second, we have made subtle changes to the resumma-
tion procedure concerning the quasiparticle renormalization
Z , based on both formal and physical grounds. On the one
hand it is clear that in the Brueckner approximation (Eq. (11)),
where the self energy is linear in the density (ΣB ∝ nt), the
dependence of the vertex Γ−1 on density is beyond the accu-
racy of the calculation, meaning one can put uq = 1, vq = 0
in (10), instead of determining them self-consistently. This
leads to a decreased influence of the hard-core ΣB (which fa-
vors the dimer state) on the Hartree-Fock self-energy Σ(k)
from (4) (which favors the Néel state). It is indeed the mutual
interplay between ΣB > 0 and Σ(k) < 0, that determines the
exact location of the QCP in the course of the Dyson’s equa-
tion iterative solution. While in the “weak-coupling” regime
K/J < 1, ΣB always dominates, in the “strong-coupling” re-
gion K/J > 2, Σ(k) starts playing a significant role, since
parametrically Σ ∝ Knt. It is physically consistent that in
the region where singlet fluctuations in the dimer background
are strong, the hard-core effect is less important, i.e. in ef-
fect the kinematic hard-core constraint is “relaxed”. We also
observe that in typical models with QCP driven by explicit
dimerization, such as the bilayer model, the described differ-
ence in approximation schemes makes a very small difference
on the location of the QCP16, since those models are always
in the “weak-coupling” regime, dominated by the hard-core
repulsion of excitations on a non-fluctuating dimer configura-
tion. The purpose of the above rather technical diversion is to
emphasize that care has been taken to take into account as ac-
curately as possible the effect of the (low-energy) two-particle
spectrum on the one-particle triplon gap.
Our results are summarized in Fig. 4 and Fig. 2 (red line) for
the gap. The critical point is shifted towards (K/J)c ≈ 2.16
(in much better agreement with QMC data), with a very
strong increase of the density towards Kc. This translates
into a decrease of the dimer order, as measured by the two
dimer order parameters that we compute from the expres-
sions: Dx = |〈S3 · S4〉 − 〈S5 · S4〉| = | − 34 + nt +
Dy = |〈S3 ·S4〉−〈S1 ·S4〉| = |− 34+nt−
Σy|. The spins are
labeled as in Fig. 1. The singlet bound state energy Es(0) also
tends towards zero at the QCP, with the corresponding binding
energy remaining quite large ǫ/J ≈ 1. All these effects point
towards a tendency of the system to restore the lattice sym-
metry, although it is certainly clear that as the critical point is
approached, our approximation scheme (low density of quasi-
particles) breaks down (dashed lines on figures). We should
point out that the sharpness of variation near Kc is not due to
divergence in any of the self-energies but is a result of rapid
cancellation at high orders (i.e. iterations in the Dyson equa-
tion). In fact cutting off our iterative procedure at finite order
gives a smooth curve, suggesting that additional classes of di-
agrams become important (although in practice their classifi-
cation is an insurmountable task). The merger of singlet and
triplet modes, which we find near the QCP, in principle reflects
a tendency towards quasiparticle fractionalization (spinon de-
confinement) and is also found in the 1D Heisenberg chain
with frustration17, where spinons are always deconfined.
Since we are now dealing with a situation where the density
is not very small nt ≈ 0.2, it is prudent to check how the next
order in the density may affect the above results. For example
at second order in the density, the self-energy Σ(k) changes
by amount δΣ(k), i.e. one has to add this contribution to the
right hand side of Eq.(12), namely:
Ak → Ak + δΣ(k), Bk → Bk + δΣ(k). (17)
We have found
δΣ(k) = −2K(Q2x − P 2x ) cos kx +
+2K(Q2y − P 2y ) cos ky −
−K(Q2xy − P 2xy) cos kx cos ky, (18)
and the following definitions are used:
Px = (1/3)
tmα〉 =
v2k cos kx,
Qx = (1/3)
〈tiαtmα〉 =
ukvk cos kx, (19)
and similarly for the other directions, for example Py =
v2k cos ky , Pxy =
v2k cos kx cos ky , etc. After includ-
ing these expressions in our numerical iterative procedure, we
have found that the QCP is shifted by a very small amount,
and the overall picture, as summarized in Fig. 4 and Fig. 2
(red line), still stands.
V. CONCLUSIONS
In conclusion, we have shown that the QCP between the
Néel and the dimer state in the model (1) is of unconven-
tional nature, in the sense that it is characterized by the pres-
ence of both triplet and singlet low-energy modes. Near the
QCP, whose location ((K/J)c ≈ 2.16) we find in fairly good
agreement with recent QMC studies, the system exhibits: (1.)
Strong rise of the triplon excitation density, due to increased
quantum fluctuations, (2.) Corresponding strong decrease
(and ultimately vanishing) of the dimer order at the QCP (3.)
Vanishing of a singlet energy scale, related to the destruction
of the dimer “columns” in Fig. 1. The above effects are all re-
lated and influence strongly one another, ultimately meaning
that the QCP reflects strong fluctuations and can not be de-
scribed in a mean-field theory framework. These results also
suggest a desire of the system to restore the lattice symmetry
at the QCP, as found in the QMC studies6.
At the same time all our improvements beyond mean-field
theory have also resulted in a very sharp transition, which ap-
pears to be first order. However in our view our approach
is not capable of addressing correctly the issue of the order
of the phase transition, basically because once we take the
strong (inter) dimer fluctuations in to account, the triplon den-
sity starts rising quickly beyond control. This is in a certain
sense natural in a situation where the system wants to restore
the lattice symmetry at the QCP and thus the ground state
acquires strong admixture of plaquette, etc. fluctuations as
the dimers begin to “disappear.” This is also manifested in the
fact that our procedure is sensitive to the number of iterations
in the Dyson equation; all presented results are for an “infi-
nite” number of iterations, so that a fixed point is reached,
but cutting off the procedure results in a smoother behavior
and a shift of the QCP, which becomes iteration dependent.
We have not previously encountered such volatile behavior in
any other spin model with a dimer to magnetic order transi-
tion. Since iterations translate into accounting of more and
more fluctuations, the sensitivity of the results seems to mean
that the situation starts spiraling out of control near the QCP,
quite likely because classes of fluctuations become important
that are not included in the dimer description, such as longer
range correlations, etc. All this suggests that the triplon quasi-
particle description breaks down near the QCP which indeed
appears natural in a model where spinon deconfinement is ex-
pected to take place at the QCP6. On the other hand, if we
put aside the arguments that our approach is not reliable near
the QCP, the natural conclusion would be that the transition is
first order.
Acknowledgments
We are grateful to A. W. Sandvik, K. S. D. Beach, S.
Sachdev, and O. P. Sushkov for numerous stimulating dis-
cussions. A.H.C.N. was supported through NSF grant DMR-
0343790; V.N.K., D.X.Y., and D.K.C. were supported by
Boston University.
1 S. Sachdev, Quantum Phase Transitions (Cambridge University
Press, Cambridge, 1999).
2 T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P.
A. Fisher, Science 303, 1490 (2004); T. Senthil, L. Balents, S.
Sachdev, A. Vishwanath, and M. P.A. Fisher, Phys. Rev. B 70,
144407 (2004).
3 S. Sachdev and N. Read, Int. J. Mod. Phys. B 5, 219 (1991).
4 R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys.
Rev. B 60, 7278 (1999); L. Capriotti, F. Becca, A. Parola, and
S. Sorella, Phys. Rev. Lett. 87, 097201 (2001); M. Mambrini,
A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422
(2006).
5 P. Henelius and A. W. Sandvik, Phys. Rev. B 62, 1102 (2000).
6 A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
7 R. G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 017203 (2008).
8 J. Lou, A. W. Sandvik, and N. Kawashima, arXiv:0908.0740.
9 L. Isaev, G. Ortiz, and J. Dukelsky, arXiv:0903.1630.
10 A. B. Kuklov, M. Matsumoto, N. V. Prokof’ev, B. V. Svistunov,
and M. Troyer, Phys. Rev. Lett. 101, 050405 (2008).
11 F.-J. Jiang, M. Nyfeler, S. Chandrasekharan, and U.-J. Wiese, J.
Stat. Mech., P02009 (2008).
12 The possibility of four-spin exchange induced dimerization has
been discussed in the context of the full ring exchange, of which
the interaction (1) is part; see e.g. A. Läuchli, J. C. Domenge,
C. Lhuillier, P. Sindzingre, and M. Troyer, Phys. Rev. Lett. 95,
137206 (2005).
13 S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
14 A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle
Systems (Dover Publications, Mineola, NY, 2003).
15 V. N. Kotov, O. P. Sushkov, Z. Weihong, and J. Oitmaa, Phys. Rev.
Lett. 80, 5790 (1998).
16 P. V. Shevchenko, A. W. Sandvik, and O. P. Sushkov, Phys. Rev.
B 61, 3475 (2000).
17 W. H. Zheng, C. J. Hamer, R. R. P. Singh, S. Trebst, and H.
Monien, Phys. Rev. B 63, 144411 (2001); ibid. 63, 144410 (2001).
http://arxiv.org/abs/0908.0740
http://arxiv.org/abs/0903.1630

  We study the S=1/2 Heisenberg antiferromagnet on a square lattice with
nearest-neighbor and plaquette four-spin exchanges (introduced by A.W. Sandvik,
Phys. Rev. Lett. {\bf 98}, 227202 (2007).)
  This model undergoes a quantum phase transition from a spontaneously
dimerized phase to N\'eel order at a critical coupling. We show that as the
critical point is approached from the dimerized side, the system exhibits
strong fluctuations in the dimer background, reflected in the presence of a
low-energy singlet mode, with a simultaneous rise in the triplet quasiparticle
density. We find that both singlet and triplet modes of high density condense
at the transition, signaling restoration of lattice symmetry. In our approach,
which goes beyond mean-field theory in terms of the triplet excitations, the
transition appears sharp; however since our method breaks down near the
critical point, we argue that we cannot make a definite conclusion regarding
the order of the transition.

Quantum Phase Transition in the Four-Spin Exchange Antiferromagnet
Valeri N. Kotov, Dao-Xin Yao, A. H. Castro Neto, and D. K. Campbell
Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, MA 02215
We study the S=1/2 Heisenberg antiferromagnet on a square lattice with nearest-neighbor and plaquette four-
spin exchanges (introduced by A.W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).) This model undergoes a
quantum phase transition from a spontaneously dimerized phase to Néel order at a critical coupling. We show
that as the critical point is approached from the dimerized side, the system exhibits strong fluctuations in the
dimer background, reflected in the presence of a low-energy singlet mode, with a simultaneous rise in the triplet
quasiparticle density. We find that both singlet and triplet modes of high density condense at the transition,
signaling restoration of lattice symmetry. In our approach, which goes beyond mean-field theory in terms of the
triplet excitations, the transition appears sharp; however since our method breaks down near the critical point,
we argue that we cannot make a definite conclusion regarding the order of the transition.
PACS numbers: 75.10.Jm, 75.30.Kz, 75.50.Ee
I. INTRODUCTION
Problems related to quantum criticality in quantum spin
systems are of both fundamental and practical importance1.
Numerous materials, such as Mott insulators, exhibit either
antiferromagnetic (Néel) order or quantum disordered (spin
gapped) ground state depending on the distribution of Heisen-
berg exchange couplings and geometry. External perturba-
tions (such as doping or frustration) can also cause quantum
transitions between these phases. Systems with spin 1/2 are
indeed the most interesting as they are the most susceptible to
such transitions. It is well understood that the quantum transi-
tion between a quantum disordered and a Néel phase is in the
O(3) universality class1, where a triplet state condenses at the
quantum critical point (QCP).
A recent exciting development in our theoretical under-
standing of QCPs originated from the proposal that if the
quantum disordered (QD) phase spontaneously breaks lattice
symmetries (e.g. is characterized by spontaneous dimer or-
der), and the transition is of second order, then exactly at
the QCP spinon deconfinement occurs, i.e. the excitations
are fractionalized2. It is assumed that the Hamiltonian itself
does not break the lattice symmetries (i.e. does not have “triv-
ial” dimer order caused by some exchanges being stronger
than the others). We use the terms “dimer order” and “va-
lence bond solid (VBS) order” interchangeably. It is expected
that the dimer order vanishes exactly at the point where Néel
order appears, i.e. there is no coexistence between the two
phases. Deconfinement thus is intimately related to disap-
pearance of VBS order; indeed if the latter persisted in the
Néel phase it would be impossible to isolate a spinon, as
“pairing” would always take place. Spontaneous VBS or-
der driven by frustration has been a common theme in quan-
tum antiferromagnetism3, although its presence and the na-
ture of criticality in specific models, such as the 2D square-
lattice frustrated Heisenberg antiferromagnet, is still some-
what controversial4. It would be particularly useful to apply
unbiased numerical approaches, such as the Quantum Monte
Carlo (QMC) method, to study frustrated spin models; how-
ever due to the fatal “sign” problem5, frustrated Heisenberg
systems are beyond the QMC reach.
In a recent study, the QMC method was applied to a four-
spin exchange quantum spin model without frustration, which
was shown to exhibit columnar dimer VBS order and a mag-
netically ordered phase with a deconfined QCP separating
them6. These conclusions were later confirmed by further
QMC studies7. Extensions of the model, which include for
example additional (six-spin) interactions, provide additional
support for a continuous QCP8. A different VBS pattern (pla-
quette order) was also proposed for the four-spin exchange
model9. At the same time, the nature of the quantum phase
transition was challenged in Refs.[10,11], where arguments
were given that the transition is in fact of (weakly) first order.
It is the objective of the present work to study the Sand-
vik model6, by approaching the quantum transition from the
dimer VBS phase. Our approach uses as a starting point a
symmetry broken state (i.e. one out of four degenerate VBS
configurations), and we thus must search for signatures that
the system attempts to restore the lattice symmetry at the QCP.
Even though full restoration is impossible within the present
framework, we find a QCP characterized by condensation of
triplet modes of high density; this is in contrast to the conven-
tional situation when the condensing particles are in the dilute
Bose gas limit. The high density itself is due to the presence
of a singlet mode that condenses at the QCP, and reflects the
strong fluctuations of the background dimer order. The above
effects lead to the vanishing of the VBS order parameter; at
the same time our method, which accounts for the strong fluc-
tuations, leads to a rather sharp phase transition. It appears
that we cannot draw a definite conclusion about the order of
the transition because in the vicinity of the QCP the triplon
density increases uncontrollably, suggesting that other states
(such a plaquette states and larger clusters) are strongly ad-
mixed into the ground state. This is generally expected in a
situation where the lattice symmetry is restored at the quan-
tum critical point.
The model under consideration is
H = J
〈a,b〉
Sa.Sb −K
a,b,c,d
(Sa · Sb)(Sc · Sd), (1)
where J > 0,K > 0, and all spins are S = 1/2. Consider
the numbers 1,2,3,4 in Fig. 1. The summation in the four-
spin term is over indexes (a, b) = (1, 2), (c, d) = (3, 4) and
http://arxiv.org/abs/0704.0114v3
FIG. 1: (Color online) Dimer pattern in the quantum disordered
(VBS) phase, K/J > (K/J)c.
(a, b) = (1, 4), (c, d) = (2, 3) on a given plaquette, and then
summation is made over all plaquettes12. The range of pa-
rameters explored in Ref.[6] is K/J ≤ 2, and the QCP is at
(K/J)c ≈ 1.85. Our coupling notation is slightly different
from the one used in Refs.[6,7]; the coupling K is related to
the parameter Q6,7 via K = Q/(1 + Q/(2J)), and the criti-
cal point in that notation is (Q/J)c ≈ 25. The dimerization
pattern is proposed to be of the “columnar” type, as shown in
Fig. 1. Four such configurations exist. We will assume a con-
figuration of this type, will show that it is stable at K/J ≫ 1,
and will then search for an instability towards the Néel state
as K/J decreases.
The rest of the paper is organized as followed. In Sec-
tion II we present results based on the mean-field approach in
terms of the dimer (triplon) operators. In Section III we extend
our treatment beyond mean-field, and even further in Section
IV, where we also take into account low-energy singlet two-
triplon excitations. Section V contains our conclusions.
II. MEAN-FIELD TREATMENT
We start by rewriting Eq. (1) in the the bond-operator
representation13, where on a dimer i, the two spins forming it
are expressed as: Sα1,2 =
tiα±t†iαsi−iǫαβγt
tiγ), and
, α = x, y, z create a singlet and triplet of states. We re-
fer to the triplet (S=1) quasiparticle, t†
, as “triplon”. The
bold indexes i, j,m, l label the dimers (see Fig. 1). Summa-
tion over repeated Greek indexes is assumed, unless indicated
otherwise.
The hard-core constraint, s†
si + t
tiα = 1, must be en-
forced on every site, which at the mean-field (MF) level can be
done by introducing a term in the Hamiltonian, −µ
(s2 +
tiα − 1). Then µ and the (condensed) singlet amplitude
s ≡ 〈si〉, are determined by the MF equations13. We obtain at
the quadratic level, in momentum representation:
tkα +
−kα + h.c.
where
Ak = J/4− µ+ s2(ξ−k +K/2) + s
4Σ(k) ,
Bk = s
+ s4Σ(k) ,
= −(J/2) coskx + (J ±K/4) cosky . (3)
The four-spin interaction from (1) acting between two
dimers (e.g. i, j in Fig. 1) contributes to the “on-site” gap and
hopping (ξ−
) via Ak, as well as to the quantum fluctuations
term Bk. The part involving four dimers has been split in a
mean-field fashion, leading to the Hartree-Fock self-energy
−Σ(k)/K = 2Σx cos kx + 2Σy cos ky +Σxy cos kx cos ky,
tmα + t
t†mα〉 , (5)
where i,m are neighboring dimers in the x (horizontal) direc-
tion (Fig. 1), and similarly for the y and the diagonal contribu-
tions. The triplon dispersion is ω(k) =
, and has
a minimum at the Néel ordering wave-vector kAF = (0, π)
(since we work on a dimerized lattice). The ground state en-
ergy is then easily computed,
EGS = E0 + 〈H2〉 , (6)
where
E0/N = −
(Js2 +Ks4) + µ(−s2 + 1) + (7)
Σ2x +Σ
, (8)
〈H2〉 =
(ω(k)−Ak) . (9)
The mean-field equations require a numerical minimization
of EGS with respect to the parameters {µ, s,Σx,Σy,Σxy}.
This amounts to the self-consistent Hartree-Fock approxima-
tion for Σ(k). The result for the triplon gap ∆ = ω(kAF ) is
presented in Fig. 2 (black curve).
The MF result (K/J)c ≈ 0.6 substantially underestimates
the location of the critical point, compared to the the QMC
calculations, where (K/J)c ≈ 1.856,7. Interestingly, if one
solves the MF equations ignoring both the hard core and the
Σ(k), one finds (K/J)c = 1. Physically, in the full MF, the
hard core contribution increases the gap (and hence the stabil-
ity of the dimer phase) while at the same time suppressing the
antiferromagnetic fluctuations (which favor the Néel state).
We also note that a recent (hierarchical) MF treatment
based on the plaquette ground state also underestimates very
strongly the QCP location ((K/J)c ≈ 19), similarly to our
result. In our view this means that both mean field approaches
are not sufficient to attack the present problem, where fluctu-
ations are apparently very strong. We choose to accept that
the numerical QMC result gives the most accurate determina-
tion of the QCP, and therefore in what follows we extend our
treatment in several directions beyond mean-field theory.
0.5 1 1.5 2 2.5
Brueckner field theory (II)
Brueckner +Singlet fluctuations (III)
Mean-Field Theory (I)
(III)
FIG. 2: (Color online) Triplon excitation gap ∆ = ω(kAF ) in vari-
ous approximations. The point ∆ → 0 corresponds to transition to
the Néel phase.
III. BEYOND MEAN-FIELD: THE DILUTE TRIPLON GAS
APPROXIMATION
A more accurate treatment of fluctuations is possible by
taking into account the hard-core constraint beyond mean-
field. One can set the singlet amplitude s = 1 in the pre-
vious formulas, but introduce an infinite on-site repulsion be-
tween the triplons, U
tβitαi, U → ∞. As long as
the triplon density (determined by the quantum fluctuations) is
low, an infinite repulsion corresponds to a finite scattering am-
plitude between excitations and can be calculated by resum-
ming ladder diagrams for the scattering vertex14. This leads
to the effective triplon-triplon vertex Γ(k, ω) which was pre-
viously calculated15:
Γ−1(k, ω)=
ω(q) + ω(k− q)− ω
u → v
ω → −ω
This vertex in turn affects the triplon dispersion via (what we
call) the Brueckner self-energy15:
ΣB(k, ω) = 4
v2qΓ(k+ q, ω − ω(q)). (11)
The corresponding parameters in the quadratic Hamiltonian
(2) in this case are
Ak = J + 2K(1− 4nt/3) + ξ−k +Σ(k) + ΣB(k, 0),
Bk = ξ
+Σ(k). (12)
The Bogolubov coefficients are defined in the usual way
= 1/2 + Ak/(2ω(k)) = 1 + v
. The various terms
in Σ(k) can be expressed through them: for example Σx =
+ vkuk) cos kx, and so on. The density of triplons
is nt = 〈t†iαtiα〉 = 3
v2k. In addition, the renormaliza-
tion of the quasiparticle residue, Z−1
= 1 − ∂ΣB(k, 0)/∂ω,
−k y y
FIG. 3: Renormalization of quantum fluctuations by resummation of
a ladder series, with (13) at the vertices.
implies the replacement uk →
Zkuk, vk →
Zkvk in
all the formulas15, and the renormalized spectrum ω(k) =
An iterative numerical evaluation of the spectrum using
the above equations, which amounts to solution of the Dyson
equation, leads to the result shown in Fig. 2 (blue curve). The
above approach appears to be well justified since the quasi-
particle density nt < 0.1. The resulting critical point is still in
the “weak-coupling” regime K/J < 1, with about 100% de-
viation from the QMC result ((K/J)c ≈ 1.85). This suggests
that the on-site triplon fluctuations are not the dominant cause
for the disagreement with the QMC results; thus we proceed
to include two-particle fluctuations (in the triplon language),
which amounts to including dimer-dimer correlations.
IV. STRONG FLUCTUATIONS IN THE SINGLET
BACKGROUND: QCP BEYOND THE DILUTE TRIPLON
GAS APPROXIMATION
It is clear that “non-perturbative” effects are responsible for
driving the QCP towards the “strong-coupling” regionK/J ∼
2. To proceed we make two improvements to the previous
low-density, weak-coupling theory.
First, we take into account fluctuations in the singlet back-
ground, i.e. the manifold on which the triplons are built and
interact. The main effect originates from the action of the
four-spin K-term from (1) on two dimers, e.g. i, j in Fig. 1.
Part of this action has led to the on-site gap 2K in (12), fa-
voring dimerization. However, a strong attraction between the
two dimers is also present, since theK-term is symmetric with
respect to the index pair exchange (1, 2)(3, 4) ↔ (1, 4)(2, 3),
leading to a “plaquettization” tendency as well. In the triplon
language this is manifested by formation of bound states of
two triplons, due to their nearest-neighbor interactions
H4,y =
〈i,j〉y,αβ
tβitαj + γ2t
tβitβj
+ γ3t
tαitβj
, (13)
γ1 = −
, γ2 = −
, γ3 = −
We also checked that on the perturbative (Hartree-Fock) level,
the effect of this term on equations (3) and (12) was negligible
(and we did not write it explicitly).
An intuitive way of taking into account the effect of two-
triplon bound states (with total spin S=0) on the one-triplon
2.1 2.2 2.3 2.4 2.5
FIG. 4: (Color online) (a.) Singlet bound state energy Es (black),
binding energy ǫ = 2∆ − Es (blue), and the triplon gap ∆ (red).
(b.) Triplon density nt. (c.) Dimer order parameters. Dashed parts
of the lines represent points corresponding to rapid growth of the
quasiparticle density.
spectrum, is to work in the “local” approximation. This means
effectively neglecting the triplon dispersion and directly eval-
uating the ladder series that renormalizes the quantum fluctu-
ation term Bk in (2), corresponding to emission of a pair of
triplons with zero total momentum. This is illustrated graphi-
cally in Fig. 3, with the result
Bk = −
cos kx +
J +K/4
1− |γ|
 cos ky +Σ(k) ,
γ ≡ γ1 + 3γ2 + γ3 = −J −
K, (14)
where γ is the effective attraction of two triplons with total
S = 0, and
∆E = 2J +
K (15)
is the energy of two (non-interacting) triplons on adjacent
sites. This calculation is justified for K/J ≫ 1 and leads
to an increase of the quantum fluctuations, and from there to
almost doubling of the triplon density nt (see Fig. 4 below).
It contributes significantly to the shift of the QCP.
We can go beyond the “local” approximation by solving
the Bethe-Salpeter equation for the bound state, formed due
to the attraction (13), and taking into account the full triplon
dispersion. The equation for the singlet bound state energy
Es(Q), corresponding to total pair momentum Q is
1 = 2γ
u4q cos
Es(Q)− ω(Q/2 + q)− ω(Q/2− q)
. (16)
Here we have, for simplicity, written only the main contribu-
tion to pairing (Eq. (13)) in the limit K/J ≫ 1, and have
neglected the on-site repulsion (which leads to slightly di-
minished pairing), as well as small pairing due to the ex-
change J from dimers in the x-direction on Fig. 1. It is
easily seen that the lowest energy corresponds to Q = 0;
we define from now on Es ≡ Es(Q = 0). The bind-
ing energy is ǫ = 2∆ − Es, where ∆ is the one-particle
gap. The bound state wave-function corresponding to Es is
|Ψ〉 =
α,i,j,qy
iqy(i−j)t
|0〉. In the “local” limit
(nearest-neighbor pairing), Ψqy =
2 cos qy .
Second, we have made subtle changes to the resumma-
tion procedure concerning the quasiparticle renormalization
Z , based on both formal and physical grounds. On the one
hand it is clear that in the Brueckner approximation (Eq. (11)),
where the self energy is linear in the density (ΣB ∝ nt), the
dependence of the vertex Γ−1 on density is beyond the accu-
racy of the calculation, meaning one can put uq = 1, vq = 0
in (10), instead of determining them self-consistently. This
leads to a decreased influence of the hard-core ΣB (which fa-
vors the dimer state) on the Hartree-Fock self-energy Σ(k)
from (4) (which favors the Néel state). It is indeed the mutual
interplay between ΣB > 0 and Σ(k) < 0, that determines the
exact location of the QCP in the course of the Dyson’s equa-
tion iterative solution. While in the “weak-coupling” regime
K/J < 1, ΣB always dominates, in the “strong-coupling” re-
gion K/J > 2, Σ(k) starts playing a significant role, since
parametrically Σ ∝ Knt. It is physically consistent that in
the region where singlet fluctuations in the dimer background
are strong, the hard-core effect is less important, i.e. in ef-
fect the kinematic hard-core constraint is “relaxed”. We also
observe that in typical models with QCP driven by explicit
dimerization, such as the bilayer model, the described differ-
ence in approximation schemes makes a very small difference
on the location of the QCP16, since those models are always
in the “weak-coupling” regime, dominated by the hard-core
repulsion of excitations on a non-fluctuating dimer configura-
tion. The purpose of the above rather technical diversion is to
emphasize that care has been taken to take into account as ac-
curately as possible the effect of the (low-energy) two-particle
spectrum on the one-particle triplon gap.
Our results are summarized in Fig. 4 and Fig. 2 (red line) for
the gap. The critical point is shifted towards (K/J)c ≈ 2.16
(in much better agreement with QMC data), with a very
strong increase of the density towards Kc. This translates
into a decrease of the dimer order, as measured by the two
dimer order parameters that we compute from the expres-
sions: Dx = |〈S3 · S4〉 − 〈S5 · S4〉| = | − 34 + nt +
Dy = |〈S3 ·S4〉−〈S1 ·S4〉| = |− 34+nt−
Σy|. The spins are
labeled as in Fig. 1. The singlet bound state energy Es(0) also
tends towards zero at the QCP, with the corresponding binding
energy remaining quite large ǫ/J ≈ 1. All these effects point
towards a tendency of the system to restore the lattice sym-
metry, although it is certainly clear that as the critical point is
approached, our approximation scheme (low density of quasi-
particles) breaks down (dashed lines on figures). We should
point out that the sharpness of variation near Kc is not due to
divergence in any of the self-energies but is a result of rapid
cancellation at high orders (i.e. iterations in the Dyson equa-
tion). In fact cutting off our iterative procedure at finite order
gives a smooth curve, suggesting that additional classes of di-
agrams become important (although in practice their classifi-
cation is an insurmountable task). The merger of singlet and
triplet modes, which we find near the QCP, in principle reflects
a tendency towards quasiparticle fractionalization (spinon de-
confinement) and is also found in the 1D Heisenberg chain
with frustration17, where spinons are always deconfined.
Since we are now dealing with a situation where the density
is not very small nt ≈ 0.2, it is prudent to check how the next
order in the density may affect the above results. For example
at second order in the density, the self-energy Σ(k) changes
by amount δΣ(k), i.e. one has to add this contribution to the
right hand side of Eq.(12), namely:
Ak → Ak + δΣ(k), Bk → Bk + δΣ(k). (17)
We have found
δΣ(k) = −2K(Q2x − P 2x ) cos kx +
+2K(Q2y − P 2y ) cos ky −
−K(Q2xy − P 2xy) cos kx cos ky, (18)
and the following definitions are used:
Px = (1/3)
tmα〉 =
v2k cos kx,
Qx = (1/3)
〈tiαtmα〉 =
ukvk cos kx, (19)
and similarly for the other directions, for example Py =
v2k cos ky , Pxy =
v2k cos kx cos ky , etc. After includ-
ing these expressions in our numerical iterative procedure, we
have found that the QCP is shifted by a very small amount,
and the overall picture, as summarized in Fig. 4 and Fig. 2
(red line), still stands.
V. CONCLUSIONS
In conclusion, we have shown that the QCP between the
Néel and the dimer state in the model (1) is of unconven-
tional nature, in the sense that it is characterized by the pres-
ence of both triplet and singlet low-energy modes. Near the
QCP, whose location ((K/J)c ≈ 2.16) we find in fairly good
agreement with recent QMC studies, the system exhibits: (1.)
Strong rise of the triplon excitation density, due to increased
quantum fluctuations, (2.) Corresponding strong decrease
(and ultimately vanishing) of the dimer order at the QCP (3.)
Vanishing of a singlet energy scale, related to the destruction
of the dimer “columns” in Fig. 1. The above effects are all re-
lated and influence strongly one another, ultimately meaning
that the QCP reflects strong fluctuations and can not be de-
scribed in a mean-field theory framework. These results also
suggest a desire of the system to restore the lattice symmetry
at the QCP, as found in the QMC studies6.
At the same time all our improvements beyond mean-field
theory have also resulted in a very sharp transition, which ap-
pears to be first order. However in our view our approach
is not capable of addressing correctly the issue of the order
of the phase transition, basically because once we take the
strong (inter) dimer fluctuations in to account, the triplon den-
sity starts rising quickly beyond control. This is in a certain
sense natural in a situation where the system wants to restore
the lattice symmetry at the QCP and thus the ground state
acquires strong admixture of plaquette, etc. fluctuations as
the dimers begin to “disappear.” This is also manifested in the
fact that our procedure is sensitive to the number of iterations
in the Dyson equation; all presented results are for an “infi-
nite” number of iterations, so that a fixed point is reached,
but cutting off the procedure results in a smoother behavior
and a shift of the QCP, which becomes iteration dependent.
We have not previously encountered such volatile behavior in
any other spin model with a dimer to magnetic order transi-
tion. Since iterations translate into accounting of more and
more fluctuations, the sensitivity of the results seems to mean
that the situation starts spiraling out of control near the QCP,
quite likely because classes of fluctuations become important
that are not included in the dimer description, such as longer
range correlations, etc. All this suggests that the triplon quasi-
particle description breaks down near the QCP which indeed
appears natural in a model where spinon deconfinement is ex-
pected to take place at the QCP6. On the other hand, if we
put aside the arguments that our approach is not reliable near
the QCP, the natural conclusion would be that the transition is
first order.
Acknowledgments
We are grateful to A. W. Sandvik, K. S. D. Beach, S.
Sachdev, and O. P. Sushkov for numerous stimulating dis-
cussions. A.H.C.N. was supported through NSF grant DMR-
0343790; V.N.K., D.X.Y., and D.K.C. were supported by
Boston University.
1 S. Sachdev, Quantum Phase Transitions (Cambridge University
Press, Cambridge, 1999).
2 T. Senthil, A. Vishwanath, L. Balents, S. Sachdev, and M. P.
A. Fisher, Science 303, 1490 (2004); T. Senthil, L. Balents, S.
Sachdev, A. Vishwanath, and M. P.A. Fisher, Phys. Rev. B 70,
144407 (2004).
3 S. Sachdev and N. Read, Int. J. Mod. Phys. B 5, 219 (1991).
4 R. R. P. Singh, Z. Weihong, C. J. Hamer, and J. Oitmaa, Phys.
Rev. B 60, 7278 (1999); L. Capriotti, F. Becca, A. Parola, and
S. Sorella, Phys. Rev. Lett. 87, 097201 (2001); M. Mambrini,
A. Läuchli, D. Poilblanc, and F. Mila, Phys. Rev. B 74, 144422
(2006).
5 P. Henelius and A. W. Sandvik, Phys. Rev. B 62, 1102 (2000).
6 A. W. Sandvik, Phys. Rev. Lett. 98, 227202 (2007).
7 R. G. Melko and R. K. Kaul, Phys. Rev. Lett. 100, 017203 (2008).
8 J. Lou, A. W. Sandvik, and N. Kawashima, arXiv:0908.0740.
9 L. Isaev, G. Ortiz, and J. Dukelsky, arXiv:0903.1630.
10 A. B. Kuklov, M. Matsumoto, N. V. Prokof’ev, B. V. Svistunov,
and M. Troyer, Phys. Rev. Lett. 101, 050405 (2008).
11 F.-J. Jiang, M. Nyfeler, S. Chandrasekharan, and U.-J. Wiese, J.
Stat. Mech., P02009 (2008).
12 The possibility of four-spin exchange induced dimerization has
been discussed in the context of the full ring exchange, of which
the interaction (1) is part; see e.g. A. Läuchli, J. C. Domenge,
C. Lhuillier, P. Sindzingre, and M. Troyer, Phys. Rev. Lett. 95,
137206 (2005).
13 S. Sachdev and R. N. Bhatt, Phys. Rev. B 41, 9323 (1990).
14 A. L. Fetter and J. D. Walecka, Quantum Theory of Many-Particle
Systems (Dover Publications, Mineola, NY, 2003).
15 V. N. Kotov, O. P. Sushkov, Z. Weihong, and J. Oitmaa, Phys. Rev.
Lett. 80, 5790 (1998).
16 P. V. Shevchenko, A. W. Sandvik, and O. P. Sushkov, Phys. Rev.
B 61, 3475 (2000).
17 W. H. Zheng, C. J. Hamer, R. R. P. Singh, S. Trebst, and H.
Monien, Phys. Rev. B 63, 144411 (2001); ibid. 63, 144410 (2001).
http://arxiv.org/abs/0908.0740
http://arxiv.org/abs/0903.1630
