Thermodynamics of a two-dimensional frustrated spin-1/2 Heisenberg ferromagnet
aa r X i v : . [ c ond - m a t . s t r- e l ] M a y Thermodynamics of a two-dimensional frustrated spin- Heisenberg ferromagnet
M. H¨artel and J. Richter
Institut f¨ur Theoretische Physik, Otto-von-Guericke-Universit¨at Magdeburg, D-39016 Magdeburg, Germany
D. Ihle
Institut f¨ur Theoretische Physik, Universit¨at Leipzig, D-04109 Leipzig, Germany
S.-L. Drechsler
Leibniz-Institut f¨ur Festk¨orper- und Werkstoffforschung Dresden, D-01171 Dresden, Germany (Dated: November 12, 2018)Using the spin-rotation-invariant Green’s function method we calculate the thermodynamic quan-tities (correlation functions h S S R i , uniform static spin susceptibility χ , correlation length ξ , andspecific heat C V ) of the two-dimensional spin-1 / J - J Heisenberg ferromagnet for J < J c ≈ . | J | , where J c is the critical frustrating antiferromagnetic next-nearest neighbor coupling atwhich the ferromagnetic ground state gives way for a ground-state phase with zero magnetization.Examining the low-temperature behavior of χ and ξ , in the limit T → χ ∝ exp( b/T ) and ξ ∝ exp( b/ T ), respectively. We find a linear decrease inthe coefficient b with increasing frustration according to b = − π ( J + 2 J ), i.e., the exponentialdivergence of χ and ξ is present up to J c . Furthermore, we find an additional low-temperaturemaximum in the specific heat when approaching the critical point, J → J c . PACS numbers:
I. INTRODUCTION
Low-dimensional quantum magnets have attractedmuch attention during the last decades.
They are pre-destined to study the influence of strong thermal andquantum fluctuations. Much attention has been paidto the theoretical investigation of the two-dimensional(2D) spin-1 / J - J quantum Heisenberg antiferromag-net which may serve as a canonical model to studythe interplay of frustration effects and quantum fluc-tuations (see, e.g., Ref. 3). An additional motivationto study this model comes from the experimental side. Very recently, several quasi-2D magnetic materials witha ferromagnetic nearest-neighbor (NN) coupling J < J > VO(PO ) , (CuCl)LaNb O , SrZnVO(PO ) , and BaCdVO(PO ) . Thequite large frustrating J drives these materials out ofthe ferromagnetic phase. The experimental findings havestimulated several theoretical studies of the ground-stateand thermodynamic properties of the J - J model with J < J > It was found thatthe ferromagnetic ground state for the spin-1 / J = J c ≈ . | J | . Note thatfor the classical model (spin s → ∞ ) the correspondingtransition point is at J = 0 . | J | .On the other hand, some materials considered as2D spin-1/2 ferromagnets, such as K CuF , Cs CuF ,Cs AgF , La BaCuO , and Rb CrCl , might havealso a weak frustrating NNN interaction J < J c . Incontrast to the previous investigations of the 2D J - J model with ferrromagnetic J which have consid-ered predominantly the case of strong frustration, in the present paper we will focus on the region of weak frus-tration J < J c .Although for J < J c the ground state remains fer-romagnetic, the frustrating J may influence the ther-modynamics substantially, in particular near a zero-temperature transition. This has been demonstratedfor the one-dimensional (1D) frustrated ferromagnet inRef. 28, where a change in the low-temperature behav-ior of the susceptibility and the correlation length as wellas an additional low-temperature maximum in the spe-cific heat have been found when approaching the zero-temperature critical point from the ferromagnetic side.The Hamiltonian of the system considered in this paperis given by H = J X h i,j i S i S j + J X [ i,j ] S i S j ; J < J ≥ , (1)where ( S i ) = 3 /
4, and h i, j i denotes NN and [ i, j ]denotes NNN bonds. The unfrustrated 2D ferromag-net ( J = 0) has been widely investigated, e.g., bythe modified spin-wave theory, renormalization groupapproaches, the quantum Monte Carlo method, and by a spin-rotation-invariant second-order Green’sfunction method (RGM). However, in presence offrustration ( J >
0) the choice of appropriate methodsfor studying temperature dependent quantities of thissystem is restricted due to frustration and dimension-ality. For instance, the powerful quantum Monte Carlomethod is not applicable due to the minus sign problem,whereas finite-temperature density matrix renormaliza-tion group studies are restricted to 1D systems.So far, only a few theoretical papers deal with the ther-modynamics of the model for J >
0, see, e.g., Refs. 13,16and 17, where finite lattices of N = 16 and N = 20sites are considered. As we will see below, for these lat-tice sizes the finite-size effects at low temperatures arelarge. Therefore, methods studying infinite systems arehighly desirable. Among others, the Green’s functionmethod is a powerful tool to study magnetic systemsat finite temperatures, see e.g. Refs. 40–42 and refer-ences therein. A specific variant appropriate for frus-tated systems is the RGM. In the present paper, we usethe RGM to investigate thermodynamic properties of thefrustrated model (1), focusing on J ≤ J c ≈ . | J | . TheRGM has been applied successfully to low-dimensional(frustrated) infinite quantum spin systems. In particular, for the 1D spin-1 / Since the RGM can be formulated also forfinite systems, we use full exact diagonalization (ED) offinite square lattices of N = 16 and N = 20 sites tocompare the RGM with ED data for finite lattices.The paper is organized as follows: In Sec. II the RGMapplied to model (1) is presented. The results for thethermodynamic quantities are discussed in Sec. III. Fi-nally, a summary is given in Sec. IV. II. ROTATION-INVARIANT GREEN’SFUNCTION METHOD
To calculate the thermodynamic quantities of themodel (1), we use the RGM for determin-ing the transverse dynamic spin susceptibility χ + − q ( ω ) = −hh S + q ; S −− q ii ω , where hh . . . ; . . . ii ω denotes the two-timecommutator Green’s function. Taking the equations ofmotion up to the second step and supposing spin rota-tional symmetry, i.e., h S zi i = 0, we obtain ω hh S + q ; S −− q ii ω = M q + hh− ¨ S + q ; S −− q ii ω , (2) M q = h (cid:2)(cid:2) S + q , H (cid:3) , S − q (cid:3) i , − ¨ S + q = (cid:2)(cid:2) S + q , H (cid:3) , H (cid:3) . (3)The moment M q is given by the exact expression M q = − X k =1 , J k C ,k − (cid:16) − γ ( k ) q (cid:17) , (4)where C n,m ≡ C R = h S +0 S − R i = 2 h S z S z R i , R = n e x + m e y , γ (1) q = (cos q x + cos q y ) /
2, and γ (2) q = cos q x cos q y .The second derivative − ¨ S + q is approximated in the spiritof the schemes employed in Refs. 28,35–39,43–45. Thatis, in − ¨ S + i we adopt the decoupling S + i S + j S − k = α i,k h S + i S − k i S + j + α j,k h S + j S − k i S + i . (5)Following the arguments of Ref. 28 we put α i,k = α inthe whole temperature region. We obtain − ¨ S + q = ω q S + q and χ + − q ( ω ) = −hh S + q ; S −− q ii ω = M q ω q − ω , (6) with ω q = 2 X k,l (=1 , J k J l (cid:16) − γ ( k ) q (cid:17) × h K k,l + 8 αC ,k − (cid:16) − γ ( l ) q (cid:17)i , (7)where K , = 1 + 2 α (2 C , + C , − C , ), K , = 1 +2 α (2 C , + C , − C , ), K , = 4 α ( C , − C , ), and K , = 4 α ( C , + C , − C , ). From the Green’s func-tion (6) the correlation functions C R = N P q C q e i qR of arbitrary range R are determined by the spectraltheorem, C q = h S + q S − q i = M q ω q [1 + 2 n ( ω q )] , (8)where n ( ω q ) = (cid:0) e ω q /T − (cid:1) − is the Bose function. Tak-ing the on-site correlator C R = and using the oper-ator identity S + i S − i = + S zi we get the sum rule C = N P q C q = . The uniform static spin suscep-tibility χ = lim q → χ q , where χ q = χ q ( ω = 0) and χ q ( ω ) = χ + − q ( ω ), is given by χ = − X k =1 , kJ k C ,k − , ∆ = X k,l =1 , kJ k J l K k,l . (9)The correlation length may be calculated by expanding χ q around q = 0, χ q = χ/ (cid:2) ξ (cid:0) q x + q y (cid:1)(cid:3) . Wefind ξ = 2 α ( J C , + 2 J J ( C , + C , ) + 4 J C , )∆ . (10)The ferromagnetic long-range order, occurring at T = 0only, is reflected by a non-vanishing quantity C (see, e.g.,Ref. 36) according to C R (0) = N P q =0 C q (0) e i qR + C ,where C , typically called condensation term, is con-nected with the magnetization m by m = 3 C/
2. Dueto the exact result C R = (0) = , C q (0) = M q (0)2 ω q (0) mustbe independent of q . Examining Eqs. (4) and (7), thisrequires K k,l (0) = 0 which leads to α (0) = and C q (0) = . Hence, in case of a ferromagnetic groundstate we find at zero temperature ω q = 2 ρ s q ( | q | << ρ s = − ( J + 2 J ) (where ρ s is the spin-stiffness) and C R (0) = δ R , + . Note that the sum rule C = isfulfilled. In Eqs. (9) and (10), we have ∆(0) = 0 so that χ and ξ diverge as T → J c (see Sec. III A), we determine the ground-state corre-lation functions and the ground-state energy per site, E = 3( J C , + J C , ), for J > J c . Thereby, we takeinto consideration the possible existence of a quantumcollinear phase described by the magnetic ordering vector Q = ( π,
0) or (0 , π ) and proceed along the lines indicatedin Ref. 44.Since we will compare the RGM data with ED data forfinite lattices, see Sec. III, we have to adopt the RGM tofinite N . In this case the quantity C is not related to themagnetization and stays nonzero in the whole tempera-ture region. In Ref. 45 it was shown that C = 2 T χ/N .To evaluate the thermodynamic quantities, a coupledsystem of six (for finite systems: seven) nonlinear alge-braic self-consistency equations, including the sum rule C = , has to be solved numerically to determine thecorrelators C , , C , , C , , C , , C , , ( C ) and the ver-tex parameter α . To solve this system we use Broyden’smethod, which yields the solutions with a relative errorof about 10 − on the average. The momentum integralsare done by Gaussian integration. To find the numericalsolution of the RGM equations for T >
0, we start athigh temperatures and decrease T in small steps. Belowa certain (low) temperature T ( J ) no solutions of theRGM equations (except at T = 0) could be found, sincethe quantity ∆( T, J ) in Eqs. (9) and (10) becomes ex-ponentially small which leads to numerical instabilities.Presenting our results in the next section we put J = − III. RESULTSA. Phase transition at T = 0 In a first step we use the RGM as described above todetermine the transition point J c , where the ferromag-netic ground state gives way for a ground-state phasewith zero magnetization. For that we solve the RGM −0.3−0.2−0.1 0 0.1 0.2 0.3 0 0.5 1 1.5 2 2.5 3 < S S R > J R=1,0R=2,0R=1,1R=1,2R=2,2−0.35−0.3−0.25 0.4 0.45 0.5J E FIG. 1: Ground-state spin-spin correlation functions h S S R i in dependence on frustration J . (Note that for J < J c ≈ .
44 all h S S R i coincide.) Inset: Ground-state energy E ver-sus J (solid line: RGM data, dashed line: extrapolated RGMdata to fix the transition point, see also text). equations at T = 0 (i) starting from J = 0 with increas-ing J and (ii) starting from large J ≫ J . In case (i) the RGM equations at T = 0 can be solveduntil the classical transition at J = 0 .
5. In case (ii) we find solutions of the RGM ground-state equations downto J = 0 .
46. This value lies certainly above the transi-tion point J c . To fix J c we may extrapolate the ground-state energy obtained for case (ii). The crossing pointof both energies at J = 0 .
44 can be considered as theRGM estimate of J c , see the inset of Fig. 1. This valueis in reasonable agreement with values for J c reported inother papers. The behavior of the spin-spin correlation functions at T = 0 near the transition point J c is illustrated in Fig. 1.For J & J c there is a noticeable variation in the cor-relation functions. With increasing J , at J &
1, thespin-spin correlation functions approach the correspond-ing values of the limiting case J ≫
1. This behavioris in qualitative agreement with data from Lanczos EDand coupled-cluster method. At the critical point thecorrelation functions jump to the exact results of the fer-romagnetic phase, i.e., h S S R i = 0 .
25. Together withthe kink in the energy, this indicates a first-order tran-sition. These results confirm the findings in previouspapers.
B. Thermodynamic properties
Now we investigate the thermodynamic properties ofthe model for J < J c ≈ .
44, where the ground stateis ferromagnetic. First we consider the temperature de- −0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0 0.5 1 1.5 2 < S S , ( , ) > T NNNNN
FIG. 2: NN and NNN spin-spin correlation functions for J =0, 0 .
15, and 0 .
3, from top to bottom, calculated by RGM for N → ∞ (lines) and ED for N = 20 (filled symbols). pendence of the spin-spin correlation functions. The NNand NNN correlators are shown in Fig. 2 for variousvalues of J . For comparison we show also ED data for N = 20 sites. The RGM results agree qualitatively withthe ED data. With increasing frustration, the correla-tion functions decrease more rapidly with temperature.Interestingly, for large frustration the NNN correlationfunction changes its sign at a certain temperature, e.g.,for J = 0 .
3, at T ≈ .
25 (RGM data). The faster decayof the spin-spin correlators due to frustration is relatedto the decrease in the correlation length ξ with increasing J , see the discussion below. c T 0 40 80 120 0 0.1 0.2 0.3 0.4 0.5 c J T=0.3
FIG. 3: Uniform static spin susceptibility χ for J =0 , . , . , . , .
4, from right to left. The inset shows thesusceptibility for T = 0 . J . The solid line represents an exponential fit, χ ( T = 0 . ≈ . − . J ), of the data points (filled squares). x T 0 4 8 12 0 0.1 0.2 0.3 0.4 0.5 x J T=0.3
FIG. 4: Correlation length ξ for J = 0 , . , . , . , . T = 0 . J . The solid linerepresents an exponential fit, ξ ( T = 0 . ≈ . − . J ),of the data points (filled squares). The uniform static spin susceptibility χ and the cor-relation length ξ depicted in Figs. 3 and 4 show a sim-ilar behavior. Due to the ferromagnetic ground state,both quantities diverge at T = 0 exponentially, see be-low. With increasing frustration J the rapid increasein both quantities is shifted to lower temperatures. Asshown in the insets of Figs. 3 and 4, at a certain fixed temperature both χ and ξ decrease rapidly with J . c T (a) 0 10 20 30 40 50 0 0.2 0.4 0.6 0.8 1 c T (b)
FIG. 5: Uniform static spin susceptibility χ for (a) J = 0 . J = 0 . N = 16 , , , , thermodynamic limit(lines, RGM) and N = 16 ,
20 (filled symbols, ED), from leftto right. Note that for J = 0 . N = 400 and N → ∞ almost coincide. Since the susceptibility (together with the correlationlength) is an important quantity to analyze the criticalproperties for T →
0, we first test the quality of ourRGM approach in more detail by comparing the RGMdata for N = 16 and N = 20 with ED data for the samesystem sizes, see Fig. 5. It is obvious, that the RGMand the ED data for χ are in excellent agreement. Thisis consistent with our previous finding that for the 1DHeisenberg ferromagnet the RGM data for χ ( T → ξ ( T →
0) are in perfect agreement with exact Betheansatz results. Moreover, the RGM calculations of χ forfinite N allow to estimate the magnitude of the finite-sizeeffects. As can be clearly seen in Fig. 5, the finite-sizeeffects at T . . T . .
4) for J = 0 . J = 0 .
3) be-come very large. Since the correlation length ξ becomessmaller with increasing J (see the inset in Fig. 4), thefinite-size effects become less pronounced with growingfrustration. Nevertheless, based on our data we arguethat the finite systems of N = 16 and N = 20, accessibleby exact full diagonalization, are not representative forthe thermodynamic limit, see also Refs. 16 and 17. Notethat for the 1D case, the ED data for N = 20 are ingood agreement with RGM data for N → ∞ . Consis-tent to that observation, in two dimensions the system of N = 400 sites with the same linear extension is alreadyclose to the thermodynamic limit, see Fig. 5.Next we use the RGM for N → ∞ to investigatethe critical behavior of χ and ξ for T → J = 0 exact Betheresults are available, we have only approximate re-sults for the 2D model. From low-temperature expan-sions of the susceptibility and the correlation length forthe model with J = 0 using renormalization groupapproaches, the modified spin-wave threory, andthe RGM it is known that for T → χ ∝ T s exp ( b/T ) and the correlation lengthas ξ ∝ T σ exp ( β/T ). While different leading exponents s and σ for the (less important) preexponential factorwere obtained by different methods, the exponential di-vergence is obtained by all authors. However, dif-ferent values for the coefficients b and β were reported,namely, b ( J = 0) = π/ b ( J = 0) = π us-ing modified spin-wave theory or renormalization groupapproach, and b ( J = 0) = 0 . π using a differentversion of the renormalization-group approach. In allthese papers it was found that the coefficients inthe exponents fulfill the relation β = b/
2. We mentionfurther, that early numerical studies based on quantumMonte-Carlo calculations (using, however, χ and ξ dataonly for quite large temperatures) give b ( J = 0) ≈ . ≈ . π (Ref. 33) and β ( J = 0) = 0 . π (Ref. 32). Letus finally argue, that the results for the coefficient b ob-tained by modified spin-wave theory as well as renor-malization group approach and by the RGM seemto be most reliable, since these methods are well-tested.Moreover, for the 1D spin-1 / as well as the modifiedspin-wave theory reproduce the exact Bethe-ansatz re-sults for the low-temperature behavior of χ and ξ . Hence,it is to some extent surprising, that there is a differenceby a factor of 2 in the coefficient b between the value ob-tained by modified spin-wave theory (as well as renor-malization group approach ) and the RGM for the 2Dunfrustrated ferromagnet. This discrepancy is a knownbut unresolved problem. However, we believe thatour general conclusions concerning the exponential diver-gence of the frustrated model, see the discussion below,are not affected by this problem.Next we use the RGM results to determine the co-efficients b and β as functions of J . We assume thatthe general low- T behavior of these quantities is pre-served for 0 < J < J c , cf. Ref. 28, and fit our nu-merical RGM data for χ and ξ at low temperatures us-ing the ansatzes χ = (cid:0) a T − + a + a T (cid:1) exp (cid:0) bT (cid:1) and ξ = p ( α T − + α + α T ) exp (cid:16) βT (cid:17) . Note that the lead-ing power T − in the preexponential functions was de-rived for the unfrustated case with the RGM in Ref. 37.We consider b , a , a , and a as well as β , α , α , and α as independent fit parameters. Recall that we can calcu-late RGM data only for temperatures down to a certain T , where T (e.g., T = 0 .
161 for J = 0 and T = 0 . J = 0 .
4) is the lowest temperature where the sys-tem of RGM equations converges (see Sec. II). Thus,in addition to the leading power T − in the preexpo-nential function it is reasonable to consider higher orderterms to achieve an optimal fit of the RGM data. Forthe fit we use 500 equidistant data points in the interval T . . . T + T cut , where T cut is set to 0 .
05. The fit of thenumerical RGM data reproduces the analytical results ofRef. 37 for J = 0, i.e., b ( J = 0) = 2 β ( J = 0) = π/ J = 0, we nowconsider the frustrated model, where (to the best of ourknowledge) no other results are available. From our nu-merical data for χ and ξ we determine the J dependenceof the coefficients b and β . We find that the numericaldata for b and β obtained by the fitting procedure de-scribed above are very well described by a linear decreasein both parameters with increasing frustration, b = 2 β = − π J + 2 J ) . (11)Obviously, both parameters would be zero at the classi-cal transition point J = 0 .
5, but they are still finite atthe transition point J c ≈ .
44 of the quantum model.Hence, the exponential divergence is present in the fullparameter range J ≤ J c where the ground state is fer-romagnetic. We emphasize that this result is contraryto the behavior observed for the 1D frustrated spin-1 / We mention furtherthat the leading coefficients a and α of the preexponen-tial functions, see the expressions for χ and ξ given above,vanish at the transition point J c , whereas the next coeffi-cients a and α remain finite at J c . Hence, the preexpo-nential temperature dependence is changed approachingthe zero-temperature transition.Let us mention, that the linear decrease in the coef-ficients b and β , see Eq. (11), found by fitting the low-temperature behavior of χ and ξ is the same as that ob-tained analytically for the zero-temperature spin-stiffness ρ s , see Sec. II. This relation between ρ s and the diver-gence of the correlation length and the susceptibility is inaccordance with general arguments concerning thelow-temperature physics of low-dimensional Heisenbergferromagnets.Another interesting quantity is the specific heat C V shown in Fig. 6. For J = 0 the specific heat ex-hibits a typical broad maximum at about T = 0 . C V T 0.1 0.15 0 0.4 0.8T C V FIG. 6: Specific heat calculated by RGM for J =0 , . , . , . , .
35, and 0 .
4, from top to bottom. The insetshows the specific heat for J = 0 . .
35 in steps of 0 . C V ( T ) curve changes for large frustration. For J & . T < . T > . J c . Then many low-lyingmultiplets appear above the fully polarized ferromagneticground-state multiplet. Hence, the appearance of theadditional low-temperature maximum in C V ( T ) can beattributed to a subtle interplay between all of these low-lying states. Note that a similar behavior was found forthe 1D case. IV. SUMMARY
In this paper we investigated the influence of the frus-trating NNN coupling J < J c on the thermodynamicproperties of the 2D spin-1/2 Heisenberg ferromagnet us-ing the spin-rotation-invariant Green’s function method(RGM). The RGM estimate for the critical frustration J c , where the ferromagnetic ground state breaks down,is J c ≈ . | J | .We tested the method by comparing RGM data forthe spin-spin correlation functions and the uniform sus-ceptibility χ calculated for finite lattices of N = 16 and N = 20 sites with corresponding exact-diagonalizationdata for the same lattice sizes and found a good agree-ment between both methods. However, the comparisonof RGM data for finite lattices and for N → ∞ indicatesstrong finite-size effects at low temperatures. This leadsto the conclusion that the finite systems of N = 16 and N = 20 are not representative for the low-temperaturethermodynamics of large systems.As it is known from the Mermin-Wagner theorem the thermal fluctuations are strong enough to suppressmagnetic LRO for the Heisenberg ferromagnet in dimen-sion D < χ and the correlation length ξ in the ferromagnetic ground-state region J < J c exhibits the exponential divergences χ ∝ exp (cid:16) π ( | J |− J )2 T (cid:17) and ξ ∝ exp (cid:16) π ( | J |− J )4 T (cid:17) . Al-though the numerator in the exponent decreases withgrowing J the exponential divergence perpetuates in thewhole parameter region J ≤ J c .The specific heat calculated within the RGM exhibitsa double-maximum structure for J & .
34 analogous tothe 1D model.To test our theoretical predictions it would be inter-esting to reconsider experimentally quasi-2D spin-1/2ferromagnets, such as K CuF , Cs CuF , Cs AgF ,La BaCuO , and Rb CrCl with respect to apossible frustration and its influence on low-temperaturethermodynamics.
Acknowledgments:
This work was supported by theDFG (project RI615/16-1 and DR269/3-1). For the exactdiagonalization J. Schulenburg’s spinpack was used. J.R.thanks N.B. Ivanov for useful discussions. Frustrated Spin Systems , H. T. Diep, Ed. (World Scientific,Singapore, 2004). Quantum Magnetism , U. Schollw¨ock, J. Richter,D. J. J. Farnell, R. F. Bishop, Eds. (Lecture Notesin Physics, ) (Springer, Berlin, 2004). P. Chandra and B. Doucot, Phys. Rev. B , 9335(1988); H.J. Schulz and T.A.L. Ziman, Europhys.Lett. , 355 (1992); H.J. Schulz, T.A.L. Ziman, andD. Poilblanc J. Phys. I , 675 (1996); J. Richter,N.B. Ivanov, and K. Retzlaff, Europhys. Lett. , 545(1994); A. E. Trumper, L. O. Manuel, C. J. Gazza, and H. A. Ceccatto, Phys. Rev. Lett. , 2216 (1997); L. Capri-otti, F. Becca, A. Parola, and S. Sorella, Phys. Rev. Lett. , 097201 (2001); R.R.P. Singh, Weihong Zheng , J. Oit-maa, O.P. Sushkov, and C.J. Hamer, Phys. Rev. Lett. , 017201 (2003); J. Sirker, Z. Weihong, O. P. Sushkov,and J. Oitmaa, Phys. Rev. B , 184420 (2006); R. Dar-radi, O. Derzhko, R. Zinke, J. Schulenburg, S. E. Kr¨uger,and J. Richter, Phys. Rev. B , 214415 (2008); V. Murg,F. Verstraete, and J. I. Cirac, Phys. Rev. B , 195119(2009), J. Richter and J. Schulenburg, Eur. Phys. J. B ,117 (2010). R. Melzi, P. Carretta, A. Lascialfari, M. Mambrini,M. Troyer, P. Millet, and F. Mila, Phys. Rev. Lett. ,1318 (2000); R. Melzi, S. Aldrovandi, F. Tedoldi, P. Car-retta, P. Millet, and F. Mila, Phys. Rev. B , 024409(2001); H. Rosner, R.R.P. Singh, W.H. Zheng, J. Oitmaa,and W.E. Pickett, Phys. Rev. B E. E. Kaul, H. Rosner, N. Shannon, R.V. Shpanchenko,and C. Geibel, J. Magn. Magn. Mater. , 922(2004). M. Skoulatos, J.P. Goff, N. Shannon, E.E. Kaul, C. Geibel,A.P. Murani, M. Enderle, and A.R. Wildes J. Magn. Magn.Mater. , 1257 (2007). P. Carretta, M. Filibian, R. Nath, C. Geibel, and P. J. C.King, Phys. Rev. B , 224432 (2009). M. Skoulatos, J.P. Goff, C. Geibel, E.E. Kaul, R. Nath, N.Shannon, B. Schmidt, A.P. Murani, P.P. Deen, M. Enderle,and A.R. Wildes, Europhys. Lett. , 57005 (2009). H. Kageyama, T. Kitano, N. Oba, M. Nishi, S. Nagai,K. Hirota, L. Viciu, J.B. Wiley, J. Yasuda, Y. Baba,Y. Ajiro, and K. Yoshimura, J. Phys. Soc. Jpn. , 1702(2005). A.A. Tsirlin and H. Rosner, Phys. Rev. B A.A. Tsirlin, B. Schmidt, Y. Skourski, R. Nath, C. Geibel,and H. Rosner, Phys. Rev. B R. Nath, A.A. Tsirlin, H. Rosner, and C. Geibel, Phys.Rev. B N. Shannon, B. Schmidt, K. Penc, and P. Thalmeier, Eur.Phys. J. B , 599 (2004). N. Shannon, T. Momoi, and P. Sindzingre, Phys. Rev. Lett. , 027213 (2006). P. Sindzingre, N. Shannon and T. Momoi, J. Magn. Magn.Mat. , 1340 (2007). B. Schmidt, N. Shannon, and P. Thalmeier, J. Phys. Cond.Mat. , 145211 (2007). B. Schmidt, N. Shannon, and P. Thalmeier, J. Magn.Magn. Mater. , 1231 (2007). J.R. Viana and J.R. de Sousa, Phys. Rev. B , 052403(2007). P. Sindzingre, L. Seabra, N. Shannon, and T. Momoi, J.Phys.: Conf. Series , 012048 (2009). R. Shindou and T. Momoi, Phys. Rev. B , 064410(2009). J. Richter, R. Darradi, J. Schulenburg, D.J.J. Farnell, andH. Rosner, arXiv:1002.2299 (2010). D. V. Dmitriev, V. Ya. Krivnov, and A. A. Ovchinnikov,Phys. Rev. B , 3620 (1997). S. Feldkemper, W. Weber, J. Schulenburg, and J. Richter,Phys. Rev. B , 313 (1995). S. Feldkemper and W. Weber, Phys. Rev. B , 7755(1998). H. Manaka, T. Koide, T. Shidara, and I. Yamada, Phys.Rev. B , 184412 (2003). S.E.McLain, D.A.Tennant, J.F.C.Turner, T. Barnes, M.R.Dolgos, Th. Proffen , B.C. Sales, and R.I Bewley,arXiv:cond-mat/0509194. D. Kasinathan, A.B. Kyker, and D.J. Singh Phys. Rev. B , 214420 (2006). M. H¨artel, J. Richter, D. Ihle, and S.-L. Drechsler, Phys.Rev. B , 174412 (2008); J. Richter, M. H¨artel, D.Ihle, and S.-L. Drechsler, J.Phys.:Conf.Series , 012064(2009). M. Takahashi, Prog. Theo. Phys. Supp. , 233 (1986);M. Takahashi, Phys. Rev. Lett.
168 (1987). P. Kopietz and S. Chakravarty, Phys. Rev. B N. Karchev, Phys. Rev. B E. Manousakis and R. Salvador, Phys. Rev. B , 575(1989). Y.C. Chen, H.H. Chen, and F. Lee, Phys. Rev. B , 11082(1991). P. Henelius, A.W. Sandvik, C. Timm, and S. M. GirvinPhys. Rev. B , 364 (2000). I. Juh´asz Junger, D. Ihle, L. Bogacz, and W. Janke, Phys.Rev. B , 174411 (2008). J. Kondo and K. Yamaji, Prog. Theor. Phys. , 807(1972); H. Shimahara and S. Takada, J. Phys. Soc. Jpn. , 2394 (1991); S. Winterfeldt and D. Ihle, Phys. Rev. B , 5535 (1997). F. Suzuki, N. Shibata, and C. Ishii, J. Phys. Soc. Jpn. ,1539 (1994). I. Junger, D. Ihle, J. Richter, and A. Kl¨umper, Phys. Rev.B , 104419 (2004). T.N. Antsygina, M.I. Poltavskaya, I.I. Poltavsky, andK.A. Chishko, Phys. Rev. B , 024407 (2008). W. Gasser, E. Heiner and K. Elk,
Greensche Funktionenin Festk¨orper- und Vielteilchenphysik (Wiley, Berlin 2001). Yu. A. Izyumov, N.I. Chaschin, and V. Yu. Yushankhai,Phys. Rev. B , 214425 (2002). P. Froebrich and P.J. Kuntz, Physics Reports , 223(2006). W. Yu and S. Feng, Eur. Phys. J. B , 265 (2000); B.H.Bernhard, B. Canals, and C. Lacroix, Phys. Rev. B ,104424 (2002); D. Schmalfuß, J. Richter, and D. Ihle, Phys.Rev. B , 184412 (2004); D. Schmalfuß, J. Richter, andD. Ihle, Phys. Rev. B , 224405 (2005). L. Siurakshina, D. Ihle, and R. Hayn, Phys. Rev. B ,104406 (2001); D. Schmalfuß, R. Darradi, J. Richter,J. Schulenburg, and D. Ihle, Phys. Rev. Lett. , 157201(2006). I. Juh´asz Junger, D. Ihle, and J. Richter, Phys. Rev. B ,064454 (2005). W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flan-nery,
Numerical Recipes in C++, The Art of Scien-tific Computing (Cambridge University Press, Cambridge,2007). M. Yamada and M. Takahashi, J. Phys. Soc. Jpn. , 2024(1986); M. Yamada, J. Phys. Soc. Jpn. , 848 (1990). P. Kopietz and G. Castilla, Phys. Rev. B N.B. Ivanov, Condens. Matter Phys. (L’viv) , 435 (2009)(see also arXiv:0909.2182). N. Mermin and H. Wagner, Phys. Rev. Lett.17