Wang-Landau study of the 3D Ising model with bond disorder
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r EPJ manuscript No. (will be inserted by the editor)
Wang-Landau study of the 3D Ising model with bond disorder
P.E. Theodorakis , , and N.G. Fytas Faculty of Physics, University of Vienna, Boltzmanngasse 5, 1090 Vienna, Austria Institute for Theoretical Physics and Center for Computational Materials Science, Vienna University of Technology, Haupt-straße 8-10, 1040 Vienna, Austria Vienna Computational Materials Laboratory, Sensengasse 8/12, 1090 Vienna, Austria Department of Materials Science, University of Patras, 26504 Patras, GreeceReceived: date / Revised version: date
Abstract.
We implement a two-stage approach of the Wang-Landau algorithm to investigate the criticalproperties of the 3D Ising model with quenched bond randomness. In particular, we consider the case wheredisorder couples to the nearest-neighbor ferromagnetic interaction, in terms of a bimodal distribution ofstrong versus weak bonds. Our simulations are carried out for large ensembles of disorder realizations andlattices with linear sizes L in the range L = 8 −
64. We apply well-established finite-size scaling techniquesand concepts from the scaling theory of disordered systems to describe the nature of the phase transitionof the disordered model, departing gradually from the fixed point of the pure system. Our analysis (basedon the determination of the critical exponents) shows that the 3D random-bond Ising model belongs tothe same universality class with the site- and bond-dilution models, providing a single universality classfor the 3D Ising model with these three types of quenched uncorrelated disorder.
PACS.
PACS. 05.50+q Lattice theory and statistics (Ising, Potts. etc.) – 64.60.De Statistical mechanicsof model systems – 75.10.Nr Spin-glass and other random models
Understanding the role of impurities on the nature ofphase transitions is of great importance, both from ex-perimental and theoretical point of view [1]. First-orderphase transitions are known to be significantly softenedunder the presence of quenched randomness [2,3,4,5,6,7,8],whereas continuous transitions may have their exponentsaltered under random fields or random bonds [9,10]. Thereare some very useful phenomenological arguments andsome, perturbative in nature, theoretical results, pertain-ing to the occurrence and nature of phase transitions un-der the presence of quenched randomness [4,5,11,12]. His-torically, the most celebrated criterion is that suggested byHarris [9]. This criterion relates directly the persistence,under random bonds, of the non random behavior to thespecific heat exponent α p of the pure system. Accordingto this criterion, if α p >
0, then disorder will be rele-vant, i.e., under the effect of the disorder, the system willreach a new critical behavior. Otherwise, if α p <
0, disor-der is irrelevant and the critical behavior will not change.Pure systems with a zero specific heat exponent ( α p = 0)are marginal cases of the Harris criterion and their study,upon the introduction of disorder, has been of particularinterest [13]. The paradigmatic model of the marginal case a e-mail:[email protected] b e-mail:[email protected] is, of course, the general random 2D Ising model, whichhas been extensively debated [14].Respectively, the 3D Ising model with quenched ran-domness - which is a clear case in terms of the Harris cri-terion having a positive specific heat exponent in its pureversion - has also been extensively studied using MonteCarlo (MC) simulations [15,16,17,18,19,20,21,22,23,24,25]and field theoretical renormalization group approaches [26,27,28].Especially, the diluted model can be treated in the low-dilution regime by analytical perturbative renormaliza-tion group methods [29,30,31], where a new fixed pointindependent of the dilution has been found, yet for thestrong dilution regime only MC results remain valid. Al-though the first numerical studies of the model suggesteda continuous variation of the critical exponents along thecritical line, it soon became clear that the concentration-dependent critical exponents found in MC simulations arethe effective ones characterizing the approach to the asymp-totic regime [17,18,19]. Note, here, that a crucial problemof the new critical exponents obtained in these studies isthat the ratios β/ν and γ/ν occurring in finite-size scal-ing (FSS) analysis are almost identical for the disorderedand pure models. In fact, for the pure 3D Ising model, ac-curate values are [32]: ν = 0 . β/ν = 0 . γ/ν = 1 . α = 0 . ν , β/ν , and γ/ν ) have been P.E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder given by the numerical works of Ballesteros et al. [21] andBerche et al. [25] are: (0 . . . . . . We start this Section with the definition of the model.In the following we consider the 3D bond-disorder Isingmodel, whose Hamiltonian with uncorrelated quenchedrandom interactions reads H = − X h ij i J ij S i S j , (1) where the spin variables S i take on the values − , +1, h ij i indicates summation over all nearest-neighbor pairsof sites, and the ferromagnetic interactions J ij > P ( J ij ) = 12 [ δ ( J ij − J ) + δ ( J ij − J )] . (2)In equation (2) we set J + J = 2, J > J , and r = J /J reflects the strength of the bond randomness. Addition-ally, in the following, we fix as usual 2 k B / ( J + J ) = 1,to set the temperature scale ( k B = 1 also for simplicity).The values of the disorder strength r considered through-out this paper are the following: r = 0 . / .
25, 0 . / . . / . E , E )for the entropic sampling of each realization have beenpresented in references [41,42]. To estimate the appropri-ate subspace from a chosen pseudocritical temperatureone should be careful to account for the shift behaviorof other important pseudocritical temperatures and ex-tend the subspace appropriately from both low- and high-energy sides (as also discussed in reference [41]) in orderto achieve an accurate estimation of all finite-size anoma-lies. Of course, taking the union of the corresponding sub-spaces, ensures accuracy for the temperature region of allstudied pseudocritical temperatures.The up to date version of our implementation uses acombination of several stages of the WL process. First, wecarry out a starting (or preliminary) multi-range (multi-R) stage, in a very wide energy subspace. This prelimi-nary stage is performed up to a certain level of the WLrandom walk. The WL refinement is G ( E ) → f · G ( E ),where G ( E ) is the density of states (DOS) and we fol-low the usual modification factor adjustment f j +1 = p f j .E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder 3 * q q L = 32 ; r = 0.5 / 1.5
Fig. 1.
Disorder distribution of the susceptibility maxima of alattice with linear size L = 32 for disorder strength r = 0 . / . and f = e [38,39]. The preliminary stage may consist ofthe levels : j = 1 , . . . , j = 18 and to improve accuracythe process may be repeated several times. However, inrepeating the preliminary process and in order to be effi-cient, we use only the levels j = 13 , . . . ,
18 after the firstattempt, using as starting DOS the one obtained in thefirst random walk at the level j = 12. From our experi-ence, this practice is almost equivalent to simulating thesame number of independent WL random walks. Also inour recent studies we have found out that is much moreefficient and accurate to loosen up the originally appliedvery strict flatness criteria [38]. Thus, a variable flatnessprocess starting at the first levels with a very loose flat-ness criteria and assuming at the level j = 18 the originalstrict flatness criteria is nowadays used. After the abovedescribed preliminary multi-R stage, in the wide energysubspace, one can proceed in a safe identification of theappropriate energy subspace using one or more alterna-tives outlined in reference [39].The process continues in two further stages (two-stageprocess), using now mainly high iteration levels, wherethe modification factor is very close to unity and there isnot any significant violation of the detailed balance con-dition during the WL process. These two stages are suit-able for the accumulation of histogram data (for instance,energy-magnetization histograms), which can be used foran accurate entropic calculation of non-thermal thermo-dynamic parameters, such as the order parameter and itssusceptibility [39]. In the first (high-level) stage, we followagain a repeated several times (typically ∼ −
10) multi-R WL approach, carried out now only in the restrictedenergy subspace. The WL levels may be now chosen as j = 18 , ,
20 and as an appropriate starting DOS for thecorresponding starting level the average DOS of the pre-liminary stage at the starting level may be used. Finally,the second (high-level) stage is applied in the refinementWL levels j = j i , . . . , j i + 3 (typically j i = 21), wherewe usually test both an one-range (one-R) or a multi- L = 8Q = 2000[ ] av r = 0.75 / 1.25 r = 0.5 / 1.5 r = 0.25 / 1.75 [ C ] a v ( T ) ; [] a v ( T ) T [C] av r = 0.75 / 1.25 r = 0.5 / 1.5 r = 0.25 / 1.75 Fig. 2. (Color online) Disorder-averaged specific heat [ C ] av (lower curves) and magnetic susceptibility [ χ ] av (upper curves)as a function of the temperature T for a lattice with linearsize L = 8. The ensemble of random realizations is Q = 2000for all cases of the disorder strength shown here. Notice theexpected shift of the pseudocritical temperatures to the leftwith increasing the disorder strength ( r → R approach with large energy intervals. In the case ofthe one-R approach we have found very convenient andmore accurate to follow the Belardinelli-Pereyra [45] ad-justment of the WL modification factor according to therule ln f ∼ t − . Finally, it should be also noted that byapplying in our scheme a separate accumulation of his-togram data in the starting multi-R stage (in the wideenergy subspace) offers the opportunity to inspect the be-havior of all basic thermodynamic functions in an alsowide temperature range and not only in the neighborhoodof the finite-size anomalies.Using this scheme we performed extensive simulationsfor several lattice sizes in the range L = 8 −
64, overlarge ensembles { , · · · , q, · · · , Q } of random realizations( Q = 1000 − Z asfollows: R Z = ([ Z ] av − [ Z ] av ) / [ Z ] av . Figure 1 presentsevidence that the above number of random realizations issufficient in order to obtain the true average behavior andnot a typical one. In particular, we plot in this figure (fora lattice size L = 32 and disorder strength r = 0 . / . χ ∗ q (where the subscript q denotes the random realization)and the corresponding running average, i.e., a series ofaverages of different subsets of the full data set - each ofwhich is the average of the corresponding subset of a largerset of data points, over the samples for the simulated en-semble of Q = 1447 disorder realizations. A first striking P.E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder
10 20 30 40 50 60 704.304.354.404.454.504.554.60
Z = C Z = Z = <|M|> / K Z = ln
Shift behavior of several pseudocritical temperaturesdefined in the text for r = 0 . / .
25. The error bars reflect thesample-to-sample fluctuations. observation from this figure is the existence of very largevariance of the values of χ ∗ q , indicating the violation ofself-averaging for this quantity [22,46]. This figure illus-trates that the simulated number of random realizationsis sufficient in order to probe correctly the average behav-ior of the system, since already for Q ≈
300 the averagevalue of χ ∗ q is stable.Closely related to the above issue of self-averaging indisordered systems is the manner of averaging over the dis-order [47]. This non-trivial process may be performed intwo distinct ways when identifying the finite-size anoma-lies, such as the peaks of the magnetic susceptibility. Thefirst way corresponds to the average over disorder real-izations ([ . . . ] av ) and then taking the maxima ([ . . . ] ∗ av ), ortaking the maxima in each individual realization first, andthen taking the average ([ . . . ∗ ] av ). In the present paper wepresent our FSS analysis using mainly the first approachof averaging, although we should note that also the sec-ond gave comparable results. As an example, we show inFigure 2 the curves of the specific heat and magnetic sus-ceptibility for a lattice with linear size L = 8 averaged over Q = 2000 random realizations. One may retrieve the loca-tion of the pseudocritical point by taking the maximum ofthese curves. Commenting on the statistical errors of ourWL scheme, we found that the statistical errors of ourscheme on the observed average behavior were of smallmagnitude (of the order of the symbol sizes) and thus areneglected in the figures. On the other hand, for the case[ . . . ∗ ] av the error bars shown reflect the sample-to-samplefluctuations. We present in this Section our numerical results and FSSanalysis for the 3D random-bond Ising model. In Figures 3- 5 we illustrate in the main panels the shift behaviorof several pseudocritical temperatures, i.e., we take the
10 20 30 40 50 60 704.14.24.34.44.54.6 r = 0.5 / 1.5
10 20 30 40 50 60 700.0000.0020.0040.0060.008 = 0.688(15) ([ T * ] a v ) L T c = 4.4219(49) = 0.6843(67) [ T * Z ] a v L Z = C Z = Z = <|M|> / K Z = ln
Same as Figure 3, but for r = 0 . / .
5. The inset showsthe FSS of the sample-to-sample fluctuations of the pseudo-critical temperature of the magnetic susceptibility. average over the pseudocritical temperatures of the in-dividual samples. In all cases the error bars reflect thesample-to-sample fluctuations. The considered pseudocrit-ical temperatures correspond to the peaks of the followingsix quantities: specific heat C , magnetic susceptibility χ ,derivative of the absolute order parameter with respect toinverse temperature K = 1 /T [48] ∂ h| M |i /∂K = h| M | H i − h| M |ih H i (3)and logarithmic derivatives of the first ( n = 1), second( n = 2), and fourth ( n = 4) powers of the order parameterwith respect to inverse temperature [48] ∂ ln h M n i /∂K = h M n H i / h M n i − h H i . (4)Fitting our data for the whole lattice range to the expectedpower-law behavior[ T ∗ Z ] av = T c + bL − /ν , (5)where Z stands for the different thermodynamic quantitiesmentioned above, we estimate the critical temperaturesas a function of the disorder strength and the respectivecritical exponents of the correlation length. In particular,for the case of very weak disorder (case r = 0 . / . T c = 4 . ν = 0 . ν islarger than the corresponding value of the pure system,it is still far away from the expected value of the randomsystem (as obtained in references [21,25,34] and Figures 4and 5). This is due to finite-size effects that dominatein the regime of very weak randomness, hindering the ap-proach of the system to the asymptotic limit. On the otherhand, the simultaneous fittings of the form (5) shown inthe main panels of Figures 4 and 5 for stronger valuesof the disorder strength r reveal estimates for the corre-lation length exponent ν around the value 0 . . . .E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder 5
10 20 30 40 50 60 704.04.14.24.34.44.5 r = 0.25 / 1.75
10 20 30 40 50 60 700.0000.0040.0080.0120.016 = 0.679(21) ([ T * ] a v ) L Z = C Z = Z = <|M|> / K Z = ln
Same as Figure 4 but for the case r = 0 . / . given by Ballesteros et al. [21] and Berche et al. [25] forthe site- and bond-diluted models. Additionally, the esti-mated values for the critical temperatures, in all Figures 3,4, and 5, follow the usual reduction trend with increasingdisorder strength ( r → T ∗ Z of the disordered system are distributed with a width δ [ T ∗ Z ] av , that scales, in the case of a new random fixedpoint, with the system size as δ ([ T ∗ Z ] av ) ∼ L − /ν . (6)In the insets of Figures 4 and 5 we plot these sample-to-sample fluctuations of the pseudocritical temperatureof the magnetic susceptibility. The solid lines show, inboth cases, a very good power-law fitting giving the values0 . . ν , which is alsoin very good agreement with the values denoted above andobtained via the classical shift behavior and the accurateestimates in the literature [21,25].We move now to investigate the magnetic propertiesof the model. In Figure 6 we provide estimates for themagnetic exponent ratio γ/ν of the 3D random-bond Isingmodel, by presenting the FSS behavior of the maxima ofthe disorder-averaged magnetic susceptibility [ χ ] ∗ av for allthe values of the disorder strength considered. The solid,dotted, and dashed lines show a simultaneous fitting ofthe form [ χ ] ∗ av ∼ L γ/ν , (7)using the L = 16 −
64 lattice-range spectrum, providingan estimate 1 . γ/ν . We re-mind the reader that the respective fitting attempts of thenumerical data for each case separately gave the estimates γ/ν = 1 . . . r = 0 . / .
10 20 30 40 50 60 70050100150200250 R [ * ] a v r = 0.75 / 1.25 r = 0.5 / 1.5 r = 0.25 / 1.75 / = 1.964(9) [] * a v L r = 0.75 / 1.25 r = 0.5 / 1.5 r = 0.25 / 1.75 Fig. 6.
FSS of the maxima of the disorder-averaged magneticsusceptibility for all the three values of the disorder strengthconsidered. The solid, dotted, and dashed lines show a simul-taneous fitting attempt of the form [ χ ] ∗ av ∼ L γ/ν for the range L = 16 −
64. The inset shows the corresponding limiting be-havior of the ratio R [ χ ∗ ] av as defined in the text. . / .
5, and 0 . / .
75, respectively. All these sets of val-ues are in excellent agreement with the accurate estimates1 . . γ/ν of the site- and bond-diluted models, re-spectively. In the corresponding inset of Figure 6 we tryto further quantify the behavior of the sample-to-samplefluctuations of the model by plotting the noise to signalratio R Z - already introduced in Section 2 - as a func-tion of the inverse lattice size, for Z = [ χ ∗ ] av . Clearly, forthe present model the limiting value of R [ χ ∗ ] av is nonzero,indicating, as expected also for marginal disordered sys-tems [22], a strong violation of self-averaging of the mag-netic properties of the 3D Ising model with bond disor-der, a property that intensifies with increasing disorderstrength ( r → R [ χ ∗ ] av has been observed alsoin several other 3D and 2D random models, some of thembeing the bond- [25] and site-diluted [21] versions of thepresent 3D Ising model and the random-bond versionsof the square Blume-Capel [43] and triangular Ising [44]model. As a final remark, we stress that the asymptoticbehavior of these R -ratios is expected, from theoreticalrenormalization-group arguments, to follow specific powerlaw functions and reach limiting values independent of thevalue of the disorder strength [46]. However, such an anal-ysis is by far much more demanding, in terms of numer-ical data, as it has already been discussed by Berche etal. [25], and, in any case, goes beyond the scope of thepresent work.Respectively, in Figure 7 we plot the disorder-averagedmagnetization at the estimated critical temperatures ofthe disordered model, as a function of the lattice size L .The solid, dotted, and dashed lines show a simultaneousfitting following the well-known relation[ M ] av ( T = T c ) ∼ L − β/ν , (8) P.E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder
10 20 30 40 50 60 700.100.150.200.250.30
20 30 40 50 60 704080120160200 [ < | M | > / K ] * a v L / = 0.512(6) / = 0.518(9) / = 0.514(7) / = 0.517(4) [ M ] a v ( T = T c ) L r = 0.75 / 1.25 r = 0.5 / 1.5 r = 0.25 / 1.75 Fig. 7.
FSS of the disorder-averaged order parameter at the es-timated critical temperatures T c ( r ) for r = 0 . / .
25, 0 . / . . / .
75. The inset illustrates the FSS of the maxima ofthe disorder-averaged inverse-temperature derivative of the ab-solute order parameter. In both main panel and inset, latticeswith linear sizes L = 16 −
64 have been used in the fittings. giving an estimate of the order of 0 . β/ν . Let us note here that, as in the case of themagnetic susceptibility in Figure 6, the respective separatefitting attempts of the numerical data gave comparable tothe 0 .
517 estimates in the regime 0 . − . ∂ h| M |i /∂K . Thus, inthe corresponding inset of Figure 7 we plot the data for ∂ h| M |i /∂K averaged over disorder as a function of L , ina double logarithmic scale. Their maxima are expected toscale with the system size as [48]( ∂ h| M |i /∂K ) ∗ ∼ L (1 − β ) /ν . (9)As in the main panel, the three lines shown in the insetcorrespond to linear fittings for the three values of the dis-order strength, verifying the estimates of the main panelof the figure for β/ν , being also in very good agreementwith accurate estimates 0 . . ν of thecorrelation length and the magnetic exponent ratios γ/ν and β/ν given in this Section, and summarized in Table 1,are in excellent agreement with accurate estimates in theliterature for the respective site- [21] and bond-dilutedmodels [25], reinforcing the scenario of a single distinctiveuniversality class in the 3D Ising model with quencheduncorrelated disorder, indicating the presence of a newrandom fixed point. In summary, we reported in the present paper the effectsinduced by the presence of quenched bond randomness on
Table 1.
Summary of critical exponents for the pure and dis-ordered 3D Ising model. The values of references [32], [21], and[25] have been used for the pure, site-, and bond-diluted ver-sions of the model.Disorder Distribution ν γ/ν β/ν
None (Pure model) 0.6304(13) 1.966(3) 0.517(3)Site dilution 0.6837(53) 1.963(5) 0.519(3)Bond dilution 0.68(2) 1.97(2) 0.515(5)Bond randomness 0.685(7) 1.964(9) 0.517(4) the critical behavior of the 3D Ising spin model embeddedin the simple cubic lattice, by implementing an efficienttwo-stage entropic simulation scheme based on the Wang-Landau algorithm. Our numerical approach enabled us tosimulate large ensembles of disorder realizations of themodel for very large lattice sizes (up to N = 64 ) andseveral disorder strengths, and thus, obtain, via finite-sizescaling techniques, accurate estimates of all the criticalexponents describing the phase transition of the model.Particular interest was paid to the sample-to-sample fluc-tuations of the pseudocritical temperatures of the modeland their scaling behavior that was used as a successfulalternative approach to criticality.The main conclusion of this study can be synopsizedas follows: the critical behavior of the 3D Ising model withquenched uncorrelated disorder is controlled by a new ran-dom fixed point, independently of the way randomness isimplemented in the system. It should be acknowledged atthis point that, Ballesteros et al. [21] and Berche et al. [25]were the first scientific groups that strongly supported theabove view for the present model. However, we believethat the results brought forward in this paper, combinedwith the existing ones for the site- and bond-diluted mod-els, complete the picture, at least for this specific case ofthe 3D disordered Ising model, contributing further to theunderstanding of the concept of universality in disorderedspin models.Closing, we stress here that it would be interesting tostudy further the universality aspects of even more com-plicated spin models, in terms of different disorder dis-tributions and couplings. One interesting case would beto consider field-disorder coupling to the local order pa-rameter, and one prominent candidate for this case is the3D random-field Ising model [49,50]. Of course for thistype of model one may be particularly careful, since inthis case there exist the well-known hyperscaling violationand new scaling concepts should be taken under seriousaccount [49,50,51]. Several attempts towards this direc-tion of understanding universality in random systems arecurrently under consideration on both numerical and the-oretical grounds. References Spin glasses and random fields , edited by A.P. Young(World Scientific, Singapore, 1998).E. Theodorakis and N.G. Fytas: Wang-Landau study of the 3D Ising model with bond disorder 72. M. Aizenman, J. Wehr, Phys. Rev. Lett. , 2503 (1989)3. M. Aizenman, J. Wehr, Phys. Rev. Lett. , 1311(E) (1990)4. K. Hui, A.N. Berker, Phys. Rev. Lett. , 2507 (1989)5. K. Hui, A.N. Berker, Phys. Rev. Lett. , 2433(E) (1989)6. S. Chen, A.M. Ferrenberg, D.P. Landau, Phys. Rev. Lett. , 1213 (1992)7. J. Cardy, J.L. Jacobsen, Phys. Rev. Lett. , 4063 (1997)8. C. Chatelain, B. Berche, Phys. Rev. Lett. , 1670 (1998)9. A.B. Harris, J. Phys. C , 1671 (1974)10. J.T. Chayes, L. Chayes, D.S. Fisher, T. Spencer, Phys.Rev. Lett. , 2999 (1986)11. V. Dotsenko, M. Picco, P. Pujol, Nucl. Phys. B , 701(1995)12. J.L. Jacobsen, J. Cardy, Nucl. Phys. B , 701 (1998)13. G. Mazzeo, R. K¨uhn, Phys. Rev. E , 3823 (1999)14. A. Gordillo-Guerrero, R. Kenna, J.J. Ruiz Lorenzo, AIPConf. Proc. , 42 (2009)15. D.P. Landau, Phys. Rev. B , 2450 (1980)16. D. Chowdhury, D. Stauffer, J. Stat. Phys. , 203 (1986)17. H.-O. Heuer, Europhys. Lett. , 551 (1990)18. H.-O. Heuer, Phys. Rev. B , 6476 (1990)19. H.-O. Heuer, J. Phys. A , L333 (1993)20. M. Hennecke, Phys. Rev. B , 6271 (1993)21. H.G. Ballesteros, L. A. Fern´andez, V. Mart´ın-Mayor, A.Mu˜noz Sudupe, G. Parisi, J. J. Ruiz-Lorenzo, Phys. Rev. B , 2740 (1998)22. S. Wiseman, E. Domany, Phys. Rev. Lett. , 22 (1998)23. P. Calabrese, V. Mart´ın-Mayor, A. Pelissetto, E. Vicari,Phys. Rev. E , 036136 (2003)24. M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, E. Vicari,J. Stat. Mech.: Theory Exp. 2007, P0201625. P.E. Berche, C. Chatelain, B. Berche, W. Janke, Eur. Phys.J. B , 463 (2004)26. R. Folk, Y. Holovatch, T. Yavorskii, Phys. Rev. B ,15114 (2000)27. D.V. Pakhnin, A.I. Sokolov, Phys. Rev. B , 15130 (2000)28. A. Pelissetto, E. Vicari, Phys. Rev. B , 6393 (2000)29. K.E. Newman, E.K. Riedel, Phys. Rev. B , 264 (1982)30. G. Jug, Phys. Rev. B , 609 (1983)31. I.O. Mayer, J. Phys. A , 2815 (1989)32. R. Guida, J. Zinn-Justin, J. Phys. A , 8103 (1998)33. A.K. Murtazaev, A.B. Babaev, J. Magn. Magn. Mater. , 2630 (2009)34. N.G. Fytas, P.E. Theodorakis, Phys. Rev. E , 062101(2010)35. M.E.J. Newman, G.T.Barkema, Monte Carlo Methods inStatistical Physics (Clarendon Press, Oxford, 1999)36. F. Wang, D.P. Landau, Phys. Rev. Lett. , 2050 (2001)37. F. Wang, D.P. Landau, Phys. Rev. E , 056101 (2001)38. A. Malakis, A. Peratzakis, N.G. Fytas, Phys. Rev. E ,066128 (2004)39. A. Malakis, S.S. Martinos, I.A. Hadjiagapiou, N.G. Fytas,P. Kalozoumis, Phys. Rev. E , 066120 (2005)40. A. Malakis, N.G. Fytas, Phys. Rev. E , 016109 (2006)41. N.G. Fytas, A. Malakis, K. Eftaxias, J. Stat. Mech.: The-ory Exp. (2008) P0301542. N.G. Fytas, A. Malakis, I.A. Hadjiagapiou, J. Stat. Mech.:Theory Exp. (2008) P1100943. A. Malakis, A.N. Berker, I.A. Hadjiagapiou, N.G. Fytas,T. Papakonstantinou, Phys. Rev. E , 041113 (2010)44. N.G. Fytas, A. Malakis, Phys. Rev. E , 041109 (2010)45. R.E. Belardinelli, V.D. Pereyra, Phys. Rev. E , 046701(2007) 46. A. Aharony, A.B. Harris, Phys. Rev. Lett. , 3700 (1996)47. K. Bernardet, F. P´azm´andi, G.G. Batrouni, Phys. Rev.Lett. , 4477 (2000)48. A.M. Ferrenberg, D.P. Landau, Phys. Rev. B , 5081(1991)49. R.L.C. Vink, T. Fischer, K. Binder, Phys. Rev. E ,051134 (2010)50. N. Sourlas, Comput. Phys. Commun. , 183 (1999)51. N.G. Fytas, A. Malakis, Eur. Phys. J. B79