Table Of ContentExotic disordered phases in the quantum J J model on the honeycomb lattice
1 2
−
Hao Zhang1,∗ and C. A. Lamas2
1Institute for Solid State Physics, University of Tokyo, Kashiwa, Chiba 277-8581, Japan
2Laboratoire de Physique Th´eorique, IRSAMC, CNRS and Universit´e de Toulouse, UPS, F-31062 Toulouse, France
We studytheground-state phasediagram of the frustrated quantum J J Heisenberg antifer-
1 2
−
romagnet on the honeycomb lattice using a mean field approach in terms of the Schwinger boson
representation of the spin operators. We present results for the ground-state energy, local magne-
tization, energy gap and spin-spin correlations. The system shows magnetic long range order for
0 J /J . 0.2075 (N´eel) and 0.398 . J /J 0.5 (spiral). In the intermediate region, we find
2 1 2 1
≤ ≤
3 two magnetically disordered phases: a gapped spin liquid phase which shows short-range N´eel cor-
1 relations (0.2075 .J /J .0.3732), and a lattice nematic phase (0.3732 .J /J .0.398), which
2 1 2 1
0 is magnetically disordered but breaks lattice rotational symmetry. The errors in the values of the
2 phase boundaries which are implicit in the number of significant figures quoted, correspond purely
tothe error in theextrapolation of our finite-size results to thethermodynamic limit.
n
a
J PACSnumbers: 75.10.Kt,75.10.Jm,75.50.Ee
1
1
I. INTRODUCTION
]
l
e Thetwo-dimensional(2D)Heisenbergmodelonbipar- e
1
- titelatticeshasbeenintensivelystudiedinthelastyears.
r
t Inthe unfrustratedcase,the classicalgroundstate isob-
s
. tained when all the spins in one sublattice are pointing J
at in a given direction whereas in the other sublattice the 2 e2
m spins are pointing in the opposite direction. However,
- in the quantum case this state is not the real ground
d state, in fact this is not an eigenstate of the Hamilto- J
n 1
nian. Thequantumgroundstateisexactlyknowninone
a
o 2
dimension1, but no exactresults for the two dimensional
c a
[ antiferromagnet are known, even for simple lattices as 3
the square lattice. However, several experimental and a1
1
numerical studies suggested that the ground state is in
v
factthe spinSU(2)symmetrybrokenN´eeltype state. In
8
contrast, when we include frustration in the system, for
2
4 example by including second neighbor interactions, the
FIG. 1. (Color online) The honeycomb lattice
2 ground state may become much more complicated.
with J and J couplings considered in this paper.
. In the quantum case, the ground state energy is lower 1 2
1 The lattice sites with different colors belong to differ-
0 thantheclassicalvalue,duetothequantumfluctuations. ent sublattices. The primitive translation vectors of
3 The effects of these fluctuations vary depending on the the direct lattice are e = √3/2,3/2 ,e = √3/2, 3/2 .
1 2
1 dimension, the spin quantum number, the presence of a =(0, 1), a = √3/2,1/2 and a = √3/2,1/2 −corre-
v: frustrating interactions and the coordination number of sp1ond t−o then2earest n(cid:2)eighb(cid:0)orbonds3. (cid:1) − (cid:0) (cid:1)(cid:3)
(cid:0) (cid:1) (cid:0) (cid:1)
i the lattice. One can ask what the quantum fluctuations
X
are when the coordination number is changed. In two
r dimensions two paradigmatic examples of unfrustrated the lowestallowedina 2Dsystem, quantumfluctuations
a
systemsarethesquarelattice,withcoordinationnumber couldbeexpectedtobestrongerthanthoseonthesquare
z = 4, and the honeycomb lattice with z = 3. Previous latticeandmaydestroytheantiferromagneticorder10–13.
results2,3 have shown that the staggered magnetization The study offrustratedquantum magnetsonthe hon-
is smaller in the z = 3 case. This behavior is in accord eycomb lattice has also experimental motivations14–20.
with the tendency towards a less classical behavior for One of the most exciting experimental progresses is one
systems of lower coordination number. kindofbismuthoxynitrate,Bi Mn O (NO ),whichwas
3 4 12 3
The inclusion of frustration in 2D quantum antiferro- obtained by Smirnova et al.14. In this compound the
magnets is expected to enhance the effect of quantum Mn4+ ions form a S = 3/2 honeycomb lattice with-
spin fluctuations and hence suppress magnetic order4. out any distortion. The magnetic susceptibility data in-
This idea has motivated many researchers to look for dicates two-dimensional magnetism. Despite the large
its realization5–9. A special scenario to check this is the AF Weiss constant of -257K, no long-range ordering
frustrated Heisenberg model on the honeycomb lattice. was observed down to 0.4K, which suggests a nonmag-
Due to the small coordination number (z = 3) which is netic ground state14–17. The substitution of Mn4+ in
2
a)
N´eel Spiral
0 1 0.5
6
b)
N´eel GSL VBC Spiral
0 0.2075 0.3732 0.398 0.5
FIG. 2. (Color online) Phase diagram as a function of the
frustration J /J . a) Classical phase diagram. b) Quantum
2 1
phasediagram corresponding toS = 1 obtained bymeansof
2
SBMFT.
Bi Mn O (NO ) by V4+ may lead to the realization of
3 4 12 3 FIG. 3. (Color online) Sketch of the staggered dimer VBC
theS =1/2Heisenbergmodelonthehoneycomblattice.
state which breaks the lattice rotational symmetry but pre-
serves thelattice translational symmetry.
The analysis of the honeycomb lattice from a more
general point of view has gained lately a lot of interest
II. MODEL AND OVERVIEW OF THE PHASE
bothcomingfromgraphene-relatedissues21andfromthe
DIAGRAM
possiblespin-liquidphasefoundintheHubbardmodelin
such geometry22–29. Due to these reasons,recently there
TheJ J Heisenbergmodelonthehoneycomblattice
ishugetheoreticalinterestinfrustratedHeisenbergmod- 1− 2
is given by
els on the honeycomb lattice, in which frustration is in-
corporatedby secondnearestneighborscouplings23,30–37
and maybe also third nearest neighbors couplings38–45. H =J1 Sˆx·Sˆy+J2 Sˆx·Sˆy, (1)
hXxyi1 hXxyi2
where Sˆ is the spin operator on site x and xy in-
Motivated by previous results, in this paper we study x h in
dicates sum over the n-th neighbors (see Fig. 1). In
the spin-1/2Heisenberg modelonthe honeycomblattice
thispaperweareinterestedintheantiferromagneticcase
with first (J ) and second (J ) neighbors couplings. Us-
1 2
(J ,J 0), and we focus on the region J /J [0,0.5].
ing a Schwinger boson mean field theory (SBMFT) we 1 2 ≥ 2 1 ∈
In the classical limit, S , the model displays dif-
find strong evidence for the existence of an intermediate
ferent zero temperature ph→as∞es48–50, see Fig. 2(a). For
disordered region where a spin gap opens and spin-spin
J /J < 1/6, the system is N´eel ordered, while for
correlations decay exponentially. This magnetically dis- 2 1
J /J > 1/6, the system shows spiral phases. For
orderedregionquantitatively agreeswell with recent nu- 2 1
merical simulation results36,37,39,45. Another key finding the quantum case, aspects of this model have been
explored previously in the literature by various ap-
of our work is the presence of two kinds of magnetically
proaches,includingspinwavetheory30,31,49,50,non-linear
disorderedphasesinthisregion. Oneisagappedspinliq-
uid (GSL)46,47 with short-rangeN´eel correlations,main- σ-modelapproach52, meanfieldtheory10,23,exactdiago-
nalization (ED)34,39,50, variational Monte Carlo (VMC)
tainingthelatticetranslationalandrotationalsymmetry.
method33,36, series expansion (SE)40, pseudofermion
The other phase is a staggereddimer valence-bond crys-
tal (VBC), which is also called lattice nematic30. This functional renormalization group (PFFRG)42 and cou-
pled cluster method (CCM)37. However, these works
phase breaks lattice rotational symmetry, but preserves
yielded conflicting physical scenarios.
lattice translational symmetry.
This model was studied by Mattsson et al.10 using
SBMFT with a mean field decoupling that considers
The restofthe paper isarrangedasfollows. InSec. II onlyantiferromagneticcorrelationsfornearestneighbors
we introduce our model and give a quick overviewof the and ferromagnetic correlations for next nearest neigh-
final phase diagram. In Sec. III the general formalism bors. ThisschemecanonlycorrectlydescribeN´eelorder.
oftheSchwingerbosonmean-fieldapproachispresented. MorerecentlyWang23studiedthismodelwithinSBMFT
In Sec. IV, using the solutions of mean field equations, including antiferromagnetic correlations for both near-
we discuss the phase diagram, especially the magneti- est and next nearest neighbors. Unfortunately, The au-
cally disordered region. We close with a summary and thor did not give the phase diagram for different val-
discussion in Sec. V. ues of J /J . Actually, for frustrated models we can
2 1
3
0.4
0.15 0.015 0.08
0.3 HLMQ00..000150 HLMQ000...000246 æ æ æ æ æ
0.000
0.2 0.204 0.208 0.00 æ
0.10 LQ æ J2(cid:144)J1 0.398 0J.420(cid:144)J21 0.406
0.2 æ
H
p M æ à
a æ àà à
G 0.05 æææ ààà
0.1 ææ à
ææ ààààà
æ
æ à
0.00 0.0 æææææææææææææææææææàààààààààà
0.0 0.1 0.2 0.3 0.4 0.5
a b c d J (cid:144)J
2 1
0.0 0.1 0.2 0.3 0.4 0.5
J /J FIG. 5. (Color online) Local magnetization determined by
2 1
N´eel GSL Spiral Eq. (31) extrapolated to thethermodynamiclimit as a func-
tion ofthefrustration J /J . Theshaded region corresponds
2 1
VBC to the magnetically disordered phases. Insets correspond to
theregionswherethemagnetizationforN´eel(left)andSpiral
FIG. 4. (Color online) Gap in the boson dispersion extrap- (right) phases becomes zero.
olated to the thermodynamic limit as a function of the frus-
tration J /J corresponding to S = 1/2. The gapped region
2 1
corresponds to two different magnetically disordered phases: phases based on the picture of the resonating valence
oneisGSL,theotherisstaggereddimerVBC.Inset: Z3order bondstates4,60–62. Asamerit,thismethoddoesnotstart
parameter defined in Eq. (32). The onset of the VBC phase fromany magnetic long rangeorderfor the groundstate
isdeterminedbythevalueofJ2/J1 where ψ isnon-zero(red (in contrast to spin wave theory), which should emerge
| |
arrows)
naturally if the Schwinger bosons condense at some mo-
mentum vector63. At this momentum vector, the lowest
excitation spectrum of the Schwinger bosons should be
not generally exclude either ferromagnetic or antiferro-
gapless. On the other hand, If the Schwinger bosons are
magnetic correlations53 and is important to use a mean
gapped, the phase is magnetically disordered. In the fol-
field decomposition that allows to include ferromagnetic
lowing,wewillpresentindetailtherotationallyinvariant
an antiferromagnetic correlations in equal footing. An- version of SBMFT which was introduced by Ceccatto et
other point is that both of them assume the bond mean al.54–56 and we use in the following sections.
fields are independent of the directions of bonds. There-
Consider the SU(2) Heisenberg Hamiltonian on a gen-
fore, these two schemes can not describe the phases in
eral lattice:
which the lattice rotational symmetry has broken. Here
1
we study the Hamiltonian (1) in the strong quantum Hˆ = J (x y)Sˆ Sˆ , (2)
2 αβ − x+rα · y+rβ
limit(S = 1/2) using a rotationally invariant version of
xyαβ
X
this technique, which has proven successful in incorpo-
where x and y are the positions of the unit cells and
rating quantum fluctuations38,53–59.
vectorsr arethe positionsofeachatomwithinthe unit
Our main results are summarized in Fig. 2(b). The α
cell. J (x y) is the exchange interaction between the
magnetic phase diagram is divided into four regions.51 αβ −
spins located in x+r and y+r .
At small values of the frustrating coupling J /J , the α β
2 1 In whatfollows we assume that the classicalorder can
system presents a N´eel-like ground state. By increasing
be parameterized as
the frustration, we find at J /J1 0.2075 a continuous
2
transitiontoagappedspinliquidp≃hase. Whenthevalue Sˆxx+rα =Ssinϕα(x) (3)
fiofndthaecforunsttinrautoiunsgtcroaunpsiltiniogneixnctoeeadsstJa2g/gJer1ed≃d0im.3e7r32V,BwCe Sˆxy+rα =0 (4)
(lattice nematic) with broken Z3 symmetry (See Fig. 3), Sˆxz+rα =Scosϕα(x), (5)
which transforms at J2/J1 0.398 into a spiral phase. with ϕ (x)=Q x+θ , where Q is the ordering vector
≃ α α
·
andθ arethe relativeanglesbetweenthe classicalspins
α
inside each unit cell.
III. SCHWINGER BOSON MEAN-FIELD ThespinoperatorsSˆ onsitexarerepresentedbytwo
APPROACH x
bosonsˆb (σ = , )
xσ
↑ ↓
It is well known that the SBMFT provides a natural Sˆ = 1 bˆ† ~σ bˆ , bˆ = ˆbr↑ , (6)
descriptionforbothmagneticallyorderedanddisordered r 2 r· · r r ˆb
(cid:18) r↓ (cid:19)
4
where ~σ =(σ ,σ ,σ ) are the Pauli matrices. Eq. (6) is where we have defined
x y z
a faithful representation of the algebra SU(2) if we take
into account the following local constraint A∗ (x y)= Aˆ† (x,y) (12)
αβ − h αβ i
2S =ˆb†x↑ˆbx↑+ˆb†x↓ˆbx↓. (7) Bα∗β(x−y)=hBˆα†β(x,y)i (13)
The exchange term can be expressed as h(Sˆ~x+~rα·Sˆ~y+~rβ)MFi=|Bαβ(~x−~y)|2−|Aαβ(~x−~y)|2,
and denotestheexpectationvalueinthegroundstate
Sˆ Sˆ =:Bˆ† (x,y)Bˆ (x,y): Aˆ† (x,y)Aˆ (x,y), at Th=i 0. It is convenient to change variables to R =
x+rα· y+rβ αβ αβ − αβ αβ x y, and eliminating x in the sums we obtain
(8) −
fiwnheedreaAsˆα,β(x,y)andBˆα,β(x,y)areSU(2)invariantsde- HˆMF = 21 Jαβ(R)(21 Bα,β(R)ˆb†R(+α)y,σˆb(yβ,σ)
RXyαβ Xσ h
1
Aˆα,β(x,y)= 2 σˆb(xα,σ)ˆb(yβ,−)σ (9) − σAα,β(R)ˆb†R(+α)y,σˆb†y(,β−)σ+H.C.
1Xσ − |Bα,β(R)|2−|Aα,β(R)|2 . i
Bˆ (x,y)= ˆb†(α)ˆb(β), (10)
α,β 2 x,σ y,σ (cid:0) (cid:1)(cid:9)
Xσ The mean field Hamiltonian is quadratic in the boson
withσ = , . Thedoubledots(:Oˆ :)indicatethenormal operators and can be diagonalized. It is convenient to
↑ ↓ transform the operators to momentum space
ordering of operator Oˆ. This decoupling is particularly
useful in the study of magnetic systems near disordered
1
phases,becauseitallowstotreatantiferromagnetismand ˆb(α) = ˆb(α)eik·(x+rα), (14)
ferromagnetisminequalfooting53–56. Ontheotherhand, x,σ √Nc k,σ
k
thisschemehasbeentestedtoobtainquantitativelyquite whereN isthe numberofXunit cells. After somealgebra
c
accurate results which show excellent agreements with and using the symmetry properties:
ED38,54–56.
ToconstructameanfieldHamiltonianweperformthe
J (R)=J ( R)
following Hartree-Fock decoupling αβ βα −
A (R)= A ( R) (15)
αβ βα
(Sˆx+rα ·Sˆy+rβ)MF =A[B∗α∗β((xx−yy))AˆBˆαβ((xx,,yy))+H.c] Bαβ(R)=B−β∗α(−−R),
− αβ − αβ
−h(Sˆx+rα ·Sˆy+rβ)MFi, (11) we obtain the following form for the Hamiltonian
1
HˆMF = 2 γαBβ(k)ˆb†k(σα)ˆb(kβσ)+γαBβ(−k)ˆb†−(kα−)σˆb(−βk)−σ−σγαAβ(k)ˆb†k(σα)ˆb†−(kβ−)σ−σγ¯αAβ(k)ˆb(kασ)ˆb(−βk)−σ
kXαβXσ n o
Nc J (R) B (R)2 A (R)2 , (16)
αβ αβ αβ
− 2 | | −| |
Rαβ
X (cid:2) (cid:3)
where with
1
γB (k)= J (R)B (R)e−ik·(R+rα−rβ) (17)
αβ 2 αβ αβ
R
X
1
γαAβ(k)=2XR Jαβ(R)Aαβ(R)e−ik·(R+rα−rβ) (18) Hˆλ = xα λ(α) σ ˆbx†(σα)ˆb(xασ)−2S!. (21)
1 X X
γ¯A (k)= J (R)A¯ (R)e−ik·(R+rα−rβ). (19)
αβ 2 αβ αβ
R
X
Now,weimposetheconstraint(7)inaverageovereach
sublattice α by means of Lagrange multipliers λ(α) Using the symmetries (15) we can see that both kinds
ofbosons( , )give the same contributionto the Hamil-
Hˆ Hˆ +Hˆ (20) tonian. Th↑en↓, we can perform the sum over σ to obtain
MF MF λ
→
5
1
Hˆ = (γB (k)+λ(α)δ )ˆb†(α)ˆb(β)+(γB ( k)+λ(α)δ )ˆb†(α)ˆb(β) σ γA (k)ˆb†(α)ˆb†(β) +γ¯A (k)ˆb(α)ˆb(β)
MF 2 αβ αβ k↑ k↑ αβ − αβ −k↓ −k↓− αβ k↑ −k↓ αβ k↑ −k↓
kXαβn (cid:16) (cid:17)o
N
c J (R) B (R)2 A (R)2 2SN λ(α).
αβ αβ αβ c
− 2 | | −| | −
Rαβ α
X (cid:2) (cid:3) X
ItisconvenienttointroducetheNambuspinorbˆ†(k)= where N is the total number of unit cells and S is the
c
bˆ† ,bˆ where spin strength. The mean field equations (28) and (29)
k↑ −k↓ mustbe solvedina self-consistentwaytogetherwith the
(cid:16) (cid:17)
bˆ† =(ˆb†(α1),ˆb†(α2),...,ˆb†(αnc)) (22) constraints (30) on the number of bosons.
k↑ k↑ k↑ k↑ Finding numerical solutions involves finding the roots
bˆ =(ˆb†(α1),ˆb†(α2),...,ˆb†(αnc)) (23) of the coupled nonlinear equations for the parameters A
−k↓ −k↓ −k↓ −k↓
and B, plus the additional constraints to determine the
and n is the number of atoms in the unit cell. Now, we
c values of the Lagrange multipliers λ(α). We perform the
can rewrite the Hamiltonian into a compact form:
calculations for finite but very large lattices and finally
H = bˆ†(k) D(k) bˆ(k) (24) we extrapolate the results to the thermodynamic limit.
MF
· · We solve numerically for different values of the frus-
k
X
tration parameter J /J and with the values obtained
(2S+1)N λ(α) H , 2 1
− c −h MFi for the MF parameters and the Lagrange multipliers we
Xα compute the energy and the new values for the MF pa-
where the 2nc 2nc dynamical matrix D(k) is given by rameters. We repeat this self-consistent procedure until
×
theenergyandtheMFparametersconverge. Afterreach-
γB (k)+λ(α)δ γA (k)
D(k)= αβ αβ − αβ . ing convergence we can compute all physical quantities
γA (k) γB (k)+λ(α)δ
(cid:18) αβ αβ αβ(cid:19) like the energy, the excitation gap, the spin-spin correla-
To diagonalize the Hamiltonian (24) we need to per- tionandthelocalmagnetization. Duringthecalculation,
form a para-unitary transformation of the matrix D(k) itisconvenienttofixtheenergyscalebysettingthevalue
which preserves the bosonic commutation relations64. of the nearest-neighbor coupling J =1.
1
We candiagonalizethe Hamiltonianby defining the new
operators aˆ =F bˆ, where the matrix F satisfy
·
IV. RESULTS
I 0
(F†)−1 τ (F)−1 =τ , τ = 2×2 . (25)
· 3· 3 3 (cid:18) 0 −I2×2 (cid:19) In Fig. 4, we show the boson dispersion relation gap
With this transformation, the Hamiltonian reads extrapolated to the thermodynamic limit as a function
of the frustration (J /J ). In the gapped region, the
2 1
HˆMF = aˆ†k·E(k)·aˆk−(2S+1)Nc λ(α)−hHˆMFi, absence of Bose condensation indicates that the ground
k α state is magnetically disordered. This result agrees well
X X
(26) with recent ED39, VMC36 and CCM37,45 studies. In the
gaplessregion,the excitationspectrumis zeroata given
where
wave vector k∗ = Q/2, where the Boson condensation
E(k)=diag(ω1(k),...,ωnc(k),ω1(k),...,ωnc(k)).(27) occurs. Thisischaracteristicofthemagneticallyordered
phases. Thestructureofthesephasescanbe understood
In terms of the original bosonic operators, the mean
through the spin-spin correlation function (SSCF) and
field parameters are
the excitation spectrum. Some typical examples for dif-
A (R)= 1 eik(R+rα−rβ) ˆb(α)ˆb(β) ferent phases will be shown later.
αβ 2N h k↑ −k↓i Topindowntheprecisephaseboundariesbetweenthe
c
Xk n magnetically ordered and disordered phases, we intro-
e−ik(R+rα−rβ) ˆb(α) ˆb(β) (28) duce the localmagnetizationM(Q) as anorder parame-
− h −k↓ k↑i
1 o ter,whichis obtainedfromthe longdistance behaviorof
Bαβ(R)=2N eik(R+rα−rβ)hˆbk†(↑β)ˆb(kα↑)i the spin-spin correlation function (SSCF)54,55:
c
Xk n
− e−ik(R+rα−rβ)hˆb†−(kβ↓)ˆb(−αk)↓i (29) |x−lyim|→∞hSx·Syi≈M2(Q)cos[Q·(x−y)]. (31)
o
andtheconstraintinthenumberofbosonscanbewritten In Fig. 5, we show the local magnetization for J /J
2 1
in the momentum space as ∈
[0,0.5]. For J /J = 0, the local magnetization is
2 1
hˆb†k(↑α)ˆb(kα↑)i+hˆb†−(kα↓)ˆb(−αk)↓i =2SNc, (30) tMhe(Qse)co=n0d.2o4r1d7e6r,spwihnicwhaviescianlceuxlacetilolenntreasuglrteeomf0e.n2t41w8i.t6h5
Xk n o
6
et al.36, as well as 0.385 0.010 by Bishop et al.37. In
-0.80 ±
thisregion,theSSCFshowsdifferentpropertiesindiffer-
entdirections,however,itexhibitslongrangeorderinall
-0.85
directions. Thegaplesspointsoftheexcitationspectrum
-0.90 movecontinuouslyinsidethefirstBrillouinzoneasJ2/J1
changes. Thisresultscorrespondtoaspiralphase. Inthe
c-0.95 classical version (S ) of the model (See Fig. 2(a)),
N → ∞
/s for J2/J1 > 1/6 there remains a line-type degeneracy in
Eg-1.00 whichthespiralwavenumberisnotdetermineduniquely
and is allowed on a ring in the Brillouin zone.49,50 Our
-1.05 results suggest that the classical degeneracy is lifted in
thequantumversion,wheresomespiralwavevectorsare
-1.10
favored by quantum fluctuations from the manifold of
classicallydegeneratespiralwavevectors. Thisspiralor-
-1.15
0.0 0.1 0.2 0.3 0.4 0.5 derbydisorderselectionwasalreadyseenbyusingaspin
waveapproachbyMulderet al.,30 andwehaverecovered
J /J
2 1 this selection with a different approach.
FIG. 6. (Color online) Ground-state energy per unit cell The most interesting part of the phase diagram is the
extrapolatedtothethermodynamiclimit asafunctionofthe intermediate region which has no classical counterpart.
frustrationJ /J . Theregionsofthefourdifferentphasesare
2 1 In this region, the nonmagnetic ground state retains
indicated using thesame colors that are used in Fig. 2.
SU(2) spin rotational symmetry and the lattice trans-
lational symmetry, However, it may break the Z direc-
3
tionalsymmetry ofthe lattice. FollowingMulder et al.30
This value is significantly reduced by quantum fluctua-
weintroducetheZ directionalsymmetrybreakingorder
tionscomparedwiththeclassicalvalue0.5. Thequantum 3
parameter ψ where
Monte Carlo (QMC) result66 is 0.2677(6), which is con- | |
siderably largerthan ours. For the unfrustrated case, all ψ = S (r) S (r) +ω S (r) S (r+e )
A B A B 1
h · i h · i
themeanfieldapproachesarequiteinaccuratecompared +ω2 S (r) S (r e ) . (32)
A B 2
with much more controlled techniques like QMC. The h · − i
difference in the M(Q) values of about 10%, provides, Here A, B correspond to the two different sublattices,
in the absence of any other quantitative evidence for the r denotes the unit cell position, and ω = exp(i2π/3).
accuracyofthemethodasappliedtothismodel,anindi- Equivalently,Okumuraetal.32 definem =ε a +ε a +
3 1 1 2 2
cationoftheaccuracyofthemethodandofalltheresults ε a ,whereε (µ=1,2,3)arebondenergiescorrespond-
3 3 µ
quoted that depend on the order parameters, including ing to the three nearest neighbor bonds a (µ =1,2,3).
µ
the phaseboundaries. However,the meanfieldapproach It is trivial to see ψ = m . This order parameter
3
| | | |
is still very useful to study gapped phases in frustrated is zero when the spin correlations along the three di-
systems. On one hand it is well known that for frus- rections are equal. We find that ψ keeps zero when
| |
trated systems QMC presents the famous sign problem. J /J . 0.3732; it becomes non-zero continuously at
2 1
On the other hand, the study of quantities like energy J /J1 0.3732 as shown in Fig. 4.51 Therefore, in the
2
≃
gap requires the study of big sizes clusters and the use region 0.2075 . J /J . 0.3732, the ground state pre-
2 1
of exact diagonalization for small size clusters makes it serves the Z lattice rotational symmetry. The SSCF
3
very difficult to extrapolate the results. shows short range antiferromagnetic correlations in all
As J /J increases, the local magnetization decreases. directions, and the minimum of the excitation spectrum
2 1
It vanishes continuously at J /J 0.2075, as shown in remains pinned at the Γ point. Namely, the system re-
2 1
Fig. 5.51 This value is inexcellenta≃greementwith recent mains to be a GSL. The appearance of the GSL agrees
numerical results, such as 0.2 by Mezzacapo et al.36 us- withrecenttwodifferentVMCstudies.33,36 Intheregion
ingVMCwithanentangled-plaquettevariationalansatz, 0.3732.J /J .0.398,theZ latticerotationalsymme-
2 1 3
as well as 0.207 0.003 by Bishop et al.37 using CCM. tryhasbroken. Wefindthatthevaluesofthemeanfields
±
The shift of N´eel boundary compared with the classical A and B: A(B) = A(B) = A(B) ; the bond ener-
a2 a3 6 a1
estimate 1/6 is due to quantum fluctuations which pre- gies have the same property: ε = ε = ε . Therefore,
2 3 1
6
fer to collinear N´eel rather than spiral phases in some the system should belong to the staggered dimer VBC
cases.39 In this region, the SSCF is antiferromagnetic in (lattice nematic). To further analyze this region, one
all directions, and the Boson condensation happens at needtocalculatethedimer-dimercorrelations. However,
the Γ pointof the firstBrillouinzone: k∗ =(0,0),which it is out ofthe scope ofthe presentpaper. The existence
corresponds to the ordering vector Q=(0,0). As J /J ofthestaggereddimerVBCisinagreementwitharecent
2 1
decreases from 0.5, the local magnetization M(Q) de- ED study,34 a bond operator mean field study,30 and a
creases. It vanishes continuously at J /J 0.398, as VMC study.33
2 1
shown in Fig. 5.51 This value is also in good≃agreement The errors in the values of the phase boundaries that
with recent numerical results, such as 0.4 by Mezzacapo are implicit here in the number of significant figures
7
0 5 10 15 20 0 5 10 15 20
0.118
> 0.102 J2/J1=0.18 > J2/J1=0.18
R R 0.059
S S
<S.0 0.051 <S.0 0.000
0.000
-0.059
-0.051 (a) (a)
0.063
0.052
S>R 0.042 J2/J1=0.36 S>R J2/J1=0.36
S.0 0.021 S.0 0.026
< <
0.000
0.000
(b) -0.026 (b)
0.069
S>R 0.046 J2/J1=0.38 S>R 0.058 J2/J1=0.38
S.0 0.023 S.0 0.029
< <
0.000
0.000
(c) -0.029 (c)
0.184
0.12
>R J2/J1=0.48 >R 0.092 J2/J1=0.48
S S
<S.0 0.06 <S.0 0.000
0.00
-0.092
-0.06 (d) (d)
0 5 10 15 20 0 5 10 15 20
R R
FIG. 7. (Color online) SSCF for a system of size N = 2 FIG. 8. (Color online) SSCF for a system of size N =
×
50 50inthezigzagdirectioncorrespondingtothefourdiffer- 2 50 50 in the armchair direction corresponding to the
× × ×
ent phases: (a)J /J = 0.18 (N´eel), (b)J /J = 0.36 (GSL), four different phases: (a)J /J = 0.18 (N´eel), (b)J /J =
2 1 2 1 2 1 2 1
(c)J /J =0.38(staggereddimerVBC),and(d)J /J =0.48 0.36 (GSL), (c)J /J = 0.38 (staggered dimer VBC), and
2 1 2 1 2 1
(spiral). (d)J /J =0.48 (spiral).
2 1
quoted, correspond purely to the error in the extrapola- the energy curve also supports that the three quantum
tionofourfinite-sizeresultstothethermodynamiclimit. phase transitions are continuous.
In no way are they intended to represent the essentially In the following we show several typical examples for
unknownerrorsimplicitinthe mean-fieldapproach,e.g., the four different phases. The SSCF along zigzag and
the10%differenceinM(Q)comparedwiththeQMCre- armchair directions for a system of 5000 sites is shown
sult in the unfrustrated limit. All the transition values in Fig. 7 and Fig. 8 for J /J =0.18(N´eel), 0.36(GSL),
2 1
presented in this paper correspondto mean field estima-
0.38(staggereddimerVBC)and0.48(spiral). Thecorre-
tions. In order to improve these values, it is necessary
sponding lowest excitation spectrum is shown in Fig. 9.
to studyindetailthe phasetransitionsbeyondthe mean
Although it is a finite size system, we can still see the
fieldlevel,whichis outofthe scopeofthe presentpaper.
corresponding properties for the four different phases as
In Fig. 6 we show the results for the ground state en- we have presented above. For J /J = 0.18, the SSCF
2 1
ergy per unit cell extrapolated to the thermodynamic in both of the zigzagandarmchairdirections shows long
limit. For the unfrustrated case (J = 0), E /N =- range N´eel correlations, and the lowest excitation spec-
2 gs c
1.09779,which is in excellent agreementwith the second trum becomes gapless at the Γ point (for a finite size
order spin wave calculation result of 1.0978.65 Com- system there is a small gap which disappears after the
pared with published QMC results by−Reger et al.67: extrapolation). For J /J = 0.36, the SSCF in both
2 1
1.0890(9),and more recently by Lo¨w68: 1.08909(39), of the zigzag and armchair directions shows short range
− −
it has appreciable difference, as our previous discussion N´eelcorrelations,andthe minimum ofthe lowestexcita-
of the difference in the M(Q) values. Since energy esti- tionspectrumremainsattheΓpoint,however,thereisa
matesalwayshaveanintrinsicquadraticerror,compared large gap which does not disappear after the extrapola-
toanintrinsiclinearerrorforotherproperties,evensmall tion. ForJ /J =0.38,theSSCFdoesnotshowanylong
2 1
errors in the energy can be of significance. The shape of range correlation, and the short range correlations are
8
a b
2 0.6300 2 0.4870
0.5515 0.4429
1 0.4730 1 0.3987
0.3945 0.3546
ky0 0.3160 ky0 0.3105
0.2375 0.2664
-1 0.1590 -1 0.2223
0.08050 0.1781
-2 0.002000-2 0.1340
-2 -1 0 1 2 -2 -1 0 1 2
c kx d kx
2 0.4090 2 0.4840
0.3630 0.4235
1 0.3170 1 0.3630
0.2710 0.3025
ky0 0.2250ky 0 0.2420
0.1790 0.1815
-1 0.1330 -1 0.1210
0.08700 0.06050
-2 0.04100 -2 0.000
-2 -1 0 1 2 -2 -1 0 1 2
kx kx
FIG. 9. (Color online) Momentum dependence of the lowest excitation spectrum for a system of size N = 2 50 50
× ×
corresponding to the four different phases: (a)J /J =0.18 (N´eel), (b)J /J =0.36 (GSL), (c)J /J =0.38 (staggered dimer
2 1 2 1 2 1
VBC), and (d)J /J =0.48 (spiral). The dashed hexagon denotes thefirst Brillouin zone of thelattice.
2 1
different along the zigzag or armchair directions, which parisonbetween the mean field results andthe corrected
is a indication that the lattice rotational symmetry is results was made. However, this hard calculation allows
broken. Simultaneously,theminimumofthelowestexci- only to calculate some special quantities like the ground
tationspectrum is awayfrom the Γ pointand the lattice state energy or spin stiffness. The corrections developed
rotational symmetry is clearly broken. There is also a by Trumper et al. could be extended to spiral phases70,
gap in this region which remains finite in the thermody- whichwouldallowtoinvestigate,forinstance,thepresent
namiclimit. ForJ /J =0.48,theSSCFshowsmagnetic model.
2 1
longrangecorrelationsinbothofthezigzagandarmchair
directions. Since one component of the ordering vector
Q = 0 (corresponding to k∗ = 0 in the lowest excita-
x x V. SUMMARY AND DISCUSSION
tion spectrum), the SSCF is N´eel-like along the zigzag
directions. This result agrees well with the spin wave
calculations by Mulder et al..30 In the present paper, we have investigated the quan-
tum J J Heisenberg model on the honeycomb lat-
1 2
−
Finally, we would like to talk about the next step of ticewithinarotationallyinvariantversionofSBMFT.In
our work. We have used a mean field approach based the region J /J [0,0.5], the quantum phase diagram
2 1
in the Schwinger boson representation of the spin oper- of the model displ∈ays four different regions.51 The mag-
ators. This mean field approach has the drawback of netic long range order of N´eel and spiral types is found
being defined in a constrained bosonic space, with un- for J /J .0.2075 and J /J & 0.398, respectively. For
2 1 2 1
physical configurations being allowed if this constraint the spiral region, we get the spiral order from quantum
is treated as an average restriction. This drawback can disorder selection which agrees with Mulder et al.30 us-
be in principle corrected by including local fluctuations ing spin wave theory. In the intermediate region, the
of the bosonic chemical potential.69 This correction was energygapis finite while the localmagnetizationiszero,
calculatedbyTrumperet al.56 fortheJ J squarelat- whichindicatesthepresenceofamagneticallydisordered
1 2
−
tice using collective coordinate methods, where a com- groundstate. WehaveusedtheZ directionalsymmetry
3
9
breakingorderparameter ψ definedinEq. (32)toclas- metry is broken or not in the region 0.2.J /J .0.4.
2 1
sify this part into two diff|er|ent magnetically disordered InarecentstudyusingPFFRG42 theauthorshaveob-
phases: one is a GSL which shows short-range N´eel cor- tainedthatwithinthemagneticallydisorderedregion,for
relations (J /J . 0.3732), the other is staggered dimer largerJ /J ,thereisastrongtendencyforthestaggered
2 1 2 1
VBC (lattice nematic), which breaks the Z directional dimerordering;forlowJ /J ,bothofplaquetteandstag-
3 2 1
symmetry (J /J & 0.3732). Considering the properties gered dimer responses are very weak. A further recent
2 1
of order parameters and the ground state energy, these studyusingCCM37 hasgotamorequantitativemagnet-
three quantum phase transitions seem to be continuous. icallydisorderedregion: 0.207 0.003<J /J <0.385
2 1
± ±
Aswehavementionedabove,recenttheoreticalstudies 0.010,inwhichthePVBCphasehasbeenreported. How-
of the phase diagram of the spin-1/2 J1 J2 Heisenberg ever,thegroundstatewithin0.21.J2/J1 .0.24ishard
−
modelonthehoneycomblatticehaveobtainedconflicting to be determined using this approach.
results. The central controversial point is the existence The other controversial point is the form of the mag-
and nature of magnetically disordered phases when the netic long range order when J2/J1 exceeds the magnet-
N´eelorderbecomesunstableasincreasingthefrustration ically disordered region. There are two proposals: the
J /J . Thereisagrowingconsensus30,33,34,36,37,39,42,45,50 anti-N´eelorder37 or the spiral order. It is difficult to get
2 1
that a magnetically disordered region should appear. a conclusion by ED since it is hard to treat the incom-
However, the nature of this region is still not clear with mensurate spin correlations due to small lattice sizes.39
differentapproachesgivingdifferentresults. AnearlyED Both of the recent SE40 and PFFRG42 studies have not
workbyFouetet al.50 firstclaimedthataGSLmightap- found any evidence for the existence of the anti-N´eel or-
pearintheregionJ /J 0.3 0.35,andforJ /J 0.4 der and concluded that the spiral state should be the
2 1 2 1
≈ − ≈
thesystemmightbeinfavorofthestaggereddimerVBC. stable ground state. However, both of the VMC with
A recentEDstudy by Mosadeqet al.34 has claimedthat EPVansatz36 andthe CCM37 studies supportthe oppo-
a plaquette valence bond crystal (PVBC) might exist site proposal. Since we areinterestedin the exoticdisor-
in the region 0.2 < J /J < 0.3, and a phase transi- dered phases in the magnetically disordered region and
2 1
tion from PVBC to the staggered dimer VBC exists at focusonJ2/J1 [0,0.5],wecannotexcludethepossibil-
∈
a point of the region 0.35 J2/J1 0.4. However, a ity that the anti-N´eel order state exists for J2/J1 >0.5.
morerecentEDworkbyAlb≤uquerque≤et al.39,whichhas Due to the existence of strong quantum fluctuations
treated larger system sizes, has been unable to discrim- and frustration, the spin-1/2 J1 J2 Heisenberg model
−
inate whether this magnetically disordered region cor- on the honeycomb lattice is a challenging model which
responds to PVBC with a small order parameter or a needs further investigation especially for the nature of
GSL. It is possible that the PVBC may just come from the intermediate phase. Unbiased numerical simulations
the finite size effects.36 ForlargerJ /J , ithas been also are still needed, such as the density matrix renormaliza-
2 1
hard to discriminate the staggereddimer VBC with spi- tion group (DMRG) method.71–73 Recently, DMRG has
ral phases, since ED is especially difficult to treat the been applied to spin-1/2 Kagome Heisenberg model74,75
incommensurate behavior of spin correlationsdue to the and square J1 J2 Heisenberg model76, and obtained
−
small lattice sizes. GSLs as the ground state. Since quantum fluctuations
There aretwo recentstudies of this model using VMC are expected to be stronger on the honeycomb lattice
with different variational wave functions. Clark et al.33 than those on the square lattice, it would be very inter-
haveusedHuse-Elserstatesandresonatingvalencebond estingtoapplyDMRGtothespin-1/2J1 J2Heisenberg
−
(RVB) states, and claimed that a GSL appears in the model on the honeycomb lattice.
region 0.08 J /J 0.3; a dimerized state which
2 1
≤ ≤
breakslatticerotationalsymmetryforJ /J &0.3. How-
2 1
ever, a more recentwork by Mezzacapo et al.36 using an ACKNOWLEDGMENTS
entangled-plaquette variational (EPV) ansatz have ob-
tained lower energy estimates, and claimed that in the We are especially grateful to Hirokazu Tsunetsugu for
magnetically disordered region 0.2 . J /J . 0.4, the his suggestions of this project and numerous enlighten-
2 1
PVBC order parameter vanishes in the thermodynamic ingdiscussionsfornumericalcalculations. We wouldlike
limit. Therefore, the PVBC may just come from the to thank Peng Li for fruitful discussions. Hao Zhang
finite size effects. Since the Z directional symmetry is supported by Japanese Government Scholarship from
3
breakingorderparameterhasnotbeenconsideredinthis MEXT of Japan. C. A. Lamas is partially supported by
paper,it is stillnot clearthat the lattice rotationalsym- CONICET (PIP 1691) and ANPCyT (PICT 1426).
∗ Corresponding author: [email protected] 3 AndersW. Sandvik.Phys.Rev.B 56, 11678 (1997)
1 Bethe HA.Z Phys.74, 205 (1931). 4 P. W. Anderson, Mater. Res. Soc. Bull. 8, 153 (1973);
2 Weihong Zheng and C. J. Hamer, Phys. Rev. B 47, 7961 P.FazekasandP.W.Anderson,Phil.Mag.30,423(1974);
(1993). P. W. Anderson,Science 235, 1196 (1987).
10
5 MaxA.MetlitskiandS.Sachdev,Phys.Rev.B77,054411 39 A.F.Albuquerque,D.Schwandt,B.Het´enyi,S.Capponi,
(2008). M.Mambrini,andA.M.L¨auchli,Phys.Rev.B84,024406
6 R.K.Kaul,M.A.Metlitski, S.SachdevandC.Xu,Phys. (2011).
Rev.B 78, 045110 (2008). 40 J. Oitmaa and R. R. P. Singh, Phys. Rev. B 84, 094424
7 L. Wang and A. W. Sandvik, Phys. Rev. B 81, 054417 (2011).
(2010). 41 D.J.J. Farnell, R.F.Bishop, P.H.Y.Li, J. Richter,and
8 R. Moessner , S.L. Sondhi and P. Chandra, Phys. Rev. B C. E. Campbell, Phys.Rev.B 84, 012403 (2011).
64, 144416 (2001). 42 J. Reuther,D.A. Abanin,and R.Thomale, Phys.Rev. B
9 A. Ralko, M. Mambrini and D. Poilblanc, Phys. Rev. B 84, 014417 (2011).
80, 184427 (2009). 43 P.H.Y.Li, R.F.Bishop, D.J.J. Farnell, J. Richter,and
10 A. Mattsson, P. Fr¨ojdh and T. Einarsson, Phys. Rev. B C. E. Campbell, Phys.Rev.B 85, 085115 (2012).
49, 3997 (1994). 44 R. F. Bishop and P. H. Y. Li, Phys. Rev. B 85, 155135
11 K. Takano Phys.Rev.B 74, 140402 (2006). (2012).
12 M. Hermele, Phys.Rev. B 76, 035125 (2007). 45 P.H.Y.Li,R.F.Bishop,D.J.J.Farnell,andC.E.Camp-
13 R. Kumar, D. Kumar and B. Kumar Phys. Rev. B 80, bell, Phys.Rev.B 86, 144404 (2012).
214428 (2009). 46 X. G. Wen , Phys.Rev.B 44, 2664 (1991).
14 O.Smirnova,M.Azuma,N.Kumada,Y.Kusano,M.Mat- 47 C.LhuillierandG.Misguich,inIntroduction toFrustrated
suda, Y.Shimakawa, T. Takei, Y. Yonesaki,and N. Kino- Magnetism, Eds. C. Lacroix, P. Mendels, and F. Mila,
mura, J. Am.Chem. Soc., 131, 8313 (2009). (Springer-Verlag, Berlin Heidelberg, 2011).
15 S.Okuboet al,J.Phys.: Conf.Series200,022042 (2010). 48 S. Katsura, T. Ide, and T. Morita, J. Stat. Phys. 42, 381
16 M.Matsuda,M.Azuma,M.Tokunaga,Y.Shimakawaand (1986)
N. KumadaPhys. Rev.Lett. 105, 187201 (2010). 49 E. Rastelli, A. Tassi, and L. Reatto, Physica B 97, 1
17 M.Azumaetal,J.Phys.: Conf.Series320,012005(2011). (1979).
18 Magnetic Properties of Layered Transition Metal Com- 50 J. B. Fouet, P. Sindzingre and C. Lhuillier, Eur. Phys. J.
pounds, Ed.L. J. DeJongh, Kluwer, Dordrecht (1990). B 20, 241 (2001).
19 A. Moller et al, Phys.Rev. B 78, 024420 (2008). 51 Theerrorsinthevaluesofthephaseboundarieswhichare
20 A.A. Tsirlin, O. Janson and H. Rosner, Phys. Rev. B 82, implicitinthenumberofsignificantfiguresquoted,onlyre-
144416 (2010). flecttheerrorintheextrapolationofourfinite-sizeresults
21 A.H. Castro Neto, F. Guinea, N.M.R. Peres, K.S. tothethermodynamiclimit.Theyarenotintendedtorep-
Novoselov, and A.K. Geim, Rev. Mod. Phys. 81, 109 resenttheessentially unknownerrorsimplicit inthemean
(2009). fieldapproach,e.g.,thedifferencesinthelocalmagnetiza-
22 Z.Y. Meng, T.C. Lang, S. Wessel, F.F. Assaad, and A. tion and the ground state energy per unit cell compared
Muramatsu, Nature464, 847 (2010). with theQMC results in the unfrustrated limit.
23 F. Wang, Phys. Rev.B 82, 024419 (2010). 52 T. Einarsson and H. Johannesson, Phys. Rev. B 43, 5867
24 Y.-M. Lu and Y.Ran, Phys.Rev.B 84, 024420 (2011). (1991).
25 H.Y. Yang and K.P. Schmidt, Europhys. Lett. 94, 17004 53 R.FlintandP.Coleman, Phys.Rev.B79,014424(2009).
(2011). 54 H. A. Ceccatto, C. J. Gazza and A. E. Trumper, Phys.
26 A.VaeziandX.G.Wen,arXiv:1010.5744v1[cond-mat.str- Rev.B 47, 12329 (1993).
el] (2010). 55 C.J.GazzaandH.A.Ceccatto,J.Phys.: Condens.Matter
27 A. Vaezi, M. Mashkoori, and M. Hosseini, Phys. Rev. B 5, L135 (1993).
85, 195126 (2012). 56 A.E.Trumper,L.O.Manuel,C.J.GazzaandH.A.Cec-
28 M.-T. Tran and K.-S. Kim, Phys. Rev. B 83, 125416 catto, Phys.Rev.Lett. 78, 2216 (1997).
(2011). 57 A.Mezio,C.N.Sposetti,L.O.ManuelandA.E.Trumper,
29 S.Sorella, Y.Otsuka,andS.Yunoki,ScientificReports2, Europhys.Lett. 94, 47001 (2011).
992 (2012). 58 H. Feldner, D.C. Cabra, and G. L. Rossini, Phys.Rev. B
30 A. Mulder, R. Ganesh, L. Capriotti and A. Paramekanti, 84, 214406 (2011).
Phys. Rev.B 81, 214419 (2010). 59 L. Messio, B. Bernu, and C. Lhuillier, Phys. Rev. Lett.
31 R. Ganesh, D.N. Sheng, Y.-J. Kim and A. Paramekanti, 108, 207204 (2012).
Phys. Rev.B 83, 144414 (2011). 60 D. P. Arovas and A. Auerbach, Phys. Rev. B 38, 316
32 S.Okumura,H.Kawamura, T.OkuboandY.Motome, J. (1988); A. Auerbach and D. P. Arovas, Phys. Rev. Lett.
Phys. Soc. Jpn. 79, 114705 (2010). 61, 617 (1988).
33 B.K.Clark,D.A.AbaninandS.L.SondhiPhys.Rev.Lett. 61 A. Auerbach, Interacting Electrons And Quantum Mag-
107, 087204 (2011). netism (Springer-Verlag, New York,1994).
34 H.Mosadeq,F.Shahbazi,andS.A.Jafari, J.Phys.: Con- 62 A. Auerbach and D. P. Arovas, in Introduction to Frus-
dens. Matter 23, 226006 (2011). trated Magnetism, Eds. C. Lacroix, P. Mendels, and F.
35 H.D. Rosales, D.C. Cabra, C. A.Lamas, P.Pujol, M. E. Mila, (Springer-Verlag, Berlin Heidelberg, 2011).
Zhitomirsky, arXiv:1208.2416 (2012) 63 J. E. Hirsch and S. Tang, Phys. Rev.B 39, 2850 (1989)
36 F. Mezzacapo and M. Boninsegni, Phys. Rev. B 85, 64 J. H. P. Colpa, Physica A 93, 327 (1978).
060402(R) (2012). 65 Weihong Zheng, J. Oitmaa, and C. J. Hamer, Phys. Rev.
37 R.F.Bishop,P.H.Y.Li,D.J.J.Farnell,andC.E.Camp- B 44, 11869 (1991).
bell J. Phys.: Condens. Matter 24, 236002 (2012). 66 Eduardo V. Castro, N. M. R. Peres, K. S. D. Beach, and
38 D.C.Cabra, C.A.Lamas, andH.D.Rosales, Phys.Rev.B AndersW. Sandvik,Phys.Rev.B 73, 054422 (2006).
83, 094506 (2011). 67 J. D. Reger, J. A.Riera, and A.P.Young, J. Phys.: Con-
dens. Matter 1, 1855 (1989).
