Emergent universal critical behavior of the 2D $N$-color Ashkin-Teller model in the presence of correlated disorder

We study the critical behavior of the 2D $N$-color Ashkin-Teller model in the presence of random bond disorder whose correlations decays with the distance $r$ as a power-law $r^{-a}$. We consider the case when the spins of different colors sitting at the same site are coupled by the same bond and map this problem onto the 2D system of $N/2$ flavors of interacting Dirac fermions in the presence of correlated disorder. Using renormalization group we show that for $N=2$, a"weakly universal"scaling behavior at the continuous transition becomes universal with new critical exponents. For $N>2$, the first-order phase transition is rounded by the correlated disorder and turns into a continuous one.


Introduction
Two-dimensional (2D) systems are of particular interest in studying phase transitions in condensed matter. On the one hand, this interest is constantly growing due to the progress in experimental techniques of producing and studying low-dimensional materials like graphene [1], 2D crystals [2], and ultrathin ferromagnetic films [3]. On the other hand, the scaling behavior of many models is much easier to analyze in two dimensions than in three dimensions, since the 2D conformal invariance much stronger restricts possible scenarios, at least in the absence of disorder [4]. Moreover, some of 2D models allow for an exact solution. The well-known examples include the 2D Ising model [5], the 2D ice-model (6vertex model) [6] and the Baxter model (symmetric 8-vertex model) [7]. The latter model is related to the so-called Ashkin-Teller model [8], which was introduced to describe cooperative phenomena in quaternary alloys [9], along a selfdual line. Both models can be described in terms of two Ising models coupled through their energy densities (i.e., by four-spin interaction): (1.1) The main feature of the model (1.1) is a continuous transition with the so-called "weak universality" [10]: contrary to the usual critical universality, the critical exponents of model (1.1) continuously depend on the coupling constant J 4 . In particular, the correlation length exponent is [7] 1 ν pure = 4 π arctan e 2J 4 /T c , (1.2) where the critical temperature of the pure system T c is given by 2J /T c = ln(1+ 2). The heat capacity exponent behaves as α pure ∼ J 4 , and thus, changes sign with J 4 . However, expressing the singular behavior of thermodynamic quantities near the critical point in terms of the inverse correlation length rather than of the reduced temperature, one finds that the critical exponents rescaled in this way are universal [11].
Generalizing the model (1.1) to an arbitrary N , one arrives at the so-called N -color Ashkin-Teller model [12]: 3) in which N Ising models are coupled pairwise through interaction J 4 . It reduces to the usual Ashkin-Teller (Baxter) model for N = 2. For N > 2, the transition properties of model (1.3) drastically depend on the sign of J 4 : the system undergoes a continuous phase transition for J 4 < 0, while for J 4 > 0, the transition is of first-order [12]. The effects of quenched disorder on phase transitions have been a hot topic of research for several decades (for review see e.g., [13][14][15]). For example, it is well-known that the critical exponents of a system undergoing a continuous phase transition may be modified by uncorrelated quenched impurities coupled to the local energy density (the so-called random-bond or random site disorder). In this case, the relevance of disorder can be predicted using the Harris criterion [16]: if the heat capacity exponent of the corresponding pure system is positive, α pure = 2−dν pure > 0, then the presence of weak uncorrelated disorder leads to a new critical behavior. Here, ν pure is the correlation length exponent of the pure system. According to the Harris criterion, the 2D Ising model corresponds to a marginal situation since the heat capacity exhibits only a logarithmic divergence in the vicinity of critical point, i.e., α pure = 0 [5]. Explicit calculations performed for the disordered 2D Ising model in [17][18][19][20][21] reveal that uncorrelated disorder modifies the logarithmic divergence to a double logarithmic behavior, while other power-law scaling laws acquire universal logarithmic corrections. The situation is quite different when the quenched disorder is correlated. According to the generalized Harris criterion [22], the Gaussian disorder, whose variance decays as a power law r −a , modifies the critical behavior of the pure system for a < d provided that it satisfies the inequality ν pure < 2/a. For a > d, the usual Harris criterion is recovered and the condition is replaced by ν pure < 2/d (see also [23,24]). For the 2D Ising model with long-range correlated disorder, this has been explicitly shown in [25] by mapping to the 2D Dirac fermions.
The effect of uncorrelated disorder on the continuous phase transition with "weak universality" exhibited by the model (1.1) has been considered in [26][27][28][29][30]. Since the heat capacity exponent α pure of the pure Baxter model is positive for J 4 > 0, one can expect that the critical behavior is modified by uncorrelated disorder for J 4 > 0. The renormalization group (RG) picture obtained using a mapping to fermions suggests [26,27] that for J 4 > 0, the "weakly universal" critical behavior of the Baxter model changes to that of the disordered 2D Ising model, e.g., the heat capacity exhibits a double logarithmic singularity. The numerical simulations of [28][29][30] support the relevance of disorder but quantitatively they are less conclusive. For instance, the numerical simulations of the random bond Ashkin-Teller model seem to require additional efforts due to large sample to sample fluctuations and are rather in favor of a logarithmic than a double logarithmic behavior of the heat capacity [29,30]. For J 4 < 0, the exponent α pure is also negative so that the Harris criterion naively suggests that the critical behavior is unaffected by uncorrelated disorder. The critical behavior deduced from the RG picture is, however, different due to a new vertex generated by the RG flow. For instance, the correlation length behavior becomes universal while the heat capacity remains finite. However, the precise behavior depends on the initial disorder, in particular, whether the disorder seen by the both coupled Ising models in equation (1.1) is correlated or not [26,27]. It seems rather striking that adding a weak short-range correlated (SR) or uncorrelated quenched disorder to the 2D N -color Ahkin-Teller (1.3) with N > 2 results in emergent critical behavior [31][32][33][34][35][36][37].
Indeed, the pure model exhibits a fluctuation-driven first-order transition characterized by runaway of the RG flow which is reversed by even weak uncorrelated disorder. The 2D three-color Ashkin-Teller model has been studied by means of large-scale Monte Carlo simulations in [38] and [39]. While the first early work excludes the possibility of continuous transition with universal scaling behavior, the second recent paper demonstrates that the first-order phase transition is rounded by the disorder and turns into a continuous one. The resulting transition seems to be in the disordered 2D Ising universality class. This agrees with perturbative RG predictions of [31][32][33][34][35][36][37].
In the present paper we study the effects of long-range correlated (LR) disorder with power-law decay

13603-2
of correlations on the phase transitions in the Baxter and N -color Ashkin-Teler models. The article is organized as follows. In section 2, we consider the formulation of the problem in terms of Dirac (complex) fermions following [25,40,41] and restrict our consideration to the case of the same disorder for all fermion flavors [i.e., the disorder potentials seen by different Ising components in the models (1.1) and (1.3) are completely correlated]. In section 2, we introduce a fermion representation for both models. In section 3 we briefly describe the renormalization scheme we use for both models while in section 4 we present the one-loop RG functions and the derived scaling behavior. Finally, we conclude in section 5.

Models and their fermion representation
In the vicinity of a critical point, the long-distance properties of the Baxter model (1.1) can be described using two Majorana (real) fermionic fields χ 1 and χ 2 with the action [27]: The corresponding action reads [40] Generalization of the action (2.2) to the N -color Ashkin-Teler model with even N is straightforward [42]: In the presence of random bond disorder which is completely correlated between different Ising colors (i.e., different fermion flavors), the mass can be written as m(r ) = m 0 + δm(r ), where δm(r ) is the local disorder strength. We assume that the disorder strength is a random Gaussian with the mean value δm(r ) = 0 and variance δm(r )δm(0) = g (r ), (2.4) where g (r ) ∼ δ(r ) and g (r ) ∼ r −a are for SR and LR disorder, respectively. We shall use dimensional regularization so that we have to generalize the problem to arbitrary d. To that end, we replace the Pauli matrices by the Clifford algebra represented by the matrices γ i , i = 1, . . . , d satisfying the corresponding anticommutation relations. To average the free energy of (2.3) over different disorder configurations, we employ the replica trick [43]. Introducing n replicas of the system (2.3) and averaging over Gaussian disorder distribution we arrive at the replicated effective action The properties of the original system with quenched disorder can be then obtained by taking the limit n → 0. It is convenient to fix the normalization of the disorder correlator (2.4) in Fourier space as where u 0 and v 0 are bare coupling constants corresponding to the SR and LR parts of disorder correlator, respectively.

RG description
We study the long-distance properties of (2.5) within the standard approach of field-theoretical RG technique [44]. Applying it one can calculate the correlation functions for the action (2.5) perturbatively in λ 0 , u 0 and v 0 . The integrals entering this perturbation series turn out to be UV divergent in the dimension we are interested in (d = 2). To make the theory finite we use the dimensional regularization [45,46] and compute all integrals in d = 2 − ε. Inspired by the works [22,41] we perform a double expansion in ε = 2 − d and δ = 2 − a so that all divergences are transformed into poles in ε and δ while the ratio ε/δ remains finite. We define the renormalized fields ψ i ,ψ i , mass m, and dimensionless coupling constants λ, u and v in such a way that all poles can be hidden in the renormalization factors Z ψ , Z m , Z λ , Z u and Z v leaving finite the correlation functions computed with the renormalized action where we have introduced a renormalization scale µ and the shortcut notation k := d d k, so that k 1 ,k 2 ,k 3 stands for the corresponding triple integral. Since the renormalized action is obtained from the bare one by the fields rescaling the bare and renormalized parameters are related by  To calculate the renormalization constants it suffices to renormalize the two-point vertex function Γ (2) and the four-point vertex function Γ (4) . We impose that they are finite at m = µ and find the renormalization constants using a minimal subtraction scheme [45,46]. To that end, it is convenient to split the four-point function in the clean Γ λ , SR disorder Γ u and LR disorder Γ v parts: Renormalization constants are determined from the condition that Γ (4) λ (0; m = µ), Γ (4) u (0; m = µ) and Γ (4) v (0; m = µ) are finite.
Since the bare vertex function does not depend on the renormalization scale µ, the renormalized vertex function satisfies the RG equation where we have introduced the RG functions Here, the subscript "0" stands for derivatives at fixed λ 0 , u 0 , v 0 and m 0 . The critical behavior, if present, should be controlled by a stable fixed point (FP) of the β-functions, which is defined as Stability of a given FP can be determined from the eigenvalues of the stability matrix The FP is stable provided that all the eigenvalues calculated at the FP (3.12) have negative real parts. Using coordinates of the stable FP we can calculate the critical exponent. For example, the correlation length exponent ν is given by (3.14) The heat capacity exponent is given by the hyperscaling relation, which in two dimensions reads as α = 2(1 − ν).

One-loop RG flow and the critical behavior
Applying the renormalization procedure described in the previous section we obtain the β-functions (4.3) and the γ function to one-loop order in the replica limit n → 0. We find that the flow equations corresponding to the βfunctions at the corresponding FP are summarized in table 1. Note that in two dimensions, i.e., for d = 2 (ε = 0), the first three FPs merge and coincide with the FP G which describes the critical behavior of N uncoupled clean Ising models up to logarithmic corrections. Let us analyze the RG flow for both the Baxter model (N = 2) and the N -color Ashkin-Teller model with N > 2 in three regimes: (i) Pure system (u 0 = v 0 = 0). For N = 2, the β λ -function (4.1) vanishes in d = 2 so that the model has a line of FPs parameterized by λ. This is not surprising since the model in this limit coincides with the O(2) Gross-Neveu model or the massive Thirring model. The β-function of the latter model is equal to zero identically leading to nonuniversal critical exponents [47]. This is consistent with the "weak universality" picture given by the exact solution of the Baxter model [7]. Indeed, according to equations (3.14) and (3.15), the correlation length and the singular part of the heat capacity behave as where λ 0 is the initial value of the dimensionless coupling constant λ. Expanding the exponent (1.2) in small J 4 we find the relation between parameters of the continuous and the lattice models: λ 0 ≈ 2J 4 /(πT c ). For N > 2, the RG flow given by the β λ -function (4.1) depends on the sign of the initial value of λ 0 : for λ 0 0, it flows to zero. Solving the flow equations (3.16)-(3.17) (see the Appendix for details) we find λ(l) ≈ −1/[2(N − 2)l] for l ≫ 1 and arrive at (see also [26]) For λ 0 > 0, the λ-flow equation (3.16) exhibits a runway, i.e., the coupling constant λ leaves the region in which the perturbative calculations are valid. As a result, the continuous (within a mean-field approximation) transition is driven by fluctuations to [32,33]: first order phase transition. (4.7) The above results may be contrasted with those for the pure 2D Ising model: ξ ∼ τ −1 and C sing ∼ ln τ −1 .
(iii) LR correlated disorder (u 0 0, v 0 0). Typical RG flows for the 2D Baxter (N = 2) and the Ncolor Ashkin-Teller models (N = 4) are shown for δ = ±1 in figures 1 and 2, respectively. For δ < 0, the LR disorder is irrelevant while the SR disorder is only marginally irrelevant. Thus, one can neglect the contribution from the LR disorder and the scaling behavior is given by equation  This result was already obtained for the 2D Ising model with LR correlated disorder in [25] which supports the conjecture that the exact correlation length exponent is ν LR = 2/a. The corresponding heat capacity exponent is α LR = 2 − a.
Let us briefly discuss the validity of the extended Harris criterion for the 2D Baxter model with LR correlated disorder. According to the extended Harris criterion, the LR correlated disorder is relevant provided that the correlation length exponent of the pure system satisfies the inequality a < 2/ν pure , i.e., δ > −4λ 0 . Although the extended Harris criterion correctly predicts the relevance of the LR correlated disorder for δ > 0, as we found above the critical behavior is in fact modified for any δ < 0. In the last case, FP LR is unstable, and the asymptotic critical behavior is described by equations (4.8) and (4.9), which corresponds to the 2D Baxter model with SR disorder rather than to the critical behavior of the pure 2D Baxter model. Thus, in the case of the 2D Baxter model, the extended Harris criterion is violated by correlated disorder in the same way as the usual Harris criterion is violated by uncorrelated disorder [27].

Conclusion
We have studied the effect of LR correlated disorder on the 2D Baxter and N -color Ashkin-Teller models. The clean 2D Baxter model exhibits a "weak universal" critical behavior with the critical exponents depending on microscopic parameters, while in the clean N -color Ashkin-Teller model, fluctuations drive the system from the second order to the first order phase transition. Using the mapping to the 2D interacting Dirac fermions in the presence of LR correlated random mass disorder and dimension regularization with double expansion in ε = 2 − d and δ = 2 − a we obtain the RG flow equations to one-loop approximation. Their analysis in d = 2 shows that (i) for a > 2 (δ < 0), the critical behavior is controlled by the Gaussian FP that gives the critical exponents of the clean 2D Ising model (up to logarithmic corrections); (ii) for a < 2 (δ > 0), the only stable FP is the LR FP (λ * = 0, u * = 0, v * = δ/4). It describes the rounding of the weak universality in the Baxter model and the first order phase transition in the N -color Ashkin-Teller model by correlated disorder. This leads to a new emergent critical behavior which is in the same universality class as the 2D Ising model with LR correlated disorder [25]. For instance, we argue that the exact values of the correlation length and heat capacity exponents are ν LR = 2/a and α LR = 2 − a, respectively. Since quantum systems can be mapped onto classical systems in d + 1 dimensions, it would be interesting if these results could be generalized to the first order phase transitions in 1D quantum systems [48,49].
Thus, the asymptotic behavior for large l is given by Substituting equation (7)