Ising model on a square lattice with second-neighbor and third-neighbor interactions
F.A. Kassan-Ogly, A.K. Murtazaev, A.K. Zhuravlev, M.K. Ramazanov, A.I. Proshkin
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J u l Ising model on a square lattice with second-neighbor and third-neighbor interactions
F.A. Kassan-Ogly a , A.K. Murtazaev b , A.K. Zhuravlev a , M.K. Ramazanov b , A.I. Proshkin a a Institute of Metal Physics, Ural Division, Russian Academy of Sciences, ul. S.Kovalevskoi 18, Ekaterinburg, 620990, Russia b Institute of Physics, DSC, Russian Academy of Sciences, ul. Yaragskogo 94, Makhachkala, Daghestan, 367003, Russia
Abstract
We studied the phase transitions and magnetic properties of the Ising model on a square lattice by the replica MonteCarlo method and by the method of transfer-matrix, the maximum eigenvalue of which was found by Lanczos method.The competing exchange interactions between nearest neighbors J , second J , third neighbors J and an externalmagnetic field were taken into account. We found the frustration points and expressions for the frustration fields, atcrossing of which cardinal changes of magnetic structures (translational invariance changes discontinuously) take place.A comparative analysis with 1D Ising model was performed and it was shown that the behavior of magnetic propertiesof the 1D model and the 2D model with J and J interactions reveals detailed similarity only distinguishing in scalesof magnetic field and temperature. Keywords:
Ising model, square lattice, competing interactions
PACS:
1. Introduction
It has been over 70 years since the publication of On-sager celebrated solution [1] of the Ising model on a squarelattice with nearest-neighbor (nn) interactions. This so-lution has served as a cornerstone for modern theoriesand as a testing ground for many approximate theoreti-cal approaches. Since that time, there have been very fewother systems for which exact solutions have been found.Among them the exact solutions of the Ising model onother 2D lattices: a triangular [2], a honeycomb [3], anda kagome lattice [4] have been found, but again only withnearest-neighbor (nn) interactions. The slight alterationsof the original model such as the addition of next-nearest-neighbor (nnn) interactions or an external magnetic fieldare no longer exactly soluble by presently available the-oretical approaches. Naturally without an exact solutiondifferent approximate methods such as mean-field approx-imations, series expansions [5–8], calculating the inter-face free energy [9], real-space renormalization-group tech-niques [10–13], finite-size scaling of the transfer matrix[14, 15], the Fisher zeros of the partition function [16, 17],and various kinds of Monte Carlo simulations [18–27] weredeveloped up to recent years.The overwhelming majority of papers on a square latticein the Ising model are devoted to multifarious topical prob-lems such as determining and elucidating the phase dia-grams of magnetic structures, the order and the universal-ity class of the phase transitions; finding the multicritical
Email address: [email protected] (F.A. Kassan-Ogly) points and critical lines that separate the different phases;determining the frustration lines and points that result inthe appearance of Quantum Phase Transitions; calculat-ing the critical exponents, etc. Nevertheless, certain issueshave remained controversial either due to different theoret-ical approaches or, more often, due to insufficient sizes ofa lattice subjected to numerical calculations.However, despite the significant amount of effort madein many years, a number of related problems still remainessentially untouched, in particular, the behavior of mag-netization, as a function of an external magnetic field and(or) the temperature.The aim of the present paper is just the study of mag-netic properties of the Ising model on a square lattice incomparison with those on a one-dimensional (linear chain)lattice.The Ising model on a N = L × L square lattice (or on alinear chain with L sites) is described by the Hamiltonian: H = J X nn S i S j + J X nnn S i S k + J X nnnn S i S l − H X i S i , (1)where all the exchange interactions J , J and J are posi-tive (antiferromagnetic); nn, nnn, and nnnn stand, respec-tively, for all next-neighbor pairs, second-neighbor pairs,and third-neighbor pairs; S i = ± H is an external mag-netic field. Here, all fields, exchange interactions, and tem-peratures are given in units of J , unless otherwise stated.We have kept J = 1 in all the calculations reported inthis work. The exchange interactions on a square latticeare shown in Fig. 1.On a square lattice, we will numerically calculate themagnetization, entropy, and heat capacity as functions of Preprint submitted to Journal of Magnetism and Magnetic Materials March 12, 2018 J J J J J Figure 1: The nearest-, second-, and third-neighbor interactions ona square lattice in the Ising model. the temperature and an external magnetic field accordingto the conventional formulas: F = − TN ln Z M = − ∂F∂H (2) S = − ∂F∂T C = − T ∂ F∂T , (3)where Z is the partition function, F is the free energy, and N is the number of sites on a square lattice. On a one-dimensional (linear chain) lattice, we will perform calcula-tions using the exact analytical solution for the maximumeigenvalue λ max of the Kramers-Wannier transfer-matrix,obtained in [28]. In this case, the magnetization and en-tropy are expressed solely in terms of λ max as follows: M = Tλ max ∂λ max ∂H , (4) S = ln λ max + Tλ max ∂λ max ∂T , (5) C = 2 Tλ max ∂λ max ∂T + T ∂∂T λ max ∂λ max ∂T . (6)The maximum transfer-matrix eigenvalue λ max is givenby λ max = p a − b + 4 y − a s(cid:18) a − b + 4 y (cid:19) − y − c − ya p a − b + 4 y , (7) h M Figure 2: Magnetizations of 1D chain and 2D square lattice. Isingmodel with only nearest-neighbor interactions. where y = rp Q − q r − p Q − q b , Q = p
27 + q ,p = − b ac − d, q = − b
27 + bac − bd − a d − c ,a = − (cid:18) − J − J T (cid:19) cosh HT ,b = − (cid:18) − J T (cid:19) sinh 2 J T ,c = − (cid:18) J − J T (cid:19) sinh 2 J T cosh HT ,d = − J T . (8)
2. Nearest-neighbor interaction
Let us first consider the simplest case, when J and J are equal to zero in the Hamiltonian (1). Figure 2 showsthe magnetizations as functions of rescaled magnetic field h = H/z , where z is the number of nearest neighbors on1D chain ( z = 2) or on 2D square lattice ( z = 4). The mag-netizations are calculated at T = 0 .
05 on 1D chain (solidline) and at T = 0 . H fr.1D = 2 J and H fr.2D = 4 J , and that can be expresseduniformly as H fr = zJ . It should be noted that on such ascale the magnetizations are almost indistinguishable fromone another.In Fig. 3, we show plots of the magnetization on a lin-ear chain and square lattice as a function of temperature2 T M H a L ææææææææææææææææææææææææææææææææææææææææææææàà à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à à ììììììììììììììììììììììììììììììììììììììììììì T M H b L
15 0.5459
Figure 3: Temperature dependence of magnetizations of 1D chain (a) and 2D square lattice (b) at: H = H fr (middle curves), H = 0 . H fr (lower curves), and H = 1 . H fr (upper curves). at different values of magnetic field: at frustration fields,above, and below them. It is again seen that the magneti-zation behavior is very similar. At T → H > H fr themagnetization approaches unity, at H < H fr it vanishes,and at H = H fr it tends to constant values, although dif-ferent.Our calculations of the entropy have shown that at T → H = H fr on both the linear chain and squarelattice it tends to non-zero value, which is shown in Fig. 4.The non-zero entropy value at T = 0 indicates vanishingof the transition temperature at these fields (see, for ex-ample, [29]) and gives evidence that these fields are reallythe frustration points. The calculations of heat capac-ity at and nearby the frustration fields also confirm thisassertion. The suppression of the transition temperatureon a square lattice at H fr.1D = 4 J was also corroboratedby M¨uller-Hartmann and Zittartz in [9], where they putforward a considerably simple conjecture for the relationbetween the transition temperature and an external mag-netic field cosh (cid:18) HT c (cid:19) = sinh (cid:18) JT c (cid:19) . (9)The relation (9) had been more or less successfully con-firmed by subsequent Monte-Carlo simulations [13, 20] andreal space renormalization group method [12]. ln2 ln 1 + T S Figure 4: Entropy of Ising model with only nearest-neighbor interac-tions as a function of temperature at frustration fields on 1D chain(solid line) and 2D square lattice (dashed line). h M Figure 5: Magnetizations of Ising model with nearest-neighbor andsecond-neighbor interactions on a linear chain at: R = 0 (solidcurve), R = 0 . R = 0 . R = 0 . T = 0 .
3. Nearest-neighbor and second-neighbor interac-tions
Having seen in the previous section many similarities inthe behavior of magnetic properties between 1D chain and2D square lattice let us consider the second case, when J = 0 and J = 0 in the Hamiltonian with a hope to findnew similar features. The majority of papers are particu-larly devoted to this case. However, despite of a plethora ofarticles, the calculations of the magnetization dependenceon temperature and magnetic field are virtually lacking.In Figs. 5 and 6, we show plots of the magnetizationon a linear chain and a square lattice as a function ofrescaled magnetic field h = H/z at different values of theratio between second-neighbor J and nearest-neighbor J interactions ( R = J /J = 0 , . , . . R = 1 and at T = 0 .
273 in [20] (see Fig. 7).Analyzing Figs. 5 and 6, we can notice some similaritiesin the magnetization behavior on both the 1D chain and2D square lattice, namely, the appearance of two frustra-tion fields (lower and upper), a plateau at intermediatevalues of a magnetic field, and an antiferromagnetic typeof magnetization increasing at very low fields (magneticsusceptibility vanishes at H → / / h M Figure 6: Magnetization of Ising model with nearest-neighbor andsecond-neighbor interactions on a square lattice at: R = 0 (solidcurve), R = 0 . R = 0 . R =0 . T = 0 . H M Figure 7: Magnetization of Ising model with nearest-neighbor andsecond-neighbor interactions on a square lattice at R = 1 .
0, [20]. .5 1.0 1.5 2.0 2.5 T c H Figure 8: Phase diagram for Ising square lattice with antiferromag-netic nearest-neighbor and second-neighbor interactions in a mag-netic field for R = 1 .
0, [29]. R = 1 / H → H lowfr.1D = ( zJ − zJ , at 0 < R < . zJ − z J , at R > . H uppfr.1D = zJ + zJ . (10) H lowfr.2D = ( zJ − zJ , at 0 < R < . zJ , at R > . H uppfr.2D = zJ + zJ . (11)In the Ising model with antiferromagnetic nearest- andsecond-neighbor interactions on a square lattice, the frus-tration fields (11) were for the first time determined byBinder and Landau in [20], in which they attempted to cal-culate the field dependence of the transition temperaturefor various values of the ratio R , using the mean-field ap-proximation. In the subsequent paper [29] after elaboratedMonte Carlo simulations, the authors obtained the correctfield dependence of the transition temperature. The phasediagram (or the transition curve) for the Ising square lat-tice with antiferromagnetic nearest- and second-neighborinteractions in a magnetic field for R = 1 is shown in Fig. 8.The calculations of the entropy as a function of tem-perature at the frustration point of the interactions ratio R = 1 / T → ln2 ln 1 + S S T S Figure 9: Temperature dependences of entropies in Ising model withnearest- and second-neighbor interactions on 1D chain (solid line)and square lattice (open circles) at the frustration point R = 0 . . Figure 10: The N´eel structure. a L H b L Figure 11: Superantiferromagnetic structure. the entropy on a linear chain tends to a constant valueln √ that is the logarithm of golden ratio. At the sametime, on a square lattice at J = 0 and J = 0 the en-tropy vanishes, nevertheless the transition temperature issuppressed to zero, which follows from the heat capacitycalculations.Let us now discuss the ground state magnetic structuresthat can be found from energy arguments, as was first doneby Fan and Wu [30]. When the nearest-neighbor interac-tion is strong enough J > J ( R < /
2) it determinesthe magnetic ordering. In this case, the energy per site isequal to E N´eel = − J + 2 J , and the magnetic structureis known as antiferromagnetic or the N´eel-ordered statedepicted in Fig. 10. Hereafter solid and open circles rep-resent the spin states S i = ± J > J / R > /
2) the energy per site onlydepends on J and equals E SAF = − J . The magneticstructure, shown in Fig. 11, was names by Fan and Wu[30] ”superantiferromagnetic”. In many other subsequentpapers it is referred to as SAF, ”striped” or ”layered”structure. It looks like vertical chains are ordered in anti-ferromagnetic manner (Fig. 11a) or like horizontal chainsare ordered antiferromagnetically (Fig. 11b) with keepingtranslational invariance (as if some ferromagnetic interac-tion is present) in the direction along one side of a squareand with doubled period in the other. But this impressionis erroneous, and the structure should be perceived as an-tiferromagnetic ordering (determined by strong J ) alongthe both diagonals.The N´eel and superantiferromagnetic phases are sepa-rated by the frustration point at J = J / R = 1 / Figure 12: Magnetic structure in intermediate field (between thelower and upper frustration fields) in a model with nearest-neighborand second-neighbor interactions. h M Figure 13: Magnetizations of Ising model with nearest-neighbor andthird-neighbor interactions on a square lattice at: R ′ = 0 (solidcurve), R ′ = 0 . R ′ = 0 . R ′ = 0 . T = 0 . which the transition temperature is suppressed to zero.Since at J = J / E N´eel ) andSAF ( E SAF ) structures are equal, a single translation ofany horizontal chain in Fig. 11a or any vertical chain inFig. 11b costs no energy. Thus, at the frustration pointthe ground state is composed of 4 · L configurations andhence degenerate on the order of 2 L (but not on the orderof 2 N ) so that the entropy vanishes at T = 0.The ground state structure in an external magnetic fieldbelow the lower frustration field (11) coincides either withthe N´eel (at R < / R > /
4. Nearest-neighbor and third-neighbor interac-tions
Let us now consider the magnetic properties of the Isingmodel on a square lattice in the third case, when J = 0and J = 0 in the Hamiltonian (1). In this case, our calcu-lations of the magnetization on a square lattice as a func-tion of rescaled magnetic field h = H/z at four values ofthe ratio between third-neighbor J and nearest-neighbor ln2 ln 1 + S S T S Figure 14: Temperature dependence of entropy in Ising model withnearest-neighbor and second-neighbor interactions on 1D chain (solidline) and with nearest-neighbor and third-neighbor interactions on asquare lattice (dashed line) at frustration points R = 0 . J interactions ( R ′ = J /J = 0 , . , . .
7) are plot-ted in Fig. 13. All of them are calculated at T = 0 .
1. Acomparison of Figs. 13 and 5 shows a striking similarity be-tween the magnetizations on a linear chain and a squarelattice. Keeping in mind different scales of temperatureand magnetic field the matched curves are almost indis-tinguishable. It should be emphasized that the obtainedfrustration fields on a square lattice, when J = 0 and J = 0, exactly duplicate those on a linear chain (10).Figure 14 shows the calculation results of entropy de-pendence on temperature at the frustration points R = R ′ = 1 / J = 0 and J = 0. The entropies on both a lin-ear chain and a square lattice display similar tending tononzero values at T → J = 0 and J = 0 (comparewith Fig. 9).In the case, when J = 0 and J = 0 and at strong J > J ( R ′ < /
2) the magnetic structure with energyequal to E = − J + 2 J again is the N´eel-ordered one(Fig. 10). When the third-neighbor interaction is strong J > / J ( R ′ > /
2) the energy per site only dependson J and is equal to E = − J . Figure 15 shows two ob-tained structures of different symmetry. Notwithstandingthis difference, they both have the same energy and arecomposed in a similar way, namely, as alternate sequenceof vertical and (or) horizontal chains of + + − − + + −− type. We here have the translational invariance withquadruple period along both sides of a square. Thesestructures, the N´eel and quadruple one, are separated by7 a L H b L Figure 15: Quadruple structures in Ising model with nearest-neighbor and third-neighbor interactions on a square lattice at R ′ > . the frustration point at J = 1 / J ( R ′ = 1 / T → R ′ = 1 / R ′ < / R ′ > / − + + − type along bothsides of a square, similar to the 1D lattice.
5. Conclusions
In this paper we have considered the Ising model withcompeting nearest-neighbor, second-neighbor, and third-neighbor interactions (all interactions are antiferromag-netic) on a square lattice in two alternative versions: first, J = 0, J = 0, J = 0 and second, J = 0, J = 0, J = 0. The main goal of this paper was to find and in-vestigate the magnetic properties of these two models, andto compare them to the properties of 1D chain with thenearest-neighbor and second-neighbor interactions J = 0and J = 0.Magnetic properties of both versions of the model havemuch in common. The frustration points in the absence Figure 16: Magnetic structure with triple translational periods atintermediate field (between the lower and upper frustration fields) ina model with nearest-neighbor and third-neighbor interactions.
8f magnetic field coincide R = R ′ = 1 /
2. The magne-tization curves have plateaus and two frustration fields,and the upper fields also coincide. In all the frustrationfields the entropy does not vanish at T →
0. In the caseof strong nearest-neighbor interaction and below the lowerfrustration field the magnetic structure is the same (N´eelantiferromagnetic). At the frustration points and at all thefrustration fields the transition temperature is suppressedto zero.However, substantial dissimilarities in the magneticproperties do not enable the two versions of the model tobe referred to as similar. In the absence of magnetic field,at the frustration point R = 1 / T →
0, while it tends to a non-zero inthe second version at identical frustration point R ′ = 1 / / / R = R ′ = 1 / R > / R ′ > /
2) the magnetic structureshave utterly diverse translational symmetry, SAF in thefirst version and quadruple in the second. The intermedi-ate structures in between the lower and upper frustrationfield also have quite different symmetry.A comparison between the second version of the model(that can equally be called a model with nearest-neighborinteraction and second-neighbor interaction along sides ofa square) and the 1D linear chain gives radically differentresult. Apart from previously established common featuresmany new are revealed. The heights of plateaus are thesame, namely, 1 / R ′ = 1 / R ′ and magnetic field) thesecond version along either side of a square coincides withthat in a linear chain. All the frustration fields coincidebeing expressed in the rescaled form h = H/z . A completeagreement between the magnetizations (Figs. 5 and 11) isthe most striking similarity of the second version of 2Dmodel and a linear chain.We may ultimately conclude that magnetic propertiesof the linear chain with nearest-neighbor and second-neighbor interactions and the 2D model with nearest-neighbor and third-neighbor interactions are alike at al-most every aspect provided that the temperature and mag-netic field are rescaled.We predict that at appropriate choice of a model themagnetic properties in 1D, 2D, and 3D lattices will besimilar. In particular, we do believe that the future nu-merical calculations of magnetizations on the Ising simplecubic lattice with nearest-neighbor interaction and second-neighbor interaction along all three cube edges should re-produce the 1D magnetizations from Fig. 5.
Acknowledgments
This work was partially supported by projects no.12-I-2-2020 of Ural Division RAS, no.12-P-2-1041 of PresidiumRAS, by Russian Foundation for Basic Research (projectsno.12-02-96504-r yug a, and no.13-02-00220).
References [1] L. Onsager, Phys. Rev. 65 (1944) 117.[2] G.H. Wannier, Phys. Rev. 79 (1950) 357.[3] R.M.F. Houtappel, Physica (Amsterdam) 16 (1950) 425.[4] K. Kanˆo, S. Naya, Prog. Theor. Phys. 10 (1953) 158.[5] N.W. Dalton and D.W. Wood, J. Math. Phys. 10 (1969) 1271.[6] C. Rapaport and C. Domb, J. Phys. C 4 (1971) 2684.[7] J. Oitmaa, J. Phys. A 14 (1981) 1159.[8] S.-F. Lee and K.-Y. Lin, Chin. J. of Phys. 34 (1996) 1261.[9] E. M¨uller-Hartmann and J. Zittartz, Z. Physik B 27 (1977) 261.[10] B. Nienhuis and M. Nauenberg, Phys. Rev. B 13 (1976) 2021.[11] K.R. Subbaswamy and G.D. Mahan, Phys. Rev. Lett. 37 (1976)642.[12] B. Schuh, Z. Phys. 31 (1978) 55.[13] W. Kinzel, Phys. Rev. B 19 (1979) 4584.[14] S.L.A. de Queiroz, Phys. Rev. E 84 (2011) 031132.[15] M.P. Nightingale and H.W.J. Blote, Physica A 251 (1998) 211.[16] J.L. Monroe. and Kim S. Phys. Rev. E 76 (2007) 021123.[17] S.-Y. Kim, Phys. Rev. Lett. 93 (2004) 130604.[18] A.K. Murtazaev, M.K. Ramazanov, and M.K. Badiev, LowTemp. Phys. 37 (2011) 1001.[19] D.P. Landau, Phys. Rev. B 21 (1980) 1285.[20] K. Binder and D.P. Landau, Phys. Rev. B 21 (1980) 1941.[21] D.P. Landau and K. Binder, Phys. Rev. B 31 (1985) 5946.[22] F.A. Kassan-Ogly, B.N. Filippov, V.V. Menshenin, A.K. Mur-tazaev, M.K. Ramazanov, M.K. Badiev, Solid State Phenom.168-169 (2011) 435.[23] F.A. Kassan-Ogly, B.N. Filippov, A.K. Murtazaev, M.K. Ra-mazanov, M.K. Badiev, JMMM. 324 (2012) 3418.[24] A. Malakis, P. Kalozoumis, and N. Tyraskis, Eur. Phys. J. B 50(2006) 6367.[25] A. Kalz, A. Honecker, S. Fuchs and T. Pruschke, Eur. Phys. J.B 65 (2008) 533.[26] A. Kalz, A. Honecker, S. Fuchs and T. Pruschke, Journal ofPhysics: Conference Series 145 (2009) 012051.[27] A.K. Murtazaev, M.K. Ramazanov, F.A. Kassan-Ogly,M.K. Badiev, JETP, 144 (2013) 1091.[28] F.A. Kassan-Ogly, Phase Transitions, 74 (2001) 353.[29] J. Yin and D.P. Landau, Phys. Rev. E 80 (2009) 051117.[30] C. Fan and F.Y. Wu, Phys. Rev. 179 (1969) 560.[1] L. Onsager, Phys. Rev. 65 (1944) 117.[2] G.H. Wannier, Phys. Rev. 79 (1950) 357.[3] R.M.F. Houtappel, Physica (Amsterdam) 16 (1950) 425.[4] K. Kanˆo, S. Naya, Prog. Theor. Phys. 10 (1953) 158.[5] N.W. Dalton and D.W. Wood, J. Math. Phys. 10 (1969) 1271.[6] C. Rapaport and C. Domb, J. Phys. C 4 (1971) 2684.[7] J. Oitmaa, J. Phys. A 14 (1981) 1159.[8] S.-F. Lee and K.-Y. Lin, Chin. J. of Phys. 34 (1996) 1261.[9] E. M¨uller-Hartmann and J. Zittartz, Z. Physik B 27 (1977) 261.[10] B. Nienhuis and M. Nauenberg, Phys. Rev. B 13 (1976) 2021.[11] K.R. Subbaswamy and G.D. Mahan, Phys. Rev. Lett. 37 (1976)642.[12] B. Schuh, Z. Phys. 31 (1978) 55.[13] W. Kinzel, Phys. Rev. B 19 (1979) 4584.[14] S.L.A. de Queiroz, Phys. Rev. E 84 (2011) 031132.[15] M.P. Nightingale and H.W.J. Blote, Physica A 251 (1998) 211.[16] J.L. Monroe. and Kim S. Phys. Rev. E 76 (2007) 021123.[17] S.-Y. Kim, Phys. Rev. Lett. 93 (2004) 130604.[18] A.K. Murtazaev, M.K. Ramazanov, and M.K. Badiev, LowTemp. Phys. 37 (2011) 1001.[19] D.P. Landau, Phys. Rev. B 21 (1980) 1285.[20] K. Binder and D.P. Landau, Phys. Rev. B 21 (1980) 1941.[21] D.P. Landau and K. Binder, Phys. Rev. B 31 (1985) 5946.[22] F.A. Kassan-Ogly, B.N. Filippov, V.V. Menshenin, A.K. Mur-tazaev, M.K. Ramazanov, M.K. Badiev, Solid State Phenom.168-169 (2011) 435.[23] F.A. Kassan-Ogly, B.N. Filippov, A.K. Murtazaev, M.K. Ra-mazanov, M.K. Badiev, JMMM. 324 (2012) 3418.[24] A. Malakis, P. Kalozoumis, and N. Tyraskis, Eur. Phys. J. B 50(2006) 6367.[25] A. Kalz, A. Honecker, S. Fuchs and T. Pruschke, Eur. Phys. J.B 65 (2008) 533.[26] A. Kalz, A. Honecker, S. Fuchs and T. Pruschke, Journal ofPhysics: Conference Series 145 (2009) 012051.[27] A.K. Murtazaev, M.K. Ramazanov, F.A. Kassan-Ogly,M.K. Badiev, JETP, 144 (2013) 1091.[28] F.A. Kassan-Ogly, Phase Transitions, 74 (2001) 353.[29] J. Yin and D.P. Landau, Phys. Rev. E 80 (2009) 051117.[30] C. Fan and F.Y. Wu, Phys. Rev. 179 (1969) 560.