Transfer-matrix renormalization group study of the spin ladders with cyclic four-spin interactions

0.3 MB·English

by  H. T. Lu

Transfer-matrix renormalization group study of the spin ladders
with cyclic four-spin interactions
H. T. Lu,1,2 L. Q. Sun,2 Shaojin Qin,2 Y. J. Wang3
1School of Physics, Peking University, Beijing 100871, China
2Institute of Theoretical Physics and Interdisciplinary Center of Theoretical Studies,
Chinese Academy of Sciences, Beijing 100080, China
3Department of Physics, Beijing Normal University, Beijing 100875, China
7
0
The temperature dependenceof thespecific heat and spin susceptibility of the spin ladders with
0
2 cyclic four-spin interactions in the rung-singlet phase is explored by making use of the transfer-
matrix renormalization group method. The values of spin gap are extracted from the specific heat
n
and susceptibility, respectively. It is found that for different relative strength between interchain
a
andintrachaininteractions,thespingapisapproximatelylinearwiththecyclicfour-spininteraction
J
in the region far away from the critical point. Furthermore, we show that the dispersion for the
1 one-triplet magnon branch can be obtained bynumerically fitting on the partition function.
1

PACSnumbers: 75.10.Jm,75.40.Cx,75.40.Mg
]
l
e
r- I. INTRODUCTION S1,i−1 Jleg S1,i S1,i+1 S1,i+2
t
s
.
at Inthepastfewyears,cyclicspinexchangesinquantum J J J J
m antiferromagnetic magnets have attracted considerable rung cyc cyc cyc
attention. Fromvariousexperiments,ithasbeenrealized
-
d thatfour-spincyclicexchangesplayanimportantrolein
S J S S S
n understanding the spin dynamics of insulating cuprate 2,i−1 leg 2,i 2,i+1 2,i+2
o materials, such as the two-leg spin ladders Cu O in the
2 3
c compounds SrCu O ,1 (Sr, Ca, La) Cu O ,2,3,4 and FIG. 1: Schematic structure of a two-leg ladder with cyclic
[ 2 3 14 24 41
1 mthueltCipuleOs2pipnlaenxechiannLgeas2Cinu3OH4e.5(,6ReSfi.m7)ilaanrdmWecihganneirsmcryos-f (tfioonuro-sfptihne) finouterr-aspctinionte.rmTsheasroguivnednainrroEwqs. s(t1a)n,dwfhoirchthisenaoct-
v tal,8whichwasfoundtobelarge,hasalsobeenproposed. completely identical to theusual cyclic permutations.
3
4 Spin ladders can be viewed as intermediates between
2 one-dimensional spin chains and two-dimensional spin
1 lattices. They are particular of interest in the under- Takhtajan-Babujian critical point.17 This statement has
0 standing of the high-temperature superconductivity in been confirmed by numerical simulations.18,19
7 two dimensions and the Haldane’s conjecture9,10 in spin The Hamiltonian describing an isotropic S =1/2 spin
0
chains. The two-leg spin ladder is the simplest sys- ladder with additional cyclic (four-spin) interaction is
/
t tem on which cyclic exchanges can be materialized. In represented by
a
m the absence of cyclic exchanges, with nonzero interchain
couplings, the S = 1/2 spin ladder lies in a Haldane H = J S S +J [S S +S S ]
- rung 1,i 2,i leg 1,i 1,i+1 2,i 2,i+1
d phase, which has a short-range resonating-valence-bond { · · ·
i
n (RVB) ground state with a spin gap, and superconduc- +2XJ [(S S )(S S )
o tivity emerges upon doping.11,12 Contrary to the case cyc 1,i· 1,i+1 2,i· 2,i+1
v:c fiinneadSba=ck1/t2octhhaein,wethlle-dsepfiinneodnsqiunastihpearltaidcldeesr(amreascsoivne- +((SS11,,ii·SS22,,ii)+(1S)1(,Si+2,1i·SS21,,ii++11))] , (1)
− · · }
i magnons) in the neighborhood of q = π, where a spec-
X
tral gap exists.13 It has been shown that four-spin cyclic where the indices 1 and 2 distinguish upper and lower
r
a exchangescanmakethespingapdecreasesrapidly.2 Fur- legs,and i labels rungs; Jrung andJleg are the interchain
thermore, tuning of the four-spin interactions can drive and intrachain couplings, respectively; Jcyc denotes the
the systeminto differentphases.14,15 At the nearestcrit- strengthofthe cyclic(four-spin)exchange. Themodelis
ical point, the system becomes gapless to cross from the schematically shown in Fig. 1. In the rest of the paper,
RVB rung-singlet phase to a staggered dimer phase. By for convenience, we call Eq. (1) as type-I Hamiltonian.
a field theory using the Majorana Fermion representa- TheaboveHamiltoniandistinguishesfromanotherone
tion,NersesyanandTsvelik16referredtothecriticalphe- where the ring exchange term is directly represented by
nomenonasthelevelk =2SU(2)Wess-Zumino-Novikov- a cyclic permutation operator Pijkl in the plaquette:
Witten model or SU(2) c = 3 conformal field theory.
2
In other words, the phase transition is of the second- K P +P−1 1 .
order type and belongs to the universality class of the ijkl ijkl − 4
(cid:18) (cid:19)
2
Wecallthisformastype-IIHamiltonian. Thecyclicper- As proposed by experiments, for the insulator cuprate
mutation operator can be expressed by spin-1 operators laddermaterials,1,2,3,4,25 bothJ andJ arepositive
2 rung leg
including bilinear frustrating and biquadratic terms.2,20 and close in magnitude. The values of K/Jrung range in
So type-II Hamiltonian can be re-written as a region 0.025 K/Jrung 0.075, which roughly corre-
≤ ≤
sponds to the parameter region0.05 J /J 0.15
cyc rung
≤ ≤
H = (J +2K)S S for the type-I Hamiltonian and the system is still in the
rung 1,i 2,i
{ · rung-singletphase. Forthisreason,inournumericalsim-
i
+(XJleg+K)(S1,i S1,i+1+S2,i S2,i+1) ulations, we confine the parameter xcyc ∈ [0,0.1], and x
· · is taken three different values: 0.5, 1.0 and 1.2, where
+K(S S +S S )
1,i 2,i+1 2,i 1,i+1
· ·
+4K[(S S )(S S ) x := J /J , (3)
1,i 1,i+1 2,i 2,i+1 leg rung
· ·
+(S S )(S S ) x := J /J . (4)
1,i· 2,i 1,i+1· 2,i+1 cyc cyc rung
(S S )(S S )] . (2)
1,i 2,i+1 2,i 1,i+1
− · · } We take x=0.5 in orderto make a comparisonwith the
typical rung-singlet phase.
There have been extensive theoretical studies on the
type-IIHamiltonian. Bythedensitymatrixrenormaliza- Thenumericalalgorithmweemploytostudythe ther-
tiongroup(DMRG)21 andexactdiagonalizationmethod, modynamical properties of the present model is the
the energy spectra about one-magnon branch and two- TMRG method,26,27,28 which is a powerful numerical
magnon continuum in the rung-singlet phase have been tool for studying the thermodynamic properties of one-
obtained.4,10 The phase diagram for J = J was dimensional quantum systems. One is able to evalu-
rung leg
studiedindetailed14andaglobalpictureforgeneralcases ate nearly all thermodynamic quantities by the max-
has been proposed.15 The phase transitions,mainly con- imum eigenvalue and the corresponding left and right
cerning the case from the rung-singlet to the dimerized eigenvectors of the transfer matrix. We will show later
phase, together with the critical exponent were stud- that the magnitude of the gap and the low-lying excita-
ied by field theory, perturbative and series expansion tions can also be extracted from the obtained data. In
method.18,19,20,22 our numerical calculations, τ = 0.05, the error caused
by the Trotter-Suzuki decomposition is less than 10−3.
Ontheotherhand,thetype-IHamiltoniancatchesless
During the TMRG iterations, 80 100 states are re-
attention. And it is widely believed that the physics
tained and the truncation error is −less than 10−4 down
should be the same for the type-I and type-II Hamil-
to k T 0.01J .
tonians, although there is a lack of detailed compari- B rung
∼
sions. Theoretically, the ring exchange interaction natu-
rallyarisesasaneffectivelow-energymodelforthemulti-
band Hubbard model near the insulating filling and in II. SPECIFIC HEAT AND SPIN
the strong coupling limit. It has been shown that the SUSCEPTIBILITY
four-spin terms are most significant ones.23 In order to
highlight the effect of this four-spin term, in this paper, Figure 2 shows the TMRG results on the temperature
weemploythetype-IHamiltonianforournumericalsim- dependence of the specific heat coefficient C/T for x =
ulations. 0.5, 1.0, 1.2 with x = 0, 0.02, 0.04, 0.06, 0.08 and
cyc
As mentioned above, the zero-temperature behaviors 0.1. Besides the exponential decrease as T 0, which
→
ofthespinladderswithcyclicexchangeinteractionshave indicates the existence of a gap, the general behavior of
been analyzedin detail, especially for the type-II Hamil- the curves for a given x is as follows: As x increased,
cyc
tonian. A study on the thermodynamic properties of the slope of (C/T) shifts to lower temperatures. The
max
the systems has not been carried out extensively, ex- reasonis obvious: thatwith the increase ofthe four-spin
cept for the study by a high-temperature series expan- cyclic exchange, the spin gap and the whole dispersion
sion method.24 These quantities can be directly com- decrease.2,22 Another interesting feature is that at x =
pared with the experimental observations and therefore 1.0 and 1.2, a shoulder gets more and more distinct at
are of particular use. Detailed and reliable numericalre- lowtemperatureswiththe increaseofx . Weattribute
cyc
sults on these quantities may offer an expedient way to this shoulder to the tendency of deconfinement of bound
determine the model parameters for the compounds we spinons from the magnons, as the crossover occurs from
areinterestedin. Inthispaper,weshowthatbyusingthe the strong coupling regime (J /J =x 1) to the
leg rung
≪
transfer-matrixrenormalizationgroup(TMRG)method, weakcouplingone(J /J =x 1)withtheextreme
leg rung
≫
we can not only obtain the temperature dependence of limiting case of a finite C/T value at T =0.
the specific heat and susceptibility of the spin ladders The corresponding results for spin susceptibility χ(T)
down to the low-temperature regime and with high ac- arepresentedinFig.3. Theexponentialdecreaseofχ(T)
curacy, but also extract extra information from the data as T 0, and the shift of the peak to lower tempera-
→
forthe magnitudeofthe gapandtheone-tripletmagnon tures with the increase of x are consistent with the
cyc
excitations,whichcanbecomparedwiththeexperiments above behaviors of C/T. Our numerical result on the
directly. susceptibility well agrees with that of Ref. 24 obtained
3
1.5 0
x=0.5
x=0.5
)
1 3/2
T -10
C
(
0.5 n
l
-20
0
0.2 0.4 0.6 0.8
T 2) -10
C/ 1 x=1.0 1/T
ng χn( -20
Jru 0.5 l x=0.5
-30
0 10 20 30
0
0.2 0.4 0.6 0.8 1/T
1 x=1.2
FIG. 4: 1/T dependence of ln CT3/2 and ln χT1/2 (in
“ ” “ ”
0.5 theuniteof Jrung) forx=0.5 with xcyc=0, 0.02, 0.04, 0.06,
0.08,and0.1inascendingorderalongthedirectionofarrows.
0
0 0.2 0.4 0.6 0.8 1
T/J
rung by high-temperature series expansion and exact diago-
nalization.
FIG. 2: (Color online) The specific heat coefficient C/T for
The x dependence of the gap is a central quantity,
type-I Hamiltonian for three values of x with xcyc = 0, 0.02, cyc
since basically it determines the relevancy of the cyclic
0.04,0.06,0.08,and0.1,inascendingorderalongthedirection
of arrows. exchange term. In the vicinity of the critical point xccyc,
generally,
∆ x xc η .
∼| cyc− cyc |
0.2
Forthetype-IIHamiltonianinbothstrongandweakcou-
x=0.5
pling limits, by the perturbative cluster expansion and
0.1 Pad´e approximation in the strong coupling regime, to-
gether with the bosonization technique in the weak cou-
pling limit, Mu¨ller et al.20 showedthatthe criticalexpo-
0 nent η for the gap is approximately equal to 1; i.e., the
0.5 1 1.5
χrung 0.1 x=1.0 gFbayoprHtiishjieliineJteleaagrll.=y19dJHerpeureneng,dcweanestew,oinltlhsxiesceyrctehsnuaetlatrfowtrhatesycparlesi-toIicHoablatmpaoiinlitneotd-.
J
nian,intheregionconsiderablyfarawayfromthecritical
0 point, xcyc (0,0.1), the linear behavior of the gap vs.
0.5 1 1.5 ∈
x still holds in almost all the parameter range.
cyc
x=1.2 For the spin-1 ladder in the rung-singlet phase, sim-
0.1 2
ilar to the case in a S = 1 chain, we can assume that
the low-lying excitation spectrum ǫ(k) for k near π is
0 approximately given by
0.5 1 1.5 2
T/J
rung ǫ(k) = υ2(k π)2+∆2
−
υ2
FIG. 3: Temperature dependence of the susceptibility χ for p∆+ (k π)2+ k π 3 , (5)
≈ 2∆ − O | − |
type-I Hamiltonian for three values of x with xcyc = 0, 0.02, (cid:0) (cid:1)
0.04,0.06,0.08,and0.1inascendingorderalongthedirection
of arrows. where υ is the spin wave velocity and ∆ is the spectrum
gap. The low temperature behaviors of the specific heat
andsusceptibilityperspincanbeexpressed,respectively,
as29
4
0.8
x = 0.5 (C)
3/2 x = 0.5 (χ)
3∆ ∆
C(T) e−∆/T, (6) 0.6 x = 1.0 (C)
≈ υ√8π T x = 1.0 (χ)
(cid:18) (cid:19)
1 ∆
χ(T) e−∆/T, (7) ∆ 0.4
≈ υ 2πT
r
for T ∆. In this situation, ln CT3/2 and ln χT1/2 0.2
≪
shouldbelinearto1/T withaslope ∆. Figure4shows
the 1/T dependence of ln CT3(cid:0)/2 −and(cid:1) ln χT(cid:0)1/2 (in(cid:1)
0
units of Jrung) for x = 0.5 with xcyc = 0, 0.02, 0.04, 0 0.02 0.04 0.06 0.08 0.1
(cid:0) (cid:1) (cid:0) (cid:1)
0.06, 0.08, and 0.1. We observe a well-behaved linearity x
at 1/T > 10. For x = 1.0 and 1.2, we have the simi- cyc
lar results. The values of the gap at various x and x
cyc
can thus be extracted by fitting the data in the interval FIG.5: Thexcyc dependenceof thespin gap for x=0.5 and
1.0 obtained by fitting the Eqs. (6) (signed as C) and (7)
β = 1/T (10,20). The results for x = 0.5 and x= 1.0
∈ (given as χ), respectively.
arepresentedinFig. 5. Thegapsobtainedfromthe spe-
cificheatandsusceptibility arewellcoincidentwitheach
other.
Tobothverifythelow-temperaturescalingbehaviorof
1
thesystemandcheckthevalidityofEq. (7),inFig.6,the
rescaled susceptibility χ/χ versus T/∆ for x = 0.5,
max
1.0 at various x is presented. For a fixed value of
cyc
x, the curves converge to a single one at low tempera-
0.5
tures,demonstratingthe dominance of∆in determining
the susceptibility at low temperatures. Furthermore, it
also suggests that the values of the gap we obtained by ax x=0.5
m
this approach are quantitatively reliable. This situation
χ 0
is quite similar to the case of the S = 1 antiferromag- /
χ
neticHeisenbergchainwithsingle-ionanisotropy.30 Here
χ serves as a common factor which incorporates the
max
unknown effect on the prefactor in Eq. (7).
0.5
From Fig. 5, we observe that the gap is approxi-
mately linearly dependent on x for fixed x in the pa-
cyc
rameter space we studied. We also notice the differ- x=1.0
ent descending rates of ∆ with xcyc for various x. At 0
x = 0, ∆(x = 0.5) > ∆(x = 1.0) with the difference 0 2 4 6
cyc
0.15. This indicates that the gap decreases as the in- T/∆
∼
trachaincoupling J increases,whichis consistentwith
leg
the perturbative picture from the dimer limit of uncou-
FIG. 6: χ/χmax versus T/∆ for x = 0.5, 1.0 with xcyc = 0,
pled rungs(strongcoupling regime). Onthe otherhand,
0.02, 0.04, 0.06, 0.08, and 0.1 in ascending order along the
at x = 0.1, the difference of the gaps becomes very
cyc direction of arrows.
small. We expectforhighervalueofx areverseofthe
cyc
gaps with ∆(x=0.5)<∆(x=1.0) should occur.
For x = 1.0 and 1.2 (not shown in Fig. 5), the values
of the gap obtained by the present method well coincide than 10−2 at most. On the other hand, in obtaining the
withthenumericalresultsfromthedensity-matrixrenor- fitting formulas (6) and (7), the low-temperature limit
malizationgroup(DMRG)studyforthetype-IIHamilto- and quadratic approximation (See Eq. (5)) have been
nian for the corresponding values of the cyclic exchange adopted. For the specific heat, the expression including
coupling.4,11 For instance, when x = 1.0 and xcyc = 0, the higher order contributions of T/∆ is given explicitly
the gap obtained here by TMRG is 0.51, as compared as31
with the DMRG result of 0.504.11 ∼
To come to a close,we estimates the errorsin the pro-
3/2
3∆ ∆
cessofgapfitting. Ontheonehand,asmentionedinthe C(T)=
previous section, combining the errors from the Trotter- υ√8π T
(cid:18) (cid:19)
Suzuki decomposition with the TMRG truncation, we 2
T 3 T
estimate that the error caused by the TMRG algorithm 1+ + e−∆/T. (8)
in the fitting interval(in our workT (0.05,0.1))is less ×" ∆ 4(cid:18)∆(cid:19) #
∈
5
The uncertainty in Eq. (6) is roughly estimated as xcyc µ0 µ1 µ2 µ3 e0 ∆ ∆C
T¯ 3 T¯ 2 0 1.58 0.59 -0.38 0.09 -1.156 0.53 0.51
δ∆ T¯ln 1+ + , (9) 0.02 1.57 0.60 -0.39 0.10 -1.137 0.47 0.46
∼ " ∆ 4(cid:18)∆(cid:19) # 0.04 1.55 0.62 -0.41 0.12 -1.119 0.42 0.40
where T¯ represents some value of temperature T lying 0.06 1.54 0.63 -0.42 0.13 -1.101 0.37 0.35
in the fitting interval. From this, in our case the uncer- 0.08 1.53 0.64 -0.43 0.15 -1.084 0.31 0.30
taintyofgapduetothefittingformulascanbeestimated 0.1 1.52 0.65 -0.44 0.15 -1.068 0.27 0.26
approximately as δ∆ 0.01 0.02. Based on these con-
≈ −
siderations, we conclude that the error in the gap fitting TABLE I: The partition function fitting results by Eq. (10)
is about 10−2, which is not exceeding the symbol size in forthex=1.0 case. ∆C inthelast columnistheenergygap
Fig. 5. obtained byEq. (6),as plotted in Fig. 5.
III. PARTITION FUNCTION AND THE
DISPERSION RELATIONS
wouldbediluteatlowtemperatures. Takingintoaccount
of the form of the effective partition function for single
In TMRG calculations, the partition function
particledispersionZ ,weexpectitwouldmatchtheac-
Z (β) is the original and the most accurate numer- eff
TMRG curateTMRGresultsatbothhighandlowtemperatures
ical quantity. In previous section we obtained the spe-
by including higher order harmonics in the dispersion:31
cific heat and spin susceptibility, which can also be got
fromthe derivativesofthe partitionfunctiondirectly. In ǫ(k)= µ +µ cos(k)+µ cos(2k)+µ cos(3k). (11)
0 1 2 3
the following we will discuss the single-particle state in
this system and show that its dispersion, as well as the The energy gap is given by ∆ = µ0 µ1 +µ2 µ3 at
− −
ground-state energy and the magnitude of the gap, can k = π. Starting with a set of trial values for µ’s given
be obtained from the analyses of the partition function. by perturbative calculation,20 we minimize the squared
In the rung singlet phase, the elementary excitations 2-normof the residualin free energycalculatedfrom the
are gapped. The energy gap and dispersion for the sin- trial dispersion:
gleparticlestatecanbeanalyzedperturbatively.20When
W(µ ,µ ,µ ,µ ,e )=
J = J = 0, the system is consisted of an array of 0 1 2 3 0
leg cyc
rungs with a vacuum state being the product of spin N
singlet state on each rung. Above this vacuum state, [lnZeff(βi) lnZTMRG(βi)]2, (12)
−
a hard core boson of spin 1 is the only possible excited i=1
X
stateoneachrung. Thesystemhasaflatbanddispersion
where β are smoothly distributed in the interval (5,15)
ǫ(k) = J for single particle excitations. Due to the i
rung withaspacing∆β =0.05. Typically,weshowthe fitting
hardcorescatteringbetweenthebosons,thestatisticsof
resultsforthex=1.0caseinTableI. Thecorresponding
thesystemisessentiallyFermi-Diractype,32whichallows
dispersions are plotted in Fig. 7. In our final numerical
us to write down the partition function for the system: results, the values of W in Eq. (12) are of order 10−6,
1 π and the mean deviation at each data point is less than
lnZeff(β)= ln 1+3e−βǫ(k) dk βe0, (10) 10−4.
2π Z−π (cid:16) (cid:17) − From Table I, we find that µ0 and µ1 monotonically
where e is the ground-state energy per rung in the decreases and increases with x , respectively. This is
0 cyc
thermodynamic limit and the prefactor 3 before e−βǫ(k) consistent with the perturbative results.20 On the other
comesfromthethreeinternalstatesofthesingleparticle hand, the third order coefficient µ increases with x .
3 cyc
excitations. At x = 0.1, it reaches 0.15, about one-third of µ .
cyc 2
Intherungsingletphase,switchinguponthecouplings This indicates that the approximation up to the third
J and J will induce a dispersion4,10,20 for the single order in the dispersion is necessary for x = 1.0 case.
leg cyc
particle state, and the interaction between the particles. On the contrary, for x=0.5, µ 1, µ 0.5 0.6, the
0 1
∼ ∼ −
Classified by the internal symmetries, the two-particle absolutevalueofbothµ andµ arelessthan0.1,andthe
2 3
interaction is either spin-channel independent or spin- dispersion can be satisfactorily approximated by single-
channel dependent. A bound state may also appear in component harmonics with only µ . This fact reflects
1
certain spin-channels. This has been observed in a pre- that with the increase of J and J , higher orders in
leg cyc
vious study,4 in which it is shown that the bound state the harmonic expansion are needed to describe the one-
is alsomomentum dependent. Therefore,the interaction particle excitation in the spectrum.
isverycomplicatedanddependsonthetotalmomentum The dispersion relations shownin Fig. 7 should corre-
of the scattering particles. spond to the one-tripletmagnonbranchin the spectrum
Sincewearedealingwithasystemwithmassiveparti- of the ladder. We compare the results with the previous
clesatzerochemicalpotential,thedensityoftheparticles studiesbytheDMRG.4,10Wefindtheyhavequitesimilar
6
bytheTMRGmethod. Thetemperaturedependencesof
2 x=1.0
the specific heat and the susceptibility of the ladder in
the rung-singlet phase are analyzed numerically. Based
1.5 on the above results, we extracted the values of the spin
gapandfoundthatitmanifestsapproximatelylinearbe-
)
ε(k 1 haviorwith the cyclic interactioninthe parameterspace
we considered. We also showed that the dispersion for
low-lying excitations can be obtained from the partition
0.5 function very effectively by combining with the mean-
field or perturbative studies. Comparing with the high
0 temperature series expansion and the exact diagonaliza-
0 0.2 0.4 0.6 0.8 1 tion approaches,24 we believe that the results presented
k/π inthisstudygivemoredirectandaccuratephysicalquan-
tities in the thermodynamic limit. Further experimental
FIG. 7: (Color online) Dispersion (in units of Jrung) for x = results on specific heat and susceptibility may allow one
1.0withxcyc =0,0.02,0.04,0.06,0.08,and0.1indescending to determine the various parameters for the two-leg spin
order along thedirection of arrows. ladder materials via this approach. In addition, the dis-
persion relations extracted from the free energy may be
used to be compared with the experimental data of the
inelastic neutron scattering spectra.
shapes and are quantitatively agreement in the vicinity
Acknowledgments
of k=π. The deviation when approaching k =0 can be
attributed partially to the two-triplet bound state aris-
ing from the magnon-magnon interaction that we have
ignored.4
The energy gaps obtained in this way are presentedin
Table I, and are compared with their counterparts from
the specific heat data. We see that they are rather coin-
cident. In addition, we have obtained the ground state
energy (e in Eq. (10)), which is compatible with the
0
results from the zero-temperature DMRG scheme.
IV. SUMMARY This workwassupportedbythe NationalNaturalSci-
ence Foundation of China. The numerical work of this
Inthis paper,weexploredthe thermodynamicproper- projectwasperformedontheHP-SC45Sigma-Xparallel
ties of the spin ladders with cyclic four-spin interactions computer of ITP and ICTS, CAS.
1 A. G¨oßling, U. Kuhlmann, C. Thomsen, A. L¨offert, C. Phys. Rev.Lett. 86, 870 (2001).
Gross, nd W. Assmus, Phys. Rev.B 67, 052403 (2003). 9 F.D.M.Haldane,Phys.Lett.93A,464(1983);Phys.Rev.
2 S.Brehmer,H.-J.Mikeska,M.Mu¨ller, N.Nagaosa, andS. Lett. 50, 1153 (1983); J. Appl.Phys. 57, 3359 (1985).
Uchida, Phys.Rev.B 60, 329 (1999). 10 G. Roux, S. R. White, S. Capponi, A. L¨auchli, and D.
3 M. Matsuda, K. Katsumata, R. S. Eccleston, S. Brehmer, Poilblanc, Phys. Rev.B 72, 014523 (2005).
andH-J.Mikeska,Phys.Rev.B62,8903(2000); M.Mat- 11 S.R.White,R.M.Noack,andD.J.Scalapino,Phys.Rev.
suda,K.Katsumata,R.S.Eccleston,S.Brehmer,andH-J. Lett. 73, 886 (1994).
Mikeska, J. Appl.Phys. 87, 6271 (2000). 12 E. Dagotto and T. M. Rice, Science 271,618 (1996), and
4 T. S. Nunner, P. Brune, T. Kopp, M. Windt, and M. references therein.
Gru¨ninger, Phys. Rev.B 66, R180404 (2002). 13 D.G. Shelton,A.A. Nersesyan,and A.M. Tsvelik, Phys.
5 Y. Honda,Y. Kuramoto, and T. Watanabe, Phys. Rev.B Rev.B 53, 8521 (1996).
47, 11329 (1993). 14 A. L¨auchli, G. Schmid, and M. Troyer, Phys. Rev. B 67,
6 R.Coldea, S. M.Hayden,G. Aeppli,T. G. Perring, C. D. R100409 (2003).
Frost,T.E.Mason,S.-W.Cheong,andZ.Fisk,Phys.Rev. 15 V. Gritsev, B. Normand, and D. Baeriswyl, Phys. Rev. B
Lett. 86, 5377 (2001). 69, 094431 (2004).
7 G. Misguich, B. Bernu, C. Lhuillier, and C. Waldtmann, 16 A. A. Nersesyan and A. M. Tsvelik, Phys. Rev. Lett. 78,
Phys. Rev.Lett. 81, 1098 (1998). 3939 (1997).
8 S. Chakravarty et al. Philos. Mag. B 79, 859 (1999); M. 17 P.P.Kulish,N.Y.Reshetikhin,andE.K.Sklyanin,Lett.
Katano and D. S. Hirashima, Phys. Rev. B 62, 2573 Math. Phys. 5, 393 (1981); L. A. Takhtajan, Phys. Lett.
(2000); B. Bernu, Ladir. Cˆandido, and D. M. Ceperley, 87A, 479 (1982); H. M. Babujian, ibid. 90A, 479 (1982).
7
18 K. Hijii and K. Nomura,Phys. Rev.B 65, 104413 (2002). 26 R.J.Bursill,T.Xiang,andG.A.Gehring,J.Phys.: Con-
19 K. Hijii, Shaojin Qin, and K. Nomura, Phys. Rev. B 68, dens. Matter 8, L583 (1996).
134403 (2003). 27 X. Wangand T. Xiang, Phys. Rev.B 56, 5061 (1997).
20 M.Mu¨ller,T.Vekua,andH.-J.Mikeska,Phys.Rev.B66, 28 T. Xiang and X. Wang, in Density-Matrix Renormaliza-
134423 (2002). tion: A New Numerical Method in Physics, edited by I.
21 S.R.White,Phys.Rev.Lett.69,2863(1992); Phys.Rev. Peschel, X. Wang, M. Kaulke, and K. Hallberg (Springer,
B 48, 10345 (1993). New York,1999), pp. 149-172.
22 K.R.Schmidt,H.Monien, andG.S.Uhrig,Phys.Rev.B 29 T. Xiang, Phys.Rev.B 58, 9142 (1998).
67, 184413 (2003). 30 D.Coombes,T.Xiang,andG.A.Gehring,J.Phys.: Con-
23 E. Mu¨ller-Hartmann and A. A. Reischl, Eur. Phys. J. B dense. Matter 10, L159 (1998).
28, 173 (2002). 31 M. Troyer, H. Tsunetsugu, and D. Wu¨rtz, Phys. Rev. B
24 A.Bu¨hler, U.L¨ow, K.P.Schmidt,andG.S.Uhrig,Phys. 50, 13515 (1994).
Rev.B 67, 134428 (2003). 32 Jizhong Lou, Shaojin Qin, Tai-Kai Ng, Zhaobin Su, and
25 M. Windt, M. Gru¨ninger, T. Nunner, C. Knetter, K. P. Ian Affleck, Phys. Rev. B 62, 3786 (2000); K. Okunishi,
Schmidt, G. S. Uhrig, T. Kopp, A. Freimuth, U. Ammer- ibid. 60, 4043 (1999).
ahl,B.Bu¨chne,andA.Revcolevschi,Phys.Rev.Lett. 87,
127002 (2001).