Universality from disorder in the random-bond Blume-Capel model
N.G. Fytas, J. Zierenberg, P.E. Theodorakis, M. Weigel, W. Janke, A. Malakis
UUniversality from disorder in the random–bond Blume–Capel model
N. G. Fytas , J. Zierenberg , , , P. E. Theodorakis , M. Weigel , W. Janke , and A. Malakis , Applied Mathematics Research Centre, Coventry University, Coventry, CV1 5FB, United Kingdom Institut f¨ur Theoretische Physik, Universit¨at Leipzig, Postfach 100 920, 04009 Leipzig, Germany Max Planck Institute for Dynamics and Self-Organization, 37077 G¨ottingen, Germany Bernstein Center for Computational Neuroscience, 37077 G¨ottingen, Germany Institute of Physics, Polish Academy of Sciences,Al. Lotnik´ow 32/46, 02-668, Warsaw, Poland and Department of Physics, Section of Solid State Physics,University of Athens, Panepistimiopolis, GR 15784 Zografou, Greece (Dated: October 8, 2018)Using high-precision Monte Carlo simulations and finite-size scaling we study the effect ofquenched disorder in the exchange couplings on the Blume-Capel model on the square lattice.The first-order transition for large crystal-field coupling is softened to become continuous, with adivergent correlation length. An analysis of the scaling of the correlation length as well as thesusceptibility and specific heat reveals that it belongs to the universality class of the Ising modelwith additional logarithmic corrections observed for the Ising model itself if coupled to weak dis-order. While the leading scaling behavior in the disordered system is therefore identical betweenthe second-order and first-order segments of the phase diagram of the pure model, the finite-sizescaling in the ex-first-order regime is affected by strong transient effects with a crossover lengthscale L ∗ ≈
32 for the chosen parameters.
PACS numbers: 75.10.Nr, 05.50.+q, 64.60.Cn, 75.10.Hk
The effect of random disorder on phase transitionsis one of the basic problems in condensed-matterphysics [1]. Examples include quantum Ising magnetssuch as LiHo x Y − x F x [2, 3], nematic liquid crystals inporous media [4], noise in high-temperature superconduc-tors [5] and the anomalous Hall effect [6]. Understand-ing random disorder in classical, equilibrium systems is acrucial step towards solving the more involved problemsin quantum systems [7], for example many-body local-ization with programmable random disorder [8], and innon-equilibrium phase transitions [9].The case of weak disorder coupled to the energy den-sity of systems with continuous transitions is rather wellunderstood: Uncorrelated disorder is relevant and leadsto new critical exponents if the specific-heat exponent α of the pure system is positive, a rule known as the Har-ris criterion [10]. If long-range correlations in the disor-der are present, this rule can be generalized leading tointeresting ramifications [11–16]. These effects, and inparticular the marginal case of a vanishing specific-heatexponent as present in the two-dimensional Ising model,are intriguing and have attracted a large research effortover the past decades [17–24].The situation is less clear for systems undergoing first-order phase transitions that are much more common innature The observation that formally ν = 1 /D and α = 1for such systems in D dimensions suggests that disorderis always relevant in this case, and the general observa-tion is that it indeed softens transitions to become con-tinuous [25]. Such a rounding of discontinuities has beenrigorously established for systems in two dimensions [26],but is believed to be more general – a view that is sup- ported by a mapping of the problem onto the random-field model [27–29]. This general picture is commonlyaccepted, and similar phenomena are recently studied inquantum systems [30–32] and for non-equilibrium phasetransitions [33, 34]. Still, a number of important ques-tions have not been answered in full generality: Is a finitestrength of disorder required to soften a first-order tran-sition? Is there a divergent correlation length? What isthe universality class of the resulting continuous transi-tion [29, 35, 36]?While a softening must occur for arbitrarily small dis-order strength in two dimensions [26–28], the situation isless clear in three dimensions [37, 38], but in both casesone finds divergent correlation lengths. The question ofthe universality class of softened transitions is perhapsthe most intriguing one. This has been studied in somedetail for the random-bond q -state Potts model [39–41].It turns out to be difficult to determine the exponentswith sufficient precision to arrive at decisive statements,but the most likely situation appears to be that ν ≈ q , while the magnetic ratio β/ν changeswith q , a scenario that has recently also found additionalsupport in perturbation theory [42].A fertile testing ground for predictions relating to thebehavior of first-order transitions under the influence ofquenched disorder is the Blume-Capel model [43, 44]. Ithas been used to describe the prime nuclear fuel ura-nium dioxide [43], Mott insulators [45, 46], He– Hemixtures [48, 49] and more general multi-component flu-ids [50], as well as potentially the hardest piezomagnetknown [47]. The pure system features a tricritical pointseparating second-order and first-order lines of transi- a r X i v : . [ c ond - m a t . d i s - nn ] M a r . . T (0,0) = 2 ferromagneti (F)paramagneti (P) (∆ t , T t )2 nd order st order T = 1 . T = 0 . T ∆ FIG. 1: (Color online) Phase diagram of the pure two-dimensional Blume-Capel model [51], showing the ferromag-netic ( F ) and paramagnetic ( P ) phases that are separated bya continuous transition for small ∆ (solid line) and a first-order transition for large ∆ (dotted line). The line segmentsmeet at a tricritical point, as indicated by the black rhombus.The horizontal arrows indicate the paths of crossing the phaseboundary implemented in the simulations of the present work. tions [51]. There are many open questions concerningthe behavior of this model in presence of quenched dis-order, be it of random-bond type as considered here orin the form of random (crystal) fields [52, 53]. In par-ticular, conflicting results have been found for the uni-versality class of the ex-first-order segment of the transi-tion line [54], and some authors have favored a scenariothat contradicts universality [55, 56]. In the following, wepresent the results of high-statistics Monte Carlo simula-tions that demonstrate that the transitions in the second-order and the first-order segments of the transition lineof the pure system are in the same universality class aftercoupling to the disorder, and this class is that of the two-dimensional (random) Ising model. Hence any doubtsabout the universality of critical behavior in this systemare dispelled.We study the spin-1 or Blume-Capel model [43, 44]with Hamiltonian H = − (cid:88) (cid:104) xy (cid:105) J xy σ x σ y + ∆ (cid:88) x σ x = E J + ∆ E ∆ = E, (1)where the spin variables σ x ∈ {− , , +1 } live on a squarelattice with periodic boundaries and (cid:104) xy (cid:105) indicates sum-mation over nearest neighbors. The couplings J xy aredrawn from a bimodal distribution P ( J xy ) = 12 [ δ ( J xy − J ) + δ ( J xy − J )] , (2)where following Refs. [55, 56] we choose J + J = 2and J > J >
0, so that r = J /J defines the disor-der strength. The crystal field ∆ controls the density ofvacancies, i.e. , sites with σ x = 0. The pure model has been studied extensively (for a review see Ref. [51]), thephase diagram in the (∆, T )–plane is shown in Fig. 1:For small ∆ there is a line of continuous transitions be-tween the ferromagnetic and paramagnetic phases thatcrosses the ∆ = 0 axis at T ≈ .
693 [56]. For large ∆,on the other hand, the transition becomes discontinuousand it meets the T = 0 line at ∆ = zJ/ z = 4 is the coordination number (here we set J = 1,and also k B = 1, to fix the temperature scale). The twoline segments meet in a tricritical point estimated to beat (∆ t ≈ . , T t ≈ . α = 0 there, the Harris criterion is inconclusive, but ex-plicit studies of the Ising model indicate that the singu-larity is only logarithmically modified [17, 59, 60]. Thefirst-order transition gets stronger as ∆ is approachedand, in fact, the interface tension increases linearly withdecreasing temperature [58]. According to the rigorousresult of Aizenman and Wehr [26], the transitions mustsoften under the presence of even arbitrarily weak disor-der, and we expect a second-order transition to emergein this regime too.As the phase boundary in the first-order regime is al-most vertical, it is most convenient to cross it by varyingthe crystal field ∆ while keeping the temperature con-stant. To this end we used a previously developed imple-mentation of the multicanonical method [61, 62] appliedonly to the crystal-field energy E ∆ of Eq. (1) [63]. Themethod iteratively yields a flat histogram along E ∆ byreplacing the canonical Monte Carlo weights exp( − βE )by exp( − βE J ) W ( E ∆ ) and adapting W ( E ∆ ). Our cal-culations are implemented in a parallel fashion followingthe scheme discussed in Refs. [63–65]. This procedureallows us to directly study the probability distribution of E ∆ . The corresponding histogram for T = T = 0 . R = 256 re-alizations of the random couplings for r = 0 .
6, is shownin Fig. 2. For small system sizes there is a clear double-peak structure, characteristic of a first-order phase tran-sition. However, with increasing system size the distribu-tion changes, exhibiting only a single, symmetric peak,clearly illustrating the second-order nature of the transi-tion in the limit L → ∞ . In fact, the inset shows that thefraction of disorder samples with a double peak quicklydecays to zero for increasing L , with R /R ≈ L ≥ L ∗ ≈
32. This is clear evidence that bond disorderwith r = 0 . T = 0 .
574 into a disorder-induced continuousone, yet, with a crossover behavior for small system sizes.To reveal the universality class of the continuous tran-sition resulting from the softening by disorder, we usedan additional array of canonical Monte Carlo simulations,employing a combination of a Wolff single-cluster up-date [66] of the ± . . . .
01 0 20 40 60 80 100 [ P ( E ∆ ) ] L E ∆ /L R p e a k s / R L FIG. 2: (Color online) Probability distribution of crystal fields E ∆ of the random-bond Blume-Capel model at T = 0 .
574 andwith disorder strength r = 0 .
6. The data are averaged over R = 256 disorder samples. With increasing system size, thedouble peak expected for a first-order transition changes to asingle broad peak typical of a continuous transition. The insetshows the fraction of disorder samples exhibiting a doublepeak. the two temperature points indicated by the arrows inFig. 1, the case T = 0 .
574 crossing the phase bound-ary in the first-order regime, and the choice T = 1 . L ∈ { , , , , , , , , , , } for disorderstrength r = 0 .
6. The ensemble sizes, R , of disorder real-izations used are as follows: R = 5 × for L = 8 − R = 3 × for L = 48 −
96, and R = 1 × for L >
96. Error bars were computed from the sample-to-sample fluctuations.We first discuss the ratio of correlation length and sys-tem size, ξ/L . This is known to be universal for a givenchoice of boundary conditions and aspect ratio. For Isingspins on a square lattice with periodic boundary condi-tions as L → ∞ it approaches the value [70] (cid:18) ξL (cid:19) ∞ = 0 .
905 048 829 2(4) . (3)The behavior of the pure, square-lattice Blume-Capelmodel in the second-order regime is found to be per-fectly consistent with Eq. (3) [51]. To determine ξ/L ,we use the second-moment definition of the correlationlength ξ [71, 72]. From the Fourier transform of thespin field, ˆ σ ( k ) = (cid:80) x σ x exp( i kx ), we determine F = (cid:10) | ˆ σ (2 π/L, | + | ˆ σ (0 , π/L ) | (cid:11) / ξ ≡
12 sin( π/L ) (cid:114) (cid:104) M (cid:105) F − , (4)where M = (cid:80) x σ x . To estimate the limiting value of ξ/L we relied on the quotients method [73–75]: The crystal-field value where ξ L /ξ L = 2, i.e. , where the curves of ξ/L .
02 0 .
04 0 .
06 0 .
08 0 . . Ising ( ξ / L ) ∗ /LT = 1 . ; r = 1 (pure) T = 1 . ; r = 0 . T = 0 . ; r = 0 . FIG. 3: (Color online) Finite-size scaling of the correlation-length ratios at their crossing points, ( ξ/L ) ∗ , for the pureand random model and the two temperatures considered inthis work. Results are shown for the following pairs ( L , 2 L )of system sizes: (8 , , , , , , , , , /L . for the sizes L and 2 L cross, defines the finite-size pseudo-critical points ∆ cross . Let us denote the value of ξ/L atthese crossing points as ( ξ/L ) ∗ . In Fig. 3 we show resultsof ( ξ/L ) ∗ for three cases, namely the pure and randommodel at T = 1 .
398 and the random model at T = 0 . T = 1 . r = 0 .
6, with the results for ( ξ/L ) ∗ practi-cally falling on top of the data for the pure system. Forstronger disorder r → L → ∞ . For the pure Blume-Capelmodel at the same temperature, it was previously foundthat ( ξ/L ) ∞ = 0 . /L for L ≥
12 (as shown by the red dashed line) yields (cid:18) ξL (cid:19) T =1 . ∞ , random = 0 . , (5)with goodness-of-fit parameter Q ≈ .
3. This is clearlyconsistent with the Ising value (3). An additional analy-sis of the scaling behavior of the magnetic susceptibilityand specific heat (not shown) is also consistent with Isinguniversality, in line with previous analyses [55, 56].We now turn to the temperature point T = 0 . T = 0 . c = 1 . ∼ L − − − − ∆ c r o ss L ∆ c r o ss − ∆ c L T = 0 . χ ∗ ∼ L γ/ν γ/ν = 1 . C ∗ ∼ ln[ln( L )]] χ ∗ L C ∗ L (a)(b)FIG. 4: (Color online) Finite-size scaling in the ex-first-orderregime of the Blume-Capel model. (a): Shift behavior of thepseudo-critical points ∆ cross estimated at the ( L , 2 L ) cross-ings of the ratio ξ/L shown in Fig. 3. (b): Scaling of the mag-netic susceptibility χ ∗ = χ (∆ cross ) (main panel) and specificheat C ∗ = C (∆ cross ) (inset) evaluated at the pseudo-criticalpoints for the smaller size of the pairs ( L , 2 L ) considered. very strong there, leading to huge and non-monotonousscaling corrections. For smaller lattice sizes, the ratios( ξ/L ) ∗ do not show any tendency of converging to theuniversal Ising value until L ≈
32, when ( ξ/L ) ∗ attains aminimum. Only for larger lattices the correlation lengthratios start to approach the Ising limit approximatelylinearly in 1 /L . Taking lattice sizes L ≥ L ∗ ≈
32 intoaccount, a linear fit in 1 /L (as shown by the blue dashedline) yields (cid:18) ξL (cid:19) T =0 . ∞ , random = 0 . , (6)with Q ≈ .
3. The limit is again fully consistent with theIsing value. Note that the point L ≈
32 of the minimumcorresponds to the crossover length scale L ∗ determinedalready as the size where the first-order nature of thetransition disappears for the chosen disorder strength r =0 .
6, see Fig. 2. While the extrapolated value (6) of the correlation-length ratio ( ξ/L ) ∗ is strong evidence for Ising behavior,universality classes are characterized by the entirety oftheir critical exponents and universal amplitude ratios.We therefore also considered the scaling of the pseudo-critical points ∆ cross , as well as the magnetic suscepti-bility χ and the specific heat C [76], both evaluated at∆ cross . We first considered the scaling for the temper-ature T = 0 .
574 in the first-order regime of the puresystem. For large system sizes, the crossing points areexpected to scale as∆ cross = ∆ c + bL − /ν , (7)where ν = 1 for the two-dimensional Ising universalityclass. Our data for the pseudo-critical points are shownin Fig. 4(a), and we again observe strong scaling cor-rections with a pronounced turnaround in the behavioraround the crossover length scale L ∗ ≈
32. As the insetillustrates, however, the behavior for L ≥
32 is in perfectagreement with the inversely linear behavior expectedfrom Eq. (7) with ν = 1. In fact, a fit of the form (7) for L ≥
32 with Q ≈ . c = 1 . ν = 1 . χ and thespecific heat C evaluated at the pseudo-critical points∆ cross at T = 0 .
574 are shown in Fig. 4(b). Follow-ing the above analysis, we exclude small system sizes L ≤
24. For the magnetic susceptibility, a fit of theform χ ∗ ∼ L γ/ν yields γ/ν = 1 . Q ≈ .
9, fullycompatible within error bars to the Ising value 1 .
75. Thespecific heat, shown in the inset of Fig. 4(b), is well de-scribed by a double logarithm C ∗ ∼ ln [ln ( L )] as pre-dicted by Ref. [59], the corresponding fit quality being Q ≈ .
9. Similarly strong corrections to scaling in sus-ceptibility data have also been reported for the dilutedIsing model [21]. An analogous analysis of our data atthe higher temperature T = 1 .
398 in the second-orderregime also yields values compatible to the Ising behav-ior, but without the strong scaling corrections observedfor T = 0 . He– He mixtures in porousmedia as well as impurities in uranium dioxide. Finally,when replacing the random bonds by random fields theBlume-Capel model might hold an answer to the intrigu-ing question of whether first-order transitions can surviverandomness if it couples to the order parameter insteadof to the energy density [29, 52].N.G.F is grateful to V. Mart´ın-Mayor for his moti-vating comments that triggered this work. The projectwas financially supported by the Deutsch-Franz¨osischeHochschule (DFH-UFA) through the Doctoral College“ L ” under Grant No. CDFA-02-07 as well as bythe EU FP7 IRSES network DIONICOS under con-tract No. PIRSES-GA-2013-612707. J.Z. received fi-nancial support from the German Ministry of Educa-tion and Research (BMBF) via the Bernstein Centerfor Computational Neuroscience (BCCN) G¨ottingen un-der Grant No. 01GQ1005B. This research has been sup-ported by the National Science Centre, Poland, undergrant No. 2015/19/P/ST3/03541. This project has re-ceived funding from the European Union’s Horizon 2020research and innovation programme under the MarieSk(cid:32)lodowska–Curie grant agreement No. 665778. This re-search was supported in part by PLGrid Infrastructure. [1] A.P. Young, ed., Spin Glasses and Random Fields (WorldScientific, Singapore, 1997).[2] S. M. A. Tabei, M. J. P. Gingras, Y.-J. Kao, P. Stasiak,and J.-Y. Fortin, Phys. Rev. Lett. , 237203 (2006).[3] D. M. Silevitch, D. Bitko, J. Brooke, S. Ghosh, G. Aeppli,and T. F. Rosenbaum, Nature , 567 (2007).[4] A. Maritan, M. Cieplak, T. Bellini, and J. R. Banavar,Phys. Rev. Lett. , 4113 (1994).[5] E. W. Carlson, K. A. Dahmen, E. Fradkin, and S. A.Kivelson, Phys. Rev. Lett. , 097003 (2006)[6] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, andN. P. Ong, Rev. Mod. Phys. , 1539 (2010).[7] T. Vojta and J. A. Hoyos, Phys. Rev. Lett. , 075702(2014).[8] J. Smith, A. Lee, P. Richerme, B. Neyenhuis, P. W. Hess, P. Hauke, M. Heyl, D. A. Huse, and C. Monroe, Nat.Phys. , 907 (2016).[9] H. Barghathi and T. Vojta, Phys. Rev. Lett. , 170603(2012).[10] A.B. Harris, J. Phys. C , 1671 (1974).[11] A. Weinrib and B.I. Halperin, Phys. Rev. B , 413(1983).[12] J.M. Luck, Europhys. Lett. , 359 (1993).[13] W. Janke and M. Weigel, Phys. Rev. B , 144208(2004).[14] H. Barghathi and T. Vojta, Phys. Rev. Lett. , 120602(2014).[15] N. Fricke, J. Zierenberg, M. Marenz, F.P. Spitzner, V.Blavatska, and W. Janke, Condens. Matter Phys. ,13004 (2017).[16] J. Zierenberg, N. Fricke, M. Marenz, F.P. Spitzner, V.Blavatska, and W. Janke, Phys. Rev. E 96, 062125(2017).[17] V.S. Dotsenko and V.S. Dotsenko, Adv. Phys. , 129(1983).[18] B.N. Shalaev, Sov. Phys. Solid State , 1811 (1984);Phys. Rep. , 129 (1994).[19] R. Shankar, Phys. Rev. Lett. , 2466 (1987); ibid. ,2390 (1988).[20] A.W.W. Ludwig, Phys. Rev. Lett. , 2388 (1988); Nucl.Phys B , 639 (1990).[21] J.-S. Wang, W. Selke, V.S. Dotsenko, and V.B. Andre-ichenko, Physica A , 221 (1990).[22] M. Hasenbusch, F.P. Toldin, A. Pelissetto, and E. Vicari,Phys. Rev. E , 011110 (2008).[23] R. Kenna and J.J. Ruiz-Lorenzo, Phys. Rev. E ,031134 (2008).[24] V. Dotsenko, Y. Holovatch, M. Dudka, and M. Weigel,Phys. Rev. E , 032118 (2017).[25] J.L. Cardy, Physica A , 215 (1999).[26] M. Aizenman and J. Wehr, Phys. Rev. Lett. , 2503(1989).[27] K. Hui and A.N. Berker, Phys. Rev. Lett. , 2507(1989).[28] A.N. Berker, Physica A , 72 (1993).[29] J.L. Cardy and J.L. Jacobsen, Phys. Rev. Lett. , 4063(1997).[30] P. Goswami, D. Schwab, and S. Chakravarty, Phys. Rev.Lett. , 015703 (2008).[31] R.L. Greenblatt, M. Aizenman, and J.L. Lebowitz, Phys.Rev. Lett. , 197201 (2009); Physica A , 2902(2010).[32] F. Hrahsheh, J.A. Hoyos, and T. Vojta, Phys. Rev. B ,214204 (2012).[33] T. Vojta, J. Phys. A: Math. Gen. , R143 (2006).[34] P. Villa Mart´ın, J.A. Bonachela, and M.A. Mu˜noz, Phys.Rev. E , 012145 (2014).[35] A. Bellafard, H.G. Katzgraber, M. Troyer, and S.Chakravarty, Phys. Rev. Lett. , 155701 (2012).[36] Q. Zhu, X. Wan, R. Narayanan, J.A. Hoyos, and T. Vo-jta, Phys. Rev. B , 224201 (2015).[37] C. Chatelain, B. Berche, W. Janke, and P.E. Berche,Phys. Rev. E , 036120 (2001).[38] C. Chatelain, B. Berche, W. Janke, and P.-E. Berche,Nucl. Phys. B , 275 (2005).[39] S. Chen, A.M. Ferrenberg, and D.P. Landau, Phys. Rev.Lett. , 1213 (1992).[40] M. Picco, Phys. Rev. Lett. , 2998 (1997).[41] B. Berche and C. Chatelain, in Order, Disorder And Crit- icality , edited by Y. Holovatch (World Scientific, Singa-pore, 2004), p. 147.[42] G. Delfino, Phys. Rev. Lett. , 250601 (2017); G.Delfino and E. Tartaglia, J. Stat. Mech. (2017) 123303.[43] M. Blume, Phys. Rev. , 517 (1966).[44] H.W. Capel, Physica (Utr.) , 966 (1966); , 295(1967); , 423 (1967).[45] K. N. Kudin, G. E. Scuseria, and R. L. Martin, Phys.Rev. Lett. , 266402 (2002).[46] N. Lanat`a, Y. Yao, X. Deng, V. Dobrosavljevi´c, and G.Kotliar, Phys. Rev. Lett. , 126401 (2017).[47] M. Jaime, et al. , Nat. Commun. , 99 (2017).[48] M. Blume, V. J. Emery, and R. B. Griffiths, Phys. Rev.A , 1071 (1971).[49] I.D. Lawrie and S. Sarbach, in Phase Transitions andCritical Phenomena , edited by C. Domb, J.L. Lebowitz.,Vol. 9 (Academic Press, London, 1984).[50] N.B. Wilding, Phys. Rev. E , 926 (1996).[51] J. Zierenberg, N.G. Fytas, M. Weigel, W. Janke, and A.Malakis, Eur. Phys. J. Special Topics , 789 (2017).[52] Sumedha and N.K. Jana, J. Phys. A: Math. Theor. ,015003 (2017).[53] P.V. Santos, F.A. da Costa, and J.M. de Ara´ujo, J. Magn.Magn. Mater. , 737 (2018).[54] P.E. Theodorakis and N.G. Fytas, Phys. Rev. E ,011140 (2012).[55] A. Malakis, A.N. Berker, I.A. Hadjiagapiou, and N.G.Fytas, Phys. Rev. E , 011125 (2009).[56] A. Malakis, A.N. Berker, I.A. Hadjiagapiou, N.G. Fy-tas, and T. Papakonstantinou, Phys. Rev. E , 041113(2010).[57] W. Kwak, J. Jeong, J. Lee, and D.-H. Kim, Phys. Rev.E , 022134 (2015).[58] M. Jung and D.-H. Kim, Eur. Phys. J. B , 245 (2017).[59] V.S. Dotsenko and V.S. Dotsenko, Sov. Phys. JETP Lett. , 37 (1981).[60] W. Selke, L.N. Shchur, and O.A. Vasilyev, Physica A , 388 (1998).[61] B.A. Berg and T. Neuhaus, Phys. Lett. B , 249(1991); Phys. Rev. Lett. , 9 (1992).[62] W. Janke, Int. J. Mod. Phys. C , 1137 (1992); PhysicaA , 164 (1998).[63] J. Zierenberg, N.G. Fytas, and W. Janke, Phys. Rev. E , 032126 (2015).[64] J. Zierenberg, M. Marenz, and W. Janke, Comput. Phys.Commun. , 1155 (2013).[65] J. Gross, J. Zierenberg, M. Weigel, and W. Janke, Com-put. Phys. Commun. , 387 (2018).[66] U. Wolff, Phys. Rev. Lett. , 361 (1989).[67] H.W.J. Bl¨ote, E. Luijten, and J.R. Heringa, J. Phys. A:Math. Gen. , 6289 (1995).[68] M. Hasenbusch Phys. Rev. B , 174433 (2010).[69] A. Malakis, A.N. Berker, N.G. Fytas, and T. Papakon-stantinou, Phys. Rev. E , 061106 (2012).[70] J. Salas, and A.D. Sokal, J. Stat. Phys. , 551 (2000).[71] F. Cooper, B. Freedman, and D. Preston, Nucl. Phys. B , 210 (1982).[72] H.G. Ballesteros, L.A. Fern´andez, V. Mart´ın-Mayor, A.Mu˜noz Sudupe, G. Parisi, and J.J. Ruiz-Lorenzo, J.Phys. A: Math. Gen. , 1 (1999).[73] M.P. Nightingale, Physica (Amsterdam) , 561(1976).[74] H.G. Ballesteros, L.A. Fern´andez, V. Mart´ın-Mayor, andA. Mu˜noz-Sudupe, Phys. Lett. B , 207 (1996). [75] N.G. Fytas and V. Mart´ın-Mayor, Phys. Rev. Lett. ,227201 (2013); Phys. Rev. E , 063308 (2016); N.G.Fytas, V. Mart´ın-Mayor, M. Picco, and N. Sourlas, Phys.Rev. Lett. , 227201 (2016); Phys. Rev. E , 042117(2017).[76] Following Refs. [51, 63] we have used the following def-initions for the magnetic susceptibility and the specificheat: χ = β (cid:0) (cid:104) M (cid:105) − (cid:104)| M |(cid:105) (cid:1) /L and C ≡ ∂ (cid:104) E J (cid:105) ∂ ∆ 1 L = − β ( (cid:104) E J E ∆ (cid:105) − (cid:104) E J (cid:105)(cid:104) E ∆ (cid:105) ) /L , respectively.[77] F.D.A. Aar˜ao Reis, S.L. de Queiroz, and R.R. dos Santos,Phys. Rev. B , R9616 (1996).[78] H.-O. Heuer, Europhys. Lett. , 503, (1991).[79] A.L. Talapov and L.N. Shchur, Europhys. Lett.27