Table Of Content

Universality Class of Thermally Diluted Ising Systems at Criticality Manuel I. Marqués and Julio A. Gonzalo Depto. de Física de Materiales C-IV Universidad Atónoma de Madrid 0 Cantoblanco 28049 Madrid Spain 0 email: [email protected] 0 2 JorgeÍn˜iguez n Depto. de Física Aplicada II a Universidad del Pais Vasco J Apdo 644 48080 Bilbao Spain 6 (25-1-2000) 2 TheuniversalityclassofthermallydilutedIsingsystems,inwhichtherealizationofthedisposition ] ofmagneticatomsandvacanciesistakenfromthelocaldistributionofspinsinthepureoriginalIsing h model at criticality, is investigated by finite size scaling techniques using the Monte Carlo method. c We find that the critical temperature, the critical exponents and therefore the universality class e of these thermally diluted Ising systems depart markedly from the ones of short range correlated m disordered systems. Our results agree fairly well with theoretical predictions previously made by - Weinrib and Halperin for systems with long range correlated disorder. t a t s . t a m I. Introduction particular, long range correlation (LRC) has been found in X-ray and Neutron Critical Scattering experiments in - d During last decades the systems with quenched ran- systems undergoing magnetic and structural phase tran- n domness have been intensively studied [1]. The Harris sitions[13,14]. Thiseffecthasbeenmodeledbyassuming o criterion [2] predicts that weak dilution does not change a spatial distribution of critical temperatures obeying a c [ the character of the critical behavior near second order power law g(x) ∼ x−a for large separations x [15]. In phase transitions for systems of dimension d with spe- generalthesesystemsbehaveinawayverydifferentfrom 1 cific heat exponent lower than zero in the pure case (the RDIS, since systems with randomly distributed impuri- v 8 so called P systems), αpure < 0 =⇒ νpure > 2/d, being ties may be considered as the limit case of short range 8 ν the correlation length critical exponent. This crite- correlated(SRC) distributions of the vacancies. The ba- 3 rion has been supported by renormalization group (RG) sic approach to the critical phenomena of LRC systems 1 [3,4,5], and scaling analyses [6]. The effect of strong di- wasestablishedbyWeinribandHalperin[16]almosttwo 0 lution was studied by Chayes et al. [7]. For α > 0 decadesago. TheyfoundthattheHarriscriterioncanbe 0 pure (the so called R systems), the Ising 3D case for exam- extendedforthesecases,showingthatfora<d, the dis- 0 t/ ple, the systemfixed pointflowsfroma pure (undiluted) order is irrelevant if aνpure−2>0, and that in the case a fixed point towards a new stable fixed point [3,4,5,6,7,8] of relevant disorder a new universality class (and a new m at which α < 0. Recently Ballesteros et al. have fixed point) with correlation length exponent ν = 2/a, random - usedtheMonteCarloapproachtostudythedilutedIsing and a specific heat exponent α=2(a−d)/a appears. In d systems in two [9], three [10] and four dimensions [11]. contrast,ifa>d,the usualHarriscriterionforSRCsys- n The existence of a new universality class for the random tems is recovered. LRC disorder has been studied also o c diluted Ising system (RDIS), different from that of the bythe MonteCarloapproach[17]. Inthiscaseacorrela- : pure Isingmodelandindependent ofthe averagedensity tion function g(x) = x−a with a = 2 (defects consisting v i of occupied spin states (p), is proved, using an infinite inrandomlyorientedlines ofmagneticvacanciesinside a X volumeextrapolationtechnique[10]baseduponthelead- threedimensionalIsingsystem)confirmedthetheoretical r ingcorrectiontoscaling. Thecriticalexponentsobtained predictions of Weinrib and Halperin. a this way could be compared with the experimental crit- In the present paper we will study Ising three dimen- ical exponents for a random disposition of vacancies in sional systems where the long range correlated dilution diluted magnetic systems [12]. has been introduced as a thermal order-disorder distri- In all cases previously mentioned frozen disorder was bution of vacancies in equilibrium, governed by a char- alwaysproducedinarandomway,thatis,vacancieswere acteristic ordering temperature (θ), in a way similar to distributed throughout the lattice randomly. Real sys- the thermal disordering in a binary alloy of magnetic tems, however, can be realized with other kinds of dis- (spins) and non-magnetic atoms (vacancies). [See pre- order, where the vacancy locations are correlated. In vious work [18] and also the application of this kind of 1 disorder in percolation problems [19]]. Making more ex- III) and by using the effective-exponent approach (Sec- plicit the analogy with order-disorder in alloys we may tion IV). A summary of the main results and concluding distinguish clearly between: remarks are given in Section V. i) thermally diluted Ising system realizations (TDIS), in which the quenched randomness is produced II.Transitiontemperatureandself-averagingof by considering a ferromagnetic Ising system at (θ = thermally diluted systems T3D): after thermalization, the spins of the dominant c type(concentrationc≥0.5)aretakenasthelocationsof For a hypercubic sample of linear dimension L and the magnetic atoms, and the rest are taken as the mag- numberofsitesN =Ld,anyobservablesingularproperty netic vacancies. The structure of the realization is fixed X has different values for the different random realiza- thereafterfor alltemperatures atwhichthe magneticin- tions of the disorder,correspondingto the same dilution teractions are subsequently investigated. probability p (grand-canonical constraint). This means ii) random diluted Ising system realizations (RDIS) that X behaves as a stochastic variable with average X (equivalent to θ >> T3D), also fixed thereafter for all (in the following, the bar indicates average over subse- c temperatures at which the spin interactions are inves- quent realizations of the dilution and the brackets indi- tigated. Of course, we chose for comparison a vacancy cate MC average). The variance would then be (∆X)2, probability p=0.5, resulting in c≈0.5. and the normalized square width, correspondingly: iii)asathirdpossibility,notconsideredhere,onecould investigate what might be called an antiferromagnet- ically diluted Ising system (with θ << T3D, which N would lead to a disposition of non-magnetic atoms (vacancies) strictly alternating with spins, with c=0.5. So,ifthisorderingtemperatureθ determiningthepar- R =(∆X)2/X2 (1) X ticular realization is high enough, the equilibrium thermal disorder will be very similar to the random (short range correlated) disorder of the usual previous inves- tigations. On the other hand, if θ happens to coin- cide with the characteristic magnetic critical temperature (T3D = 4.511617) of the undiluted system, we will Asystemissaidtoexhibitself-averaging(SA)ifR →0 c X have vacancies in randomly located points, but with a asL→∞. Ifthesystemisawayfromcriticality,L>>ξ long range correlated distribution. [Note that the situa- (being ξ the correlation length). The central limit the- tiondiffersmarkedlyfromthatofpreviouslystudiedLRC orem indicates that strong SA must be expected in this systems,inwhichlinesorplanesofvacancieswerecon- case. However,the behavior ofa ferromagnetat critical- sidered]. The correlationofour TDISis givenby a value ity (with ξ >>L)is notsoobvious. This pointhas been a=2−η , where η is the correlationfunction ex- studied recently for short rangecorrelatedquenched dis- pure pure ponentforthepuresystem.Sinced=3andη =0.03 order. AharonyandHarris(AH),usingarenormalization pure for the three dimensional Ising system, we have a long groupanalysisind=4−εdimensions,provedtheexpec- range correlated disorder with a = 1.97 < 3 = d. So we tationof a rigorousabsence of self-averagingin critically are in the case where LRC disorder is relevant and we random ferromagnets [20]. More recently, Monte Carlo should detect a change of universality class with respect simulationswereusedto investigatethe self-averagingin to the SRC case (following Weinrib and Halperin we ex- critically disordered magnetic systems [10,21,22]. The pect for the thermally diluted Ising system ν ≈ 1 and absenceofself-averagingwasconfirmed. Thenormalized α ≈ −1). Details about the construction of these ther- square width R is an universal quantity affected just X mally diluted Ising systems (TDIS) can be found in Ref. by correction to scaling terms. LRC diluted systems are [18]. In the present work we study the critical behav- expected to have different critical exponents and differ- iorandtheuniversityclassofthreedimensionalTDISat entnormalizedsquarewidthswithrespecttothoseofthe criticalityusingtheMonteCarloapproach. Wewillcom- usual randomly disordered systems studied previously. pare our results with the critical behavior of the RDIS. We perform Monte Carlo calculations of the magne- The structure of the paper is as follows: In Section tization and the susceptibility (χ = ((cid:10)M2(cid:11)−hMi2)/T) II we study the dependence of the critical temperature per spin at different temperatures for different realiza- (and of the self-averaging at criticality) with the size of tions of TDIS, and for randomly diluted systems with the system for both TDIS and RDIS. Once the critical p = 0.5 (restricting to c > 0.5) using in both cases point is determined, we investigatewhether or notTDIS the Wolff [23] single cluster algorithm [24] with peri- andRDISbelongtothesameuniversalityclass. Inorder odic boundary conditions, on lattices of different sizes toproceed,westudythecriticalbehaviorofbothkindsof L = 10,20,40,60,80,100. Results for susceptibility vs. systemsbyapplyingfinitesizescalingtechniques(Section temperature are shown in Fig. 1. 2 has been determined by means of Monte Carlo data for L=20 2.00 the r anRdomDcasIeSby Ballesteros et al. [10]. They found a value ν = 0.683. On the other hand the result 40 random by Weinrib and Halperin [16] indicates that the criti- ( n =cal0ex.p6on8en3t e)xpected for the thermal case should be ν = 2/a = 1.015. Fig. 2 represents the depen- thermal dence of the critical temperature with respect to the length of the systems for random and thermal dilutions. 1.95 T (L) c 0 -1/ n 1.90 ~ L L=40 1.85 200 10 100 c 3.45 FIG. 1. Susceptibility χ vs. temperature T for random TDIS dilution realizations (p = 0.5) (black) and critical thermal dilution0realizations (white). The size of the systems under consideration is given byL=20,40,80. ( n =1.015) Note how due to the existence of randomness the sus- 3.40 ceptibilitypointsdonotcollapseintoasinglecurve,since each realization has a different value of the critical temperature and of the concentration. This is even more L=80 RDIS T (cLlear)for small values of L and for the critically thermal FIG.2. Semi log representation of critical temperature Tc case (white points). Fig. 1 indicates that the critical c1000 vs. lenght L for random dilution (p = 0.5) (black) and for temperature of the TDIS clearly differs from that of the TtDhermIaSl -di1lut/io nn (white). The continuous line is the fit ob- ~ L RDIS,andalsothanthe effectofthedilutiononthelack tainedusingν =0.6837(shortrangecorrelatedrandomexpo- 3.35 of self-averaging seems to be stronger in the thermally nent) and dashed line the fit obtained using ν = 1.015 (long diluted Ising case. The critical temperature may be ob- range correlated exponent,a=1.97) in eithercase. tained at the point where the normalized square width forthesusceptibility,R ,reachesitsmaximum. Adiffer- In both cases a fit to Eq. 2 has been performed for χ ent value of the critical temperature does not imply, of both values ν = νrandom = 0.683 (continuous line) and course, a different universality class. However from the ν = νthermal = 1.015 (dashed line). Note how the ther- dependenceofthecriticaltemperaturewiththelengthof mal data fit better Eq. 2 for Tc(L) using νthermal, indi- th3e sy.s03tem0we expect to detect a change in universality catingapossiblechangeinuniversalityclasswithrespect between TDIS and RDIS following the scaling relation: totherandomcase. [Ifwefitthedataleavingallparam- eters free, we find ν ≈ .0.7 and ν ≈ 1.2, random thermal 1T.(0L)=T 1(∞).+5AL−1/ν2.0 2(2).5which3are.ve0ry near3the.ex5pected r4esu.lt0s]. The e4xtr.ap5o- c c lated values of critical te1mpe0rat0ures for infinite systems being ν the critical exponent associated with the spe- obtainedthiswayareTrandom(∞)=1.845±0.003(close c cific system’s correlation length. This critical exponent to thLe values previously obtained by Ballesteros et al. T 3 [10]) and Tthermal(∞) = 3.269±0.002 (clearly different cording to Weinrib and Halperin [16] the analysis would c 00f1ro..m550T (∞) for the SRC case). Incidentally ν can be even more complicated, due to the fact that the long c thermal be compared with ν for the observed sharp component range correlated disordered systems present complex os- inneutronscatteringlineshapes,whichisaround1.3for cillating corrections to scaling. Tb [14]. This point deserves more careful analysis and RDIS will be taken up elsewhere. III. Critical behavior and exponents of ther- Once the critical temperatures are known we can perT- DmaIllySdiluted systems form simulations for the magnetization TandDthIeSsus cep- tibility at criticality for several realizations of thermal The dispersion in concentration and magnetization at andrandomdilutedsystemsinanefforttodeterminethe criticality between the different realizations is shown at value of R , the normalizedsquarewidth for the suscep- a glance in scattered plots as in Ref. [18]. Each point 00..44 χ tibility. We wentupto500realizationsforL=10,20,40 in Fig. 4 represents a single realization with magnetiza- and up to 200 realizations for L = 60,80. Results are tion at criticality (M) and concentration (c). Note that shown in Fig.3. in both cases (TDIS and RDIS), the dispersion on the magnetization and on the concentration decreases with L, but this is more clearly shown in the thermal case. Fig. 4 shows clearly the difference in behavior between 1 the random and the thermal cases, at least up to the RDIS values of L considered. 00..33 R c M 00..22 FIG. 3. Log-log plot of the normalized square width for Ballesteros 0.1 the susceptibility at criticallity Rχ vs. lenght L for random et al. (R.10) dilution (p=0.5) (black) and thermal dilution (white). The arrowrepresents the (concentrationindependent) R value obtainedby Ballesteroset al [10]. The straight χ continuousline representsthe averagevalue obtainedfor 00ran..d11om dilution data, and the straight dashed line gives L=20 the average value obtained for thermal dilution data. FIG. 4. Scattered plot of magnetization M at the criti- Note that the TDIS presents a lack of self-averaging cal temperature vs. concentration c for the realizations con- aroundoneorderofmagnitudelargerthantheRDIS.We L=40 sidered of the random (p = 0.5) and the thermal dilutions. havealreadypresentedasimilaranalysisforbothkindsof (L=20,40,60,80). dilution [18], but only at the critical temperature char- L=60 acteristic of the random system. Ourresultsarenot From Fig. 4 we can extract averaged values for the precise enough to specify accurately the evolution of the magnetizationandthe inversesusceptibility, M andχ−1 normalized square width as a function of L governed by (withχ=(cid:10)M2(cid:11)atthecrit icLal=poi8nt)0fortheTDIS.Both correctionstoscalingterms. Howevertheaverageweob- averaged values are expected to fit the following scaling 000.t0ai..n001for Rrandom =0.155 is close to the value previously χ laws at criticality: reported[10]bymeansofinfinite volumeextrapolations. For the therm1al c0ase we obtain Rthermal = 1.19, about 1M(0L)0∝L−β/ν (3) χ 0.50 0.55 0.60 0.65 0.70 0.75 one order of magnitude larger than for random dilution. In this case an evolution of Rthermal vs. L given by cor- L χ−1(L)∝L−γ/ν (4) χ rection to scaling terms may be also expected, but ac- c (concentration) Considering (1/ν) = a/2 = 0.985, we can ob- thermal 4 tain from our data the values of β and γ (c(L)−0.5)∝L−(β/ν)3D (6) thermal thermal (see f.i. the fitting for the inverse of the susceptibility in Fig. 5). We get βRDIS= 0.56 ± 0.05 and where (β/ν) ≃T0.5D2 giIvSes the values corresponding to thermal 3D γ =1.91±0.06,very close to the predicted values the pure three dimensional Ising case, because in criti- thermal byWeinribandHalperin[16(]g. U/ sinn g=th2e.s0ca7lin(g5re)l)ation: cally thermally di(lug t e/d nIs i=ng1sy.s8te8m(s,5v)a)cancies are distributed with the same long range correlation spin α=2−2β−γ (5) distribution function as in the pure case [18]. The fit to the average value of the concentrations shown in Fig. 4 we obtain the following specific heat critical exponent: give,forthethermalcase,avalue(β/ν) ≃0.55±0.08. 3D αthermal =−1±0.1. WeinribandHalperingiveforLRC This implies also a clear difference between RDIS and systems [16] α=−1, in good agreementwith our result. TDIS, since Ising systems with vacancies randomly dis- 0A.n0ana1logo2usanalysishas been performedbut using the tributed are not expected to follow an scaling behav- dispersioninmagnetizationandinversesusceptibilityin- ior with (β/ν) . [Fitting to an scaling law the results 3D stead of the averaged values. The results are similar. for RDIS gives an exponent around 1.4, which implies a Thesamestudyhasbeenmadefortherandomcase(also much faster convergence to c=0.5]. shown in Fig.5). IV. Effective exponents and possible crossover for finite systems The difference between the universality class of RDIS and TDIS can be detected also by means of the effective critical exponents. In the case of the magnetization the effective critical exponent is defined by: -1 c β ≡∂log(M)/∂log(t) (t=Ti−T) (7) eff c -g/n withTci the criticaltemperature ofthe particular real- ~ L ization (i) characterized by-ag/nmaximum of the susceptibility (χ = (cid:10)M2(cid:11)−~hMLi2). For L → ∞ and t → 0, β = β. Finite size effects force the effective critical eff exponent to drop to zero before the critical value is at- 0.006 tained. However, since β = 0.3546 (Ref. [10] ), random and β is expected to be around 0.5 (Ref. [16] ) thermal , effective critical exponents (for the thermal case) may rise to values greater than 0.3546 and lower than 0.5, before finite size effects appear, indicating the difference of universality class between both kinds of systems. In the RDIS the effective critical exponent is expected to be always lower than 0.3546. Monte Carlo simulations of magnetization vs. temperature have been performed for randomly (p = 0.5) and thermally diluted systems, with L = 80. In Fig. 6 we show the results for β eff vs. log(t) for two samples of the TDIS and the RDIS type respectively. The random effective critical exponent is always under 0.35 and it seems to tend towards this value for large enough L, as expected, but for thermal systems the behavior is completely different: In the FIG. 5. Average value of the inverse of the susceptibility χ−1(L)vs. thelenghtofthesystemLfortherandomdilution figurethe value ofthe criticalthermaleffective exponent (p=0.5) case (black) and the thermal dilution case (white). is between 0.35 and 0.5. An analogous investigation has Continuous line and dashed line indicate fits to the random been done for different values of L and different realiza- 0.000 and the thermal cases, respectively. tions. The same effect has been found for lower L values. However, we may note that in TDIS the effective The criticalexp0onents o2bta0inedar4ein0agreem6en0twith80critical 0exponent2ar0rives a4t th0e max6imu0m in a8ver0y dif- those of Ballesteros et al. [10] within our error bars. ferent way depending on the realization. The reason is The average concentration for the thermal system is twofold: (1)thedifferentdispositionofthevacanciesin also expected to show a scaling law behaLvior given by: each particular realization, andL(2) the large differences 5 in concentration for the thermal case (the rise of β the size of the fixed vacancies and spin clusters, the size eff towardstheexpecteddiluteduniversalityvalueshouldbe of the system itself and the thermal spin fluctuations faster for c closer to 0.5). It may be noted that as the might all play a role. The possibility that for L −→ ∞ b size(L)ofthes=amp0lei.nc5rea0sesthepossibilityoflocalin- and c = 0.5 the apparent crossover could disappear all homogeneities in the TDIS realizations increases, giving together should not be excluded, entailing no crossover 0ris.e5to such phenomena as pseudo double peaks in the from critical thermal to random universality class. susceptibility, reduced values for the overall critical exponents,etc. However,itis importantto remarkthanin V. Concluding remarks all realizations investigatedthe β values of the TDIS, eff (before finite size effects take over), have been clearly Tosummarize,anewwaytoproducedilutedIsingsys- superior to the β for the RDIS. tems with a long range correlated distribution of vacancies has been analyzed by means of Monte Carlo calculations and finite size scaling techniques. We find an universalityclassdifferentfromthatreportedfordiluted 0.4 Ising systems with short range correlated disorder. Our systems may be included in the universality class pre- b =0.35 dicted by Weinrib and Halperin for a ≈ 2. This kind of thermal disorder had been already applied in percolation problems, but as far as we know it had never been applied to magnetic systems, in which long range correlated disorder has been previously produced mostly by random distribution of lines or planes of vacancies. The present dilution procedure may be applicable to systems 0.3 where the long range correlated disorder is not due to dislocations,preferentiallinesorplanesofvacancies,but tosystemswherethevacancies(points)arecriticallydis- f tributed in clusters as in the case of order-disorder sys- f e tems. RDIS We thank H.E. Stanley for encourament and we are b gratefulto H.G. Ballesteros,L.A. Fernańdez, V. Mart´ın- Mayor,andA.MunõzSudupeforhelpfulcorrespondence. (c=0.506) We thank P.A.Serenaforhelpful comments andforgen- 0.2 erous access to his computing facilities. We acknowl- edgefinancialsupportfromDGCyTthroughgrantPB96- 0037 and from the Basque Regional Government (J.I.). TDIS (c=0.587) 0.1 [1] For a review, see R.B. Stinchcombe, in Phase Transitions FIG.6. Effective critical magnetization exponent βeff vs. and Critical Phenomena, edited by C. Domb and J.L. log(t),with t=Tci−T,for thethermal dilution case (white) Lebowitz (Academic, NewYork, 1983), Vol. 7. and the random dilution (p=0.5) case (black). The realiza- [2] A.B. Harris, J. Phys.C 7, 1671 (1974). tions considered correspond to L=80. [3] T.C.LubenskyandA.B.Harris,inMagnetism andMag- netic Materials, edited by C. D. Graham, G. H. Lander InFig. 6theeffectivecriticalexponentforthethermal andJ.J.Rhyne,AIPConf.Proc. No.24(AIP,NewYork, caseseemstoproduceacrossovertowardsβ =0.35asthe 1975), p.99; A. B. Harris and T. C. Lubensky,Phys. Rev. temperature approaches the critical point, before finite Lett. 33, 1540 (1974). 0siz.e0effectsfinallyappear. Inprinciplethismightbeinflu- [4] T. C. Lubensky,Phys.Rev. B11, 3573 (1975). [5] G.GrinsteinandA.Luther,Phys.Rev.B13,1329(1976). enced by the indetermination in the measurement of the part-icu2lar.c5ritical tem-pe2ratu.r0e of the re-ali1zati.on5. How- -[61] A..A0harony,in-Ph0ase.T5ransitionsan0dC.rit0icalPhenomena, editedbyC.DombandM.S.Green(Academic,NewYork, ever, this crossover may point out that different length 1976), Vol. 6, p. 357. scales might be important at different distance from the [7] J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, critical point. In principle characteristic lengths suchlaos g P(hyts.)Rev. Lett.57, 2999 (1986). 6 [8] W. Kinzel and E. Domany, Phys. Rev.B 23, 3421 (1981); M. Gehring and C. F. Majkrzak, Phys. Rev. B 49, 11967 D. Andelman and A. N. Berker, Phys. Rev. B 29, 2630 (1994). (1984). [15] M.Altarelli,M.D.Nuń˜ez-RegueiroandM.Papoular,Phys. [9] H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Rev.Lett. 74, 3840 (1995). MunõzSudupe,G.Parisi andJ.J.Ruiz-Lorenzo,J.Phys. [16] A.WeinribandB.I.Halperin,Phys.Rev.B27,413(1983). A 30, 8379 (1997). [17] H.G. Ballesteros and G. Parisi, cond-mat/9903230. [10] H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. [18] M. I. Marqués and J. A. Gonzalo, Phys. Rev. E 60, 2394 Munõz Sudupe, G. Parisi and J. J. Ruiz-Lorenzo, Phys. (1999). Rev.B 58, 2740 (1998). [19] W. Klein, H. E. Stanley, P.J. Reynolds and A. Coniglio, [11] H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Phys.Rev. Lett.41, 1145 (1978). Munõz Sudupe, G. Parisi and J. J. Ruiz-Lorenzo, Nucl. [20] A. Aharony and A. B. Harris, Phys. Rev. Lett. 77, 3700 Phys. B 512[FS], 681 (1998). (1996). [12] R. Folk, Yu. Holovatch and T. Yavors’kii, cond- [21] S. Wiseman and E. Domany, Phys. Rev. Lett. 81, 22 mat/9909121. (1998). [13] T. R. Thurston, G. Helgesen, D. Gibbs, J. P. Hill, B. D. [22] A.Aharony,A.B.HarrisandS.Wiseman,Phys.Rev.Lett Gaulin and G. Shirane, Phys. Rev. Lett. 70, 3151 (1993); 81, 252 (1998). T. R. Thurson, G. Helgesen, J. P. Hill, D. Gibbs, B. D. [23] U.Wolff, Phys. Rev.Lett. 62, 361 (1989). Gaulin andP. J.Simpson, Phys.Rev.B 49, 15730 (1994). [24] S.WangandR.H.SwendsenPhysica(Amsterdam)167A, [14] P.M.Gehring,K.Hirota,C.F.MajkrzakandG.Shirane, 565 (1990). Phys.Rev.Lett.71,1087(1993);K.Hirota,G.Shirane,P. 7