Stability of the Griffiths phase in the 2D Potts model with correlated disorder

A Griffiths phase has recently been observed by Monte Carlo simulations in the 2D $q$-state Potts model with strongly correlated quenched random couplings. In particular, the magnetic susceptibility was shown to diverge algebraically with the lattice size in a broad range of temperatures. However, only relatively small lattice sizes could be considered so one can wonder whether this Griffiths phase will not shrink and collapse into a single point, the critical point, as the lattice size is increased to much larger values. In this paper, the 2D eight-state Potts model is numerically studied for four different disorder correlations. It is shown that the Griffiths phase cannot be explained as a simple spreading of local transition temperatures caused by disorder fluctuations. As a consequence, the vanishing of the latter in the thermodynamic limit does not necessarily imply the collapse of the Griffiths phase into a single point. In contrast, the width of the Griffiths phase is controlled by the disorder strength. However, for disorder correlations decaying slower than $1/r$, no cross-over to a more usual critical behavior could be observed as this strength is tuned to weaker values.


Introduction
The influence of disorder on phase transitions and critical phenomena has attracted a considerable interest in the last decades. In the absence of frustration, it is now well established that a first-order phase transition is smoothed by the introduction of randomness and can be made continuous at large enough disorder strength [1]. In 2D, an infinitesimal disorder is sufficient to remove any discontinuity [2][3][4]. When the pure system undergoes already a continuous phase transition, the Harris criterion predicts that the universality class of the pure model will be affected by disorder if the specific heat diverges, i.e. if the critical exponent α is positive [5]. In this context, the q-state Potts model has been a useful toy model, because it displays a rich phase diagram involving two lines of respectively first and second-order phase transition. Along the latter, the universality class depends on the number of states q. On the practical side, efficient Monte Carlo and transfer matrix algorithms exist for this model and Conformal Invariance can be used in 2D in combination with Renormalisation Group (RG).
In comparison, correlated disorder was much less studied. Nevertheless in some experimental situations, impurities cannot be considered as uncorrelated. This is in particular the case when they carry an electric charge or a magnetic moment and are coupled via an electromagnetic interaction. On the theoretical side, Weinrib and Halperin studied the φ 4 model with a random mass and showed that a new RG fixed point, distinct from the random and the pure ones, emerges in the phase diagram when the correlations of this mass decay algebraically [6]. For a sufficiently slow decay of these disorder correlations, the new fixed point becomes stable. Denoting a the exponent of the algebraic decay of disorder correlations, the perturbation is relevant when a < d if the correlation length exponent ν of the pure model satisfies the inequality ν < 2/a. At the new fixed point, correlated disorder is marginally irrelevant, which implies that that ν = 2/a [7]. The magnetic exponent η remains small compared to ǫ = 4 − d. Even though still controversial, these predictions were confirmed by Monte Carlo simulations of the 3D Ising model with a = 2 [8,9].
We recently studied by large-scale Monte Carlo simulations the influence of correlated couplings on the 2D Potts model [10,11]. Like in the absence of disorder correlation, the first-order phase transition, occurring for the pure model when q > 4, was shown to be smoothed and replaced by a continuous transition. However, the new universality class was shown to be q-independent, a feature shared by the strong-disorder fixed point of the q-state Potts model with a layered McCoy-Wu-like disorder. This result is remarkably different from the continuous increase of the magnetic scaling dimension x σ observed for the Potts model with an uncorrelated disorder. More intriguing is the fact that the phase diagram displays a Griffiths phase, as in the McCoy-Wu model, where the magnetic susceptibility diverges with the lattice size. Interestingly, such a phase has been predicted by Weinrib and Halperin, but only above the upper critical dimension d c = 4. Finally, the hyperscaling relation γ ν = d − 2 β ν was observed to be broken in the Griffiths phase, as a result of large disorder fluctuations.
However, these observations were made for finite systems so one cannot exclude the possibility that, at much larger lattice sizes, the Griffiths phase collapse into a single point, the critical point, where the hyperscaling relation would be restored. Moreover, the estimate of the correlation length exponent ν is incompatible with Weinrib-Halperin exact result ν = 2/a, which substantiates the idea that the lattice sizes considered could be too small and that a cross-over would be observed at much larger lattice sizes. On the other hand, no significant evolution of the Griffiths phase could be observed in the range of lattice sizes studied [11]. Moreover, the conspiracy of two amplitudes that leads to the violation of the hyperscaling relation is well verified and no sign of deviation at large lattice sizes is observed.
Since larger lattice sizes are not accessible by Monte Carlo simulations, we turn our attention in this work to larger exponents a of the disorder correlations. The fact that Weinrib-Halperin predictions were confirmed by Monte Carlo simulations of the 3D Ising model with a = 2 could indicate that finite-size effects are weaker for larger values of a. In refs [10,11], only small values of a were considered because disorder configurations were generated by simulating an auxiliary Ashkin-Teller model on a selfdual line where its critical exponents are known exactly. The polarisation density was then used to construct the couplings J i j of the Potts model. Disorder correlations correspond therefore to the polarisationpolarisation correlations of the auxiliary Ashkin-Teller model. When moving along the self-dual line, only exponents in the range a ∈ [1/4; 3/4] can be obtained.
In this work, we present new data for disorder correlation exponents a = 1/3 and 2/3 obtained by using an auxiliary Ashkin-Teller model. In order to investigate the possible existence of a cross-over towards the Weinrib-Halperin fixed point, we considered also the values a ≃ 1.036 and a = 2 obtained using the 3D and 4D Ising models as auxiliary models to generate the disorder configurations. In the first section, details about the models and the Monte Carlo simulations are given. In the second section, the behaviour of the magnetic susceptibilityχ is discussed. As already observed in Ref. [11],χ diverge algebraically with the lattice size in a broad interval of temperatures, identified as a Griffiths phase, when a is sufficiently small. A simple explanation of this phenomena is to assume that disorder fluctuations induce a spreading of local transition temperatures. Because these fluctuations vanish as L −a/2 , this would imply that a single peak would be recovered at large lattice sizes. Moreover, the smaller the exponent a and the larger the lattice sizes needed to observe a single peak. In the second section, numerical evidence is given that disorder fluctuations are not sufficient to explain the observed Griffiths phase, and therefore, that the latter phase cannot be expected to collapse as L −a/2 . In the third section, the possibility of a cross-over controlled by the amplitude of disorder correlations is considered. These amplitudes are compared for the different disorder correlations considered and, then different disorder strengths are studied. Finally, conclusions follow.

Models and simulation
The classical q-state Potts model is the lattice spin model defined by the Hamiltonian [12,13] where the spin σ i takes q possible values and is located on the i -th node of the lattice. The sum extends over all pairs (i , j ) of nearest neighbours on the lattice. In the following, the Potts model will be considered on the square lattice. As mentioned in the introduction, the phase transition is continuous for q ≤ 4 and of first-order when q > 4. We will restrict ourselves to the case q = 8, i.e. a point in the regime of first-order transition. Disorder is now introduced as bond-dependent random exchange couplings J i j .
The Hamiltonian becomes The spatial correlations between these couplings is assumed to decay algebraically with an exponent a at large distance: For convenience, we will restrict ourselves in the following to a binary coupling distribution, i.e. J i j = J 1 or J 2 . The presence of disorder correlations does not affect the self-duality condition of the random Potts model. Imposing J 1 and J 2 to be equiprobable and self-dual of each other, the self-dual line is given by the condition [14] e βJ 1 − 1 e βJ 2 − 1 = q.
The coupling configurations are generated by independent Monte Carlo simulations of two auxiliary models: the Ising and Ashkin-Teller models. The former is defined by the Hamiltonian and is equivalent to the q = 2 Potts model. It is well known that this model undergoes a second-order phase transition in any dimension d > 1. We considered hypercubic lattices of dimension d = 3 and d = 4. A few thousand spin configurations are generated at the critical point, corresponding to βJ I M ≃ 0.221655 for d = 3 [15] and βJ I M ≃ 0.149694 for d = 4 [16]. For each spin configuration, a two-dimensional section is cut and random couplings for the 2D Potts model are constructed as for each pair (i , j ) of nearest neighbours in the 2D section. Note that, at any site i , two couplings, in two different directions, are identical. By construction, disorder correlation functions J i j J kl − J i j J kl decay as the spin-spin correlation functions of the auxiliary Ising model. Therefore, the decay is algebraic at large distances with an exponent a = 2β/ν ≃ 1.036 for the 3D Ising model [15] and a = 2 for the 4D Ising model. Note that, in the second case, the exponent a is equal to the dimension d = 2 of the Potts model.
Therefore, according to Weinrib and Halperin, disorder correlations are expected to be irrelevant and the system falls into the same universality class as the Potts model with independent random couplings. The second auxiliary model is the 2D Ashkin-Teller model defined by the Hamiltonian [17,18] and corresponding to two Ising models coupled by their energy densities. On the square lattice, the model is self-dual along the line of the phase diagram given by e −2K AT = sinh 2J AT . Thanks to a mapping onto the eight-vertex model, the critical exponents are known exactly along this line. The random couplings for the Potts model are constructed from the polarisation density as The disorder correlations therefore decay as the polarisation-polarisation correlation functions of the auxiliary Ashkin-Teller model. In this work, we considered two points on the self-dual line of the Ashkin-Teller model (y = 0.50 and y = 1.25 in the language of the eight-vertex model) corresponding to exponents a = 1/3 and a = 2/3.
The above-described spin models were simulated using Monte Carlo cluster algorithms to reduce the critical slowing-down. For the Ising and Potts models, the Swendsen-Wang algorithm was employed [19]. The Ashkin-Teller was simulated using a cluster algorithm introduced by Wiseman and Domany [20,21]. ?????-3

Griffiths phase and disorder fluctuations
The magnetic susceptibility χ of a finite system undergoing a continuous phase transition in the thermodynamic limit is expected to display a peak whose maximum diverges with the lattice size L as L γ/ν . The location of this maximum goes towards the critical temperature T c in the limit of an infinite system.
A very different situation was observed in the 2D Potts model with strongly correlated disorder [11]. As can be seen on figure 1, two peaks are present for a = 1/3 and 2/3. The data show an algebraic increase of the average magnetic susceptibility for all temperatures between these two peaks. For this reason, this region was conjectured to be a Griffiths phase, similar to the one observed in the McCoy-Wu model. The absence of any evolution of the location of the two peaks was reported in the case a = 0.4. In contrast, figure 1 shows a slow evolution in the case a = 2/3. Since only lattice sizes up to L = 128 were studied, the possibility of a collapse of the Griffiths phase into a single point in the thermodynamic limit cannot be excluded. Moreover, such a collapse is even more clearly seen on figure 1 for a ≃ 1.036. Two peaks are still visible but they tend to come closer when the lattice size is increased. It seems natural in this case to assume that the two peaks will merge into a single one at larger lattice sizes. For disorder correlations with a faster decay a = 2, only one peak is observed (Fig. 1) and its location tends towards the critical value β c = 1, expected from self-duality arguments. As mentioned in the introduction, it may be assumed that the width of the Griffiths phase is due to large disorder fluctuations. It seems indeed natural to assume that the first peak is caused by the ferromagnetic ordering of large clusters with a high concentration of weak bonds J 2 while the second one corresponds to clusters of strong bonds J 1 . Such large clusters are more probable when disorder correlations decay slowly. In the following, disorder fluctuations will be compared for the different values ?????-4 of a considered. To be more specific, consider the general case of a lattice model with an energy density denoted ǫ i j = ǫ(σ i , σ j ) on the edge between the spins on sites i and j . The weak disorder limit of the partition function can be calculated using the replica trick: Introducing the interaction energy ǫ α i j = ǫ(σ α i , σ α j ) between the two spins σ α i and σ α j of the α-th replica, the partition function of n replicas reads The first contribution of disorder to the partition function involves the correlations J i j J kl − J i j J kl , and is obviously a function of a. In order to characterise the disorder strength by a scalar, we considered the sum of these correlations, which also corresponds to the fluctuations of the couplings:  The variance of the average coupling ∆J 2 is plotted on figure 2 versus the lattice size in the four cases a = 1/3, 2/3, 1.036 and 2. Note that in the last two cases (a ≃ 1.036 and a = 2), only the magnetisation ?????-5 in the two-dimensional section that was used to construct the exchange couplings J i j is considered. As expected, an algebraic decay with an exponent compatible with a/2 is observed. On figure 1, the collapse of the two peaks of the magnetic susceptibility is observed for a ≃ 1.036 for lattice sizes L ∼ O(10 2 ) when a ≃ 1.036. According to figure 2, this corresponds to disorder fluctuations of order ∆J 2 ≃ 0.06. For a = 1/3 and 2/3, none of the lattice sizes that were considered correspond to so small disorder fluctuations. Indeed, ∆J 2 = 0.159(4) when a = 1/3 for the largest lattice size L = 128 and ∆J 2 = 0.090(3) when a = 2/3. This strengthens the idea that the collapse will be observed for larger lattice sizes for a = 1/3 and 2/3. Using the scaling law (3.4), one can even predict these sizes to be of the order of L * ≃ 128(0.06/0.16) −2/a ≃ 44, 000 for a = 1/3 and L * ≃ 128(0.06/0.09) −2/a ≃ 432 for a = 2/3. On the other hand, disorder fluctuations are small for a = 2 (∆J 2 = 0.0618(4) already for L = 24), smaller than for a ≃ 1.036 with lattice sizes L ∼ O(10 2 ).
These results do not depend on the quantity used to measure disorder fluctuations. The scaling law (3.4) suggests to use the order parameter, polarisation |σ i τ i | or magnetisation |σ i |, of the auxiliary models as an alternative measure of the fluctuations of the couplings. This quantity will be denoted ∆J 1 in the following. ∆J 1 and ∆J 2 give essentially the same information and, as can be seen on figure 2, take sensibly the same value, but ∆J 1 presents the advantage of being more stable numerically. More surprising is the fact that the same conclusions can be drawn from the second contribution of disorder to the partition function (3.2). Expanding further, the next term will involve the connected four-point correlation function J i J j J k J l c of disorder. This quantity was assumed to be irrelevant by Weinrib and Halperin. We considered the fourth-order cumulant As can be seen on figure 2, no qualitative difference between the three quantities ∆J 1 , ∆J 2 and ∆J 4 is observed.
However, there are small differences between the four cases a = 1/3, 2/3, 1.036 and a that cannot be explained only in terms of disorder fluctuations. The value of ∆J for the largest lattice size L = 128 at a = 2/3 is close to the one estimated for L = 48 at a ≃ 1.036. Therefore, the average susceptibility should look qualitatively the same for a = 2/3 at L = 128 and for a ≃ 1.036 at L = 48. It is not clear that it is indeed the case on figure 2. Moreover, a nice collapse of the magnetic susceptibility is observed at large β for a = 1/3 and 2/3 while it is not case for a ≃ 1.036 and 2. Stronger statements might be formulated by comparing thermodynamic quantities displaying universal properties. The natural candidate is the 4th-order Binder cumulant whose value at the intersection of two curves with respect to temperature is expected to be universal in the thermodynamic limit. A notable difference between the different values of a is that the crossing points occur for inverse temperatures β well below β c = 1 when a = 1/3 and 2/3 and very close to β c = 1 when a ≃ 1.036 and 2. Unfortunately, the error bars are large and do not allow to be conclusive.
Another quantity displaying universal properties is the ratio [22] that measures the sample-to-sample fluctuations of magnetisation. Outside of a critical point, all disorder realisations are expected to lead to the same average magnetisation 〈m〉 in the thermodynamic limit. Therefore, the ratio R m vanishes as L → +∞ and magnetisation is said to be self-averaging. This is no longer true at a fixed point where disorder is relevant. In this case, R m goes towards a finite value in the thermodynamic limit. This limit is expected to be a universal quantity [23]. Numerical data for this ratio R m are plotted on figure 3. Two distinct behaviours are observed. For a = 1/3 and a = 2/3, R m displays a peak in the paramagnetic phase (small β = 1/k B T ), followed by a broad shouldering. The latter extends over a range of temperatures which roughly corresponds to the range between the two peaks of the average magnetic susceptibility (see figure 1). Interestingly, the estimates of R m at any temperature in this shouldering are compatible, within error bars, for all lattices sizes L ∈ [32; 128]. Unless a sudden decay of R m occurs at much larger lattices sizes, we are led to the conclusion that magnetisation is a non-selfaveraging quantity in the whole range of temperatures between the two peaks of the susceptibility. This conclusion is consistent with the assumption of the existence of a Griffiths phase. On the other hand, for a ≃ 1.036 and a = 2, the peak in the paramagnetic phase is softer and is not followed by a shouldering but by a monotonous decay. More interesting is the fact that the curves corresponding to different lattice sizes cross each other at a single point, close to the self-dual point β c = 1. This is consistent with the existence of a unique critical point at β c = 1. Would it be possible that, in the case a = 2/3, the shouldering disappears at large lattice sizes to be replaced by a monotonous decay with a single crossing point for different lattice sizes? If the coupling fluctuations ∆J provides a measure of the width of the Griffiths phase as discussed above, it should also determine the range of temperatures around β c for which R m is finite and sizeindependent. Then the ratio R m should look similar for a = 2/3 at L = 128 and for a ≃ 1.036 at L = 48. This is definitely not the case on figure 3. Therefore, the Griffiths phase is not solely the consequence of disorder fluctuations and there is no reason to expect the Griffiths phase to collapse into a single point as L −a/2 .

Griffiths phase and disorder strength
All data presented in the previous section correspond to a disorder strength r = J 1 /J 2 = 7.5. Because the two peaks of the susceptibility were interpreted as the ordering of macroscopic clusters with a majority of strong, or weak, couplings, the disorder strength r controls the width of the Griffiths phase. One can therefore wonder whether disorder is not too strong in the cases a = 1/3 and 2/3 which implies that where the two amplitudes v and w are irrelevant scaling fields at the long-range random fixed point. When simulating a finite system with an amplitude w much larger (or much weaker) than the value w * taken at the fixed point, the critical behaviour may be affected by strong scaling corrections. Indeed, in the neighbourhood of the Weinrib-Halperin random fixed point, the free energy density can be assumed to scale under a dilatation with a scale factor b as: where γ/ν = 2y h − d. The scaling function F involves a cross-over length ℓ ∼ (w − w * ) −1/y w associated to disorder. The dominant finite-size scaling behaviour χ ∼ L γ/ν will be hidden by scaling corrections if L ≪ ℓ.
In the previous section, the fluctuations of the average coupling have been compared for different exponents a. To compare now the amplitudes w of disorder correlations, note that integrating out disorder correlations leads on one hand to while on the other hand, the same quantity is equal to (∆J 2 ) 2 according to equation (3.3). The amplitude w can therefore be recomputed as w ≃ (2 − a)(∆J 2 ) 2 L a . This estimate is plotted on figure 4. Note that the amplitude w is not plotted for a = 2 because the definition is inappropriate in this case (the integration of the correlations involves a logarithm) and leads to w = 0. As can be seen on figure 4, the amplitude w does not evolve monotonously with a. This should not be a surprise because the couplings have been generated from different auxiliary models. The amplitude w for a ≃ 1.036 lies in between the amplitudes for a = 1/3 and a = 2/3. Therefore, the Griffiths phase and the small ν exponents reported in [11] for a = 1/3 and a = 2/3 cannot be explained as the result of strong scaling corrections. Indeed, if one assume that the amplitude w is close to w * in the case a ≃ 1.036, which would explain why the collapse of the two peaks ofχ is observed for reachable lattice sizes, one can conceive that the Griffiths phase is the result of too strong disorder correlations, i.e. w > w * , in the case a = 1/3. However, it is hard to understand how weak correlations, i.e. w < w * , would lead to a similar result for a = 2/3. The explanation in terms of a cross-over is therefore not supported by the numerical data. However, as predicted by Weinrib and Halperin for the φ 4 model, the amplitude w * at the fixed point can be a function of a. In the following, the effect of a change of the disorder strength, and therefore of w − w * , is studied in the case a = 2/3. Since F depends on the scaling variable L y w (w − w * ), tuning the amplitude w to come closer to w * is expected to be equivalent to increasing the lattice size L. Therefore, the cross-over, if any, should be observed when w is decreased. On figure 5, the magnetic susceptibility is plotted for two different disorder strengths r = J 1 /J 2 = 3.5 and r = 2. Comparing with figure 1 where the case r = 7.5 was plotted, it is clear that the width of the Griffiths phase is directly proportional to the disorder strength r , and therefore the amplitude w. The two peaks come closer as the disorder strength is reduced. However, even for r = 2, the magnetic susceptibility is still very different from what is observed on figure 1 for a ≃ 1.036 and r = 7.5. One can doubt that the two curves will look similar for an even weaker disorder in the case a = 2/3. Very probably, the two peaks of the magnetic susceptibility for ?????-9 a = 2/3 will collapse but only when approaching r = 1, i.e. the pure model.
The ratio R m = 〈m〉 2 /〈m〉 2 − 1 provides a stronger evidence that what is observed for a ≃ 1.036 is not what should be expected for a = 2/3 at weaker disorder. As discussed when the figure 1 was commented, the signature of the Griffiths phase is a size-independent ratio R m over a finite range of temperatures. In contrast, for a ≃ 1.036, the ratios R m computed at two different lattice sizes display a single crossing point at a temperature evolving towards the self-dual point β c = 1. As observed on figure 5, the diminution of the disorder strength induces a reduction of the width of the Griffiths phase, and, as expected, of the range of temperatures where the ratio R m appears to be size-independent. However, the ratio R m looks surprisingly similar, up to a temperature rescaling, at different disorder strengths. Even at a disorder strength r = 2, the behaviour is still very different from what is observed for a ≃ 1.036.

Conclusions
New Monte Carlo simulations of the 2D 8-state Potts model with a disorder involving algebraically decaying correlations C (r ) ∼ r −a with a = 1/3, 2/3, 1.036 and 2 are presented. While the analysis of the magnetic susceptibility does not allow to exclude the possibility of a collapse of the Griffiths phase into a single critical point, the study of the self-averaging ratio R m = 〈m〉 2 /〈m〉 2 −1 allows to be more conclusive.
Two different behaviours are indeed observed for a = 1/3 and 2/3 on one hand and a ≃ 1.036 and 2 on the other hand. The first case is compatible with the assumption of the existence of a Griffiths phase while in the second case, the signature of a single critical point is observed. These difference cannot be explained by larger disorder fluctuations in the first case. Moreover, if the width of the Griffiths phase depends on the disorder strength for a = 2/3, no cross-over towards the Weinrib-Halperin critical behaviour is observed at weak disorder. These numerical results call for a theoretical understanding of the precise mechanism behind this Griffiths phase. We hope that a theoretician will find them sufficiently surprising to get interested into this problem.