arXiv:0704.1441v3  [cond-mat.str-el]  6 Aug 2007
Heisenberg antiferromagnet with anisotropic exchange on the kagomé lattice:
Description of the magnetic properties of volborthite
T. Yavors’kii,1 W. Apel,2 and H.-U. Everts3
Department of Physics and Astronomy, University of Waterloo,
200 University Avenue W, Waterloo, N2L 3G1, Canada.
Physikalisch-Technische Bundesanstalt, Bundesallee 100, D-38116 Braunschweig, Germany.
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, D-30167 Hannover, Germany.
(Dated: September 12, 2021)
We study the properties of the Heisenberg antiferromagnet with spatially anisotropic nearest-
neighbour exchange couplings on the kagomé net, i.e., with coupling J in one lattice direction and
couplings J ′ along the other two directions. For J/J ′ >∼ 1, this model is believed to describe the
magnetic properties of the mineral volborthite. In the classical limit, it exhibits two kinds of ground
states: a ferrimagnetic state for J/J ′ < 1/2 and a large manifold of canted spin states for J/J ′ > 1/2.
To include quantum effects self-consistently, we investigate the Sp(N ) symmetric generalisation of
the original SU(2) symmetric model in the large-N limit. In addition to the dependence on the
anisotropy, the Sp(N ) symmetric model depends on a parameter κ that measures the importance of
quantum effects. Our numerical calculations reveal that in the κ-J/J ′ plane, the system shows a rich
phase diagram containing a ferrimagnetic phase, an incommensurate phase, and a decoupled chain
phase, the latter two with short- and long-range order. We corroborate these results by showing that
the boundaries between the various phases and several other features of the Sp(N ) phase diagram
can be determined by analytical calculations. Finally, the application of a block-spin perturbation
expansion to the trimerised version of the original spin-1/2 model leads us to suggest that in the
limit of strong anisotropy, J/J ′ ≫ 1, the ground state of the original model is a collinearly ordered
antiferromagnet, which is separated from the incommensurate state by a quantum phase transition.
PACS numbers: 75.10.Jm,75.30.Kz,75.50.Ee
I. INTRODUCTION
In the ongoing search for novel states of condensed mat-
ter, frustrated antiferromagnets have played a key role
(for a recent review, see Ref. 1). Among the many sub-
stances that have been investigated experimentally and
the numerous spin models that have been studied theo-
retically, those in which the magnetic ions occupy the ver-
tices of corner-sharing frustrating entities have attracted
particular attention in this context. The best known ex-
amples are the kagomé antiferromagnet (KAF), consist-
ing of corner sharing triangles, and the pyrochlore anti-
ferromagnet, consisting of corner sharing tetrahedra (see
Fig. 1).
The main distinction between the KAF, the pyrochlore
antiferromagnet and other frustrated and unfrustrated
AAAA AAAA
AAAAAAAA
(a) (b)
FIG. 1: kagomé lattice (a), pyrochlore lattice (b)
magnets is the large ground-state degeneracy of the for-
mer: classical Heisenberg antiferromagnets with nearest-
neighbour interactions on corner-sharing lattices have a
large ground-state degeneracy, which in the above two ex-
amples even leads to a finite ground-state entropy (see,
e.g., Ref. 2 and references therein). Quantum effects
may lift this degeneracy, and, indeed, in numerical stud-
ies of small cells of the spin 1
KAF, an exponentially
large number of very low-lying quantum states has been
observed3,4. It has been suggested that this abundance
of low-lying states can be understood in a description of
the low-energy physics of the quantumKAF as spin liquid
consisting of nearest-neighbour spin singlets5,6. However,
a complete picture of the ground state and of the excita-
tions of the KAF is still missing. Further theoretical, but
also experimental studies with emphasis on the quantum
properties of the KAF are therefore highly desirable. In
this last respect, the mineral volborthite is a very promis-
ing candidate. It has been the subject of several recent
experimental investigations7–10. The magnetic lattice of
this natural antiferromagnet consists of the S = 1/2 spins
of Cu2+ ions that are located on the vertices of well sep-
arated planar kagomé-like nets. A monoclinic distortion
of the lattice leads to a slight difference between the ex-
change couplings along one lattice direction (J) and the
two other directions (J ′)(see Fig. 2). Since neither signs
of long-range order nor signs of a spin-gapped singlet
ground-state were found in experiments on volborthite,
the substance seems to be a good candidate for the ob-
servation of the low-energy features that are thought to
http://arxiv.org/abs/0704.1441v3
PSfrag replacements
[i, j]
〈i, k
FIG. 2: (Color online) Anisotropic kagomé model. The cou-
pling J ′ and the nearest neighbour distance will be set equal
to unity in the calculations. δ1(‖ êx), δ2, and δ3 are the three
primitive lattice vectors of the kagomé net.
be typical for kagomé type antiferromagnets1.
Whether and to what extent the different exchange cou-
plings along different lattice directions of the kagomé net
of volborthite influence the low-energy physics of the sys-
tem is presently unknown. In the present paper, we study
this question on the basis of the model Hamiltonian
HAKAF = J
[i,j]
SiSj + J
〈k,i〉
SkSi . (1)
The symbols [i, j] and 〈k, i〉 denote, respectively, bonds
between nearest-neighbour sites on the horizontal chains
(a, b) and bonds between the middle sites (c) and the
sites a, b, see Fig. 2. Since the physics of this model de-
pends only on the ratio J/J ′ of the exchange constants,
we set J ′ = 1 in the sequel. We will consider the spatially
anisotropic kagomé antiferromagnet (AKAF), Eq. (1), in
the full range of J , 0 < J < ∞ since this is of theoretical
interest: one expects to see quantum phase transitions
as J is increased. It is of particular interest to find out
whether there is a transition from two-dimensional mag-
netic states to a set of decoupled chains with free spins
on the axes between the chains for large values of J .
The paper is organised as follows. In Sec. II, we consider
the model (1) in the classical limit. At this level, we find
no sign of a transition from the two-dimensional magnet
to a set of decoupled chains as J increases to infinity.
Nonetheless, the ground-state degeneracy, as well as the
spin wave spectrum are found to change qualitatively
as the anisotropy of the model varies. In Sec. III, we
consider a generalisation of the SU(2) symmetric model
(1) to the Sp(N ) symmetric version11,12 and describe
its properties in the large-N limit, where a mean-field
treatment of the model is adequate. We obtain a de-
tailed description of how possible ground states of the
model depend on the coupling J and on the spin length
S. A fairly rich phase diagram with a ferrimagnetic
phase for small J , long-range ordered and short-ranged
incommensurate phases for intermediate values of J , and
a decoupled-chain phase for large J emerges. Parts of
these results have been published previously, see Ref. 13.
In Sec. IV, we devise trial quantum ground states of the
original S = 1/2 model. We chose the states such that
they are exact eigenstates of HAKAF, if the couplings on
the upward pointing triangles of Fig. 2 are switched off,
and we then treat these couplings perturbatively. In the
limit J −→ ∞ this yields an effective Hamiltonian for
the spins on the c sites which represents an anisotropic
triangular antiferromagnet. The conclusions of Starykh
and Balents14 about the ground state of this effective
model lead us to conjecture the existence of a quantum
phase transition in the AKAF for large J . In Sec. V, we
summarise and discuss our results. In two Appendices,
we present technical details of the counting procedure for
the classical ground-states, and of the Ginzburg-Landau
type procedure that allows us to determine the bound-
aries in the phase diagram analytically.
II. CLASSICAL AND SEMICLASSICAL
ASPECTS
Similar to other isotropic spin models on lattices with
triangular elementary cells, the classical ground states of
HAKAF, Eq. (1), are spin configurations, which satisfy the
condition that for each elementary triangular plaquette
of the lattice, Fig. 2, the energy is minimal.
For J = 0, this yields a ferrimagnetic state with the chain
spins aligned in one direction and the middle spins point-
ing in the opposite direction, so that the total magneti-
sation is M = N▽S (N▽: number of downward pointing
triangles, N▽ = Ns/3 where Ns is the number of sites of
the system). We illustrate this situation in Fig. 3. Ac-
cording to the Lieb-Mattis theorem, the exact quantum
ground state (GS) of the model HAKAF also has total
spin Stot = N▽ S for J = 0, see Ref. 15, i.e., for J = 0,
the quantum GS is ferrimagnetic too. By continuity, one
expects the quantum GS to remain ferrimagnetic for suf-
ficiently small finite J . This will be confirmed by our
considerations of the large-N limit of the Sp(N ) version
of our model (see the analytical and numerical work in
Sects. III, III C and Appendix B) and by the block spin
perturbation approach (Sec. IV). Classically, the ferri-
FIG. 3: (Color online) Ferrimagnetic state for J = 0,
i.e., when there is no coupling between chain spins, cf. Fig. 2.
ωα(q)
FIG. 4: (Color online) Spin-wave frequencies ωα(q), α =
1, 2, 3 for J = 0.4; the contour at the top of the plot marks
half the Brillouin zone.
magnetic state remains stable up to J = 1/2. The ex-
citation spectrum of the ferrimagnetic state obtained in
linear spin-wave (LSW) approximation is shown in Fig. 4.
The analytic expressions for these three frequency sur-
faces ωα(q), α = 1, 2, 3, are obtained as solutions of a
third-order secular equation and are too lengthy to be
presented here. However, one can easily assure oneself
that the dispersion of the gapless mode is quadratic at
the origin. Thus, one has the typical mode structure of
a ferrimagnet here with one ferromagnetic mode and two
optical modes, see, e.g., Ref. 16. As J increases towards
1/2, the ferromagnetic frequency surface looses its disper-
sion and turns into a plane of zero modes, one zero mode
for each wave vector in the magnetic Brillouin zone (BZ),
at J = 1/2. The gap of the lower optical mode closes at
this value of J in the centre of the BZ and the dispersion
of this mode becomes linear for small wave vectors as for
an antiferromagnetic spin-wave mode.
At J = 1/2, the classical GS configuration changes from
the unique ferrimagnetic state to an ensemble of degener-
ate canted coplanar states. These states are characterised
by two variables: the angle θ, which the middle spin of
a given triangular plaquette forms with the two chain
spins of the same plaquette (see Fig. 5), and the two
valued chirality χ = ±1, which denotes the direction in
which the spins turn as one moves around the plaquette
in the mathematically positive sense.
For J ≥ 1/2, the requirement that the energy of any of
the elementary triangular plaquettes of the lattice Fig. 2
be minimal is θ = arccos(−1/(2J)), (θ > 0). The differ-
ent degenerate canted states arise from different possibil-
FIG. 5: (Color online) Canted spins of the AKAF at J > 1/2.
ωα(q)
FIG. 6: (Color online) Same as Fig.4 for J = 0.6.
ities to assign positive or negative chiralities to the pla-
quettes of the lattice. We show in the Appendix A that
for the general case of θ 6= 2π/3 (J 6= 1), the number of
spin configurations, NanisoGS does not grow exponentially
with the number of sites. Rather, NanisoGS < 2
where α < 3. This implies that the ground-state en-
tropy per spin of the classical AKAF vanishes in the
thermodynamic limit. In this respect, the anisotropic
model differs qualitatively from the isotropic KAF in
the classical limit, which has an extensive entropy per
spin. In the limit J → 1, the anisotropic model ap-
proaches the isotropic KAF. Hence, one expects that for
the anisotropic model there is an extensive number of
low-lying excited states that become degenerate with the
GS in the isotropic limit.
As in the case of the isotropic KAF, the spin-wave Hamil-
tonian is in linear order independent of the particular
classical GS that has been chosen as the starting point
of the expansion, Ref. 17. This implies that lowest-order
quantum fluctuation do not select one or a group of clas-
sical GSs as true GSs, i.e., the possible ordering effects
of quantum fluctuations are not captured by the linear
spin-wave (LSW) approximation. Figs. 6, 7 show the
spin-wave frequency surfaces for J = 0.6 and for J = 3. It
is easy to show analytically that, as is illustrated in these
figures, the plane of zero frequency modes persists for all
values of J greater than 1/2. The surfaces for J < 1/2
and for J > 1/2 join smoothly at J = 1/2. Thus, in the
ωα(q)
FIG. 7: (Color online) Same as Fig.4 for J = 3.
LSW approximation, the transition from the ferrimag-
netically ordered state to the canted spin states appears
to be of second order. For J ≫ 1, the nonzero frequen-
cies gradually loose their dispersion perpendicular to the
strong-J direction and take the shape of the spin-wave
spectrum of antiferromagnetic chains parallel to this di-
rection. However, no sign of a further transition from
the canted spin states to a set of decoupled spin chains is
found in this semiclassical approach. In the next section,
we will consider the symplectic Sp(N ) generalisation of
the antiferromagnetic model HAKAF in the large-N limit.
This approach, which was first proposed by Read and
Sachdev, Refs. 11,12, as a method to study frustrated
antiferromagnets, has the benefit of including the or-
dering effects of quantum fluctuations self-consistently.
It is of particular interest for spin models with two or
more competing exchange couplings in the different lat-
tice directions or over different lattice distances such as
the present model, the J1-J2-J3 model
11, the Shastry-
Sutherland antiferromagnet18 and the anisotropic trian-
gular antiferromagnet19. For these models, it has pro-
vided an unbiased selection of possible GSs that may or
may not be ordered depending on the value of a param-
eter κ, which is connected with the spin length S (see
below).
III. MEAN FIELD Sp(N ) APPROACH
A. Brief review of the method
For a general antiferromagnetic Heisenberg model with a
positive interaction matrix Jij ,
Jij Si · Sj , (2)
the Sp(N ) generalisation reads
HSp(N ) = −
2N (J
jβ)(Jγδb
j) . (3)
Here,
. . .
is the 2N × 2N generalisation of the 2× 2 antisymmetric
tensor
, (5)
and bαi with α = 1, . . . , 2N are the Sp(N ) boson anni-
hilation operators. (Here and in the sequel, we closely
follow the notation of Ref. 12; in particular, summation
over repeated upper and lower indices is implied.) Thus,
J αβb†iαb
jβ is the generalisation of the creation operator
jβ for a singlet on the bond (i, j). For the special
case N = 1, one finds
(J αβb†iαb
jβ)(Jγδb
j) = −2Si · Sj + nbinbj/2 + δijnbi ,
where
nbi = b
i (7)
is the boson number operator at site i and where
Si = b
i /2 (8)
is the usual SU(2) spin operator at site i. (τ are the
Pauli matrices). Then, if one imposes the constraint
that the number of bosons is the same for all lattice sites,
nbi ≡ nb, the HamiltonianHSp(1) is the familiar SU(2) in-
variant antiferromagnetic Heisenberg Hamiltonian (plus
some constants) with nb = 2S.
In the subsequent exposition, we shall consider a Hamil-
tonian of the form (3) in the large-N limit following the
strategy of Refs. 11,12. Depending on the values of the
couplings Jij and of κ, the GS of HSp(N ) may either
break the global Sp(N ) symmetry and exhibit LRO or
it may be Sp(N ) symmetric with only SRO. Breaking of
the Sp(N ) symmetry will happen through condensation,
i.e., by macroscopic occupation of one of the Bose fields
bα. To allow for this, we introduce the parametrisation
bmσi =
b̃m̃σi
with α = (mσ), m = 1, · · · ,N , m̃ = 2, · · · ,N and
σ =↑, ↓. The field xσi is proportional to the condensate
amplitude, 〈bmσi 〉 =
N δm1 xσi . Aiming at a mean field
treatment of the HamiltonianHSp(N ), which becomes ex-
act in the large N limit, we decouple the quartic part by
the Hubbard-Stratonovich technique with complex fields
Qij = −Qji and with Lagrange multipliers λi that en-
force the local constraints (7). The variables Qij which
are defined on nearest neighbour bonds of the lattice are
expectation values of the bond singlet creation opera-
tors in the GS, Qij = 〈
σσ′ ε
j mσ′〉 and are to
be determined self-consistently from the mean field type
Hamiltonian
HMF =
Jij |Qij |2 −
Qijεσσ′
Nxσi xσ
b̃m̃σi b̃
+ h.c.
N|xσi |2 +
im̃σ b̃
i − nb
The variational ground state energy, EMF, of HMF in the
large-N limit is obtained by diagonalising the bosonic
part of HMF, integrating over the 2(N − 1)Ns bosonic
fields b̃m̃σi in the action associated with HMF. One ob-
tains:
Jij |Qij |2 −
Qijεσσ′x
j + h.c.
ωµ(k;Q, λ) +
|xσi |2 − 1− κ
. (11)
Here ωµ(k;Q, λ) are the positive eigenvalues of the
bosonic part of HMF, and κ = nb/N is kept fixed in the
limiting procedure11,12. The parameter κ is a measure
for the importance of quantum fluctuations: by vary-
ing κ from small to large values, one drives the system
from the regime dominated by quantum fluctuations to
the classical regime, i.e., from the disordered into the or-
dered region. Finally, the GS is obtained by finding the
saddlepoint of EMF in the space of the variables Qij and
xσi subject to the constraints
∂EMF(Q, λ)/∂λi = 0 . (12)
In addition to the GS itself, the spin-spin correlation
function Gij = 〈Si ·Sj〉 in the GS is an important piece of
information. In particular, by considering its behaviour
in the limit |i − j| −→ ∞, one can distinguish between
LRO and SRO. According to Sachdev12, to obtain Gij in
the Sp(N ) symmetric approach, the SU(2) invariant ex-
pression Si ·Sj must be replaced by the Sp(N ) invariant
expression
4N 2 (b
j − J αγJβδb
j). (13)
Within the mean field approach, Gij can then be calcu-
lated straightforwardly.
B. The anisotropic kagomé antiferromagnet
1. Choice of mean-field variables
We wish to apply the procedure described above to the
AKAF represented by the Hamiltonian (1). To render
the problem of finding the eigenvalues ωµ in Eq. (10)
and of optimising EMF tractable, we have to restrict the
number of variables Qij and λi. We do so by demand-
ing that the mean field Hamiltonian HMF for the spinon
operators b(†) is symmetric under transformations of the
projective symmetry group (PSG) that is related to the
symmetry group of the spin Hamiltonian HAKAF (Eq. 1)
(see Ref. 20). We include two translations, a rotation
by π and a mirror axis orthogonal to the preferred direc-
tion of the exchange constants (J). Thus generalising the
treatment of Wang and Vishwanath to our model, we find
eight mean-field states with different symmetries. Seven
of them have flux in the sense of Ref. 21 in various cells
of the lattice. Following the arguments in Ref. 21, we ex-
clude all flux-carrying states and end up with the solution
(cf. Fig. 8) P1,2,3 = Q1,2,3, Q3 = Q2, and λb = λa.
PSfrag replacements Q1
FIG. 8: (Color online) Arrangement of mean field parameters:
Q1 ≡ Qab, Q2 ≡ Qbc, and Q3 ≡ Qca denote the intra triangle
bonds, P1 ≡ Qba′ , P2 ≡ Qcb′ and P3 ≡ Qac′ denote the inter
triangle bonds. λa, λb, and λc are the Lagrange multipliers
needed to implement the constraints on the sites a, b, and c.
In order to check the flux-argument in Ref. 21, we have
explicitely studied the solution P1,2,3 = −Q1,2,3 and
found that it is always of higher energy (For J = 1, this
agrees with the result of Ref. 12).
Thus, the expression Eq. (11) can now be cast into the
= J |Q1|2 + 2 |Q2|2 − (2λa + λc)(κ+ 1)
ωµ(k)
1 + |xµ(k)|2
, (14)
where the condensate is written in diagonalized form and
ωµ(k) are the three positive solutions of
det D̂(ω) = 0 . (15)
Here,
D̂(ω) =
Λ̂− ωÎ Q̂
Q̂† Λ̂+ ωÎ
, (16)
Λ̂= diag(λa, λc, λa) , (17)
0 Q̃2(k) −JQ̃1(−k)
−Q̃2(−k) 0 Q̃3(k)
JQ̃1(k) −Q̃3(−k) 0
 , (18)
and Q̃a(k) =
eiδak/2 − e−iδak/2
, a = 1, 2, 3,
δ1,2,3, see Fig. 2.
2. Technical details of the numerical extremalisation
Determination of the ground state of the AKAF in the
considered approximation has been reduced to minimiza-
tion of the Eq. (14) with respect to two variables Q1 and
Q2, subject to the Lagrange constraints with respect to
two parameters λa and λc. Being apparently trivial, the
optimization procedure turns out to be quite involved
technically.
First, we find it crucial to consider at least two differ-
ent chemical potentials. Other than for the spatially
isotropic KAF, J = 1, we were not able to find a non-
trivial solution if we used a single λ, λa = λb = λc. If λa
and λc are different, [Λ̂, Q̂] 6= 0, the Lagrange multipliers
enter the expressions for the frequencies ωµ non-trivially,
other than in the case of a global uniform chemical po-
tential (cf. Ref. 12). In turn, the Lagrange constraints
cannot be satisfied semi-analytically, and require a nu-
merical treatment. Second, we choose to work directly
in the thermodynamic limit Ns → ∞ of the model (14)
by performing a numerical self-adapting integration over
the BZ. In this limit, the singularities can be integrated,
and symmetry breaking is signalled by the appearance of
a finite value of the condensate amplitude xµ(k) at a cer-
tain wavevector k = qord, which characterises the type of
magnetic order. We mention here that the extremalisa-
tion of a mean-field energy of the type of Eq. (14) can also
be achieved by solving the pertinent stationarity condi-
tions numerically for finite systems, i.e., for finite Ns,
see e.g., Ref. 22. Then, the type of magnetic order has
to be detected by calculating the structure factor. Third,
we see that the Eq. (14) has a minimum with respect to
the physical bond parameters Q1 and Q2 only after the
elimination of the chemical potentials. In the full Q− λ
space we face an extremalization problem.
Technically, we find it convenient to use a polar coordi-
nate parametrisation for the variablesQ1, Q2 and λa, λc:
Q1 = Q cos(α), Q2 = Q sin(α) , (20)
λa = Λ sin(β), λc = Λcos(β) . (21)
We perform an optimization with respect to the variables
Q,Λ, α, β, as well as condensate densities xµ(k) in accord
with the following algorithm (J and κ are kept fixed).
i. We fix the angles α, β and the amplitude Q, and
first exploit the stationarity condition for EMF with re-
spect to Λ. It is convenient to write the corresponding
equation in the following form:
[2 sin(β) + cos(β)] (κ+ 1)
|xµ(k)|2 ∂Λωµ(k)
∂Λ ωµ(k) , (22)
where Ω = 8π2/
3 is the volume of the unit cell. One
finds that Q and Λ enter the Eq. (22) only via the ratio
ξ = Λ/Q.
The requirement that the frequencies must be positive,
ωµ(k) ≥ 0, defines a lower limit ξmin(α, β) for ξ: the
frequencies ωµ(k) are positive for ξ > ξmin(α, β); for
ξ = ξmin(α, β), the lowest mode ωµ0 vanishes at some
point(s) k0 in the BZ. When this happens, the corre-
sponding condensate density xµ0(k0) can be put non-
zero, if this is necessary to satisfy Eq. (22). It is im-
portant to note that in order to determine the actual
value of ξmin(α, β) (as well as those of Q, α and β) it
suffices to only consider Eq. (22) at xµ(k) = 0, irrespec-
tive of whether there is condensate, ωµ0(k0) = 0, or not,
ωµ(k) 6= 0 for all k, µ.
We solve the Eq. (22) for ξ numerically in two steps.
First, we determine ξmin(α, β): we decrease ξ from large
positive values until the condition ωµ0(k0) = 0 signals
that ξ = ξmin(α, β). Second, we set xµ(k) ≡ 0 and at-
tempt to satisfy Eq. (22) in the interval ξ ≥ ξmin(α, β).
To this end, we set Λ = ξQ in Eq. (14) and vary ξ to de-
termine the extremum of EMF (i.e., Eq. (22)). We find
that the extremum is a maximum. If this maximum oc-
curs for some ξ > ξmin(α, β), then Eq. (22) is satisfied
with xµ(k) = 0. If, however, EMF(α, β, ξQ,Q) decreases
monotonously as we lower ξ down to ξ = ξmin(α, β) ,
then the Eq. (22) cannot be solved with xµ(k) = 0. In
this case, a finite condensate density xµ0 (k0) 6= 0, is re-
quired, in order to “compensate” for too large a value of
the lhs. of Eq. (22). This fixes both ξ = ξmin(α, β) and
the value xµ0 (k0) (cf. sects III B and IV B of Ref. 12).
ii. Having determined the value of ξ, we notice
that the function EMF(α, β,Λ, Q) is quadratic in Q and
bounded from below, which allows an analytical deter-
mination of Q as the position of the minimum.
iii. Finally, knowing the values of Λ and Q, we pro-
ceed by a numerical extremalization of EMF with respect
to the angles. The calculations show that EMF as a func-
tion of the angle β possesses a maximum, and a minimum
as a function of the angle α after β has been eliminated.
Thus, the variational energy EMF is bounded from below
in the variables Q1 and Q2, as expected.
iv. We iterate this procedure (i)-(iii) until conver-
gence is achieved.
C. Numerical results of the Sp(N) formalism
The results of the Sp(N ) approach in the large-N limit
are summarised in the zero temperature phase diagram
of the AKAF, Fig. 9. The central part of the phase dia-
gram is occupied by the incommensurate (IC) phase with
LRO at sufficiently small 1/κ. The phase boundary that
separates the region with SRO from the region with LRO
was found by checking whether for a given pair of J and
1/κ the lowest branch of the one spinon spectrum ωµ(k)
has zeros in the BZ or not, i.e., whether there will be
condensate at one or several points in the Brillouin zone
or not. As one might expect, LRO is maximally sup-
pressed by quantum fluctuations for J = 1, which is the
case of maximal frustration.
For J = 0, the exact quantum ground state of the
AKAF is ferrimagnetic (FM) according to the Lieb-
Mattis theorem15. In this state, the expectation value
Q1 which measures the singlet weight on the horizontal
bonds vanishes. As shown in Fig. 10, our Sp(N ) calcu-
lations recover this exact result and extend it to a finite
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
IC/SRO
IC/LRO
DC/SRO
FIG. 9: (Color online) Phase diagram of HAKAF as obtained
in the Sp(N ) approach. Symbols and lines, respectively, de-
note numerical and analytical results for the phase boundaries
(see text, Subsec. III B 2 and Appendix B). Quantum fluctu-
ations increase along the vertical axis. LRO: Long Range
Order; SRO: short range order; FM: ferrimagnet; IC: incom-
mensurate phase; DC: decoupled chains. At J = 1, the results
of Ref. 12 are recovered. Incommensurate order (see Fig. 12)
occurs between the boundaries of the ferrimagnetic phase (×)
and of the decoupled chain phase (∗).
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 10: (Color online) Mean field parameter Q1 as function
of the anisotropy.
interval 0 ≤ J ≤ Jferri(κ), which narrows as 1/κ in-
creases. The parameter Q2, which measures the singlet
weight on the diagonal bonds, is independent of J in this
interval; its value decreases as 1/κ increases (see Fig. 11).
Remarkably, the FM state retains its LRO in its entire
region of existence.
As J is increased beyond Jferri(κ), Q1 increases in the
manner of an order parameter at a second order phase
transition. At the same time, the parameter Q2 begins
to decrease, and eventuallly it drops to zero at some
J = JDC(κ). Thus, the large-N approach predicts the
existence of a decoupled-chain phase in the region above
the phase boundary JDC(κ). Q2 decreases to zero con-
tinuously so that the phase transition at JDC(κ) appears
to be of second order again.
Both LRO and SRO phases may be characterised by an
ordering wave vector qord = 2kmin, where kmin is that
wave vector at which the one-spinon excitation spectrum
ωµ(k) has its minimum. The static spin structure factor
S(q) develops a peak at qord. In Fig. 12, we display
the x-component of the ordering vector qxord = q
ord(J)
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 11: (Color online) Mean field parameter Q2 as function
of the anisotropy.
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 12: (Color online) Ordering wave vector qxord as function
of the anisotropy.
ord = 0). At the kagomé point J = 1, |qxord| = 4π/3 is
independent of the value of κ. For 1/κ <∼ 3, the behaviour
of qxord as a function of J is as expected: as J increases, it
increases monotonously until the phase boundary JDC(κ)
is reached and remains constant inside the DC phase.
However, for 1/κ >∼ 3 the function qordx (J) develops a
minimum at J ≈ 1.5, which becomes more pronounced
as 1/κ increases.
In Sec. III B 2 we emphasised that contrary to previ-
ous applications of the large-N approach to spin models
on kagomé and anisotropic triangular lattices11,12,19, we
found it essential to consider two chemical potentials λa
and λc here, one for the spins on the horizontal lattice
lines (λa) and one for the middle spins (λc). We display
the values of these parameters as functions of J in Fig. 13.
We have no physical explanation for the behaviour of λa,
λc as functions of J and κ but it is gratifying to see that
λa = λc at J = 1 independent of κ in accordance with
earlier work12.
As indicated above, along with numerical study of Eq. 11,
we performed extensive analytical calculations, both to
1/κ=1
1/κ=1.4
1/κ=5
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 13: (Color online) Lagrange multipliers λa, λc (chemical
potentials) as functions of the anisotropy.
corroborate the numerics and to obtain additional insight
into the problem. Details of the analytical techniques are
presented in Appendix B. Here we state that we were
able to analytically determine Sp(N ) phase boundaries
between the SRO and LRO DC phase, between the DC
and IC phase, and between the FM and IC phase, see
Fig. 9. Moreover, our analytical calculations allowed us
to explicitly confirm the existence of LRO inside the FM
phase and immediately to the right of the FM-IC phase
boundary. Likewise, the regions with SRO and LRO in-
side and immediately to the left of IC-DC phase bound-
ary were determined analytically. This was achieved by
evaluating in these regions the Sp(N ) generalisation of
the spin-spin correlation function 〈Si,u ·Sj,v〉 of the model
defined by expression (13). (u, v = a, b, c denote the sites
of the triangular cells i and j of the model, see Fig. 8).
On the right-hand side of the FM-IC boundary and in-
side the FM phase, we find for large distances between
the cells, |rj − ri| ≫ 1,
〈Si,u · Sj,v〉 ∼ Su Sv , (23a)
where
|x3(kmin)|2
N▽ (λc + λa)
w = a, b
−λa w = c
(23b)
and u, v = a, b, c denote the sites of the triangu-
lar cells i and j of the model, see Fig. 8. Here,
|x3(kmin)|2/N▽ is the condensate density at kmin =
(−π, 0), |x3(kmin)|2/N▽ = κ, see Eq. (B8). On the
FM-IC transition line and inside the FM phase, where
Eqs. (23) are valid, the parameters λa and λc are not
independent but can be expressed in terms of a single
parameter δ, see Eqs. (B4), (B12). The sign pattern on
the right-hand side of Eq. (23) and the ordering wave
vector qord = 2kmin = (−2π, 0) are indeed the proper-
ties one expects to find for the long-distance behaviour
of the spin-spin correlation function of a ferrimagneti-
cally ordered state. Since |x3(kmin)|2/N▽ remains finite
for arbitrarily small values of κ, the mean-field Sp(N )
approach predicts that this order persists in the ex-
treme quantum limit of our model, 1/κ ≫ 1. Together
with Eqs. (23), the fact that the condensate density
|x3(kmin)|2/N▽ remains constant inside the FM region,
see Eq. (B8), implies that the magnetisation of the FM
phase remains constant up to the FM-IC phase boundary.
The same behaviour of the magnetisation of a ferrimag-
netic phase has previously been observed in an exact-
diagonalisation study of a one-dimensional kagomé-like
antiferromagnet23. At the FM-IC phase boundary, the
magnetisation becomes spatially modulated with an in-
commensurate wave vector qord = 2kmin.
On the left-hand side of the IC-DC boundary and in-
side the DC phase we find the following large distance
behaviour of the spin-spin correlation function:
〈Si,c · Sj,c〉∼
cos [2kmin (ri − rj)]
1 + q21
|x3(kmin)|2 + |x3(−kmin)|2
3 (kmin)
,(24a)
〈Si,u · Sj,v〉∼ 0 for u, v 6= c, c . (24b)
Here, q1 and λa denote the saddle-point values of these
variables obtained from Eqs. (B31), (B32). q2 is a func-
tion of q1, determined by Eq. (B36) or by Eq. (B44)
depending on whether 1/κ < 1/κs or 1/κ > 1/κs
(κs = 0.181, see Fig. 22). ω
3 (kmin) is the value of the
second-order expansion coefficient of the lowest spinon
frequency ω3(k), cf. Eqs. (B20), (B21c), at its min-
imum, and 2kmin is the ordering wave vector imme-
diately to the left on the IC-DC phase boundary and
inside the DC phase; it is determined by Eq. (B35).
|x3(kmin)|2/N▽ = |x3(−kmin)|2/N▽ are the condensate
densities at the wave vectors ±kmin. As is shown in Ap-
pendix B, ω
3 (kmin) remains finite for 1/κ > 1/κs and
hence |x3(±kmin)|2/N▽ vanishes. Thus, 〈Si,c ·Sj,c〉 ∼ 0,
i.e., there is no LRO in this region. By contrast, for
1/κ < 1/κs both, |x3(±kmin)|2/N▽ and ω̄(2)3 (kmin) van-
ish when the IC-DC phase boundary is approached from
the left. However their ratio, which determines the spin-
spin-correlation function, Eqs. (24), remains finite in this
limit according to Eq. (B42). Thus, for 1/κ < 1/κs,
Eqs. (24) show that while the chain spins Sa, Sb remain
disordered, there is long-range IC order between the mid-
dle spins Sc along the IC-DC phase boundary and inside
the DC phase for sufficiently large κ. The middle spins
occupy the sites of a triangular lattice. Remarkably, the
correlations between these spins predicted by Eqs. (24)
are compatible with the spin pattern
Sj,c = S [cos(2kminrj)êx + sin(2kminrj)êy] (25)
that would obtain if the middle spins Sc were classical
spins coupled by a classical Heisenberg model with ex-
change constant J̃ along one lattice direction and cou-
plings J̃ ′ along the other two directions with a ratio J̃ ′/J̃
such that incommensurate order with wave vector 2kmin
would be established. This persistence of long range or-
der in the DC phase of the AKAF distinguishes our result
from the result obtained by Chung et al., Ref. 19, in their
large-N Sp(N ) treatment of the anisotropic triangular
antiferromagnet: there the DC phase consists of uncor-
related linear spin chains. Qualitative considerations of
the finite-N corrections to the mean-field Sp(N ) result
led the authors of Ref. 19 to the conclusion that instead
of the DC phase there is spin-Peierls order in the large–J
region of their model. In the next section, we will present
a different approach, a block-spin perturbation theory, to
get further insight into the properties of the AKAF for
the physical spin-1/2 case.
FIG. 14: (Color online) The kagomé lattice as a triangular
lattice of downward pointing triangles. The coupling strength
is J on the horizontal bond and unity on the other two bonds.
IV. BLOCK-SPIN PERTURBATION
APPROACH
The basic idea of the block-spin perturbation theory is
to calculate the states of small clusters of a given lattice
exactly and to treat the coupling between these clusters
perturbatively. The basic building blocks of the kagomé
lattice are triangles. Thus it is natural to consider the
trimerised kagomé lattice in which the spins on the down-
ward pointing triangles are assumed to be strongly cou-
pled whereas the coupling on the bonds of the upward
pointing triangles are assumed to be weak, see Fig. 14.
(Clearly, the exchange of the roles of the upward and
the downward pointing triangles will not affect the fur-
ther development to be presented in the current section.)
The Hamiltonian for this trimerised model reads
H(J, γ) = H▽(J) + γH△(J) , 0 ≤ γ ≤ 1 , (26)
where H▽(J) (H△(J)) denote those terms in the
Hamiltonian (Eq. 1) that act on the bonds of the down-
ward (upward) pointing triangles. We will determine
approximate GSs of this trimerised model in different
ranges of J in a perturbation expansion w.r.t. γ. The
hope is that the results will provide some qualitative
insight into the GS properties of the non-trimerised
model H(J, 1) which is our original model Eq. (1). The
same strategy has previously been applied sucessfully to
frustrated spin models by several authors5,24–26.
Obviously, the GSs of the unperturbed Hamiltonian
H(J, 0) are products of GSs of the individual downward
pointing triangular plaquettes. The GSs of a single
plaquette and the corresponding energies are
1. for J < 1:
|α〉 = 1√
|↑↑↓〉 − |↓↑↑〉
|↑↓↑〉 − |↓↑↑〉
, (27a)
PSfrag replacements
|α(ᾱ)〉:
∣β(β̄)
↑ (↓)↑ (↓)
↑ (↓)
FIG. 15: Ground-states of triangular plaquettes. Heavy lines
depict singlets. The coupling strength is J on the horizontal
bond and unity on the other two bonds.
|ᾱ〉 = 1√
|↓↓↑〉 − |↑↓↓〉
|↓↑↓〉 − |↑↓↓〉
, (27b)
εα = εᾱ = −1 + J/4 ; (27c)
2. J > 1:
|β〉 = 1√
|↑↓↑〉 − |↑↑↓〉
, (28a)
|↓↑↓〉 − |↓↓↑〉
, (28b)
εβ = εβ̄ = −3/4J . (28c)
Here, the ket vectors denote the spin state of the
plaquette in the Sz basis. The three arrows inside the
|cba〉 symbol denote from left to right the spin direction
at the sites c, b and a of the plaquettes in Fig. 15. The
states |α〉 (|ᾱ〉) and |β〉 (
) have total z-spin 1/2
(−1/2). They can be depicted graphically as shown in
Fig. 15. From these plaquette states, the zeroth order
GSs of the Hamiltonian H(J, γ) will be constructed. We
treat the cases J < 1 and J > 1 separately.
1. J < 1: Since the states |α〉, |ᾱ〉 are the GSs of the
individual downward pointing plaquettes in this case, the
states
|A(M)〉 =
iǫ{M}
jǫ{N▽−M}
|ᾱj〉 , (29)
are here the zeroth order GSs of H(J, γ). The set {M}
is a subset of M out of the N▽ downward pointing tri-
angles of the 3N▽-site kagomé lattice; the subscripts i, j
denote the position of individual triangular plaquettes in
the lattice of these plaquettes which is also triangular, see
Fig. 14. The zeroth order energy eigenvalues associated
with the states |A(M)〉 do not depend on M :
= N▽(−1 + J/4) . (30)
Hence, there are in total 2N▽ degenerate zeroth order
GSs |A(M)〉. The single plaquette states |α〉, |ᾱ〉 satisfy
the conditions for the validity of the Lieb-Mattis theorem,
Ref. 15: after a canonical transformation which rotates
the spins on the sites a and b by π around the z-axis
|↑〉 → i|↑〉, |↓〉 → −i|↓〉, and which leaves the spins on the
site c fixed the coefficients of all basis states on the right
sides of Eqs. (27a, 27b) become positive (+1/
6). As a
consequence, all the GSs |A(M)〉 satisfy the conditions
for the validity of the Lieb-Mattis theorem. For J = 0
it follows from this theorem that the total magnetisation
of the exact quantum GS |Φexact〉 of the Hamiltonian
HAKAF must be an eigenstate of the total magnetisation
m̂tot =
(Szi,a + S
i,b + S
i,c) (31)
with eigenvalue mtot = N▽/2, i.e., |Φexact〉 must be a
ferrimagnetic state. By continuity, one expects this to
be the case not only for J = 0, but up to a certain fi-
nite value of J . This suggests that the state |A(M=0)〉,
c.f. Eq (29), is the appropriate zeroth order GS in this
case and that the degeneracy of the states |A(M)〉 is lifted
by the perturbation H△ in favour of the state |A(0)〉.
To confirm this, we determine the creation energy of a
flipped plaquette in first order in γ, i.e., the difference
of the energy of the state with one plaquette spin flipped
relative to the ferrimagnetic state, and the energy of the
ferrimagnetic state:
δ(1)EA(M=1) = EA(1) − EA(0). (32)
A simple calculation yields
δ(1)EA(M=1) =
γ(1− J) , (33)
i.e., to first order, |A(M=0)〉, the ferrimagnetic GS is
stable w.r.t. a flip of a single plaquette spin, as long as
J < 1.
As a further check on the stability of the state
|A(M=0)〉, we calculate the dispersion of the excitation
energy of a propagating single flipped plaquette spin. For
this purpose, we need to determine the overlap matrix el-
ements between the state with a flipped plaquette spin
at the site j and states with a flipped spin at one of the
neighbouring sites,
tj,j±δ1 = 〈ᾱj |〈αj±δ1 |γJSjSj±δ1 |αj〉|ᾱj±δ1〉 =
(34a)
tj,j±δ2,3 = 〈ᾱj |
αj±δ2,3
∣γSjSj±δ2,3 |αj〉
∣ᾱj±δ2,3
(34b)
Here, δν , ν = 1, 2, 3, are the primitive lattice vectors
of the kagomé net, see Fig. 2; they connect the sites of
the plaquette lattice. Then, by diagonalising the ensuing
transfer Hamiltonian
Htrans = γ
J ( |j + δ1〉〈j|+ |j − δ1〉〈j| )
(|j + δ2〉〈j|+ |j − δ2〉〈j|)
(|j + δ3〉〈j|+ |j − δ3〉〈j|)
,(35)
where |j〉 denotes the state with a flipped plaquette spin
at site j, we obtain for the kinetic energy of this excita-
tion:
ε(k) =
J cos(kx)− cos(
) cos(
. (36)
Adding the energy for the creation of a single flipped pla-
quette spin, we find for the total energy of the excitation
in the limit of small wave vector k
ω(k) =
− J)k2x +
k2y +O(k4)
. (37)
Obviously, the ferrimagnetic state |A(M=0)〉 becomes
unstable against a propagating flipped plaquette spin
already at J = 1/4, i.e., much earlier than suggested
by the excitation energy of a static flipped spin (see
Eq. (33)). We remark that this bound is independent of
the actual magnitude of the perturbation parameter γ
and therefore, the qualitative result may survive in the
limit γ = 1.
2. J > 1: In this region, the states
|B(M)〉 =
iǫ{M}
jǫ{N▽−M}
with eigenenergy
= N▽(−3J/4) . (39)
are the zeroth order eigenstates of H(J, γ). These states
consist of free spins on the c-sites and of spin-singlet
dimers that cover every second bond of the horizontal
chains of the lattice. We wish to answer the question of
whether the 2N▽-fold degeneracy of these states, which
results from the degrees of freedom of the free spins, is
lifted by the perturbation γH△; in other words, we want
to find out whether the middle spins remain decoupled
from the chain spins. We proceed as in case (i). We
compare in a perturbation expansion w.r.t. γ the en-
ergy of the state |B(0)〉 with the energy of |B(1)〉, i.e.
with the state with one plaquette spin flipped relative to
|B(0)〉. We denote this difference by δ(1)EB(M = 1) =
EB(1) − EB(0). Surprisingly, we find that the matrix el-
ements 〈B(M)|H△|B(M)〉 vanish for any choice of M .
There is no first order correction to the energy E
δ(1)EB(M = 1) = 0. Moreover, we observe that the
off-diagonal matrix elements 〈B′(M)|H△|B(M)〉, where
|B′(M)〉 and |B(M)〉 contain identical numbers of states
but differ in their distribution over the N▽ down-
ward pointing triangles, also vanish. This implies that, in
contrast to case (i), a flipped plaquette spin cannot hop
to a neighbouring site in a first order process. Coupling
between the spins on the c-sites occurs only in second
order in γ. It is succinctly described by an effective spin
Hamiltonian for the c-site spins, which are at the same
time total spins of the downward pointing plaquettes (see
Fig. 15):
Heff =
Szi S
Sxi S
.(40)
Here, Sαi , α = x, y, z, denote plaquette spin operators;
i is the position of a downward pointing plaquette on
the triangular lattice formed by these plaquettes. The
exchange couplings J
and J⊥
are given as second order
matrix elements of H△:
〈Bi↑, i′↑|H△|X〉 〈X |H△|Bi↑, i′↑〉
2 εB − εX
〈Bi↑, i′↓|H△|Y 〉 〈Y |H△|Bi↑, i′↓〉
2 εB − εY
(41a)
〈Bi↓, i′↑|H△|X〉 〈X |H△|Bi↑, i′↓〉
2 εB − εX
, (41b)
and i′ ≡ i+ δν . Here, the states |Biσ, i′σ′〉 are zeroth
order GSs, Eq. (38), whose spin patterns are identical
on all sites except for the sites i and i′ where the z-
components of the spins take the values σ and σ′, re-
spectively; |X〉 and |Y 〉 are excited states of H▽. Of
course, since the SU(2) symmetry of the original Hamil-
tonian H(J, γ) must be conserved in the derivation of
Heff , the expressions Eqs. (41) must yield identical re-
sults, J
≡ Jδν . Non-zero contributions to J
and J⊥
are obtained if either the same term SiSi′ of
H△ acts in both matrix elements of the numerators of
Eq. (41) (two-block contributions) or the terms SiSk,
SkSi′ act in the left and right elements, respectively,
where the plaquette geometry must be as shown in Fig. 16
(three-block contributions). In contrast to the case of the
isotropic KAF studied by Zhitomirsky26, the three-block
PSfrag replacements
FIG. 16: Configurations of ▽ blocks contributing to the in-
terblock couplings J
and J⊥δν . Double dashed lines indicate
that the same term element of H△ acts twice between the ▽
blocks at sites i and i′ (see also text).
contributions do not produce three-spin interactions in
the present case. Rather, they contribute to the ex-
change interactions J
and J⊥
of the HamiltonianHeff ,
Eq. (40). The evaluation of the expressions (41) yields
Jδ1= γ
1 + 1
1 +O(J−1)
Jδ2 = Jδ3 =
1 + 1
1 + O(J−1)
. (43)
Obviously, these results are useful for J ≫ 1. There,
Heff represents a spin 1/2 Heisenberg Hamiltonian on
the triangular lattice of the c-sites with a coupling
along the δ1 direction that is strong (O(γ2/J)) in com-
parison to the couplings in the two other directions
(O(γ2/J2)). This limiting case of the anisotropic trian-
gular Heisenberg antiferromagnet (ATHAF) has recently
been analysed by Starykh and Balents with field theo-
retical methods14. These authors find that in the limit
of strong anisotropy, K ≡ Jδ1/Jδ2 → ∞, the GS of
the model Eq. (40) is a collinearly ordered antiferromag-
net (CAF) as depicted in Fig. 17. Since the ordering
wave vector qCAF = (π, π/2) of this phase does not
evolve continuously from the ordering wave vector qIC
of the incommensurate (IC) spiral phase of the ATHAF,
FIG. 17: (Color online) Collinear antiferromagnetic state
(CAF) on the triangular lattice14.
FIG. 18: (Color online) Tentative ground-state of the
anisotropic kagomé antiferromagnet in the limit J ≫ 1. Dou-
ble lines: dimers between the spins on the end points.
(qIC = (qx(K), 0) with −3π/2 ≤ qx(K) ≤ −π for
1/2 ≤ K ≤ ∞), they conclude that the IC phase and
the CAF phase must be separated by a quantum phase
transition. For the trimerised anisotropic kagomé model,
Eq. 26, these results have the following implications:
i) While in the limit of strong anisotropy J ≫ 1 there is
long-range collinear antiferromagnetic order among the
c-site spins, the a- and b-site spins are paired in singlets,
see Fig. 18.
ii) This picture of the GS of the trimerised anisotropic
kagomé model Eq. (26) differs from the result obtained in
the Sp(N ) approach insofar as for sufficiently large κ the
Sp(N ) approach predicts long range IC order among the
c-site spins up to arbitrarily large values of J . Thus, if
the picture of a CAF phase for large anisotropy persists in
the non-trimerised limit H(J, γ = 1) of the model Eq. 26,
one would expect a quantum phase transition between
the IC phase and the CAF phase of the AKAF similarly
as for the ATHAF. In closing this section, we remark
that the calculation that led to the effective Hamiltonian
Heff , Eq. 40, i.e., to the coupling between the c spins
in the strongly anisotropic limit, shows clearly that this
coupling arises from quantum fluctuations of the a and b
spins.
V. SUMMARY AND DISCUSSION
In this work, we have studied the ground state (GS) phase
diagram of the quantum Heisenberg antiferromagnet on
the kagomé lattice with spatially anisotropic exchange
(AKAF). The model is relevant for a description of mag-
netic properties of volborthite, which is a natural reali-
sation of a spin 1/2 antiferromagnet consisting of weakly
coupled slightly distorted kagomé layers. A small mono-
clinic distortion along one of the three lattice directions
causes the exchange coupling along this direction, J , to
differ from the couplings in the other two directions, J ′,
which we set equal to unity, cf. Fig. 2. We have inves-
tigated the problem in the full range of the anisotropy,
0 ≤ J ≤ ∞, using three different approximate methods:
the classical and semiclassical approach, the mean-field
Sp(N ) approach, and a block-spin perturbation theory.
The case J = 1 is the much studied isotropic kagomé
antiferromagnet (KAF). Exact diagonalisation studies of
this model3,4 are available. Their results speak conclu-
sively in favour of a spin liquid ground state1. This view
is supported by block-spin approaches5,6. Conflicting re-
sults have been found in Refs. 27–31, where various va-
lence bond crystal (VBC) states are proposed as ground
states of the KAF. However, a recent comparison of the
exact spectrum of the 36-site sample of the KAF against
the excitation spectra allowed by the symmetries of these
states, casts doubts on their validity32.
Within the whole anisotropy range, the case J = 0 is spe-
cial, since it allows for an exact characterisation of the
quantum GS as ferrimagnetic (FM) with a total mag-
netisation of M = S Ns/3 for a system of Ns spins of
magnitude S. In the classical picture, this state corre-
sponds to a unique staggered layout of spins with a non-
zero net magnetisation of the lattice unit cell (cf. Fig. 3).
In the classical limit, the ferrimagnetic ground state sur-
vives up to J = 1/2. For J > 1/2, the “chain” spins
(i.e., spins coupled by J) begin to tilt gradually towards
the middle (remaining) spins (see Fig. 5). This allows for
a formation of a large degenerate manifold of canted spin
states. In contrast to the isotropic case J = 1, where
the degeneracy grows exponentially with the system size
Ns, its growth is weaker: 2
Ns for J 6= 1. This im-
plies that there must be an increasingly large number of
classical low energy configurations as J approaches unity.
In the linear semiclassical approximation, the spin-wave
spectrum has one zero-frequency mode for each point
of the magnetic Brillouin zone (BZ). The spectrum is
identical for the different canted states for all J > 1/2.
Thus, in this order of the semiclassical approximation,
no order-by-disorder mechanism appears that would se-
lect one particular state or a particular group of states
from the manifold of canted states as true ground states.
In the limit J → ∞, the frequency spectrum of non-zero
modes gradually takes the shape of the spectrum that
one would expect for a set of uncoupled antiferromag-
netic spin chains parallel to the strong-J direction. No
qualitative change from the set of canted spin states to
the set of decoupled chains at a finite value of J is found.
We have further explored the nature of the phases at var-
ious J exploiting the mean field (MF) Sp(N ) approach,
that incorporates the effect of quantum fluctuations not
only perturbatively, but self-consistently. The strength
of quantum fluctuations is controlled by a parameter κ,
which is the analogue of the spin value S in the origi-
nal SU(2) symmetric model. In fact, for N = 1, when
the Sp(1) symmetric model is equivalent to the SU(2)
model, κ = 2S. For general N , this last identity does
not hold, but κ is still a measure for the importance of
quantum fluctuations that are strong for κ ≪ 1 and weak
for κ ≫ 1. In the MF Sp(N ) approach, the nature of the
phases that occur can be read from the values of the
mean field parameters Q1 and Q2 and from the spec-
trum of the bosonic spinon excitations. While the mean
field parameters Q1 and Q2 (cf. Fig. 8) are the GS ex-
pectation values of singlet bond operators, the structure
of the spinon spectrum, ωµ(k;Q, λ), determines the ex-
istence or non existence of long-range order (LRO): If
the spectrum becomes gapless at some wavevector qord,
a Bose condensate will form and a modulated structure
with the wavevector 2qord will acquire LRO.
As was to be expected, the phase diagram of the AKAF
obtained by the MF Sp(N ) approach contains an incom-
mensurate (IC) phase in the vicinity of the isotropic point
J = 1 which is ordered for sufficiently large κ according
to this approach, see Fig. 9. Qualitatively, we may gauge
the value of κ against the spin length S by looking at
the line J = 1 of the phase diagram which is the location
of the Sp(N ) analogue of the isotropic SU(2) symmet-
ric kagomé model: since, as we have argued above, the
SU(2) model is disordered for S = 1/2, we may con-
clude from Fig. 9 that the value of 1/κ that corresponds
to S = 1/2 must be greater than two. Somewhat sur-
prisingly, the FM phase remains long-range ordered for
arbitrarily small κ. This may reflect the fact that in
the SU(2) version of the model, the FM phase is ordered
even for the smallest physical spin value S = 1/2. A new
feature of the phase diagram is the prediction of a decou-
pled chain (DC) phase for large enough J , which has no
classical analogue. In this phase, the chains of strongly
coupled a- and b-site spins show no magnetic order. The
c-site spins which are interspersed between these chains
and which occupy the sites of triangular sublattice are
decoupled from the chain spins. However, they may or
may not exhibit long range order among themselves de-
pending on the magnitude of κ. Remarkably, the spin-
spin correlations, whose asymptotics were obtained ana-
lytically, are compatible with the spin-spin correlations
of an anisotropic classical Heisenberg antiferromagnet on
the triangular lattice whose exchange couplings differ in
one direction from those in the other two directions.
In order to tackle the problem of the GSs of the AKAF
from a third corner, we have used a block-spin pertur-
bation theory. This method has the advantage of being
applicable directly to the spin 1/2 version of the model.
In applying this approach, one has to initially group the
spins of the model in clusters. For the kagomé lattice, it
is natural to choose the spins around either the upward
or the downward pointing triangles as clusters of strongly
coupled units and to consider the coupling between these
clusters, γ, as the small expansion parameter. Thus one
trimerises the original model (see Fig. 14) and in so do-
ing, one breaks the translational invariance of the original
model. In the zeroth order of this expansion, two regions
can be distinguished by the eigenenergies of the individ-
ual trimers: J < 1 and J > 1. For sufficiently small J ,
one recovers the FM state as the GS in first order w.r.t. γ.
For J > 1, there are no first order corrections to the en-
ergy. Following an earlier application of the block-spin
technique to the isotropic KAF26, we determine for J > 1
in second order in γ an effective HamiltonianHeff for the
block-spins which can be identified as the middle spins of
the original model and that occupy the sites of a trian-
gular lattice. Heff is a Heisenberg Hamiltonian with a
coupling Jδ1 of the order of γ
2/J along the δ1 direction
(cf. Fig. 2) and couplings Jδ2 = Jδ3 of the order of γ
along the other two directions. The calculations that lead
to these results show clearly that the couplings between
the c-spins of the AKAF are due to fluctuations of the sin-
glets between the a- and b-spins into excited states. In a
very recent field theoretical study, Starykh and Balents14
arrive at the conclusion that for Jδ1 ≫ Jδ2,δ3 , the ground
state of the anisotropic triangular antiferromagnet rep-
resented by Heff is a collinearly ordered spin state, see
Fig. 17. Then, together with the singlet dimers between
the a- and b-spins of the downward pointing triangles, the
state depicted in Fig. 18 emerges as the candidate ground
state of the AKAF in the limit of large anisotropy J ≫ 1:
while nearest neighbour spins on the strongly coupled a-b
chains form singlets and decouple magnetically from the
spins on the c sites, the latter order in a collinear antifer-
romagnetic structure. This structure cannot be obtained
by a continuous deformation of the spiral IC structure
that is predicted by the Sp(N ) approach and is believed
to prevail for sufficiently large κ in the region of moder-
ate anisotropy. As a consequence of the trimerisation, the
state depicted in Fig. 18 breaks the translational symme-
try of our original model, Eq. (1). If this state survives as
the ground state of the non-trimerised model, i.e., when
the expansion parameter γ approaches unity, then, ow-
ing to their incompatible symmetries, the spiral IC phase
and the large J phase of our model must be separated by
a quantum phase transition.
Acknowledgments
One of the authors (HUE) acknowledges a useful discus-
sion with F. Mila. The work at the University of Water-
loo was supported by the Canada Research Chair (Tier
I, Michel Gingras). We thank M. Gingras for a critical
reading of the manuscript and numerous helpful sugges-
tions.
APPENDIX A: GROUND STATE DEGENERACY
FOR GENERAL J
We first derive the constraint on the chiralities that leads
to the reduction in the number of degenerate ground
states for general J relative to the special case J = 1.
Let χ1, · · · , χ6 be the chiralities of the six triangles sur-
rounding one of the hexagons of the kagomé lattice, and
let φ1, · · · , φ6 denote the angles that define the direc-
tions of the spin vectors on the six corners of the hexagon,
see Fig. 19.
Then, as is seen in Fig. 19, the following relations between
the angle φ1, and the angles φ2 · · · , φ6 are an immediate
consequence of these definitions:
φ2=φ1 − θ χ2 , (A1a)
φ3=φ1 − θ (χ2 + χ3) , (A1b)
φ4=φ1 − θ (χ2 + χ3)− (2π−2θ) χ4 , (A1c)
φ5=φ1 − θ (χ2 + χ3 + χ5)− (2π−2θ) χ4 , (A1d)
φ6=φ1 − θ (χ2 + χ3 + χ5 + χ6)− (2π−2θ) χ4, (A1e)
and φ6 = φ1 + (2π − 2θ)χ1 . (A1f)
From the last two of these relations it follows that the
chiralities χ1, · · · , χ6 are constrained by the sum rule
χ2 + χ3 + χ5 + χ6 − 2χ1 − 2χ4 = 0 . (A2)
For the isotropic kagomé system, J = 1, θ = 2π/3, one
finds instead of the constraint (A2) the sum rule
χj = n where n = 0, 1, 2 (A3)
which is obviously less restrictive than (A2).
Next, we present the arguments that lead to the estimate
NanisoGS (N▽)
N▽ with α < 3 (A4)
for the number NanisoGS (N▽) of classical GSs of an
anisotropic kagomé AF with N▽ downward pointing tri-
angles (the number of sites is 3N▽). Any planar config-
uration of a cell of the kagomé lattice can be constructed
by decorating the successive rows of up and down point-
ing triangles with chirality values χ = ±1 starting with
the first row. We consider only square cells with
rows with
N▽ downward pointing triangles. Then,
each row consists of 2
N▽ triangles, see Fig. 20.
Obviously, there are 22
N▽ ways to decorate the first
row. Disregarding certain exceptions, which will be dis-
cussed below, one can, for a given configuration of the
first row, choose the chirality of an arbitrary triangle
of the second row to be either +1 or −1. After this
choice has been made, the constraint (A2) fixes the chi-
ralities of all the remaining triangles of the second row
PSfrag replacements
FIG. 19: Chiralities around heaxagonal plaquette
FIG. 20: Example of a chirality distribution; dark and light
shaded triangles represent positive and negative chirality, re-
spectively. Chirality configurations in boxes fix the chirality
distribution of the row above them uniquely. An empty cir-
cle inside a triangle indicates that its chirality can be chosen
freely to be positive or negative.
uniquely. Proceeding in this manner from row to row
one would generate 22
N▽ · 2
N▽ distributions of chi-
ralities over the N▽ downward pointing triangles of the
cell. For finite lattice cells, the requirement of periodic
boundary conditions imposes further constraints on the
number of possible chirality distributions in these cells,
but the effect of these constraints will become negligible
in the thermodynamic limit N▽ → ∞. However, there
is a further reduction of the number of possible chirality
distributions: For a given distribution in a row it is not
always possible to find two distributions for the successive
row which both satisfy the constraint (A2). If in a row
the lower half of a hexagon of the next row is decorated
by chiralities in the manner − + − or + − + (see boxes
in Fig. 20), then the chiralities of the next row are fixed
uniquely. This reduces the number of possible chirality
distributions. Obviously, this reduction of the number
of possible chirality distributions survives in the thermo-
dynamic limit so that the exponent in (A4) is less than
N▽, the value one would have expected without this
reduction. We have calculated the number of distribu-
tions for cells of up to N▽ = 13× 13 and have found the
value α ≃ 2.18 for the constant in the expression (A4),
see Fig. 21.
As we have mentioned above, the sum rule (A3) which
applies for the isotropic kagomé AF is less restrictive than
the sum rule (A2). Consequently, the number of chirality
distributions in the isotropic model33,
N isoGS ∼ 1.18333N▽ (A5)
is larger than in the anisotropic model. Since the tran-
sition from the anisotropic model to the isotropic model
happens through a continuous variation of the coupling
constant J , there should be a continuous transition be-
tween the numbers of GS configurations in these two
cases. Presumably, this transition implies that the den-
sity of low-energy states of the anisotropic model in-
creases exponentially with an exponent ∼
N▽ so that
for J → 1 a sufficient number of states collapses to the
GS to bring about the transition between the laws (A4)
and (A5).
 2  4  6  8  10  12  14
linear fit
minPSfrag replacements
FIG. 21: (Color online) Number of chirality distributions,
NanisoGS , of cells of up to N▽ = 13 × 13. Dotted line: min =
N▽ ln 2 (lower bound); dashed line: max = 3
N▽ ln 2
(upper bound); full line: ln(NanisoGS ) = 0.65 + 2.18
N▽ ln 2
(linear fit to the numerical results).
APPENDIX B: PHASE BOUNDARIES
The FM phase and the DC phase are chacterised by
the vanishing of the parameters Q1 and Q2, respec-
tively. Our numerical results in section III C show that
at the respective phase boundaries, Q1 and Q2 decrease
to zero like order parameters at second order phase tran-
sitions. This suggests that we expand the mean field
energy EMF, Eq. (14), w.r.t. either Q1 or Q2 in the
manner of a Landau-Ginzburg (LG) expansion and de-
termine the phase boundaries and the properties of the
FM and the DC phase from this expansion. We write
EMF /(N▽N ) = e(α)LG(Qα) where
LG(Qα) = eα + rα |Qα|
2 + gα |Qα|4 +O(|(Qα|6) . (B1)
The coefficients eα, rα and gα are functions of the vari-
ables κ and J , of the parameters λa, λc and of Q2, Q1
for α = 1, 2, respectively. The saddle point of e
LG(Qα)
w. r. t. λa, λc and Qβ, β 6= α, determines the physical
values of these parameters. For e
LG(Qα) to qualify as
a bona fide Landau-Ginzburg energy describing a second
order phase transition with Qα playing the role of an or-
der parameter, the coefficients gα have to be positive at
the saddle point. For g1, i.e., inside and on the bound-
ary of the FM phase, this follows from the numerical
result: Q1 is found to remain zero for all J ≤ JF (κ). By
contrast, we have no numerical results for J ≥ JDC(κ),
i.e., inside and on the boundary of the DC phase. There-
fore, we need to show by analytic means that g2 > 0.
1. The FM phase and the FM-IC phase boundary
Since, as we have just remarked, we know that g1 > 0,
the remaining task is to determine the coefficients e1 and
r1 of e
LG. To this end, we have to expand the mean
field energy EMF, Eq. (14), w.r.t. Q1 which amounts to
expanding the frequencies ωµ(k) w.r.t. Q1. As can be in-
ferred from the expressions (16), (18) the frequencies de-
pend on Q1 only through the combination ε
2 = J2|Q̃1|2.
Therefore, we write the expansion in the form
ωµ(k; ε) = ω
µ (k) + ε
2ω(1)µ (k) +O(ε4)
ω̄(0)µ (k) + ε̄
2ω̄(1)µ (k) +O(ε̄4)
with λ+ = (λa + λc)/2 , ω̄
µ = ω
µ /λ+ , ω̄
µ (k) =
λ+ ∂ε2 ωµ(k; ε)|ε=0 and ε̄ = ε/λ+.
Here, the introduction of the “dimensionless” quantities
µ , ε̄ looks like an unneccessary complication but it will
help to keep expressions further below simple. Setting
Q1 = 0 in the matrix D̂(ω), Eq. (16), and solving Eq. (15)
for ω we find
1 (k) = wF (k) + δ , (B3a)
2 (k) = 1− δ , (B3b)
3 (k) = wF (k) − δ . (B3c)
δ = λ−/λ+ with λ− = (λc − λa)/2 (B4)
wF (k) =
1− q̄22[sin2(s2/2) + sin2(s3/2)] (B5)
q̄2 = |Q2|/λ+ (B6)
and sa = δak , a = 2, 3 (see Fig. 8) .
From our numerical results, Fig. 13, we know that λc >
λa and hence δ > 0. Therefore, ω̄
3 (k) < ω̄
1,2(k),
and hence, if the minimum of ω̄
3 (k) vanishes at the
point kmin in the Brillouin zone, ω̄
1,2(kmin) will be fi-
nite. Thus, since condensate can only occur when one of
the frequencies ω̄
µ , µ = 1, 2, 3 vanishes there may be
a finite condensate density |x3(kmin)|2, but the densities
|x1|2 and |x2|2 will certainly be zero. With these remarks
and with the above results for ω̄
µ we find from Eq. (14)
e1/λ+ = 2λ+q̄
2 − (3− δ)(κ+ 1)
1 (k) + ω̄
2 (k) + ω̄
3 (k)
3 (kmin) |x3(kmin)|2/N▽ . (B7)
Stationarity of e1 w.r.t. λ−, λ+, and q̄
2 (which is equiv-
alent to stationarity w.r.t. λa, λc), and Q
2 requires the
following three conditions to be fulfilled:
∂e1/∂λ− = 0 :
|x3(kmin)|2 = κ ; (B8)
∂e1/∂λ+ = 0 :
2λ+q̄2
2 − 3
κ− 1 +E2(q̄2) +
wF (kmin) = 0 ;
∂e1/∂q̄
2 = 0 :
2λ+ −
[K2(q̄2)−E2(q̄2)]−
wF (kmin)
= 0 ;
(B10)
with K2(q̄2) =
ds3 wF (k)
E2(q̄2) =
ds3 wF (k) .
(B11)
According to Eq. (B8), condensate must be present in
the FM region. This requires that ω̄
3 (kmin) vanishes.
From Eq. (B3c) it is seen that kmin = (−π, 0), so that
3 (kmin) = 0, if
wF (kmin) =
1− 2q̄22 =
. (B12)
Within the FM region and on the FM-IC boundary
(i.e., for Q1 = 0) the saddle-point values of q̄2, λ+,
and λ− are then determined as functions of κ by the
Eqs. (B9), (B10) and (B12). Remarkably, within this re-
gion these quantities are independent of the value of the
exchange constant J . The solution of these equations
shows that 0 ≤ q̄2 ≤ 2/3 for 0 < κ < ∞, cf. Figs. 11, 13.
The FM-IC phase boundary is the solution of r1(κ, J) =
0 (cf. Eq. (B1), where
r1 = ∂e
LG/∂Q
(B13)
with e
LG (EMF) from Eq. (14).
We obtain
s2 + s3
Ω(1)(k)
k→kmin
s2 + s3
3 (k)
(B14)
with Ω(1)(k) = −ω̄(1)1 (k) − ω̄
2 (k)− ω̄
3 (k).
To obtain the expansion coefficients ω̄
µ (k) which appear
in the last equation, we solve Eq. (15) to first order in
the expansion w.r.t. ε̄2. We find
Ω(1)(k)=
wF (k)(wF (k) + 1− 2δ)
wF (k) + 1
4 sin2
(2wF (k) + 1− δ)
(wF (k) + 1)(1− δ) (wF (k)2 − δ2)
(B15)
k→kmin
s2 + s3
3 (k)
= − 1
1 + δ
1− δ .
(B16)
With these results Eq. (B14) can, in the thermodynamic
limit, be cast into the form
− I3(q̄2)
1 + δ
1− δ =
JF (κ)
, (B17)
where
I3(q̄2) =
ds3 2 sin
Ω(1)(k) .
(B18)
Then, with q̄2 = q̄2(κ) and λ± = λ±(κ) as obtained from
Eqs. (B9), (B10) and (B12), the condition r1 = 0 is
an equation for the FM-IC phase boundary J = JF (κ)
which yields the graph shown in Fig. 9. As we have
mentioned above, inside the FM region, i.e., for J <
JF (κ), the saddle-point values of the quantities q̄ and λ±
and hence of Q2, λa, λc and |x3(kmin)| are independent
of the exchange coupling J , i.e., they retain the values
they attain on the FM-IC phase boundary, cf. Figs. 11,
2. The DC phase and the IC-DC phase boundary
Proceeding in exact analogy to the development in the
previous subsection we now expand EMF/(N▽ N ) in
powers of |Q2|2. However, instead of working with the
variables Q1, Q2, λa, λc we work with q1, Q2, λa, q2 here,
where
J |Q1|
, (B19a)
|Q2|√
. (B19b)
The replacement of |Q1| is purely a matter of conve-
nience. By contrast, the replacement of variables Q2, λc,
which according to the numerics vanish simultaneously
as J approaches the IC-DC phase boundary, by the pair
Q2, q2 leaves us with only one vanishing variable, since,
as will be seen below, q2 remains finite throughout.
a. Expansion of e
We write
ωµ(k) = ω
µ (k) + ω
µ (k)Q
2 + ω
µ (k)Q
2 +O(Q62)
(B20)
and determine the coefficients ω
µ , n = 1, . . . 4, by solv-
ing Eq. (15) for ω iteratively. We obtain
1 (k) + ω
2 (k) = 2λa wDC(k) , ω
3 = 0 , (B21a)
1 (k) + ω
2 (k) = −
1− cos kx cos ky
wDC(k)
, (B21b)
3 (k) =
q22 λa
C(k)2 −D(k)2
. (B21c)
Here,
C(k) = 1 − q22
1 − cos kx cos ky
wDC(k)2
, (B22)
D(k) = q1 sin k
x q22
cos kx − cos ky
wDC(k)2
, (B23)
wDC(k) =
1− q21 sin2 kx . (B24)
The coefficients ω
µ (k), µ = 1, 2 , 3, will only be needed
in the determination of the coefficient g2 of the fourth
order term of e
LG(Q2) which will be discussed later. We
will first concentrate on the determination of the zeroth
order term, e2, and of the coefficient r2 of the second
order term of e
LG(Q2). Under the assumption that g2 is
positive, this will provide us with an expression for the
IC-DC phase boundary.
With the above expressions for ω
1 + ω
2 and ω
ν = 0, 1, we obtain for the coefficients of the Landau
Ginzburg energy from Eqs. (14), (B1)
e2(q1, λa) =
λ2a q
− 2λa
1 + κ− 1
wDC(k)
(B25)
r2(q1, q2, λa, |x3(kmin)|2) =
(κ + 1)
1− cos kx cos ky
wDC(k)
3 (k) + |x3(kmin)|2 ω
3 (kmin)
. (B26)
These are valid for arbitrary values of the parameters q1,
λa, q2, and |x3(kmin)|. In the next subsection, we will
calculate their saddle point values for given Q2 and thus
fix the parameters. Here, we have only allowed for the
existence of a condensate component |x3(kmin)|2. This
is justified since, as Eqs. (B20) and (B2) show, ω3 < ω1,2
for sufficiently small Q2 so that conceivably ω3(k) may
vanish at some point kmin in the the Brillouin zone, while
ω1(k) and ω2(k) remain finite at kmin, and hence a finite
condensate density |x3(kmin)|2 may occur at this point.
b. Saddle point, phase boundary
Next we need to determine the saddle point of e
LG(Q2)
in the space of the variables q1, λa, q2, and |x3(kmin)|.
First, the saddle point values of q1 and λa are obtained
as expansions in powers of Q2,
a + λ
2 +O(Q42) , (B27a)
1 + q
2 +O(Q42) , (B27b)
where λ
a , q
1 are the solutions of
∂λae2=0 , (B28a)
∂q1e2=0 . (B28b)
Since the first derivatives of e2 vanish at λa = λ
a , q1 =
1 , Eqs. (B28), we have
e2 = e
2 + e
2 +O(Q62), (B29)
r2 = r
2 + r
2 +O(Q42) (B30)
Here, e
2 and r
2 are the expressions (B25) and (B26)
with λa and q1 replaced by λ
a and q
1 . The fourth
order term of e2, Eq. (B29), and the second order term
of r2 contribute only to the fourth order term of e
which will be determined later. Therefore, we postpone
the presentation of explicit expressions for λ
a , q
1 and
the ensuing expressions for e
2 and r
2 until later. With
e2 from Eq. (B25), Eqs. (B28) yield the equations
1 )− 1, (B31)
1 ) − E(q
, (B32)
which determine the saddle point values q
1 and λ
( K and E are the elliptic integrals of the first and the
second kind.)
Next we seek the extremum of e
LG w.r.t. q2. Since e2 is
independent of q2 we, neglecting terms of order Q
2, have
0=∂q2 r
q32 λa
κ + 1 − I1(q1, q2)
C(kmin)
|x3(kmin)|2
3 (kmin)
. (B33)
I1(q1, q2) =
∫ π/2
[C(k)2 −D(k)2]1/2
(B34)
(In these expressions and in the sequel, we use an abbre-
viated notation: λa, q1 and q2 denote the zeroth order
quantities λ
a , q
1 and q
2 .) kmin is the location of the
minimum of ω
3 (k),
min = 0 ;
kxmin
,−π ≤ kxmin ≤ −
(B35)
From (B21c) and (B35) it follows that ω
3 (kmin) = 0,
q22 = (1 − q21)/2. (B36)
As a function of q2 the integral I1(q1, q2) increases
monotonously,
1 = I1(q1, 0) ≤ I1(q1, q2) ≤ I1 (q1)
for 0 ≤ q2 ≤
(1 − q21)/2 . (B37)
We have defined
I1(q1) := max
I1(q1, q2) = I1
(1 − q21)/2
(B38)
As is seen in Fig. 22, the graphs of the functions
κ = κ(q1), Eq. (B31), and of I1 = I1(q1) intersect at
q1s ≃ 0.708, κs ≃ 0.181. Therefore, in solving Eq. (B33)
for q2, two cases have to be considered separately:
i. q1 > q1s, κ > κs. In this case, a solution ex-
ists only, if the last term in parentheses in Eq. (B33) is
positive. This requires that ω
3 (kmin) = 0 because, as
has been discussed before, |x3(kmin)| and hence the ra-
tio |x3(kmin)|2/ω(2)3 (kmin) would vanish otherwise. The
 0.95
 1.05
 1.15
 1.25
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1
PSfrag replacements
q1q1s
κ(q1) + 1
I1(q1)
FIG. 22: (Color online) I1 and κ+ 1 as functions of q1.
condition ω3(kmin) = 0 implies that q
2 = (1− q21)/2, cf.
Eq. (B36). Using this result and Eq. (B32) to eliminate
q2 and λa from Eq. (B26) we find
2 = 2
1− JDC(κ)
, (B39)
where
JDC(κ)=
(κ+ 1) (3− q21)/2 + Ĩ2(q1,
(1− q21)/2)
(1 − q21)
q21 π
4[K(q1)−E(q1)]
(B40)
Ĩ2(q1, q2) =
C(k)2 −D(k)2
(B41)
is the IC-DC phase boundary for κ > κs, i.e., in the
region where the ratio |x3(kmin)|2/ω(2)3 (kmin) is finite.
According to the discussion at the end of section III C,
cf Eq. (24), this is the region where LRO prevails along
the decoupled chains, cf. Fig. 9.
In the development leading to Eq. (B40) for the phase
boundary, we have not needed the solution of Eq. (B40)
explicitly, but we note it here for completeness:
C(kmin)
|x3(kmin)|2
3 (kmin)
|x3|2
1+3q21
1−q21
=1 + κ− I1(q1) > 0. (B42)
These relations show that while |x3(kmin)| = 0, the ratio
|x3(kmin)|2/ω(2)3 (kmin) remains finite.
ii. q1 < q1s, κ < κs. In this case, we must have
I1(q1, q2)) < I1
(1 − q21)/2
, (B43)
(see Eq. (B38)). Consequently q22 < (1 − q21)/2 so that
3 (kmin) > 0 and hence no condensate can develop,
|x3|2 = 0. Then, Eq. (B33) yields the equation
I1(q1, q2) = 1 + κ (B44)
which replaces Eq. (B36) and determines q2 as a function
of q1, q2 = q2(q1). Then, proceeding as in case (i) one
finds for the IC-DC phase boundary in the region κ < κs
JDC(κ)=
(κ+ 1)(1 + q22) + Ĩ2(q1, q2)
q21 π
4 [K(q1)−E(q1)]
. (B45)
Here, q1 = q1(κ) from Eq. (B31) and q2 = q2(κ) from
Eq. (B44)(with q1 = q1(κ)).
We note here that inside the DC phase, i.e., for J >
JDC(κ), where Q2 = λc = 0, the saddle-point values
of q1 and λa/J and hence of Q1 are independent of J ,
cf. Eqs. (B31), (B32). Hence the graphs of Q1 and λa for
J < JDC and for J > JDC join smoothly at J = JDC , cf.
Figs. 10, 13. Furthermore, it follows from Eq. (B42), that
the ratio (|x3(±kmin)|2/N▽)/(λa q22ω
3 (kmin)), which
occurs in the amplitude of the spin-spin correlation func-
tion, cf. Eq. (24a), is also independent of J inside the DC
phase and retains the value that it has attained at the
IC-DC transition line.
c. Stability of the phase boundary
In deriving the phase boundary from the condition r
0 we have tacitly assumed that the coefficient g2 of the
fourth order term in the LG expansion, Eq. (B1), is pos-
itive. In the remaining part of this appendix we will
sketch the steps which lead to the conclusion that this is
indeed the case.
Expanding in the expression (B1) for e
LG the coefficients
e2 and r2 w.r.t. the second order contributions to q1 and
λa, q
1 and λ
a , cf. Eqs. (B27) we obtain
LG = e
2 + r
2 + (g2 + g
2 +O(Q62) , (B46)
where
1 + ω
2 + ω
(B47)
is the contribution to the fourth order term of e
LG that
arises from the fourth order terms of the frequencies ωµ
in the sum in Eq. (14) whereas the contribution to e
of the expansion of e2 and r2 is
∂2q1e2|0 ∂q1∂λae2|0
∂λa∂q1e2|0 ∂2λae2|0
∂q1r2|0
∂λar2|0
. (B48)
(In Eq. (B48) the notations ∂2q1e2|0 etc. indicate that af-
ter the derivatives have been taken the variables q1, λa
etc. have to be replaced by their zeroth order values q
a etc.)
The evaluation of the contribution (B47) is straight-
foward: the coefficients ω
µ , µ = 1, 2, 3, were obtained
by solving Eq. (15) for ω iteratively to fourth order. As
the explicit expressions are rather lengthy and contain
no direct information, we refrain from presenting them
here. The sum over k that is required in Eq. (B47) was
done numerically. g2 was obtained in the form
g̃2(q1) , (B49)
where g̃2(q1) is a function of q1 alone which is always
positive so that g2 > 0 throughout. Remarkably, no
explicit dependence on the coupling constant J appears
in these results.
The evaluation of g′, Eq. (B48) requires the knowledge of
explicit expressions for q
1 and λ
a . These are obtained
by expanding e2 to first order in q
1 and λ
a , inserting
the results into the expression (B1) for e
LG and requiring
that the terms of order Q22 satisfy the extremum condi-
tions w.r.t. q1 and λa:
0 = q
e2|0 + λ(2)a ∂q1∂λae2|0 + ∂q1r2|0 ,(B50a)
0 = q
1 ∂q1∂λae2|0 + λ(2)a ∂2λae2|0 + ∂λar2|0 .(B50b)
The solution of these equations reads
= −M̂
∂q1r2|0
∂λar2|0
, (B51)
M̂−1 =
∂2q1e2|0 ∂q1∂λae2|0
∂λa∂q1e2|0 ∂2λae2|0
. (B52)
Inserting these results into Eq. (B48) one finds
g′2 = −
∂q1r2|0 ∂λar2|0
∂q1r2|0
∂λar2|0
. (B53)
While the second derivatives of e2 are obtained straight-
forwardly from Eq. (B25) the derivatives ∂q1r2|0 and
∂λar2|0 have to be calculated separately for the region
q1 < q1s, where there is no condensate, |x3(kmin)|2 = 0,
and for the region q1 > q1s, where |x3(kmin)|2 > 0. Fi-
nally, the result for g′ can be cast into the form
g′2 =
xq xλ
(B54)
where
M̂ ′ =
4q21 Λ(κ+ 1− Λ)
1− q21
−Λ (2− q21)Λ− κ− 1
(B55)
Λ ≡ 1
[K(q1)−E(q1)] (K, E : elliptic integrals)
(B56)
xq = q
1(Λ− κ− 1) +
1− q21
q1∂q1 Ĩ2(q1, q2)
q2=q2(q1)
−Θ(q1 − q1s)
1− q21
(κ+ 1− I1(q1)) ,
xλ = (
+ 1)(κ+ 1)− 1
Ĩ2(q1, q2)
q2=q2(q1)
. (B57)
Here Θ is the step function; the integrals I1(q1) and
Ĩ2(q1, q2) have been defined above, cf. Eqs. (B38) and
(B41), respectively. After numerical evaluation of these
integrals, we find that g′2 = g
2(q1) is positive for all values
of q1.
1 G. Misguich and C. Lhuillier, Frustration in Two-
Dimensional Quantum Antiferromagnets (World Scien-
tific, Singapore, 2004), chap. 5.
2 R. Moessner, Can. J. Phys. 79, 1283 (2001).
3 P. Lecheminant, B. Bernu, C. Lhuillier, L. Pierre, and
P. Sindzingre, Phys. Rev. B 56, 2521 (1997).
4 C. Waldtmann, H.-U. Everts, B. Bernu, C. Lhuillier,
P. Sindzingre, P. Lecheminant, and L. Pierre, Eur. Phys.
J. B 2, 501 (1998).
5 F. Mila, Phys. Rev. Lett. 81, 2356 (1998).
6 M. Mambrini and F. Mila, Eur. Phys. J. B 17, 651 (2000).
7 Z. Hiroi, M. Hanawa, N. Kobayashi, M. Nohara, H. Takagi,
Y. Kato, and M. Takigawa, J. Phys. Soc. Jpn. 70, 3377
(2001).
8 A. Fukaya, Y. Fudamoto, I. M. Gat, T. Ito, M. I. Larkin,
A. T. Savici, Y. J. Uemura, P. P. Kyriaku, G. M. Luke,
M. T. Rovers, et al., Phys. Rev. Lett. 91, 207603 (2003).
9 F. Bert, D. Bono, P. Mendels, J.-C. Trombe, P. Millet,
A. Amato, C. Baines, and A. Hillier, J. Phys. Condens.
Matter 16, 9829 (2004).
10 F. Bert, D. Bono, P. Mendels, F. Ladieu, F. Duc, J.-C.
Trombe, and P. Millet, Phys. Rev. Lett. 95, 087203 (2005).
11 S. Sachdev and N. Read, Int. J. Mod. Phys. B 5, 219
(1991).
12 S. Sachdev, Phys. Rev. B 45, 12 377 (1992).
13 W. Apel, T. Yavors’kii, and H.-U. Everts, J. Phys. Con-
dens. Matter 19, 145255 (2007); W. Apel, T. Yavors’kii,
and H.-U. Everts, ibid., Corrigendum, in press.
14 O. A. Starykh and L. Balents, Phys. Rev. Lett. 98, 077205
(2007).
15 E. Lieb and D. Mattis, J. Math. Phys. 3, 749 (1962).
16 S. Brehmer, H.-J. Mikeska, and S. Yamamoto, J. Phys.
Condens. Matter 9, 3921 (1997).
17 A. B. Harris, C. Kallin, and A. J. Berlinsky, Phys. Rev. B
45, 2899 (1992).
18 C. H. Chung, J. B. Marston, and S. Sachdev, Phys. Rev.
B 64, 134407 (2001).
19 C. H. Chung, J.B.Marston, and R. H. McKenzie, J. Phys.
Condens. Matter 13, 5159 (2001).
20 F. Wang and A. Vishwanath, Phys. Rev. B 74, 174423
(2006).
21 O. Tchernyshyov, R. Moessner, and S. L. Sondhi, Euro-
phys. Lett. 73, 278 (2006).
22 G. Misguich, B. Bernu, and C. Lhuillier, J. Low Temp.
Phys. 110, 327 (1998).
23 C. Waldtmann, H. Kreutzmann, U. Schollwöck,
K. Maisinger, and H.-U. Everts, Phys. Rev. B 62,
9472 (2000).
24 V. Subrahmanyam, Phys. Rev. B 52, 1133 (1995).
25 C. Raghu, I. Rudra, S. Ramasesha, and D. Sen, Phys. Rev.
B 62, 9484 (2000).
26 M. E. Zhitomirsky, Phys. Rev. B 71, 214413 (2005).
27 P. Nikolic and T. Senthil, Phys. Rev. B 68, 214415 (2003).
28 J. B. Marston and C. Zeng, J. Appl. Phys. 69, 5962 (1991).
29 C. Zeng and V. Elser, Phys. Rev. B 51, 8318 (1995).
30 A. V. Syromyatnikov and S. V. Maleyev, Phys. Rev. B 66,
132408 (2002).
31 R. Budnik and A. Auerbach, Phys. Rev. Lett. 93, 187205
(2004).
32 G. Misguich and P. Sindzingre, J. Phys. Condens. Matter
19, 145202 (2007).
33 R. J. Baxter, J. Math. Phys. 11, 784 (1970).

  We study the properties of the Heisenberg antiferromagnet with spatially
anisotropic nearest-neighbour exchange couplings on the kagome net, i.e. with
coupling J in one lattice direction and couplings J' along the other two
directions. For J/J' > 1, this model is believed to describe the magnetic
properties of the mineral volborthite. In the classical limit, it exhibits two
kinds of ground states: a ferrimagnetic state for J/J' < 1/2 and a large
manifold of canted spin states for J/J' > 1/2. To include quantum effects
self-consistently, we investigate the Sp(N) symmetric generalisation of the
original SU(2) symmetric model in the large-N limit. In addition to the
dependence on the anisotropy, the Sp(N) symmetric model depends on a parameter
kappa that measures the importance of quantum effects. Our numerical
calculations reveal that in the kappa-J/J' plane, the system shows a rich phase
diagram containing a ferrimagnetic phase, an incommensurate phase, and a
decoupled chain phase, the latter two with short- and long-range order. We
corroborate these results by showing that the boundaries between the various
phases and several other features of the Sp(N) phase diagram can be determined
by analytical calculations. Finally, the application of a block-spin
perturbation expansion to the trimerised version of the original spin-1/2 model
leads us to suggest that in the limit of strong anisotropy, J/J' >> 1, the
ground state of the original model is a collinearly ordered antiferromagnet,
which is separated from the incommensurate state by a quantum phase transition.

arXiv:0704.1441v3  [cond-mat.str-el]  6 Aug 2007
Heisenberg antiferromagnet with anisotropic exchange on the kagomé lattice:
Description of the magnetic properties of volborthite
T. Yavors’kii,1 W. Apel,2 and H.-U. Everts3
Department of Physics and Astronomy, University of Waterloo,
200 University Avenue W, Waterloo, N2L 3G1, Canada.
Physikalisch-Technische Bundesanstalt, Bundesallee 100, D-38116 Braunschweig, Germany.
Institut für Theoretische Physik, Leibniz Universität Hannover, Appelstraße 2, D-30167 Hannover, Germany.
(Dated: September 12, 2021)
We study the properties of the Heisenberg antiferromagnet with spatially anisotropic nearest-
neighbour exchange couplings on the kagomé net, i.e., with coupling J in one lattice direction and
couplings J ′ along the other two directions. For J/J ′ >∼ 1, this model is believed to describe the
magnetic properties of the mineral volborthite. In the classical limit, it exhibits two kinds of ground
states: a ferrimagnetic state for J/J ′ < 1/2 and a large manifold of canted spin states for J/J ′ > 1/2.
To include quantum effects self-consistently, we investigate the Sp(N ) symmetric generalisation of
the original SU(2) symmetric model in the large-N limit. In addition to the dependence on the
anisotropy, the Sp(N ) symmetric model depends on a parameter κ that measures the importance of
quantum effects. Our numerical calculations reveal that in the κ-J/J ′ plane, the system shows a rich
phase diagram containing a ferrimagnetic phase, an incommensurate phase, and a decoupled chain
phase, the latter two with short- and long-range order. We corroborate these results by showing that
the boundaries between the various phases and several other features of the Sp(N ) phase diagram
can be determined by analytical calculations. Finally, the application of a block-spin perturbation
expansion to the trimerised version of the original spin-1/2 model leads us to suggest that in the
limit of strong anisotropy, J/J ′ ≫ 1, the ground state of the original model is a collinearly ordered
antiferromagnet, which is separated from the incommensurate state by a quantum phase transition.
PACS numbers: 75.10.Jm,75.30.Kz,75.50.Ee
I. INTRODUCTION
In the ongoing search for novel states of condensed mat-
ter, frustrated antiferromagnets have played a key role
(for a recent review, see Ref. 1). Among the many sub-
stances that have been investigated experimentally and
the numerous spin models that have been studied theo-
retically, those in which the magnetic ions occupy the ver-
tices of corner-sharing frustrating entities have attracted
particular attention in this context. The best known ex-
amples are the kagomé antiferromagnet (KAF), consist-
ing of corner sharing triangles, and the pyrochlore anti-
ferromagnet, consisting of corner sharing tetrahedra (see
Fig. 1).
The main distinction between the KAF, the pyrochlore
antiferromagnet and other frustrated and unfrustrated
AAAA AAAA
AAAAAAAA
(a) (b)
FIG. 1: kagomé lattice (a), pyrochlore lattice (b)
magnets is the large ground-state degeneracy of the for-
mer: classical Heisenberg antiferromagnets with nearest-
neighbour interactions on corner-sharing lattices have a
large ground-state degeneracy, which in the above two ex-
amples even leads to a finite ground-state entropy (see,
e.g., Ref. 2 and references therein). Quantum effects
may lift this degeneracy, and, indeed, in numerical stud-
ies of small cells of the spin 1
KAF, an exponentially
large number of very low-lying quantum states has been
observed3,4. It has been suggested that this abundance
of low-lying states can be understood in a description of
the low-energy physics of the quantumKAF as spin liquid
consisting of nearest-neighbour spin singlets5,6. However,
a complete picture of the ground state and of the excita-
tions of the KAF is still missing. Further theoretical, but
also experimental studies with emphasis on the quantum
properties of the KAF are therefore highly desirable. In
this last respect, the mineral volborthite is a very promis-
ing candidate. It has been the subject of several recent
experimental investigations7–10. The magnetic lattice of
this natural antiferromagnet consists of the S = 1/2 spins
of Cu2+ ions that are located on the vertices of well sep-
arated planar kagomé-like nets. A monoclinic distortion
of the lattice leads to a slight difference between the ex-
change couplings along one lattice direction (J) and the
two other directions (J ′)(see Fig. 2). Since neither signs
of long-range order nor signs of a spin-gapped singlet
ground-state were found in experiments on volborthite,
the substance seems to be a good candidate for the ob-
servation of the low-energy features that are thought to
http://arxiv.org/abs/0704.1441v3
PSfrag replacements
[i, j]
〈i, k
FIG. 2: (Color online) Anisotropic kagomé model. The cou-
pling J ′ and the nearest neighbour distance will be set equal
to unity in the calculations. δ1(‖ êx), δ2, and δ3 are the three
primitive lattice vectors of the kagomé net.
be typical for kagomé type antiferromagnets1.
Whether and to what extent the different exchange cou-
plings along different lattice directions of the kagomé net
of volborthite influence the low-energy physics of the sys-
tem is presently unknown. In the present paper, we study
this question on the basis of the model Hamiltonian
HAKAF = J
[i,j]
SiSj + J
〈k,i〉
SkSi . (1)
The symbols [i, j] and 〈k, i〉 denote, respectively, bonds
between nearest-neighbour sites on the horizontal chains
(a, b) and bonds between the middle sites (c) and the
sites a, b, see Fig. 2. Since the physics of this model de-
pends only on the ratio J/J ′ of the exchange constants,
we set J ′ = 1 in the sequel. We will consider the spatially
anisotropic kagomé antiferromagnet (AKAF), Eq. (1), in
the full range of J , 0 < J < ∞ since this is of theoretical
interest: one expects to see quantum phase transitions
as J is increased. It is of particular interest to find out
whether there is a transition from two-dimensional mag-
netic states to a set of decoupled chains with free spins
on the axes between the chains for large values of J .
The paper is organised as follows. In Sec. II, we consider
the model (1) in the classical limit. At this level, we find
no sign of a transition from the two-dimensional magnet
to a set of decoupled chains as J increases to infinity.
Nonetheless, the ground-state degeneracy, as well as the
spin wave spectrum are found to change qualitatively
as the anisotropy of the model varies. In Sec. III, we
consider a generalisation of the SU(2) symmetric model
(1) to the Sp(N ) symmetric version11,12 and describe
its properties in the large-N limit, where a mean-field
treatment of the model is adequate. We obtain a de-
tailed description of how possible ground states of the
model depend on the coupling J and on the spin length
S. A fairly rich phase diagram with a ferrimagnetic
phase for small J , long-range ordered and short-ranged
incommensurate phases for intermediate values of J , and
a decoupled-chain phase for large J emerges. Parts of
these results have been published previously, see Ref. 13.
In Sec. IV, we devise trial quantum ground states of the
original S = 1/2 model. We chose the states such that
they are exact eigenstates of HAKAF, if the couplings on
the upward pointing triangles of Fig. 2 are switched off,
and we then treat these couplings perturbatively. In the
limit J −→ ∞ this yields an effective Hamiltonian for
the spins on the c sites which represents an anisotropic
triangular antiferromagnet. The conclusions of Starykh
and Balents14 about the ground state of this effective
model lead us to conjecture the existence of a quantum
phase transition in the AKAF for large J . In Sec. V, we
summarise and discuss our results. In two Appendices,
we present technical details of the counting procedure for
the classical ground-states, and of the Ginzburg-Landau
type procedure that allows us to determine the bound-
aries in the phase diagram analytically.
II. CLASSICAL AND SEMICLASSICAL
ASPECTS
Similar to other isotropic spin models on lattices with
triangular elementary cells, the classical ground states of
HAKAF, Eq. (1), are spin configurations, which satisfy the
condition that for each elementary triangular plaquette
of the lattice, Fig. 2, the energy is minimal.
For J = 0, this yields a ferrimagnetic state with the chain
spins aligned in one direction and the middle spins point-
ing in the opposite direction, so that the total magneti-
sation is M = N▽S (N▽: number of downward pointing
triangles, N▽ = Ns/3 where Ns is the number of sites of
the system). We illustrate this situation in Fig. 3. Ac-
cording to the Lieb-Mattis theorem, the exact quantum
ground state (GS) of the model HAKAF also has total
spin Stot = N▽ S for J = 0, see Ref. 15, i.e., for J = 0,
the quantum GS is ferrimagnetic too. By continuity, one
expects the quantum GS to remain ferrimagnetic for suf-
ficiently small finite J . This will be confirmed by our
considerations of the large-N limit of the Sp(N ) version
of our model (see the analytical and numerical work in
Sects. III, III C and Appendix B) and by the block spin
perturbation approach (Sec. IV). Classically, the ferri-
FIG. 3: (Color online) Ferrimagnetic state for J = 0,
i.e., when there is no coupling between chain spins, cf. Fig. 2.
ωα(q)
FIG. 4: (Color online) Spin-wave frequencies ωα(q), α =
1, 2, 3 for J = 0.4; the contour at the top of the plot marks
half the Brillouin zone.
magnetic state remains stable up to J = 1/2. The ex-
citation spectrum of the ferrimagnetic state obtained in
linear spin-wave (LSW) approximation is shown in Fig. 4.
The analytic expressions for these three frequency sur-
faces ωα(q), α = 1, 2, 3, are obtained as solutions of a
third-order secular equation and are too lengthy to be
presented here. However, one can easily assure oneself
that the dispersion of the gapless mode is quadratic at
the origin. Thus, one has the typical mode structure of
a ferrimagnet here with one ferromagnetic mode and two
optical modes, see, e.g., Ref. 16. As J increases towards
1/2, the ferromagnetic frequency surface looses its disper-
sion and turns into a plane of zero modes, one zero mode
for each wave vector in the magnetic Brillouin zone (BZ),
at J = 1/2. The gap of the lower optical mode closes at
this value of J in the centre of the BZ and the dispersion
of this mode becomes linear for small wave vectors as for
an antiferromagnetic spin-wave mode.
At J = 1/2, the classical GS configuration changes from
the unique ferrimagnetic state to an ensemble of degener-
ate canted coplanar states. These states are characterised
by two variables: the angle θ, which the middle spin of
a given triangular plaquette forms with the two chain
spins of the same plaquette (see Fig. 5), and the two
valued chirality χ = ±1, which denotes the direction in
which the spins turn as one moves around the plaquette
in the mathematically positive sense.
For J ≥ 1/2, the requirement that the energy of any of
the elementary triangular plaquettes of the lattice Fig. 2
be minimal is θ = arccos(−1/(2J)), (θ > 0). The differ-
ent degenerate canted states arise from different possibil-
FIG. 5: (Color online) Canted spins of the AKAF at J > 1/2.
ωα(q)
FIG. 6: (Color online) Same as Fig.4 for J = 0.6.
ities to assign positive or negative chiralities to the pla-
quettes of the lattice. We show in the Appendix A that
for the general case of θ 6= 2π/3 (J 6= 1), the number of
spin configurations, NanisoGS does not grow exponentially
with the number of sites. Rather, NanisoGS < 2
where α < 3. This implies that the ground-state en-
tropy per spin of the classical AKAF vanishes in the
thermodynamic limit. In this respect, the anisotropic
model differs qualitatively from the isotropic KAF in
the classical limit, which has an extensive entropy per
spin. In the limit J → 1, the anisotropic model ap-
proaches the isotropic KAF. Hence, one expects that for
the anisotropic model there is an extensive number of
low-lying excited states that become degenerate with the
GS in the isotropic limit.
As in the case of the isotropic KAF, the spin-wave Hamil-
tonian is in linear order independent of the particular
classical GS that has been chosen as the starting point
of the expansion, Ref. 17. This implies that lowest-order
quantum fluctuation do not select one or a group of clas-
sical GSs as true GSs, i.e., the possible ordering effects
of quantum fluctuations are not captured by the linear
spin-wave (LSW) approximation. Figs. 6, 7 show the
spin-wave frequency surfaces for J = 0.6 and for J = 3. It
is easy to show analytically that, as is illustrated in these
figures, the plane of zero frequency modes persists for all
values of J greater than 1/2. The surfaces for J < 1/2
and for J > 1/2 join smoothly at J = 1/2. Thus, in the
ωα(q)
FIG. 7: (Color online) Same as Fig.4 for J = 3.
LSW approximation, the transition from the ferrimag-
netically ordered state to the canted spin states appears
to be of second order. For J ≫ 1, the nonzero frequen-
cies gradually loose their dispersion perpendicular to the
strong-J direction and take the shape of the spin-wave
spectrum of antiferromagnetic chains parallel to this di-
rection. However, no sign of a further transition from
the canted spin states to a set of decoupled spin chains is
found in this semiclassical approach. In the next section,
we will consider the symplectic Sp(N ) generalisation of
the antiferromagnetic model HAKAF in the large-N limit.
This approach, which was first proposed by Read and
Sachdev, Refs. 11,12, as a method to study frustrated
antiferromagnets, has the benefit of including the or-
dering effects of quantum fluctuations self-consistently.
It is of particular interest for spin models with two or
more competing exchange couplings in the different lat-
tice directions or over different lattice distances such as
the present model, the J1-J2-J3 model
11, the Shastry-
Sutherland antiferromagnet18 and the anisotropic trian-
gular antiferromagnet19. For these models, it has pro-
vided an unbiased selection of possible GSs that may or
may not be ordered depending on the value of a param-
eter κ, which is connected with the spin length S (see
below).
III. MEAN FIELD Sp(N ) APPROACH
A. Brief review of the method
For a general antiferromagnetic Heisenberg model with a
positive interaction matrix Jij ,
Jij Si · Sj , (2)
the Sp(N ) generalisation reads
HSp(N ) = −
2N (J
jβ)(Jγδb
j) . (3)
Here,
. . .
is the 2N × 2N generalisation of the 2× 2 antisymmetric
tensor
, (5)
and bαi with α = 1, . . . , 2N are the Sp(N ) boson anni-
hilation operators. (Here and in the sequel, we closely
follow the notation of Ref. 12; in particular, summation
over repeated upper and lower indices is implied.) Thus,
J αβb†iαb
jβ is the generalisation of the creation operator
jβ for a singlet on the bond (i, j). For the special
case N = 1, one finds
(J αβb†iαb
jβ)(Jγδb
j) = −2Si · Sj + nbinbj/2 + δijnbi ,
where
nbi = b
i (7)
is the boson number operator at site i and where
Si = b
i /2 (8)
is the usual SU(2) spin operator at site i. (τ are the
Pauli matrices). Then, if one imposes the constraint
that the number of bosons is the same for all lattice sites,
nbi ≡ nb, the HamiltonianHSp(1) is the familiar SU(2) in-
variant antiferromagnetic Heisenberg Hamiltonian (plus
some constants) with nb = 2S.
In the subsequent exposition, we shall consider a Hamil-
tonian of the form (3) in the large-N limit following the
strategy of Refs. 11,12. Depending on the values of the
couplings Jij and of κ, the GS of HSp(N ) may either
break the global Sp(N ) symmetry and exhibit LRO or
it may be Sp(N ) symmetric with only SRO. Breaking of
the Sp(N ) symmetry will happen through condensation,
i.e., by macroscopic occupation of one of the Bose fields
bα. To allow for this, we introduce the parametrisation
bmσi =
b̃m̃σi
with α = (mσ), m = 1, · · · ,N , m̃ = 2, · · · ,N and
σ =↑, ↓. The field xσi is proportional to the condensate
amplitude, 〈bmσi 〉 =
N δm1 xσi . Aiming at a mean field
treatment of the HamiltonianHSp(N ), which becomes ex-
act in the large N limit, we decouple the quartic part by
the Hubbard-Stratonovich technique with complex fields
Qij = −Qji and with Lagrange multipliers λi that en-
force the local constraints (7). The variables Qij which
are defined on nearest neighbour bonds of the lattice are
expectation values of the bond singlet creation opera-
tors in the GS, Qij = 〈
σσ′ ε
j mσ′〉 and are to
be determined self-consistently from the mean field type
Hamiltonian
HMF =
Jij |Qij |2 −
Qijεσσ′
Nxσi xσ
b̃m̃σi b̃
+ h.c.
N|xσi |2 +
im̃σ b̃
i − nb
The variational ground state energy, EMF, of HMF in the
large-N limit is obtained by diagonalising the bosonic
part of HMF, integrating over the 2(N − 1)Ns bosonic
fields b̃m̃σi in the action associated with HMF. One ob-
tains:
Jij |Qij |2 −
Qijεσσ′x
j + h.c.
ωµ(k;Q, λ) +
|xσi |2 − 1− κ
. (11)
Here ωµ(k;Q, λ) are the positive eigenvalues of the
bosonic part of HMF, and κ = nb/N is kept fixed in the
limiting procedure11,12. The parameter κ is a measure
for the importance of quantum fluctuations: by vary-
ing κ from small to large values, one drives the system
from the regime dominated by quantum fluctuations to
the classical regime, i.e., from the disordered into the or-
dered region. Finally, the GS is obtained by finding the
saddlepoint of EMF in the space of the variables Qij and
xσi subject to the constraints
∂EMF(Q, λ)/∂λi = 0 . (12)
In addition to the GS itself, the spin-spin correlation
function Gij = 〈Si ·Sj〉 in the GS is an important piece of
information. In particular, by considering its behaviour
in the limit |i − j| −→ ∞, one can distinguish between
LRO and SRO. According to Sachdev12, to obtain Gij in
the Sp(N ) symmetric approach, the SU(2) invariant ex-
pression Si ·Sj must be replaced by the Sp(N ) invariant
expression
4N 2 (b
j − J αγJβδb
j). (13)
Within the mean field approach, Gij can then be calcu-
lated straightforwardly.
B. The anisotropic kagomé antiferromagnet
1. Choice of mean-field variables
We wish to apply the procedure described above to the
AKAF represented by the Hamiltonian (1). To render
the problem of finding the eigenvalues ωµ in Eq. (10)
and of optimising EMF tractable, we have to restrict the
number of variables Qij and λi. We do so by demand-
ing that the mean field Hamiltonian HMF for the spinon
operators b(†) is symmetric under transformations of the
projective symmetry group (PSG) that is related to the
symmetry group of the spin Hamiltonian HAKAF (Eq. 1)
(see Ref. 20). We include two translations, a rotation
by π and a mirror axis orthogonal to the preferred direc-
tion of the exchange constants (J). Thus generalising the
treatment of Wang and Vishwanath to our model, we find
eight mean-field states with different symmetries. Seven
of them have flux in the sense of Ref. 21 in various cells
of the lattice. Following the arguments in Ref. 21, we ex-
clude all flux-carrying states and end up with the solution
(cf. Fig. 8) P1,2,3 = Q1,2,3, Q3 = Q2, and λb = λa.
PSfrag replacements Q1
FIG. 8: (Color online) Arrangement of mean field parameters:
Q1 ≡ Qab, Q2 ≡ Qbc, and Q3 ≡ Qca denote the intra triangle
bonds, P1 ≡ Qba′ , P2 ≡ Qcb′ and P3 ≡ Qac′ denote the inter
triangle bonds. λa, λb, and λc are the Lagrange multipliers
needed to implement the constraints on the sites a, b, and c.
In order to check the flux-argument in Ref. 21, we have
explicitely studied the solution P1,2,3 = −Q1,2,3 and
found that it is always of higher energy (For J = 1, this
agrees with the result of Ref. 12).
Thus, the expression Eq. (11) can now be cast into the
= J |Q1|2 + 2 |Q2|2 − (2λa + λc)(κ+ 1)
ωµ(k)
1 + |xµ(k)|2
, (14)
where the condensate is written in diagonalized form and
ωµ(k) are the three positive solutions of
det D̂(ω) = 0 . (15)
Here,
D̂(ω) =
Λ̂− ωÎ Q̂
Q̂† Λ̂+ ωÎ
, (16)
Λ̂= diag(λa, λc, λa) , (17)
0 Q̃2(k) −JQ̃1(−k)
−Q̃2(−k) 0 Q̃3(k)
JQ̃1(k) −Q̃3(−k) 0
 , (18)
and Q̃a(k) =
eiδak/2 − e−iδak/2
, a = 1, 2, 3,
δ1,2,3, see Fig. 2.
2. Technical details of the numerical extremalisation
Determination of the ground state of the AKAF in the
considered approximation has been reduced to minimiza-
tion of the Eq. (14) with respect to two variables Q1 and
Q2, subject to the Lagrange constraints with respect to
two parameters λa and λc. Being apparently trivial, the
optimization procedure turns out to be quite involved
technically.
First, we find it crucial to consider at least two differ-
ent chemical potentials. Other than for the spatially
isotropic KAF, J = 1, we were not able to find a non-
trivial solution if we used a single λ, λa = λb = λc. If λa
and λc are different, [Λ̂, Q̂] 6= 0, the Lagrange multipliers
enter the expressions for the frequencies ωµ non-trivially,
other than in the case of a global uniform chemical po-
tential (cf. Ref. 12). In turn, the Lagrange constraints
cannot be satisfied semi-analytically, and require a nu-
merical treatment. Second, we choose to work directly
in the thermodynamic limit Ns → ∞ of the model (14)
by performing a numerical self-adapting integration over
the BZ. In this limit, the singularities can be integrated,
and symmetry breaking is signalled by the appearance of
a finite value of the condensate amplitude xµ(k) at a cer-
tain wavevector k = qord, which characterises the type of
magnetic order. We mention here that the extremalisa-
tion of a mean-field energy of the type of Eq. (14) can also
be achieved by solving the pertinent stationarity condi-
tions numerically for finite systems, i.e., for finite Ns,
see e.g., Ref. 22. Then, the type of magnetic order has
to be detected by calculating the structure factor. Third,
we see that the Eq. (14) has a minimum with respect to
the physical bond parameters Q1 and Q2 only after the
elimination of the chemical potentials. In the full Q− λ
space we face an extremalization problem.
Technically, we find it convenient to use a polar coordi-
nate parametrisation for the variablesQ1, Q2 and λa, λc:
Q1 = Q cos(α), Q2 = Q sin(α) , (20)
λa = Λ sin(β), λc = Λcos(β) . (21)
We perform an optimization with respect to the variables
Q,Λ, α, β, as well as condensate densities xµ(k) in accord
with the following algorithm (J and κ are kept fixed).
i. We fix the angles α, β and the amplitude Q, and
first exploit the stationarity condition for EMF with re-
spect to Λ. It is convenient to write the corresponding
equation in the following form:
[2 sin(β) + cos(β)] (κ+ 1)
|xµ(k)|2 ∂Λωµ(k)
∂Λ ωµ(k) , (22)
where Ω = 8π2/
3 is the volume of the unit cell. One
finds that Q and Λ enter the Eq. (22) only via the ratio
ξ = Λ/Q.
The requirement that the frequencies must be positive,
ωµ(k) ≥ 0, defines a lower limit ξmin(α, β) for ξ: the
frequencies ωµ(k) are positive for ξ > ξmin(α, β); for
ξ = ξmin(α, β), the lowest mode ωµ0 vanishes at some
point(s) k0 in the BZ. When this happens, the corre-
sponding condensate density xµ0(k0) can be put non-
zero, if this is necessary to satisfy Eq. (22). It is im-
portant to note that in order to determine the actual
value of ξmin(α, β) (as well as those of Q, α and β) it
suffices to only consider Eq. (22) at xµ(k) = 0, irrespec-
tive of whether there is condensate, ωµ0(k0) = 0, or not,
ωµ(k) 6= 0 for all k, µ.
We solve the Eq. (22) for ξ numerically in two steps.
First, we determine ξmin(α, β): we decrease ξ from large
positive values until the condition ωµ0(k0) = 0 signals
that ξ = ξmin(α, β). Second, we set xµ(k) ≡ 0 and at-
tempt to satisfy Eq. (22) in the interval ξ ≥ ξmin(α, β).
To this end, we set Λ = ξQ in Eq. (14) and vary ξ to de-
termine the extremum of EMF (i.e., Eq. (22)). We find
that the extremum is a maximum. If this maximum oc-
curs for some ξ > ξmin(α, β), then Eq. (22) is satisfied
with xµ(k) = 0. If, however, EMF(α, β, ξQ,Q) decreases
monotonously as we lower ξ down to ξ = ξmin(α, β) ,
then the Eq. (22) cannot be solved with xµ(k) = 0. In
this case, a finite condensate density xµ0 (k0) 6= 0, is re-
quired, in order to “compensate” for too large a value of
the lhs. of Eq. (22). This fixes both ξ = ξmin(α, β) and
the value xµ0 (k0) (cf. sects III B and IV B of Ref. 12).
ii. Having determined the value of ξ, we notice
that the function EMF(α, β,Λ, Q) is quadratic in Q and
bounded from below, which allows an analytical deter-
mination of Q as the position of the minimum.
iii. Finally, knowing the values of Λ and Q, we pro-
ceed by a numerical extremalization of EMF with respect
to the angles. The calculations show that EMF as a func-
tion of the angle β possesses a maximum, and a minimum
as a function of the angle α after β has been eliminated.
Thus, the variational energy EMF is bounded from below
in the variables Q1 and Q2, as expected.
iv. We iterate this procedure (i)-(iii) until conver-
gence is achieved.
C. Numerical results of the Sp(N) formalism
The results of the Sp(N ) approach in the large-N limit
are summarised in the zero temperature phase diagram
of the AKAF, Fig. 9. The central part of the phase dia-
gram is occupied by the incommensurate (IC) phase with
LRO at sufficiently small 1/κ. The phase boundary that
separates the region with SRO from the region with LRO
was found by checking whether for a given pair of J and
1/κ the lowest branch of the one spinon spectrum ωµ(k)
has zeros in the BZ or not, i.e., whether there will be
condensate at one or several points in the Brillouin zone
or not. As one might expect, LRO is maximally sup-
pressed by quantum fluctuations for J = 1, which is the
case of maximal frustration.
For J = 0, the exact quantum ground state of the
AKAF is ferrimagnetic (FM) according to the Lieb-
Mattis theorem15. In this state, the expectation value
Q1 which measures the singlet weight on the horizontal
bonds vanishes. As shown in Fig. 10, our Sp(N ) calcu-
lations recover this exact result and extend it to a finite
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
IC/SRO
IC/LRO
DC/SRO
FIG. 9: (Color online) Phase diagram of HAKAF as obtained
in the Sp(N ) approach. Symbols and lines, respectively, de-
note numerical and analytical results for the phase boundaries
(see text, Subsec. III B 2 and Appendix B). Quantum fluctu-
ations increase along the vertical axis. LRO: Long Range
Order; SRO: short range order; FM: ferrimagnet; IC: incom-
mensurate phase; DC: decoupled chains. At J = 1, the results
of Ref. 12 are recovered. Incommensurate order (see Fig. 12)
occurs between the boundaries of the ferrimagnetic phase (×)
and of the decoupled chain phase (∗).
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 10: (Color online) Mean field parameter Q1 as function
of the anisotropy.
interval 0 ≤ J ≤ Jferri(κ), which narrows as 1/κ in-
creases. The parameter Q2, which measures the singlet
weight on the diagonal bonds, is independent of J in this
interval; its value decreases as 1/κ increases (see Fig. 11).
Remarkably, the FM state retains its LRO in its entire
region of existence.
As J is increased beyond Jferri(κ), Q1 increases in the
manner of an order parameter at a second order phase
transition. At the same time, the parameter Q2 begins
to decrease, and eventuallly it drops to zero at some
J = JDC(κ). Thus, the large-N approach predicts the
existence of a decoupled-chain phase in the region above
the phase boundary JDC(κ). Q2 decreases to zero con-
tinuously so that the phase transition at JDC(κ) appears
to be of second order again.
Both LRO and SRO phases may be characterised by an
ordering wave vector qord = 2kmin, where kmin is that
wave vector at which the one-spinon excitation spectrum
ωµ(k) has its minimum. The static spin structure factor
S(q) develops a peak at qord. In Fig. 12, we display
the x-component of the ordering vector qxord = q
ord(J)
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 11: (Color online) Mean field parameter Q2 as function
of the anisotropy.
1/κ=1
1/κ=1.4
1/κ=3
1/κ=5
1/κ=7
1/κ=10
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 12: (Color online) Ordering wave vector qxord as function
of the anisotropy.
ord = 0). At the kagomé point J = 1, |qxord| = 4π/3 is
independent of the value of κ. For 1/κ <∼ 3, the behaviour
of qxord as a function of J is as expected: as J increases, it
increases monotonously until the phase boundary JDC(κ)
is reached and remains constant inside the DC phase.
However, for 1/κ >∼ 3 the function qordx (J) develops a
minimum at J ≈ 1.5, which becomes more pronounced
as 1/κ increases.
In Sec. III B 2 we emphasised that contrary to previ-
ous applications of the large-N approach to spin models
on kagomé and anisotropic triangular lattices11,12,19, we
found it essential to consider two chemical potentials λa
and λc here, one for the spins on the horizontal lattice
lines (λa) and one for the middle spins (λc). We display
the values of these parameters as functions of J in Fig. 13.
We have no physical explanation for the behaviour of λa,
λc as functions of J and κ but it is gratifying to see that
λa = λc at J = 1 independent of κ in accordance with
earlier work12.
As indicated above, along with numerical study of Eq. 11,
we performed extensive analytical calculations, both to
1/κ=1
1/κ=1.4
1/κ=5
 0  0.2  0.4  0.6  0.8  1
PSfrag replacements
J/(J + 1)
FIG. 13: (Color online) Lagrange multipliers λa, λc (chemical
potentials) as functions of the anisotropy.
corroborate the numerics and to obtain additional insight
into the problem. Details of the analytical techniques are
presented in Appendix B. Here we state that we were
able to analytically determine Sp(N ) phase boundaries
between the SRO and LRO DC phase, between the DC
and IC phase, and between the FM and IC phase, see
Fig. 9. Moreover, our analytical calculations allowed us
to explicitly confirm the existence of LRO inside the FM
phase and immediately to the right of the FM-IC phase
boundary. Likewise, the regions with SRO and LRO in-
side and immediately to the left of IC-DC phase bound-
ary were determined analytically. This was achieved by
evaluating in these regions the Sp(N ) generalisation of
the spin-spin correlation function 〈Si,u ·Sj,v〉 of the model
defined by expression (13). (u, v = a, b, c denote the sites
of the triangular cells i and j of the model, see Fig. 8).
On the right-hand side of the FM-IC boundary and in-
side the FM phase, we find for large distances between
the cells, |rj − ri| ≫ 1,
〈Si,u · Sj,v〉 ∼ Su Sv , (23a)
where
|x3(kmin)|2
N▽ (λc + λa)
w = a, b
−λa w = c
(23b)
and u, v = a, b, c denote the sites of the triangu-
lar cells i and j of the model, see Fig. 8. Here,
|x3(kmin)|2/N▽ is the condensate density at kmin =
(−π, 0), |x3(kmin)|2/N▽ = κ, see Eq. (B8). On the
FM-IC transition line and inside the FM phase, where
Eqs. (23) are valid, the parameters λa and λc are not
independent but can be expressed in terms of a single
parameter δ, see Eqs. (B4), (B12). The sign pattern on
the right-hand side of Eq. (23) and the ordering wave
vector qord = 2kmin = (−2π, 0) are indeed the proper-
ties one expects to find for the long-distance behaviour
of the spin-spin correlation function of a ferrimagneti-
cally ordered state. Since |x3(kmin)|2/N▽ remains finite
for arbitrarily small values of κ, the mean-field Sp(N )
approach predicts that this order persists in the ex-
treme quantum limit of our model, 1/κ ≫ 1. Together
with Eqs. (23), the fact that the condensate density
|x3(kmin)|2/N▽ remains constant inside the FM region,
see Eq. (B8), implies that the magnetisation of the FM
phase remains constant up to the FM-IC phase boundary.
The same behaviour of the magnetisation of a ferrimag-
netic phase has previously been observed in an exact-
diagonalisation study of a one-dimensional kagomé-like
antiferromagnet23. At the FM-IC phase boundary, the
magnetisation becomes spatially modulated with an in-
commensurate wave vector qord = 2kmin.
On the left-hand side of the IC-DC boundary and in-
side the DC phase we find the following large distance
behaviour of the spin-spin correlation function:
〈Si,c · Sj,c〉∼
cos [2kmin (ri − rj)]
1 + q21
|x3(kmin)|2 + |x3(−kmin)|2
3 (kmin)
,(24a)
〈Si,u · Sj,v〉∼ 0 for u, v 6= c, c . (24b)
Here, q1 and λa denote the saddle-point values of these
variables obtained from Eqs. (B31), (B32). q2 is a func-
tion of q1, determined by Eq. (B36) or by Eq. (B44)
depending on whether 1/κ < 1/κs or 1/κ > 1/κs
(κs = 0.181, see Fig. 22). ω
3 (kmin) is the value of the
second-order expansion coefficient of the lowest spinon
frequency ω3(k), cf. Eqs. (B20), (B21c), at its min-
imum, and 2kmin is the ordering wave vector imme-
diately to the left on the IC-DC phase boundary and
inside the DC phase; it is determined by Eq. (B35).
|x3(kmin)|2/N▽ = |x3(−kmin)|2/N▽ are the condensate
densities at the wave vectors ±kmin. As is shown in Ap-
pendix B, ω
3 (kmin) remains finite for 1/κ > 1/κs and
hence |x3(±kmin)|2/N▽ vanishes. Thus, 〈Si,c ·Sj,c〉 ∼ 0,
i.e., there is no LRO in this region. By contrast, for
1/κ < 1/κs both, |x3(±kmin)|2/N▽ and ω̄(2)3 (kmin) van-
ish when the IC-DC phase boundary is approached from
the left. However their ratio, which determines the spin-
spin-correlation function, Eqs. (24), remains finite in this
limit according to Eq. (B42). Thus, for 1/κ < 1/κs,
Eqs. (24) show that while the chain spins Sa, Sb remain
disordered, there is long-range IC order between the mid-
dle spins Sc along the IC-DC phase boundary and inside
the DC phase for sufficiently large κ. The middle spins
occupy the sites of a triangular lattice. Remarkably, the
correlations between these spins predicted by Eqs. (24)
are compatible with the spin pattern
Sj,c = S [cos(2kminrj)êx + sin(2kminrj)êy] (25)
that would obtain if the middle spins Sc were classical
spins coupled by a classical Heisenberg model with ex-
change constant J̃ along one lattice direction and cou-
plings J̃ ′ along the other two directions with a ratio J̃ ′/J̃
such that incommensurate order with wave vector 2kmin
would be established. This persistence of long range or-
der in the DC phase of the AKAF distinguishes our result
from the result obtained by Chung et al., Ref. 19, in their
large-N Sp(N ) treatment of the anisotropic triangular
antiferromagnet: there the DC phase consists of uncor-
related linear spin chains. Qualitative considerations of
the finite-N corrections to the mean-field Sp(N ) result
led the authors of Ref. 19 to the conclusion that instead
of the DC phase there is spin-Peierls order in the large–J
region of their model. In the next section, we will present
a different approach, a block-spin perturbation theory, to
get further insight into the properties of the AKAF for
the physical spin-1/2 case.
FIG. 14: (Color online) The kagomé lattice as a triangular
lattice of downward pointing triangles. The coupling strength
is J on the horizontal bond and unity on the other two bonds.
IV. BLOCK-SPIN PERTURBATION
APPROACH
The basic idea of the block-spin perturbation theory is
to calculate the states of small clusters of a given lattice
exactly and to treat the coupling between these clusters
perturbatively. The basic building blocks of the kagomé
lattice are triangles. Thus it is natural to consider the
trimerised kagomé lattice in which the spins on the down-
ward pointing triangles are assumed to be strongly cou-
pled whereas the coupling on the bonds of the upward
pointing triangles are assumed to be weak, see Fig. 14.
(Clearly, the exchange of the roles of the upward and
the downward pointing triangles will not affect the fur-
ther development to be presented in the current section.)
The Hamiltonian for this trimerised model reads
H(J, γ) = H▽(J) + γH△(J) , 0 ≤ γ ≤ 1 , (26)
where H▽(J) (H△(J)) denote those terms in the
Hamiltonian (Eq. 1) that act on the bonds of the down-
ward (upward) pointing triangles. We will determine
approximate GSs of this trimerised model in different
ranges of J in a perturbation expansion w.r.t. γ. The
hope is that the results will provide some qualitative
insight into the GS properties of the non-trimerised
model H(J, 1) which is our original model Eq. (1). The
same strategy has previously been applied sucessfully to
frustrated spin models by several authors5,24–26.
Obviously, the GSs of the unperturbed Hamiltonian
H(J, 0) are products of GSs of the individual downward
pointing triangular plaquettes. The GSs of a single
plaquette and the corresponding energies are
1. for J < 1:
|α〉 = 1√
|↑↑↓〉 − |↓↑↑〉
|↑↓↑〉 − |↓↑↑〉
, (27a)
PSfrag replacements
|α(ᾱ)〉:
∣β(β̄)
↑ (↓)↑ (↓)
↑ (↓)
FIG. 15: Ground-states of triangular plaquettes. Heavy lines
depict singlets. The coupling strength is J on the horizontal
bond and unity on the other two bonds.
|ᾱ〉 = 1√
|↓↓↑〉 − |↑↓↓〉
|↓↑↓〉 − |↑↓↓〉
, (27b)
εα = εᾱ = −1 + J/4 ; (27c)
2. J > 1:
|β〉 = 1√
|↑↓↑〉 − |↑↑↓〉
, (28a)
|↓↑↓〉 − |↓↓↑〉
, (28b)
εβ = εβ̄ = −3/4J . (28c)
Here, the ket vectors denote the spin state of the
plaquette in the Sz basis. The three arrows inside the
|cba〉 symbol denote from left to right the spin direction
at the sites c, b and a of the plaquettes in Fig. 15. The
states |α〉 (|ᾱ〉) and |β〉 (
) have total z-spin 1/2
(−1/2). They can be depicted graphically as shown in
Fig. 15. From these plaquette states, the zeroth order
GSs of the Hamiltonian H(J, γ) will be constructed. We
treat the cases J < 1 and J > 1 separately.
1. J < 1: Since the states |α〉, |ᾱ〉 are the GSs of the
individual downward pointing plaquettes in this case, the
states
|A(M)〉 =
iǫ{M}
jǫ{N▽−M}
|ᾱj〉 , (29)
are here the zeroth order GSs of H(J, γ). The set {M}
is a subset of M out of the N▽ downward pointing tri-
angles of the 3N▽-site kagomé lattice; the subscripts i, j
denote the position of individual triangular plaquettes in
the lattice of these plaquettes which is also triangular, see
Fig. 14. The zeroth order energy eigenvalues associated
with the states |A(M)〉 do not depend on M :
= N▽(−1 + J/4) . (30)
Hence, there are in total 2N▽ degenerate zeroth order
GSs |A(M)〉. The single plaquette states |α〉, |ᾱ〉 satisfy
the conditions for the validity of the Lieb-Mattis theorem,
Ref. 15: after a canonical transformation which rotates
the spins on the sites a and b by π around the z-axis
|↑〉 → i|↑〉, |↓〉 → −i|↓〉, and which leaves the spins on the
site c fixed the coefficients of all basis states on the right
sides of Eqs. (27a, 27b) become positive (+1/
6). As a
consequence, all the GSs |A(M)〉 satisfy the conditions
for the validity of the Lieb-Mattis theorem. For J = 0
it follows from this theorem that the total magnetisation
of the exact quantum GS |Φexact〉 of the Hamiltonian
HAKAF must be an eigenstate of the total magnetisation
m̂tot =
(Szi,a + S
i,b + S
i,c) (31)
with eigenvalue mtot = N▽/2, i.e., |Φexact〉 must be a
ferrimagnetic state. By continuity, one expects this to
be the case not only for J = 0, but up to a certain fi-
nite value of J . This suggests that the state |A(M=0)〉,
c.f. Eq (29), is the appropriate zeroth order GS in this
case and that the degeneracy of the states |A(M)〉 is lifted
by the perturbation H△ in favour of the state |A(0)〉.
To confirm this, we determine the creation energy of a
flipped plaquette in first order in γ, i.e., the difference
of the energy of the state with one plaquette spin flipped
relative to the ferrimagnetic state, and the energy of the
ferrimagnetic state:
δ(1)EA(M=1) = EA(1) − EA(0). (32)
A simple calculation yields
δ(1)EA(M=1) =
γ(1− J) , (33)
i.e., to first order, |A(M=0)〉, the ferrimagnetic GS is
stable w.r.t. a flip of a single plaquette spin, as long as
J < 1.
As a further check on the stability of the state
|A(M=0)〉, we calculate the dispersion of the excitation
energy of a propagating single flipped plaquette spin. For
this purpose, we need to determine the overlap matrix el-
ements between the state with a flipped plaquette spin
at the site j and states with a flipped spin at one of the
neighbouring sites,
tj,j±δ1 = 〈ᾱj |〈αj±δ1 |γJSjSj±δ1 |αj〉|ᾱj±δ1〉 =
(34a)
tj,j±δ2,3 = 〈ᾱj |
αj±δ2,3
∣γSjSj±δ2,3 |αj〉
∣ᾱj±δ2,3
(34b)
Here, δν , ν = 1, 2, 3, are the primitive lattice vectors
of the kagomé net, see Fig. 2; they connect the sites of
the plaquette lattice. Then, by diagonalising the ensuing
transfer Hamiltonian
Htrans = γ
J ( |j + δ1〉〈j|+ |j − δ1〉〈j| )
(|j + δ2〉〈j|+ |j − δ2〉〈j|)
(|j + δ3〉〈j|+ |j − δ3〉〈j|)
,(35)
where |j〉 denotes the state with a flipped plaquette spin
at site j, we obtain for the kinetic energy of this excita-
tion:
ε(k) =
J cos(kx)− cos(
) cos(
. (36)
Adding the energy for the creation of a single flipped pla-
quette spin, we find for the total energy of the excitation
in the limit of small wave vector k
ω(k) =
− J)k2x +
k2y +O(k4)
. (37)
Obviously, the ferrimagnetic state |A(M=0)〉 becomes
unstable against a propagating flipped plaquette spin
already at J = 1/4, i.e., much earlier than suggested
by the excitation energy of a static flipped spin (see
Eq. (33)). We remark that this bound is independent of
the actual magnitude of the perturbation parameter γ
and therefore, the qualitative result may survive in the
limit γ = 1.
2. J > 1: In this region, the states
|B(M)〉 =
iǫ{M}
jǫ{N▽−M}
with eigenenergy
= N▽(−3J/4) . (39)
are the zeroth order eigenstates of H(J, γ). These states
consist of free spins on the c-sites and of spin-singlet
dimers that cover every second bond of the horizontal
chains of the lattice. We wish to answer the question of
whether the 2N▽-fold degeneracy of these states, which
results from the degrees of freedom of the free spins, is
lifted by the perturbation γH△; in other words, we want
to find out whether the middle spins remain decoupled
from the chain spins. We proceed as in case (i). We
compare in a perturbation expansion w.r.t. γ the en-
ergy of the state |B(0)〉 with the energy of |B(1)〉, i.e.
with the state with one plaquette spin flipped relative to
|B(0)〉. We denote this difference by δ(1)EB(M = 1) =
EB(1) − EB(0). Surprisingly, we find that the matrix el-
ements 〈B(M)|H△|B(M)〉 vanish for any choice of M .
There is no first order correction to the energy E
δ(1)EB(M = 1) = 0. Moreover, we observe that the
off-diagonal matrix elements 〈B′(M)|H△|B(M)〉, where
|B′(M)〉 and |B(M)〉 contain identical numbers of states
but differ in their distribution over the N▽ down-
ward pointing triangles, also vanish. This implies that, in
contrast to case (i), a flipped plaquette spin cannot hop
to a neighbouring site in a first order process. Coupling
between the spins on the c-sites occurs only in second
order in γ. It is succinctly described by an effective spin
Hamiltonian for the c-site spins, which are at the same
time total spins of the downward pointing plaquettes (see
Fig. 15):
Heff =
Szi S
Sxi S
.(40)
Here, Sαi , α = x, y, z, denote plaquette spin operators;
i is the position of a downward pointing plaquette on
the triangular lattice formed by these plaquettes. The
exchange couplings J
and J⊥
are given as second order
matrix elements of H△:
〈Bi↑, i′↑|H△|X〉 〈X |H△|Bi↑, i′↑〉
2 εB − εX
〈Bi↑, i′↓|H△|Y 〉 〈Y |H△|Bi↑, i′↓〉
2 εB − εY
(41a)
〈Bi↓, i′↑|H△|X〉 〈X |H△|Bi↑, i′↓〉
2 εB − εX
, (41b)
and i′ ≡ i+ δν . Here, the states |Biσ, i′σ′〉 are zeroth
order GSs, Eq. (38), whose spin patterns are identical
on all sites except for the sites i and i′ where the z-
components of the spins take the values σ and σ′, re-
spectively; |X〉 and |Y 〉 are excited states of H▽. Of
course, since the SU(2) symmetry of the original Hamil-
tonian H(J, γ) must be conserved in the derivation of
Heff , the expressions Eqs. (41) must yield identical re-
sults, J
≡ Jδν . Non-zero contributions to J
and J⊥
are obtained if either the same term SiSi′ of
H△ acts in both matrix elements of the numerators of
Eq. (41) (two-block contributions) or the terms SiSk,
SkSi′ act in the left and right elements, respectively,
where the plaquette geometry must be as shown in Fig. 16
(three-block contributions). In contrast to the case of the
isotropic KAF studied by Zhitomirsky26, the three-block
PSfrag replacements
FIG. 16: Configurations of ▽ blocks contributing to the in-
terblock couplings J
and J⊥δν . Double dashed lines indicate
that the same term element of H△ acts twice between the ▽
blocks at sites i and i′ (see also text).
contributions do not produce three-spin interactions in
the present case. Rather, they contribute to the ex-
change interactions J
and J⊥
of the HamiltonianHeff ,
Eq. (40). The evaluation of the expressions (41) yields
Jδ1= γ
1 + 1
1 +O(J−1)
Jδ2 = Jδ3 =
1 + 1
1 + O(J−1)
. (43)
Obviously, these results are useful for J ≫ 1. There,
Heff represents a spin 1/2 Heisenberg Hamiltonian on
the triangular lattice of the c-sites with a coupling
along the δ1 direction that is strong (O(γ2/J)) in com-
parison to the couplings in the two other directions
(O(γ2/J2)). This limiting case of the anisotropic trian-
gular Heisenberg antiferromagnet (ATHAF) has recently
been analysed by Starykh and Balents with field theo-
retical methods14. These authors find that in the limit
of strong anisotropy, K ≡ Jδ1/Jδ2 → ∞, the GS of
the model Eq. (40) is a collinearly ordered antiferromag-
net (CAF) as depicted in Fig. 17. Since the ordering
wave vector qCAF = (π, π/2) of this phase does not
evolve continuously from the ordering wave vector qIC
of the incommensurate (IC) spiral phase of the ATHAF,
FIG. 17: (Color online) Collinear antiferromagnetic state
(CAF) on the triangular lattice14.
FIG. 18: (Color online) Tentative ground-state of the
anisotropic kagomé antiferromagnet in the limit J ≫ 1. Dou-
ble lines: dimers between the spins on the end points.
(qIC = (qx(K), 0) with −3π/2 ≤ qx(K) ≤ −π for
1/2 ≤ K ≤ ∞), they conclude that the IC phase and
the CAF phase must be separated by a quantum phase
transition. For the trimerised anisotropic kagomé model,
Eq. 26, these results have the following implications:
i) While in the limit of strong anisotropy J ≫ 1 there is
long-range collinear antiferromagnetic order among the
c-site spins, the a- and b-site spins are paired in singlets,
see Fig. 18.
ii) This picture of the GS of the trimerised anisotropic
kagomé model Eq. (26) differs from the result obtained in
the Sp(N ) approach insofar as for sufficiently large κ the
Sp(N ) approach predicts long range IC order among the
c-site spins up to arbitrarily large values of J . Thus, if
the picture of a CAF phase for large anisotropy persists in
the non-trimerised limit H(J, γ = 1) of the model Eq. 26,
one would expect a quantum phase transition between
the IC phase and the CAF phase of the AKAF similarly
as for the ATHAF. In closing this section, we remark
that the calculation that led to the effective Hamiltonian
Heff , Eq. 40, i.e., to the coupling between the c spins
in the strongly anisotropic limit, shows clearly that this
coupling arises from quantum fluctuations of the a and b
spins.
V. SUMMARY AND DISCUSSION
In this work, we have studied the ground state (GS) phase
diagram of the quantum Heisenberg antiferromagnet on
the kagomé lattice with spatially anisotropic exchange
(AKAF). The model is relevant for a description of mag-
netic properties of volborthite, which is a natural reali-
sation of a spin 1/2 antiferromagnet consisting of weakly
coupled slightly distorted kagomé layers. A small mono-
clinic distortion along one of the three lattice directions
causes the exchange coupling along this direction, J , to
differ from the couplings in the other two directions, J ′,
which we set equal to unity, cf. Fig. 2. We have inves-
tigated the problem in the full range of the anisotropy,
0 ≤ J ≤ ∞, using three different approximate methods:
the classical and semiclassical approach, the mean-field
Sp(N ) approach, and a block-spin perturbation theory.
The case J = 1 is the much studied isotropic kagomé
antiferromagnet (KAF). Exact diagonalisation studies of
this model3,4 are available. Their results speak conclu-
sively in favour of a spin liquid ground state1. This view
is supported by block-spin approaches5,6. Conflicting re-
sults have been found in Refs. 27–31, where various va-
lence bond crystal (VBC) states are proposed as ground
states of the KAF. However, a recent comparison of the
exact spectrum of the 36-site sample of the KAF against
the excitation spectra allowed by the symmetries of these
states, casts doubts on their validity32.
Within the whole anisotropy range, the case J = 0 is spe-
cial, since it allows for an exact characterisation of the
quantum GS as ferrimagnetic (FM) with a total mag-
netisation of M = S Ns/3 for a system of Ns spins of
magnitude S. In the classical picture, this state corre-
sponds to a unique staggered layout of spins with a non-
zero net magnetisation of the lattice unit cell (cf. Fig. 3).
In the classical limit, the ferrimagnetic ground state sur-
vives up to J = 1/2. For J > 1/2, the “chain” spins
(i.e., spins coupled by J) begin to tilt gradually towards
the middle (remaining) spins (see Fig. 5). This allows for
a formation of a large degenerate manifold of canted spin
states. In contrast to the isotropic case J = 1, where
the degeneracy grows exponentially with the system size
Ns, its growth is weaker: 2
Ns for J 6= 1. This im-
plies that there must be an increasingly large number of
classical low energy configurations as J approaches unity.
In the linear semiclassical approximation, the spin-wave
spectrum has one zero-frequency mode for each point
of the magnetic Brillouin zone (BZ). The spectrum is
identical for the different canted states for all J > 1/2.
Thus, in this order of the semiclassical approximation,
no order-by-disorder mechanism appears that would se-
lect one particular state or a particular group of states
from the manifold of canted states as true ground states.
In the limit J → ∞, the frequency spectrum of non-zero
modes gradually takes the shape of the spectrum that
one would expect for a set of uncoupled antiferromag-
netic spin chains parallel to the strong-J direction. No
qualitative change from the set of canted spin states to
the set of decoupled chains at a finite value of J is found.
We have further explored the nature of the phases at var-
ious J exploiting the mean field (MF) Sp(N ) approach,
that incorporates the effect of quantum fluctuations not
only perturbatively, but self-consistently. The strength
of quantum fluctuations is controlled by a parameter κ,
which is the analogue of the spin value S in the origi-
nal SU(2) symmetric model. In fact, for N = 1, when
the Sp(1) symmetric model is equivalent to the SU(2)
model, κ = 2S. For general N , this last identity does
not hold, but κ is still a measure for the importance of
quantum fluctuations that are strong for κ ≪ 1 and weak
for κ ≫ 1. In the MF Sp(N ) approach, the nature of the
phases that occur can be read from the values of the
mean field parameters Q1 and Q2 and from the spec-
trum of the bosonic spinon excitations. While the mean
field parameters Q1 and Q2 (cf. Fig. 8) are the GS ex-
pectation values of singlet bond operators, the structure
of the spinon spectrum, ωµ(k;Q, λ), determines the ex-
istence or non existence of long-range order (LRO): If
the spectrum becomes gapless at some wavevector qord,
a Bose condensate will form and a modulated structure
with the wavevector 2qord will acquire LRO.
As was to be expected, the phase diagram of the AKAF
obtained by the MF Sp(N ) approach contains an incom-
mensurate (IC) phase in the vicinity of the isotropic point
J = 1 which is ordered for sufficiently large κ according
to this approach, see Fig. 9. Qualitatively, we may gauge
the value of κ against the spin length S by looking at
the line J = 1 of the phase diagram which is the location
of the Sp(N ) analogue of the isotropic SU(2) symmet-
ric kagomé model: since, as we have argued above, the
SU(2) model is disordered for S = 1/2, we may con-
clude from Fig. 9 that the value of 1/κ that corresponds
to S = 1/2 must be greater than two. Somewhat sur-
prisingly, the FM phase remains long-range ordered for
arbitrarily small κ. This may reflect the fact that in
the SU(2) version of the model, the FM phase is ordered
even for the smallest physical spin value S = 1/2. A new
feature of the phase diagram is the prediction of a decou-
pled chain (DC) phase for large enough J , which has no
classical analogue. In this phase, the chains of strongly
coupled a- and b-site spins show no magnetic order. The
c-site spins which are interspersed between these chains
and which occupy the sites of triangular sublattice are
decoupled from the chain spins. However, they may or
may not exhibit long range order among themselves de-
pending on the magnitude of κ. Remarkably, the spin-
spin correlations, whose asymptotics were obtained ana-
lytically, are compatible with the spin-spin correlations
of an anisotropic classical Heisenberg antiferromagnet on
the triangular lattice whose exchange couplings differ in
one direction from those in the other two directions.
In order to tackle the problem of the GSs of the AKAF
from a third corner, we have used a block-spin pertur-
bation theory. This method has the advantage of being
applicable directly to the spin 1/2 version of the model.
In applying this approach, one has to initially group the
spins of the model in clusters. For the kagomé lattice, it
is natural to choose the spins around either the upward
or the downward pointing triangles as clusters of strongly
coupled units and to consider the coupling between these
clusters, γ, as the small expansion parameter. Thus one
trimerises the original model (see Fig. 14) and in so do-
ing, one breaks the translational invariance of the original
model. In the zeroth order of this expansion, two regions
can be distinguished by the eigenenergies of the individ-
ual trimers: J < 1 and J > 1. For sufficiently small J ,
one recovers the FM state as the GS in first order w.r.t. γ.
For J > 1, there are no first order corrections to the en-
ergy. Following an earlier application of the block-spin
technique to the isotropic KAF26, we determine for J > 1
in second order in γ an effective HamiltonianHeff for the
block-spins which can be identified as the middle spins of
the original model and that occupy the sites of a trian-
gular lattice. Heff is a Heisenberg Hamiltonian with a
coupling Jδ1 of the order of γ
2/J along the δ1 direction
(cf. Fig. 2) and couplings Jδ2 = Jδ3 of the order of γ
along the other two directions. The calculations that lead
to these results show clearly that the couplings between
the c-spins of the AKAF are due to fluctuations of the sin-
glets between the a- and b-spins into excited states. In a
very recent field theoretical study, Starykh and Balents14
arrive at the conclusion that for Jδ1 ≫ Jδ2,δ3 , the ground
state of the anisotropic triangular antiferromagnet rep-
resented by Heff is a collinearly ordered spin state, see
Fig. 17. Then, together with the singlet dimers between
the a- and b-spins of the downward pointing triangles, the
state depicted in Fig. 18 emerges as the candidate ground
state of the AKAF in the limit of large anisotropy J ≫ 1:
while nearest neighbour spins on the strongly coupled a-b
chains form singlets and decouple magnetically from the
spins on the c sites, the latter order in a collinear antifer-
romagnetic structure. This structure cannot be obtained
by a continuous deformation of the spiral IC structure
that is predicted by the Sp(N ) approach and is believed
to prevail for sufficiently large κ in the region of moder-
ate anisotropy. As a consequence of the trimerisation, the
state depicted in Fig. 18 breaks the translational symme-
try of our original model, Eq. (1). If this state survives as
the ground state of the non-trimerised model, i.e., when
the expansion parameter γ approaches unity, then, ow-
ing to their incompatible symmetries, the spiral IC phase
and the large J phase of our model must be separated by
a quantum phase transition.
Acknowledgments
One of the authors (HUE) acknowledges a useful discus-
sion with F. Mila. The work at the University of Water-
loo was supported by the Canada Research Chair (Tier
I, Michel Gingras). We thank M. Gingras for a critical
reading of the manuscript and numerous helpful sugges-
tions.
APPENDIX A: GROUND STATE DEGENERACY
FOR GENERAL J
We first derive the constraint on the chiralities that leads
to the reduction in the number of degenerate ground
states for general J relative to the special case J = 1.
Let χ1, · · · , χ6 be the chiralities of the six triangles sur-
rounding one of the hexagons of the kagomé lattice, and
let φ1, · · · , φ6 denote the angles that define the direc-
tions of the spin vectors on the six corners of the hexagon,
see Fig. 19.
Then, as is seen in Fig. 19, the following relations between
the angle φ1, and the angles φ2 · · · , φ6 are an immediate
consequence of these definitions:
φ2=φ1 − θ χ2 , (A1a)
φ3=φ1 − θ (χ2 + χ3) , (A1b)
φ4=φ1 − θ (χ2 + χ3)− (2π−2θ) χ4 , (A1c)
φ5=φ1 − θ (χ2 + χ3 + χ5)− (2π−2θ) χ4 , (A1d)
φ6=φ1 − θ (χ2 + χ3 + χ5 + χ6)− (2π−2θ) χ4, (A1e)
and φ6 = φ1 + (2π − 2θ)χ1 . (A1f)
From the last two of these relations it follows that the
chiralities χ1, · · · , χ6 are constrained by the sum rule
χ2 + χ3 + χ5 + χ6 − 2χ1 − 2χ4 = 0 . (A2)
For the isotropic kagomé system, J = 1, θ = 2π/3, one
finds instead of the constraint (A2) the sum rule
χj = n where n = 0, 1, 2 (A3)
which is obviously less restrictive than (A2).
Next, we present the arguments that lead to the estimate
NanisoGS (N▽)
N▽ with α < 3 (A4)
for the number NanisoGS (N▽) of classical GSs of an
anisotropic kagomé AF with N▽ downward pointing tri-
angles (the number of sites is 3N▽). Any planar config-
uration of a cell of the kagomé lattice can be constructed
by decorating the successive rows of up and down point-
ing triangles with chirality values χ = ±1 starting with
the first row. We consider only square cells with
rows with
N▽ downward pointing triangles. Then,
each row consists of 2
N▽ triangles, see Fig. 20.
Obviously, there are 22
N▽ ways to decorate the first
row. Disregarding certain exceptions, which will be dis-
cussed below, one can, for a given configuration of the
first row, choose the chirality of an arbitrary triangle
of the second row to be either +1 or −1. After this
choice has been made, the constraint (A2) fixes the chi-
ralities of all the remaining triangles of the second row
PSfrag replacements
FIG. 19: Chiralities around heaxagonal plaquette
FIG. 20: Example of a chirality distribution; dark and light
shaded triangles represent positive and negative chirality, re-
spectively. Chirality configurations in boxes fix the chirality
distribution of the row above them uniquely. An empty cir-
cle inside a triangle indicates that its chirality can be chosen
freely to be positive or negative.
uniquely. Proceeding in this manner from row to row
one would generate 22
N▽ · 2
N▽ distributions of chi-
ralities over the N▽ downward pointing triangles of the
cell. For finite lattice cells, the requirement of periodic
boundary conditions imposes further constraints on the
number of possible chirality distributions in these cells,
but the effect of these constraints will become negligible
in the thermodynamic limit N▽ → ∞. However, there
is a further reduction of the number of possible chirality
distributions: For a given distribution in a row it is not
always possible to find two distributions for the successive
row which both satisfy the constraint (A2). If in a row
the lower half of a hexagon of the next row is decorated
by chiralities in the manner − + − or + − + (see boxes
in Fig. 20), then the chiralities of the next row are fixed
uniquely. This reduces the number of possible chirality
distributions. Obviously, this reduction of the number
of possible chirality distributions survives in the thermo-
dynamic limit so that the exponent in (A4) is less than
N▽, the value one would have expected without this
reduction. We have calculated the number of distribu-
tions for cells of up to N▽ = 13× 13 and have found the
value α ≃ 2.18 for the constant in the expression (A4),
see Fig. 21.
As we have mentioned above, the sum rule (A3) which
applies for the isotropic kagomé AF is less restrictive than
the sum rule (A2). Consequently, the number of chirality
distributions in the isotropic model33,
N isoGS ∼ 1.18333N▽ (A5)
is larger than in the anisotropic model. Since the tran-
sition from the anisotropic model to the isotropic model
happens through a continuous variation of the coupling
constant J , there should be a continuous transition be-
tween the numbers of GS configurations in these two
cases. Presumably, this transition implies that the den-
sity of low-energy states of the anisotropic model in-
creases exponentially with an exponent ∼
N▽ so that
for J → 1 a sufficient number of states collapses to the
GS to bring about the transition between the laws (A4)
and (A5).
 2  4  6  8  10  12  14
linear fit
minPSfrag replacements
FIG. 21: (Color online) Number of chirality distributions,
NanisoGS , of cells of up to N▽ = 13 × 13. Dotted line: min =
N▽ ln 2 (lower bound); dashed line: max = 3
N▽ ln 2
(upper bound); full line: ln(NanisoGS ) = 0.65 + 2.18
N▽ ln 2
(linear fit to the numerical results).
APPENDIX B: PHASE BOUNDARIES
The FM phase and the DC phase are chacterised by
the vanishing of the parameters Q1 and Q2, respec-
tively. Our numerical results in section III C show that
at the respective phase boundaries, Q1 and Q2 decrease
to zero like order parameters at second order phase tran-
sitions. This suggests that we expand the mean field
energy EMF, Eq. (14), w.r.t. either Q1 or Q2 in the
manner of a Landau-Ginzburg (LG) expansion and de-
termine the phase boundaries and the properties of the
FM and the DC phase from this expansion. We write
EMF /(N▽N ) = e(α)LG(Qα) where
LG(Qα) = eα + rα |Qα|
2 + gα |Qα|4 +O(|(Qα|6) . (B1)
The coefficients eα, rα and gα are functions of the vari-
ables κ and J , of the parameters λa, λc and of Q2, Q1
for α = 1, 2, respectively. The saddle point of e
LG(Qα)
w. r. t. λa, λc and Qβ, β 6= α, determines the physical
values of these parameters. For e
LG(Qα) to qualify as
a bona fide Landau-Ginzburg energy describing a second
order phase transition with Qα playing the role of an or-
der parameter, the coefficients gα have to be positive at
the saddle point. For g1, i.e., inside and on the bound-
ary of the FM phase, this follows from the numerical
result: Q1 is found to remain zero for all J ≤ JF (κ). By
contrast, we have no numerical results for J ≥ JDC(κ),
i.e., inside and on the boundary of the DC phase. There-
fore, we need to show by analytic means that g2 > 0.
1. The FM phase and the FM-IC phase boundary
Since, as we have just remarked, we know that g1 > 0,
the remaining task is to determine the coefficients e1 and
r1 of e
LG. To this end, we have to expand the mean
field energy EMF, Eq. (14), w.r.t. Q1 which amounts to
expanding the frequencies ωµ(k) w.r.t. Q1. As can be in-
ferred from the expressions (16), (18) the frequencies de-
pend on Q1 only through the combination ε
2 = J2|Q̃1|2.
Therefore, we write the expansion in the form
ωµ(k; ε) = ω
µ (k) + ε
2ω(1)µ (k) +O(ε4)
ω̄(0)µ (k) + ε̄
2ω̄(1)µ (k) +O(ε̄4)
with λ+ = (λa + λc)/2 , ω̄
µ = ω
µ /λ+ , ω̄
µ (k) =
λ+ ∂ε2 ωµ(k; ε)|ε=0 and ε̄ = ε/λ+.
Here, the introduction of the “dimensionless” quantities
µ , ε̄ looks like an unneccessary complication but it will
help to keep expressions further below simple. Setting
Q1 = 0 in the matrix D̂(ω), Eq. (16), and solving Eq. (15)
for ω we find
1 (k) = wF (k) + δ , (B3a)
2 (k) = 1− δ , (B3b)
3 (k) = wF (k) − δ . (B3c)
δ = λ−/λ+ with λ− = (λc − λa)/2 (B4)
wF (k) =
1− q̄22[sin2(s2/2) + sin2(s3/2)] (B5)
q̄2 = |Q2|/λ+ (B6)
and sa = δak , a = 2, 3 (see Fig. 8) .
From our numerical results, Fig. 13, we know that λc >
λa and hence δ > 0. Therefore, ω̄
3 (k) < ω̄
1,2(k),
and hence, if the minimum of ω̄
3 (k) vanishes at the
point kmin in the Brillouin zone, ω̄
1,2(kmin) will be fi-
nite. Thus, since condensate can only occur when one of
the frequencies ω̄
µ , µ = 1, 2, 3 vanishes there may be
a finite condensate density |x3(kmin)|2, but the densities
|x1|2 and |x2|2 will certainly be zero. With these remarks
and with the above results for ω̄
µ we find from Eq. (14)
e1/λ+ = 2λ+q̄
2 − (3− δ)(κ+ 1)
1 (k) + ω̄
2 (k) + ω̄
3 (k)
3 (kmin) |x3(kmin)|2/N▽ . (B7)
Stationarity of e1 w.r.t. λ−, λ+, and q̄
2 (which is equiv-
alent to stationarity w.r.t. λa, λc), and Q
2 requires the
following three conditions to be fulfilled:
∂e1/∂λ− = 0 :
|x3(kmin)|2 = κ ; (B8)
∂e1/∂λ+ = 0 :
2λ+q̄2
2 − 3
κ− 1 +E2(q̄2) +
wF (kmin) = 0 ;
∂e1/∂q̄
2 = 0 :
2λ+ −
[K2(q̄2)−E2(q̄2)]−
wF (kmin)
= 0 ;
(B10)
with K2(q̄2) =
ds3 wF (k)
E2(q̄2) =
ds3 wF (k) .
(B11)
According to Eq. (B8), condensate must be present in
the FM region. This requires that ω̄
3 (kmin) vanishes.
From Eq. (B3c) it is seen that kmin = (−π, 0), so that
3 (kmin) = 0, if
wF (kmin) =
1− 2q̄22 =
. (B12)
Within the FM region and on the FM-IC boundary
(i.e., for Q1 = 0) the saddle-point values of q̄2, λ+,
and λ− are then determined as functions of κ by the
Eqs. (B9), (B10) and (B12). Remarkably, within this re-
gion these quantities are independent of the value of the
exchange constant J . The solution of these equations
shows that 0 ≤ q̄2 ≤ 2/3 for 0 < κ < ∞, cf. Figs. 11, 13.
The FM-IC phase boundary is the solution of r1(κ, J) =
0 (cf. Eq. (B1), where
r1 = ∂e
LG/∂Q
(B13)
with e
LG (EMF) from Eq. (14).
We obtain
s2 + s3
Ω(1)(k)
k→kmin
s2 + s3
3 (k)
(B14)
with Ω(1)(k) = −ω̄(1)1 (k) − ω̄
2 (k)− ω̄
3 (k).
To obtain the expansion coefficients ω̄
µ (k) which appear
in the last equation, we solve Eq. (15) to first order in
the expansion w.r.t. ε̄2. We find
Ω(1)(k)=
wF (k)(wF (k) + 1− 2δ)
wF (k) + 1
4 sin2
(2wF (k) + 1− δ)
(wF (k) + 1)(1− δ) (wF (k)2 − δ2)
(B15)
k→kmin
s2 + s3
3 (k)
= − 1
1 + δ
1− δ .
(B16)
With these results Eq. (B14) can, in the thermodynamic
limit, be cast into the form
− I3(q̄2)
1 + δ
1− δ =
JF (κ)
, (B17)
where
I3(q̄2) =
ds3 2 sin
Ω(1)(k) .
(B18)
Then, with q̄2 = q̄2(κ) and λ± = λ±(κ) as obtained from
Eqs. (B9), (B10) and (B12), the condition r1 = 0 is
an equation for the FM-IC phase boundary J = JF (κ)
which yields the graph shown in Fig. 9. As we have
mentioned above, inside the FM region, i.e., for J <
JF (κ), the saddle-point values of the quantities q̄ and λ±
and hence of Q2, λa, λc and |x3(kmin)| are independent
of the exchange coupling J , i.e., they retain the values
they attain on the FM-IC phase boundary, cf. Figs. 11,
2. The DC phase and the IC-DC phase boundary
Proceeding in exact analogy to the development in the
previous subsection we now expand EMF/(N▽ N ) in
powers of |Q2|2. However, instead of working with the
variables Q1, Q2, λa, λc we work with q1, Q2, λa, q2 here,
where
J |Q1|
, (B19a)
|Q2|√
. (B19b)
The replacement of |Q1| is purely a matter of conve-
nience. By contrast, the replacement of variables Q2, λc,
which according to the numerics vanish simultaneously
as J approaches the IC-DC phase boundary, by the pair
Q2, q2 leaves us with only one vanishing variable, since,
as will be seen below, q2 remains finite throughout.
a. Expansion of e
We write
ωµ(k) = ω
µ (k) + ω
µ (k)Q
2 + ω
µ (k)Q
2 +O(Q62)
(B20)
and determine the coefficients ω
µ , n = 1, . . . 4, by solv-
ing Eq. (15) for ω iteratively. We obtain
1 (k) + ω
2 (k) = 2λa wDC(k) , ω
3 = 0 , (B21a)
1 (k) + ω
2 (k) = −
1− cos kx cos ky
wDC(k)
, (B21b)
3 (k) =
q22 λa
C(k)2 −D(k)2
. (B21c)
Here,
C(k) = 1 − q22
1 − cos kx cos ky
wDC(k)2
, (B22)
D(k) = q1 sin k
x q22
cos kx − cos ky
wDC(k)2
, (B23)
wDC(k) =
1− q21 sin2 kx . (B24)
The coefficients ω
µ (k), µ = 1, 2 , 3, will only be needed
in the determination of the coefficient g2 of the fourth
order term of e
LG(Q2) which will be discussed later. We
will first concentrate on the determination of the zeroth
order term, e2, and of the coefficient r2 of the second
order term of e
LG(Q2). Under the assumption that g2 is
positive, this will provide us with an expression for the
IC-DC phase boundary.
With the above expressions for ω
1 + ω
2 and ω
ν = 0, 1, we obtain for the coefficients of the Landau
Ginzburg energy from Eqs. (14), (B1)
e2(q1, λa) =
λ2a q
− 2λa
1 + κ− 1
wDC(k)
(B25)
r2(q1, q2, λa, |x3(kmin)|2) =
(κ + 1)
1− cos kx cos ky
wDC(k)
3 (k) + |x3(kmin)|2 ω
3 (kmin)
. (B26)
These are valid for arbitrary values of the parameters q1,
λa, q2, and |x3(kmin)|. In the next subsection, we will
calculate their saddle point values for given Q2 and thus
fix the parameters. Here, we have only allowed for the
existence of a condensate component |x3(kmin)|2. This
is justified since, as Eqs. (B20) and (B2) show, ω3 < ω1,2
for sufficiently small Q2 so that conceivably ω3(k) may
vanish at some point kmin in the the Brillouin zone, while
ω1(k) and ω2(k) remain finite at kmin, and hence a finite
condensate density |x3(kmin)|2 may occur at this point.
b. Saddle point, phase boundary
Next we need to determine the saddle point of e
LG(Q2)
in the space of the variables q1, λa, q2, and |x3(kmin)|.
First, the saddle point values of q1 and λa are obtained
as expansions in powers of Q2,
a + λ
2 +O(Q42) , (B27a)
1 + q
2 +O(Q42) , (B27b)
where λ
a , q
1 are the solutions of
∂λae2=0 , (B28a)
∂q1e2=0 . (B28b)
Since the first derivatives of e2 vanish at λa = λ
a , q1 =
1 , Eqs. (B28), we have
e2 = e
2 + e
2 +O(Q62), (B29)
r2 = r
2 + r
2 +O(Q42) (B30)
Here, e
2 and r
2 are the expressions (B25) and (B26)
with λa and q1 replaced by λ
a and q
1 . The fourth
order term of e2, Eq. (B29), and the second order term
of r2 contribute only to the fourth order term of e
which will be determined later. Therefore, we postpone
the presentation of explicit expressions for λ
a , q
1 and
the ensuing expressions for e
2 and r
2 until later. With
e2 from Eq. (B25), Eqs. (B28) yield the equations
1 )− 1, (B31)
1 ) − E(q
, (B32)
which determine the saddle point values q
1 and λ
( K and E are the elliptic integrals of the first and the
second kind.)
Next we seek the extremum of e
LG w.r.t. q2. Since e2 is
independent of q2 we, neglecting terms of order Q
2, have
0=∂q2 r
q32 λa
κ + 1 − I1(q1, q2)
C(kmin)
|x3(kmin)|2
3 (kmin)
. (B33)
I1(q1, q2) =
∫ π/2
[C(k)2 −D(k)2]1/2
(B34)
(In these expressions and in the sequel, we use an abbre-
viated notation: λa, q1 and q2 denote the zeroth order
quantities λ
a , q
1 and q
2 .) kmin is the location of the
minimum of ω
3 (k),
min = 0 ;
kxmin
,−π ≤ kxmin ≤ −
(B35)
From (B21c) and (B35) it follows that ω
3 (kmin) = 0,
q22 = (1 − q21)/2. (B36)
As a function of q2 the integral I1(q1, q2) increases
monotonously,
1 = I1(q1, 0) ≤ I1(q1, q2) ≤ I1 (q1)
for 0 ≤ q2 ≤
(1 − q21)/2 . (B37)
We have defined
I1(q1) := max
I1(q1, q2) = I1
(1 − q21)/2
(B38)
As is seen in Fig. 22, the graphs of the functions
κ = κ(q1), Eq. (B31), and of I1 = I1(q1) intersect at
q1s ≃ 0.708, κs ≃ 0.181. Therefore, in solving Eq. (B33)
for q2, two cases have to be considered separately:
i. q1 > q1s, κ > κs. In this case, a solution ex-
ists only, if the last term in parentheses in Eq. (B33) is
positive. This requires that ω
3 (kmin) = 0 because, as
has been discussed before, |x3(kmin)| and hence the ra-
tio |x3(kmin)|2/ω(2)3 (kmin) would vanish otherwise. The
 0.95
 1.05
 1.15
 1.25
 0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8  0.9  1
PSfrag replacements
q1q1s
κ(q1) + 1
I1(q1)
FIG. 22: (Color online) I1 and κ+ 1 as functions of q1.
condition ω3(kmin) = 0 implies that q
2 = (1− q21)/2, cf.
Eq. (B36). Using this result and Eq. (B32) to eliminate
q2 and λa from Eq. (B26) we find
2 = 2
1− JDC(κ)
, (B39)
where
JDC(κ)=
(κ+ 1) (3− q21)/2 + Ĩ2(q1,
(1− q21)/2)
(1 − q21)
q21 π
4[K(q1)−E(q1)]
(B40)
Ĩ2(q1, q2) =
C(k)2 −D(k)2
(B41)
is the IC-DC phase boundary for κ > κs, i.e., in the
region where the ratio |x3(kmin)|2/ω(2)3 (kmin) is finite.
According to the discussion at the end of section III C,
cf Eq. (24), this is the region where LRO prevails along
the decoupled chains, cf. Fig. 9.
In the development leading to Eq. (B40) for the phase
boundary, we have not needed the solution of Eq. (B40)
explicitly, but we note it here for completeness:
C(kmin)
|x3(kmin)|2
3 (kmin)
|x3|2
1+3q21
1−q21
=1 + κ− I1(q1) > 0. (B42)
These relations show that while |x3(kmin)| = 0, the ratio
|x3(kmin)|2/ω(2)3 (kmin) remains finite.
ii. q1 < q1s, κ < κs. In this case, we must have
I1(q1, q2)) < I1
(1 − q21)/2
, (B43)
(see Eq. (B38)). Consequently q22 < (1 − q21)/2 so that
3 (kmin) > 0 and hence no condensate can develop,
|x3|2 = 0. Then, Eq. (B33) yields the equation
I1(q1, q2) = 1 + κ (B44)
which replaces Eq. (B36) and determines q2 as a function
of q1, q2 = q2(q1). Then, proceeding as in case (i) one
finds for the IC-DC phase boundary in the region κ < κs
JDC(κ)=
(κ+ 1)(1 + q22) + Ĩ2(q1, q2)
q21 π
4 [K(q1)−E(q1)]
. (B45)
Here, q1 = q1(κ) from Eq. (B31) and q2 = q2(κ) from
Eq. (B44)(with q1 = q1(κ)).
We note here that inside the DC phase, i.e., for J >
JDC(κ), where Q2 = λc = 0, the saddle-point values
of q1 and λa/J and hence of Q1 are independent of J ,
cf. Eqs. (B31), (B32). Hence the graphs of Q1 and λa for
J < JDC and for J > JDC join smoothly at J = JDC , cf.
Figs. 10, 13. Furthermore, it follows from Eq. (B42), that
the ratio (|x3(±kmin)|2/N▽)/(λa q22ω
3 (kmin)), which
occurs in the amplitude of the spin-spin correlation func-
tion, cf. Eq. (24a), is also independent of J inside the DC
phase and retains the value that it has attained at the
IC-DC transition line.
c. Stability of the phase boundary
In deriving the phase boundary from the condition r
0 we have tacitly assumed that the coefficient g2 of the
fourth order term in the LG expansion, Eq. (B1), is pos-
itive. In the remaining part of this appendix we will
sketch the steps which lead to the conclusion that this is
indeed the case.
Expanding in the expression (B1) for e
LG the coefficients
e2 and r2 w.r.t. the second order contributions to q1 and
λa, q
1 and λ
a , cf. Eqs. (B27) we obtain
LG = e
2 + r
2 + (g2 + g
2 +O(Q62) , (B46)
where
1 + ω
2 + ω
(B47)
is the contribution to the fourth order term of e
LG that
arises from the fourth order terms of the frequencies ωµ
in the sum in Eq. (14) whereas the contribution to e
of the expansion of e2 and r2 is
∂2q1e2|0 ∂q1∂λae2|0
∂λa∂q1e2|0 ∂2λae2|0
∂q1r2|0
∂λar2|0
. (B48)
(In Eq. (B48) the notations ∂2q1e2|0 etc. indicate that af-
ter the derivatives have been taken the variables q1, λa
etc. have to be replaced by their zeroth order values q
a etc.)
The evaluation of the contribution (B47) is straight-
foward: the coefficients ω
µ , µ = 1, 2, 3, were obtained
by solving Eq. (15) for ω iteratively to fourth order. As
the explicit expressions are rather lengthy and contain
no direct information, we refrain from presenting them
here. The sum over k that is required in Eq. (B47) was
done numerically. g2 was obtained in the form
g̃2(q1) , (B49)
where g̃2(q1) is a function of q1 alone which is always
positive so that g2 > 0 throughout. Remarkably, no
explicit dependence on the coupling constant J appears
in these results.
The evaluation of g′, Eq. (B48) requires the knowledge of
explicit expressions for q
1 and λ
a . These are obtained
by expanding e2 to first order in q
1 and λ
a , inserting
the results into the expression (B1) for e
LG and requiring
that the terms of order Q22 satisfy the extremum condi-
tions w.r.t. q1 and λa:
0 = q
e2|0 + λ(2)a ∂q1∂λae2|0 + ∂q1r2|0 ,(B50a)
0 = q
1 ∂q1∂λae2|0 + λ(2)a ∂2λae2|0 + ∂λar2|0 .(B50b)
The solution of these equations reads
= −M̂
∂q1r2|0
∂λar2|0
, (B51)
M̂−1 =
∂2q1e2|0 ∂q1∂λae2|0
∂λa∂q1e2|0 ∂2λae2|0
. (B52)
Inserting these results into Eq. (B48) one finds
g′2 = −
∂q1r2|0 ∂λar2|0
∂q1r2|0
∂λar2|0
. (B53)
While the second derivatives of e2 are obtained straight-
forwardly from Eq. (B25) the derivatives ∂q1r2|0 and
∂λar2|0 have to be calculated separately for the region
q1 < q1s, where there is no condensate, |x3(kmin)|2 = 0,
and for the region q1 > q1s, where |x3(kmin)|2 > 0. Fi-
nally, the result for g′ can be cast into the form
g′2 =
xq xλ
(B54)
where
M̂ ′ =
4q21 Λ(κ+ 1− Λ)
1− q21
−Λ (2− q21)Λ− κ− 1
(B55)
Λ ≡ 1
[K(q1)−E(q1)] (K, E : elliptic integrals)
(B56)
xq = q
1(Λ− κ− 1) +
1− q21
q1∂q1 Ĩ2(q1, q2)
q2=q2(q1)
−Θ(q1 − q1s)
1− q21
(κ+ 1− I1(q1)) ,
xλ = (
+ 1)(κ+ 1)− 1
Ĩ2(q1, q2)
q2=q2(q1)
. (B57)
Here Θ is the step function; the integrals I1(q1) and
Ĩ2(q1, q2) have been defined above, cf. Eqs. (B38) and
(B41), respectively. After numerical evaluation of these
integrals, we find that g′2 = g
2(q1) is positive for all values
of q1.
1 G. Misguich and C. Lhuillier, Frustration in Two-
Dimensional Quantum Antiferromagnets (World Scien-
tific, Singapore, 2004), chap. 5.
2 R. Moessner, Can. J. Phys. 79, 1283 (2001).
3 P. Lecheminant, B. Bernu, C. Lhuillier, L. Pierre, and
P. Sindzingre, Phys. Rev. B 56, 2521 (1997).
4 C. Waldtmann, H.-U. Everts, B. Bernu, C. Lhuillier,
P. Sindzingre, P. Lecheminant, and L. Pierre, Eur. Phys.
J. B 2, 501 (1998).
5 F. Mila, Phys. Rev. Lett. 81, 2356 (1998).
6 M. Mambrini and F. Mila, Eur. Phys. J. B 17, 651 (2000).
7 Z. Hiroi, M. Hanawa, N. Kobayashi, M. Nohara, H. Takagi,
Y. Kato, and M. Takigawa, J. Phys. Soc. Jpn. 70, 3377
(2001).
8 A. Fukaya, Y. Fudamoto, I. M. Gat, T. Ito, M. I. Larkin,
A. T. Savici, Y. J. Uemura, P. P. Kyriaku, G. M. Luke,
M. T. Rovers, et al., Phys. Rev. Lett. 91, 207603 (2003).
9 F. Bert, D. Bono, P. Mendels, J.-C. Trombe, P. Millet,
A. Amato, C. Baines, and A. Hillier, J. Phys. Condens.
Matter 16, 9829 (2004).
10 F. Bert, D. Bono, P. Mendels, F. Ladieu, F. Duc, J.-C.
Trombe, and P. Millet, Phys. Rev. Lett. 95, 087203 (2005).
11 S. Sachdev and N. Read, Int. J. Mod. Phys. B 5, 219
(1991).
12 S. Sachdev, Phys. Rev. B 45, 12 377 (1992).
13 W. Apel, T. Yavors’kii, and H.-U. Everts, J. Phys. Con-
dens. Matter 19, 145255 (2007); W. Apel, T. Yavors’kii,
and H.-U. Everts, ibid., Corrigendum, in press.
14 O. A. Starykh and L. Balents, Phys. Rev. Lett. 98, 077205
(2007).
15 E. Lieb and D. Mattis, J. Math. Phys. 3, 749 (1962).
16 S. Brehmer, H.-J. Mikeska, and S. Yamamoto, J. Phys.
Condens. Matter 9, 3921 (1997).
17 A. B. Harris, C. Kallin, and A. J. Berlinsky, Phys. Rev. B
45, 2899 (1992).
18 C. H. Chung, J. B. Marston, and S. Sachdev, Phys. Rev.
B 64, 134407 (2001).
19 C. H. Chung, J.B.Marston, and R. H. McKenzie, J. Phys.
Condens. Matter 13, 5159 (2001).
20 F. Wang and A. Vishwanath, Phys. Rev. B 74, 174423
(2006).
21 O. Tchernyshyov, R. Moessner, and S. L. Sondhi, Euro-
phys. Lett. 73, 278 (2006).
22 G. Misguich, B. Bernu, and C. Lhuillier, J. Low Temp.
Phys. 110, 327 (1998).
23 C. Waldtmann, H. Kreutzmann, U. Schollwöck,
K. Maisinger, and H.-U. Everts, Phys. Rev. B 62,
9472 (2000).
24 V. Subrahmanyam, Phys. Rev. B 52, 1133 (1995).
25 C. Raghu, I. Rudra, S. Ramasesha, and D. Sen, Phys. Rev.
B 62, 9484 (2000).
26 M. E. Zhitomirsky, Phys. Rev. B 71, 214413 (2005).
27 P. Nikolic and T. Senthil, Phys. Rev. B 68, 214415 (2003).
28 J. B. Marston and C. Zeng, J. Appl. Phys. 69, 5962 (1991).
29 C. Zeng and V. Elser, Phys. Rev. B 51, 8318 (1995).
30 A. V. Syromyatnikov and S. V. Maleyev, Phys. Rev. B 66,
132408 (2002).
31 R. Budnik and A. Auerbach, Phys. Rev. Lett. 93, 187205
(2004).
32 G. Misguich and P. Sindzingre, J. Phys. Condens. Matter
19, 145202 (2007).
33 R. J. Baxter, J. Math. Phys. 11, 784 (1970).
