Multicritical Points and Crossover Mediating the Strong Violation of Universality: Wang-Landau Determinations in the Random-Bond d=2 Blume-Capel model
A. Malakis, A. Nihat Berker, I. A. Hadjiagapiou, N. G. Fytas, T. Papakonstantinou
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] J a n Multicritical Points and Crossover Mediating the Strong Violation of Universality:Wang-Landau Determinations in the Random-Bond d = 2 Blume-Capel model
A. Malakis , A. Nihat Berker , , I. A. Hadjiagapiou , N. G. Fytas , and T. Papakonstantinou Department of Physics, Section of Solid State Physics,University of Athens, Panepistimiopolis, GR 15784 Zografos, Athens, Greece Faculty of Engineering and Natural Sciences, Sabanc ı University, Orhanl ı , Tuzla 34956, Istanbul, Turkey and Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, U.S.A. (Dated: November 2, 2018)The effects of bond randomness on the phase diagram and critical behavior of the square lat-tice ferromagnetic Blume-Capel model are discussed. The system is studied in both the pure anddisordered versions by the same efficient two-stage Wang-Landau method for many values of thecrystal field, restricted here in the second-order phase transition regime of the pure model. For therandom-bond version several disorder strengths are considered. We present phase diagram points ofboth pure and random versions and for a particular disorder strength we locate the emergence of theenhancement of ferromagnetic order observed in an earlier study in the ex-first-order regime. Thecritical properties of the pure model are contrasted and compared to those of the random model.Accepting, for the weak random version, the assumption of the double logarithmic scenario for thespecific heat we attempt to estimate the range of universality between the pure and random-bondmodels. The behavior of the strong disorder regime is also discussed and a rather complex and yetnot fully understood behavior is observed. It is pointed out that this complexity is related to theground-state structure of the random-bond version.
PACS numbers: 75.10.Nr, 05.50.+q, 64.60.Cn, 75.10.Hk
I. INTRODUCTION
The effect of quenched randomness on the equilib-rium and dynamic properties of macroscopic systemsis a subject of great theoretical and practical interest.It is well known that quenched bond randomness mayproduce drastic changes on phase transitions dependingon the type of the transition [1–6]. Thus, symmetry-breaking first-order transitions are converted to second-order phase transitions by infinitesimal bond randomnessfor spatial dimensionality d = 2 [3, 4] and by bond ran-domness beyond a threshold strength in d > d ≤
2) including a large variety of types of randomnessin classical and quantum spin systems [7].Historically, the effects of disorder on phase transitionshave been studied in two extreme cases, i.e. in the lim-its of weak and strong (near the percolation point) dis-order. The first important conjecture, known today asthe Harris criterion [1], relates the value of the specificheat exponent α in a continuous transition with the ex-pected effects of uncorrelated weak disorder in ferromag-nets. According to the Harris criterion, for continuousphase transitions with a negative exponent α , the intro-duction of weak randomness is expected to be an irrel-evant field and the disordered system to remain in thesame universality class. On the other hand, the weaklydisordered system is expected to be in a different univer-sality class in the case of a pure system having a positiveexponent α . Pure systems with a zero specific heat expo- nent ( α = 0) are marginal cases of Harris criterion (sincethe criterion does not give any information) and theirstudy, upon the introduction of disorder, has been of par-ticular interest. The paradigmatic model of the marginalcase is, of course, the general random 2d Ising model(random-site, random-bond, and bond-diluted) and thismodel has been extensively investigated and debated [8–24]. Several recent studies, both analytical (renormal-ization group and conformal field theories) and numeri-cal (mainly Monte Carlo (MC) simulations) devoted tothis model, have provided very strong evidence in favorof the so-called logarithmic corrections’s scenario. Ac-cording to this, the effect of infinitesimal disorder givesrise to a marginal irrelevance of randomness and besideslogarithmic corrections, the critical exponents maintaintheir 2d Ising values. In particular, the specific heat isexpected to slowly diverge with a double-logarithmic de-pendence [11–14]. Here, we should mention that there isnot full agreement in the literature and a different sce-nario predicts a negative specific heat exponent α leadingto a saturating behavior [16], with a corresponding cor-relation length exponent ν ≥ /d [25].In general, a unitary and rigorous physical descriptionof critical phenomena in disordered systems still lacksand certainly, lacking such a description, the study offurther models for which there is a general agreementin the the behavior of the corresponding pure cases isvery important. Historically, such a suitable candidatefor testing the above predictions, that has been alsoquite extensively studied, is the general 2d q -state Pottsmodel [17, 26–30]. This model includes the Ising model( q = 2), cases of a pure system having continuous transi-tions with a positive exponent α ( q = 3 ,
4) and also thelarge q -cases ( q >
4) for which one could observe and tryto classify the above mentioned softening of first-ordertransitions in 2d models. Another similarly interestingcandidate, not yet as much studied in the random-bondversion, is the 2d Blume-Capel (BC) model [31, 32]. Wemay note here that most of the existing literature on theBC model with randomness concerns randomness appliedto the crystal field and/or spin glass exchange interac-tions [33–36]. As it is well known, the pure version of theBC model undergoes an Ising-like continuous phase tran-sition to an ordered ferromagnetic phase as the tempera-ture is lowered for crystal-field couplings less than a tri-critical value and a first-order transition for larger valuesof the crystal-field coupling. Therefore, this model pro-vides also the opportunity to study two different and veryinteresting topics of the above described effects of disor-der in critical phenomena, namely the double-logarithmicscenario for the specific heat in the regime where the 2dBC model is in the same universality class with the 2drandom-bond Ising model and also the softening of thetransition in the first-order regime. Recently the presentauthors [37] have considered this model and providedstrong numerical evidence clarifying two of the abovementioned effects induced in 2d systems by bond random-ness. By implementing a two-stage Wang-Landau (WL)approach [37–43], we presented essentially exact informa-tion on the 2d BC model under quenched bond random-ness. In this investigation, we found dramatically dif-ferent critical behaviors of the second-order phase transi-tions emerging from the first- and second-order regimes ofthe pure BC model and since, these second-order transi-tions were found to have different critical exponents, ourstudy indicated an interesting strong violation of univer-sality [37]. Namely, different sets of critical exponentson two segments of the same critical line appeared todescribe the two regimes: still-second-order and ex-first-order.In this paper, we extend our earlier work [37], by im-plementing essentially the same two-stage WL approach(Sec. II) and try to give a more complete picture by con-centrating in the weak (still-second-order) regime andsimulate the model for several disorder strengths andmany values of the crystal-field coupling. The abovestatement means that, effectively we will restrict ourstudy to moderate values of the crystal field and moder-ate values of the disorder, intending to observe the fron-tier between the weak- and the strong-disorder universal-ity classes from the disappearance of the expected 2d ran-dom Ising universality class behavior. Thus, in Sec. IIIwe will produce phase diagram points for the random-bond model but also for the pure model, reporting for thepure case a comparison with existing estimates in the lit-erature. More generally, in carrying out this project wehave also considered the pure 2d BC model for severalvalues of the crystal-field, in the second-order regime,observing its finite-size scaling (FSS) behavior. Sec. IVpresents such a comparative study between random andpure models concerning the behavior of all thermody- namic parameters used in the traditional FSS analysis ofMC data. This study enables us to observe some pecu-liarities of the pure model, due to the onset of tricritical-ity, and compare them with the corresponding behaviorof the random model. Furthermore, we try to focus, un-derstand and shed light to the extent of universality ofthe random-bond 2d BC model with the correspondingrandom-bond Ising model, for which the scenario of log-arithmic corrections seems to be the strongest option inthe current literature [24]. The prediction of the range ofsuch universality is far from trivial and the two regimes(weak and strong) have many dissimilarities which arealso reflected in the ground-state structure, as furtherdiscussed below. The attempt to estimate the range ofthe above mentioned universality is accomplished by thenovel idea which assumes the truth of the double logarith-mic scenario for the specific heat in a suitable restrictedrange. This is presented in Sec. V together with some fur-ther crucial observations concerning the behavior of thestrong disorder regime, i.e. the regime where the Isingclass universality does not apply and the system has arather complex and yet not understood behavior. Ourconclusions are summarized in Sec. VI.
II. DEFINITION OF THE MODELS AND THETWO-STAGE WANG-LANDAU APPROACHA. The pure and random-bond Blume-Capelmodels
The (pure) BC model [31, 32] is defined by the Hamil-tonian H p = − J X
16 [49]. As mentioned already in the in-troduction the phase diagram of the model consists of asegment of continuous Ising-like transitions at high tem-peratures and low values of the crystal field which endsat a tricritical point, where it is joined with a secondsegment of first-order transitions between (∆ t , T t ) and(∆ = 2 , T = 0).The model given by Eq. (1) is studied here on thesquare lattice and will be referred to as the pure BCmodel. However, our main focus, on the other hand, isthe case with bond disorder given by the bimodal distri-bution P ( J ij ) = 12 [ δ ( J ij − J ) + δ ( J ij − J )] ; (2) J + J J > J > r = J J , so that r reflects the strength of the bond randomness.The resulting quenched disordered (random-bond) ver-sion of the Hamiltonian defined in Eq. (1) reads now as H = − X
25 = 0 . B. An outline of our implementation of theWang-Landau approach
In the last few years we have used an entropic sam-pling implementation of the WL algorithm [38, 39] tostudy some simple [40, 41], but also some more complexsystems [8, 42, 43]. One basic ingredient of this imple-mentation is a suitable restriction of the energy subspacefor the implementation of the WL algorithm. This wasoriginally termed as the critical minimum energy sub-space (CrMES) restriction [40, 41] and it can be carriedout in many alternative ways, the simplest being that ofobserving the finite-size behavior of the tails of the energyprobability density function (e-pdf) of the system [41].Complications that may arise in random systems canbe easily accounted for by various simple modificationsthat take into account possible oscillations in the e-pdfand expected sample-to-sample fluctuations of individ-ual disorder realizations. In our recent papers [8, 37, 43],we have presented details of various sophisticated routesfor the identification of the appropriate energy subspace( E , E ) for the entropic sampling of each random real-ization. In estimating the appropriate subspace from achosen pseudocritical temperature one should be carefulto account for the shift behavior of other important pseu-docritical temperatures and extend the subspace appro-priately from both low- and high-energy sides in order to achieve an accurate estimation of all finite-size anoma-lies. Of course, taking the union of the correspondingsubspaces, insures accuracy for the temperature regionof all studied pseudocritical temperatures.The up to date version of our implementation uses acombination of several stages of the WL process. First,we carry out a starting (or preliminary) multi-range(multi-R) stage, in a very wide energy subspace. Thispreliminary stage is performed up to a certain level ofthe WL random walk. The WL refinement is G ( E ) → f ∗ G ( E ), where G ( E ) is the density of states (DOS)and we follow the usual modification factor adjustment f j +1 = p f j and f = e . The preliminary stage may con-sist of the levels : j = 1 , . . . , j = 18 and to improve accu-racy the process may be repeated several times. However,in repeating the preliminary process and in order to beefficient, we use only the levels j = 13 , . . . ,
18 after thefirst attempt, using as starting DOS the one obtained inthe first random walk at the level j = 12. From our expe-rience, this practice is almost equivalent of 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 [40, 41]. Thus, a variable flat-ness process starting at the first levels with a very looseflatness criteria and assuming at the level j = 18 the orig-inal strict flatness criteria is now days used. After theabove described preliminary multi-R stage, in the wideenergy subspace, one can proceed in a safe identificationof the appropriate energy subspace using one or more al-ternatives outlined in Refs. [40, 41]. In random systems,where one needs to simulate many disorder realizations,it is also possible and advisable to avoid the identifica-tion of the appropriate energy subspace separately foreach disorder realization by extrapolating from smallerlattices and/or by prediction from preliminary runs onsmall numbers of disorder realizations. In any case, theappropriate subspaces should be defined with sufficienttolerances. In our implementation we use such advanceinformation to proceed in the next stages of the entropicsampling.The process continues in two further stages (two-stageprocess), using now mainly high iteration levels, wherethe modification factor is very close to unity and thereis not any significant violation of the detailed balancecondition during the WL process. These two stages aresuitable for the accumulation of histogram data (for in-stance energy-magnetization histograms), which can beused for an accurate entropic calculation of non-thermalthermodynamic parameters, such us the order parame-ter and its susceptibility [41]. In the first (high-level)stage, we follow again a repeated several times (typically ∼ −
10) multi-R WL approach, carried out now onlyin the restricted energy subspace. The WL levels maybe now chosen as j = 18 , ,
20 and as an appropriatestarting DOS for the corresponding starting level the av-erage DOS of the preliminary stage at the starting levelmay be used. Finally, the second (high-level) stage is ap-plied in the refinement WL levels j = j i , . . . , j i + 3 (typ-ically j i = 21), where we usually test both an one-range(one-R) or a multi-R approach with large energy inter-vals. In the case of the one-R approach we have foundvery convenient and in most cases more accurate to fol-low the Belardinelli and Pereyra [50] adjustment of theWL modification factor according to the rule ln f ∼ t − .Finally, it should be also noted that by applying in ourscheme a separate accumulation of histogram data in thestarting multi-R stage (in the wide energy subspace) of-fers the opportunity to inspect the behavior of all basicthermodynamic functions in an also wide temperaturerange and not only in the neighborhood of the finite-sizeanomalies. The approximation outside the dominant en-ergy subspace is not of the same accuracy with that of therestricted dominant energy subspace but is good enoughfor the observation of the general behavior and providesalso a route of inspecting the degree of approximation.The above described numerical approach was used toestimate the properties of a large number of 100 bonddisorder realizations, for lattice sizes L = 20 −
100 for allcrystal fields and disorder strengths used in this paper,with the exception of the case ∆ = 1 . r = 0 . / . L = 20 − L = 100), and depending on the ther-modynamic parameter, the statistical errors of the WLmethod were found to be of reasonable magnitude andin some cases to be of the order of the symbol sizes, oreven smaller. This was true for both the pure version andthe individual random-bond realizations. However, sincethe nature of the present study is qualitative, not aimingto an accurate exponent estimation, the WL errors willnot be presented in our figure illustrations. We also notethat, for the random-bond version mainly the averagesover the disorder realizations, denoted as [ . . . ] av , will beconsidered in the text and their finite-size anomalies, de-noted as [ . . . ] ∗ av , will be used in our FSS attempts. Dueto very large sample-to-sample fluctuations, mean valuesof individual maxima ([ . . . ∗ ] av ) have not been used inthis study except for illustrative purposes, as in the case∆ = 1 . r = 0 . / . III. PHASE DIAGRAMS: PURE ANDRANDOM-BOND 2D BLUME-CAPEL MODELS
This Section presents, as mentioned in the introduc-tion, phase diagram points for the pure and random-bondmodels and in the case of the pure model compare thecorresponding phase diagram points with the existing es-timates in the literature. This gives also the opportunityto observe the reliability of our numerical approach. Fol-
20 40 60 80 1001.141.161.181.201.221.241.26 20 40 60 80 1000.650.660.670.680.69 (a) T c = 1.149(1) = 0.999(5) T [ Z ] * a v L Z = C Z = Z = <|M|> / K Z = ln
1. The fitting range is L = 20 − . T c = 0 .
659 (with the dotted lines indicating its errors bar-riers) and the dashed line the estimate T c = 0 .
650 given inRef. [48]. lowing the practice of our earlier paper [37], we estimatephase diagram points by fitting our data to the expectedpower-law shift behavior T = T c + bL − /ν of several pseu-docritical temperatures. The traditionally used specificheat and magnetic susceptibility peaks, as well as, thepeaks corresponding to the following logarithmic deriva-tives of the powers n = 1 , , K = 1 /T [51], ∂ ln h M n i ∂K = h M n H ih M n i − h H i , (4)and the peak corresponding to the absolute order-parameter derivative ∂ h| M |i ∂K = h| M | H i − h| M |ih H i , (5)will be implemented for a simultaneous fitting attempt.Such simultaneous fitting attempts are presented inFig. 1. In particular Fig. 1(a) presents the shift be-havior for the random-bond 2d BC model in the weak cross = 1.66(2)T cross = 1.02(3) k B T / J / J Tricritical Pure Random (r = 0.6)
FIG. 2: Phase diagrams of the pure 2d BC model and itsrandom version at the disorder strength r = 0 . / .
25. Theinset illustrates the crossing of the phase boundaries includingan approximate estimate of the crossing point. disorder case r = 0 . / . . .
95, further discussed atthe end of Sec. IV B. As noted already, the data fittedfor the random version of the model are only those ofthe pseudocritical temperatures of the peaks of the aver-aged over disorder thermodynamic parameters, as indi-cated in the figure, where the asterisk denotes the peak ofthe averaged over disorder parameter [ . . . ] ∗ av . The alter-native route of using averages of individual sample pa-rameters gives almost identical estimates. However, incases of strong lack of self-averaging, very large sample-to-sample fluctuations may be present. For simplicity,in all fitting attempts, the whole range L = 20 − r = 0 . / .
25 = 0 .
6. The tricrit-ical point shown is taken from the estimate given byBeale (∆ t , T t )=(1 . , . r = 0 . / .
25 the points of the phase dia-gram were chosen with the intension to be able to approx-imately locate the emergence of the enhancement of fer-romagnetic order observed in our earlier study [37] in theex-first-order regime at ∆ = 1 . ∼ s i = ± cross , T cross ) = (1 . , . . − . r = 0 . / . r = 0 . / .
25, and r = 0 . / .
4, together with corre-sponding phase diagram points given by Beale [48] forthe pure model. Note that, in Silva et al. [49] one canfind an analogous table for the pure 2d BC model includ-ing the phase diagram points given by Beale [48] andthe points produced in their two-parametric WL sam-pling [49]. It can be seen that there is an excellent co-incidence of our points with those of Beale [48]. Thepoints of Beale [48] are based on the very accurate phe-nomenological FSS scheme using a strip geometry [47],whereas our points are obtained via the present simul-taneous FSS analysis, based on lattices with linear sizes L = 20 − et al. [49] are based inmuch smaller lattices of the two-parameter WL sampling(linear sizes L ≤ .
95 ourphase diagram point does not agree, within errors, withthat of Beale and for this reason our estimation with arather generous error bars (shown also on the panel) hasbeen illustrated in Fig. 1(b). The rest of Table I con-tains our estimates for the random-bond versions. Onemay note from this table that for small values of ∆ forinstance ∆ = 1 . r = 0 . / . r = 0 . / .
25, and r = 0 . / . . IV. PHASE TRANSITIONS OF THE PURE ANDRANDOM-BOND 2D BLUME-CAPEL MODELSA. Strong violation of universality: Theex-first-order regime
As pointed out in the introduction, in our recent inves-tigation of the 2d random-bond BC model [37] we founddramatically different critical behaviors of the second-order phase transitions emerging from the first- andsecond-order regimes of the pure model. Namely, dif-ferent sets of critical exponents on two segments of thesame critical line appeared to describe the two regimes:the still-second-order and ex-first-order regimes. Thestudy in Ref. [37] was carried out for two values ofthe crystal-field coupling corresponding to the second-order (∆ = 1) and first-order (∆ = 1 . TABLE I: Transition temperatures of the pure and random-bond 2d BC model obtained in this paper. Second column fromreference [48], third and last entries of third and fifth columns from Ref. [37].∆ /J k B T /J
Pure RandomRef. [48] r = 0 . / . r = 0 . / . r = 0 . / .
40 1.695 1.693(3) 1.674(2)0.5 1.567 1.564(3) 1.547(2)1 1.398 1.398(2) 1.381(1)1.2 1.277(3)1.4 1.184(3)1.5 1.150 1.151(1) 1.149(1) 1.144(2) 1.131(2)1.6 1.084(1) 1.071(3)1.7 1.005(1)1.75 0.958(1) 0.960(2)1.8 0.908(1) 0.917(3)1.9 0.769(1) 0.774(2) 0.786(4)1.95 0.650 0.659(2) 0.702(3)1.975 0.574(2) 0.626(2) the pure model and for the random version the disor-der strength r = 0 . / .
25 was chosen in both cases.The strong violation of universality observed appearedto be the result of the softening of the first-order tran-sition due to bond-randomness. Specifically, it was con-cluded that the new strong-disorder universality class iswell described by a correlation length exponent in therange ν = 1 . − . γ/ν and β/ν very close to the Ising values 1 .
75 and 0 . . . r = 0 . / .
15 and r = 0 . / .
25. Thesaturation of the specific heat is very clear in both casesof the disorder strength.It is of interest to point out here that these findings forthe strong-disorder universality class appear to be fullycompatible with the classification of phase transitions indisordered systems proposed recently by Wu [55]. Ac-cording to this classification the strong-disorder transi-tion is expected to be inhomogeneous and percolativewith an expected exponent of the order ν = 1 .
34 [56].Furthermore, it has been suggested to us by Wu [57],that the strong lack of self-averaging of this transitionstems from the above properties. This violation of self- averaging, together with the strong finite-size effects,make the systematic MC approach of the strong-disorderregime very demanding, if not impractical. On the otherhand, the weak regime (or Ising universality regime) suf-fers a much weaker lack of self-averaging, by at least afactor of ∼
12 [37], and a smooth behavior is observedat moderate lattice sizes. Thus, aiming here to observe,even approximately, the extent of the involved universal-ity classes we carried out our study at moderate valuesof the crystal field and disorder, and found a behaviorquite convincing from which the frontier of the stronguniversality class can be estimated by observing the dis-appearance of the expected 2d random Ising universalityclass.
B. Pure and random-bond 2d BC model: Range ofuniversality with the 2d Ising model
Let us now proceed with the analysis of our numericaldata for the disorder strengths and crystal fields givenin Table I and observe and contrast their FSS behav-ior with that of the pure model. Starting this com-parative study with the FSS of the specific heat max-ima (using for the random-bond version at the strength r = 0 . / .
25 the corresponding quantity averaged overdisorder, i.e. [ C ] ∗ av ), we present in Fig. 4 fitting at-tempts for the same range of ∆ for the pure model[Fig. 4(a)] and the random-bond version [Fig. 4(b)]. Asindicated by the scales in the x-axis and the functionsin the corresponding panels, the expected Ising logarith-mic divergence has been assumed for the pure model C ∗ = C + C ln L , whereas the double-logarithmic di-vergence [ C ] ∗ av = C + C ln(ln L ) has been assumed forthe random version. Although, it is very difficult to ir-
10 1000.1110100 S p ec i f i c H ea t L Pure r = 0.85 / 1.15 r = 0.75 / 1.25
FIG. 3: Behavior of the random-bond 2d BC model at∆ = 1 .
975 from Ref. [37]. Illustration of the divergence ofthe specific heat of the pure model (first-order regime) andthe clear saturation of the specific heat for the random-bond(open symbols) 2d BC model for two disorder strengths in alog-log scale. refutably distinguish between a double-logarithmic diver-gence and a very weak power-law divergence, the theoret-ically well-grounded double-logarithmic scenario appliesvery well [37] and this fact can be observed also now formore values of ∆ in Fig. 4(b). There are some furtherfeatures one can observe from Fig. 4. First, for the puremodel, and in the range ∆ > .
9, we observe a sud-den change in the behavior of the specific heat peaks aswe approach the tricritical point which is apparently astrong cross-over effect. For the random version and thesame values of ∆ no such strong effects are noticeableand most probably the general softening effects of bondrandomness extends also to the expected cross-over phe-nomena between the two different universality classes ofthe random 2d BC model. From Fig. 4(b) the slopes ofthe double-logarithmic fittings appear to obey a rathersensible decreasing tendency from which we try in thenext Section to locate the frontier of strong-disorder uni-versality class.A second interesting comparison follows now in Fig. 5,where again Fig. 5(a) presents the FSS of the suscepti-bility maxima for the pure model, whereas Fig. 5(b) cor-responds to the random-bond version at r = 0 . / . γ/ν , estimated by the simple power law χ ∗ ∼ L γ/ν . These large fluctuations are a known pe-culiarity of the pure model near tricriticality [48]. Theeffective value of the exponent is now closer to the ex-pected value at the tricritical point γ/ν = 1 . γ/ν = 1 . .
95 in Fig. 5(a). Note here thatin Fig. 5(a) the simultaneous fitting is applied to the firstsix values of ∆ = 0 − .
9, whereas in Fig. 5(b) all values [ C ] * a v ln (ln L) [C] *av = C +C ln(ln L) (b) Pure Model C * = C +C ln L (a) C * ln L = 1 = 1.5 = 1.75 = 1.9 = 1.95 FIG. 4: FSS of the specific heat maxima for the pure andrandom-bond 2d BC model at the same values of ∆. (a) Puremodel: illustration of linear fittings assuming the expectedIsing logarithmic divergence. (b) Random model: illustrationof linear fittings assuming a double logarithmic divergence.Note the steady fall of the double logarithmic amplitude C . of ∆ = 0 − .
95 are used. However, even so, the com-parison of the two panels of Fig. 5, points out the muchbetter limiting behavior of the random version towardsthe expected Ising value γ/ν = 1 .
75. Therefore, we mayconvincingly suggest that, for moderate values of crystalfield and disorder, the Ising universality scenario (withpossible logarithmic corrections) is well obtained in thedisordered case.Figures 6 and 7 illustrate further alternative routes,used commonly in traditional FSS analysis, that pro-vide clear evidence to the above suggestion namely thatthe weak-disorder version belong to the 2d Ising classfor suitable moderate values of disorder and crystal-fieldcoupling. Noteworthy is the fact that the estimation ofthe critical exponents via the traditional FSS, such asthat shown in Figs. 5 - 7 for the random version, yieldsestimates very close to the expected values, i.e. the ex-ponents of the 2d Ising model. On the other hand, aspointed also earlier, for the pure model, as we approachthe tricritical point (∆ = 1 .
95) the effective exponentsremind more those expected at the tricritical point thanthose of the 2d Ising model. This is particularly true forthe correlation length’s exponent estimated in Fig. 1(b)by the simultaneous fitting of the six pseudocritical tem-
20 40 60 80 100075150225300375450 20 40 60 80 10004080120160200240 = 1.95 / = 1.44(7) / = 1.77(3) (a)
Pure Model * L = 0 = 0.5 = 1 = 1.5 = 1.75 = 1.9 (b) Random Model: r = 0.75 / 1.25 [] * a v L / = 1.74(2) FIG. 5: FSS behavior of the susceptibility maxima for thepure and random-bond 2d BC model at the same values of∆. (a) Pure model: simultaneous fitting to a simple powerlaw for the first five values of ∆ and a separate fitting closeto the tricritical point for ∆ = 1 .
95. (b) Random model:simultaneous fitting of the averaged susceptibility peaks for allvalues of ∆. Note the better behavior for the random modeland the improved estimation of the exponent ratio γ/ν . peratures. As can be seem from this figure the estimatefor the exponent ν , in the case ∆ = 1 .
95, is closer to theexpected tricritical value 40 /
77 = 0 . · · · [48], than tothe 2d Ising value ν = 1. V. MULTICRITICAL POINTS AND THESTRONG DISORDER REGIMEA. A novel estimation of multicritical points
From Fig. 4(a) one can observe the expected Ising loga-rithmic divergence of the specific heat maxima. Avoidingthe value ∆ = 1 .
95, which suffers from strong cross-overeffects, we attempted to estimate the tricritical value ofthe crystal field by fitting the decreasing logarithmic am-plitudes C to a suitable power law, as shown in Fig. 8.This may appear a naive or questionable idea, since thebehavior of specific heat data is the Achilles’ heel ofFSS analysis. Yet, Fig. 8 shows that besides the largefluctuations (errors) in logarithmic amplitudes C , onecould approximately estimate the tricritical crystal field
20 40 60 80 100050100150200250 20 40 60 80 1000.480.500.520.540.560.580.600.62 = 0.5 ; r = 0.75 / 1.25 [ l n < M n > / K ] * a v L n = 1 n = 2 n = 4 = 1.02(3) (a)(b) / = 0.127(3) = 0.5 ; r = 0.75 / 1.25 [ M ] a v ( T = T c ) L FIG. 6: (a) FSS analysis of the three logarithmic derivatives( n = 1 , ,
4) of the order parameter with respect to tem-perature for a particular crystal field ∆ = 0 . r = 0 . / .
25. (b) Traditional FSS analysis of theorder parameter at the estimated critical temperature for thesame value of crystal field and disorder strength, as in panel(a). ∆ t ≈ . r = 0 . / .
25, are better-matched to the double-logarithmic divergence than thecorresponding specific heat maxima of the pure modelto a simple logarithmic divergence. Therefore, it appearsrealistic to try to obtain the multicritical point (more pre-cisely the value of the crystal field ∆ sd where we expectthe emergence of strong-disorder regime at a given dis-order strength), where the 2d random-bond (BC) Isinguniversality class meets the strong-disorder universalityclass, by fitting the decreasing double-logarithmic am-plitudes to a similar power law. Figure 9 presents nowthese fittings and also the estimated multicritical valuesof the crystal field ∆ sd for the three disorder strengthsconsidered in this paper. The statistical errors of the cor-responding double logarithmic amplitudes C are seen tobe quite smaller, compared to the corresponding ampli-tudes of the pure model, and the estimated values ofthe values of the multicritical field, shown in the panel ofFig. 9, appear more convincing. From their general trendwe observe that as we increase the disorder strength the
20 40 60 80 1005101520253035 (1- ) / = 0.872(6) / = 0.128(6) r = 0.75 / 1.25 [ < | M | > / K ] * a v L = 0 = 0.5 = 1 FIG. 7: Simultaneous fitting to a simple power law of the av-eraged peaks corresponding to the absolute order-parameterderivative for three values of ∆ and disorder strength r = 0 . − β ) /ν and as indicated in the panel an alternative estimate for theexponent ratio β/ν , by assuming ν = 1. frontier of the strong-disorder universality class moves tolower values of the crystal field, namely: ( r = 0 . / . sd = 1 . r = 0 . / .
25, ∆ sd = 1 . r = 0 . / .
4, ∆ sd = 1 . sd approaches ∆ t , as expected. B. Strong disorder regime: The case ∆ = 1 . , r = 0 . / . and general observations As pointed out earlier, the strong lack of self-averaging,together with possible strong finite-size effects, makethe MC approach to the strong-disorder regime a verydifficult task. The self-averaging properties along thetwo segments (ex-first-order and still-second-order) of thecritical line were observed and discussed in Ref. [37]. Itwas shown in this paper that the usual finite-size mea-sure [59] of relative variance R X = V X / [ X ] av , where V X = [ X ] av − [ X ] av (and X = χ ∗ is the susceptibil-ity maxima), exhibits lack of self-averaging in both casesof marginal second-order transition and the transition inthe strong-disorder or ex-first-order regime. In particu-lar, it was shown [37] that the case studied (∆ = 1 . r = 0 . / .
25 by a fac-tor of ∼
12. Therefore, phase diagram points very closeto the frontier of the strong-disorder regime may be theworst cases to study, because, besides the very strong t = 1.96(1) C Pure Model C = const ( t - ) x FIG. 8: An approximate estimation of the tricritical value ofthe crystal field by fitting the decreasing logarithmic ampli-tudes C of the pure model at suitable values of the crystal-field coupling to a suitable power law shown in the figure. lack of self-averaging one may also expect large finite-sizeand cross-over effects. However, such cases are elucida-tory not only for observing the intrinsic difficulties butalso for giving us the opportunity to reflect on possiblelinks with basic properties of the system as for instancethe ground-state structure.Thus, Fig. 10 illustrates two important characteristicsof the case ∆ = 1 . r = 0 . / .
5. In particular,Fig. 10(a) shows the huge sample-to-sample fluctuationsin all the pseudocritical temperatures used in this pa-per. The simulation here was extended to larger lattices( L = 20 − C ] ∗ av ) and the finite-size behavior of the average of in-dividual maxima ([ C ∗ ] av ), which as shown suffer largesample-to-sample fluctuations. The behavior here resem-bles in many aspects the well known and still challeng-ing specific heat behavior of the 3d random-field Isingmodel [42]. However, a very clear tendency of [ C ] ∗ av for a saturating behavior is observed in Fig. 10(b) andwe may speculate that this saturating behavior is thecorrect asymptotic behavior for both maxima shown inFig. 10(b), although the behavior of [ C ∗ ] av will settledown only in very large lattice sizes, when the influenceof the strong finite-size effects on the individual maximawill diminish.The above illustrations open the possibility that thecase ∆ = 1 . r = 0 . / . C = const ( sd - ) x C r = 0.9 / 1.1: sd = 1.963(8) r = 0.75 / 1.25: sd = 1.955(5) r = 0.6 / 1.4: sd = 1.879(12) FIG. 9: Estimation of multicritical values of the crystal field,where the 2d random-bond (BC) Ising universality class meetsthe strong-disorder universality class. The decreasing doublelogarithmic amplitudes C of the random version has beenfitted to a power law shown and the estimates for the threedisorder strengths are shown in the panel. aspects of the ground-state structure of the 2d random-bond BC model. From Fig. 11, reproduced here fromRef. [58], one can see that approximately at this point(∆ = 1 . , r = 0 . / . s i = 0) as we increase the disorder strength. In the pres-ence of bond randomness the competition between theferromagnetic interactions with the crystal field results ina destabilization of the ferromagnetic ground state. De-pending on the realization, weak clusters exist in T = 0and their points are frozen in the s i = 0 state. This isan interesting subject, which is presently under furtherconsideration in both 2d and 3d by the present authors.The novel behavior illustrated in Fig. 10 appears nowas a consequence of the onset of the unsaturated groundstate at (∆ = 1 . , r = 0 . / .
5) which is thus related withthe critical behavior of the 2d random-bond BC model.The presented calculation of the ground states has becarried out in polynomially bounded computing time bymapping the system into a network and searching for aminimum cut by using a maximum flow algorithm (seefor instance Ref. [60]) and can be easily extended to largelattices and also to the 3d BC random-bond model.
VI. CONCLUSIONS
By carrying out an extensive two-stage Wang-Landauentropic sampling of both the pure and the random-bond 2d Blume-Capel model we have produced phasediagram points for several disorder strengths. Also forthe pure model we found an excellent coincidence ofour points with those of Beale [48]. For a particulardisorder strength ( r = 0 . / .
25) we found that, as aresult of the enhancement of ferromagnetic order, the
20 40 60 80 100 120 1401.21.31.41.51.61.71.8 20 40 60 80 100 120 1401.101.151.201.251.30 (b) = 1.5 ; r = 0.5 / 1.5 [ C * ] a v ; [ C ] * a v L [C * ] av [C] *av = 1.5 ; r = 0.5 / 1.5 T [ Z * ] a v L Z = C Z = Z = <|M|> / K Z = ln
5. (a) Behavior of the averaged pseudocriticaltemperatures, corresponding to individual maxima, with anillustration of the huge sample-to-sample fluctuations. (b)Behavior of the maxima of the averaged specific heat curves([ C ] ∗ av ) and the average of individual maxima ([ C ∗ ] av ) withtheir large sample-to-sample fluctuations. phase diagram of the random version crosses that of thepure model at approximately the point (∆ cross , T cross ) =(1 . , . . r = 0 . / . r = J / J [ M ( T = )] a v L = 30
FIG. 11: Ground-state behavior of the order parameter of the2d random-bond BC model versus r for various values of ∆averaged over 250 disorder realizations. Blume-Capel model. It has been pointed out that thebehavior of this system is very interesting and in partic-ular the strong-disorder regime may include many furtherchallenges and open problems that, at the moment, arenot fully understood. In our opinion, this is a rathercomplex subject deserving further research.
Acknowledgments