Critical behavior of the three-state random-field Potts model in three dimensions

Enormous advances have been made in the past 20 years in our understanding of the random-field Ising model, and there is now consensus on many aspects of its behavior at least in thermal equilibrium. In contrast, little is known about its generalization to the random-field Potts model which has wide-ranging applications. Here we start filling this gap with an investigation of the three-state random-field Potts model in three dimensions. Building on the success of ground-state calculations for the Ising system, we use a recently developed approximate scheme based on graph-cut methods to study the properties of the zero-temperature random fixed point of the system that determines the zero and non-zero temperature transition behavior. We find compelling evidence for a continuous phase transition. Implementing an extensive finite-size scaling (FSS) analysis, we determine the critical exponents and compare them to those of the random-field Ising model.

Introduction: Understanding the effect of quenched disorder on phase transitions is crucial for many experiments, such as magnetic systems with impurities, and technological application areas, such as quantum computers [1]. At the same time, past progress in this direction has profoundly shaped the theory of equilibrium and non-equilibrium statistical mechanics, and theoretical concepts such as replica-symmetry breaking and the cavity method have found applications even in seemingly distant fields such as gene regulation [2], neural networks [3], and the modelling of bird flocks [4].
Most of the recent focus has been on spin glasses, where competing and random interactions lead to a merely short-range ordered state, as well as on random-field systems [5]. The latter are at the heart of such diverse problems as the behavior of the quantum magnet LiHo x Y 1−x F 4 [6] and the random first-order transition scenario in structural glasses [7,8]. For such problems, destruction of order is complete even for weak fields in d = 2 dimensions [9,10], where ferromagnetic domains break up on length scales that vanish with growing strength of the random fields [11]. For continuous O(n) spins, the lower critical dimension is even elevated to d = 4. The random-field Ising model (RFIM), on the other hand, orders at non-zero temperatures already for d ≥ 3, and the transition is of second order, at least for continuous field distributions [12][13][14]. Hence the proposal of dimensional reduction [15] suggested by field theory, where the RFIM in d dimensions would be in the universality class of the d − 2-dimensional ferromagnet, does not apply in low dimensions, but is only recovered for d ≥ d c ≈ 5 [16][17][18].
Much less understood is the case of discrete spins with more than two states, i.e., the random-field Potts model (RFPM) [19,20].
The Potts model has a plethora of applications ranging from finite-temperature QCD [21], over mixed antiferromagnets [22], orientational glasses [23], to soap froths [24]. As disorder is inescapable, the RFPM is of even greater relevance for their study. Additionally, it is of profound theoretical interest since in the pure Potts model [25] the transition order can be tuned by changing the number of states q, such that there is a line q pure c (d) of tricritical points with q pure c (2) = 4 [26,27] and q pure c (3) ≈ 2.35 [28]. Since disorder tends to soften first-order transitions [10,29], one expects a shift of the line q pure c (d) to q RF c (d). Just as for the RFIM, dimensional reduction is not likely to hold in low dimensions [20]. Instead, one expects [9,11] d = 2, and hence absence of long-range order in 2d -a scenario that we recently confirmed numerically [30]. [20,31]. Even once q c (d) is known, however, one needs to ask whether for all strengths ∆ of random fields the first-order transition for d pure will be softened. This would be the case for d = 2 [10], but there is no FM order there. In d > 2, one might expect a line of tricritical points (or even two lines [32]) to appear in the (∆, T ) plane, where the transition changes from first to second order [32,33].
Very little is known about the details of this rich phase diagram (but see Ref. [34]). The purpose of this Letter is to address such issues for the physically most relevant system of the q = 3, d = 3 RFPM, which in experiments has been used to describe trigonal-to-tetragonal structural transitions in SrTiO 3 [35] when stressed along [111], and the mixed antiferromagnet Fe 1−x Co x Cl 2 [36]. Such experimental systems show continuous transitions, but the theoretical situation is unclear as for this case q RF c ≥ q pure c ≈ 2.35. Some early simulational work [37] found first-order transitions for all considered field strengths, but later studies claimed a continuous transition for intermediate fields combined with first-order behavior for small and large fields [32,33]. This scenario agreed with the prediction of Ref. [20] but contradicted Refs. [32,38], who performed a 1/q-expansion of the qstate RFPM in 3D and found a first-order transition for q ≥ 3, irrespective of the field strength. The question of a softening of the discontinuous transition hence has remained undecided.
Due to frustration and the ensuing slow relaxation, Monte Carlo methods are not very efficient in the presence of random fields. For the RFIM, much of the re-   cent progress in understanding is due to the availability of efficient combinatorial optimization methods that allow to find exact ground states (GSs) in polynomial time [12,[39][40][41]. Since the relevant renormalizationgroup fixed point is located at temperature T = 0, such ground-state calculations are also relevant for the finitetemperature transitions. Unfortunately, the same methods do not extend to the RFPM, since the ground-state problem is NP hard for q > 2 [39,42]. As we have recently shown, however, combinatorial graph-cut methods [43,44] can still be used in combination with embedding techniques to efficiently compute high-quality approximate GSs [30,45]. To further improve the accuracy, we run the GS method for n different random initial spin configurations and extrapolate the thermodynamic quantities in the limit n → ∞. The extrapolated results enable us to uncover a clear-cut picture of the phase transition.
Model and Methodology: We consider the q-state RFPM with Hamiltonian [20] where δ x,y is the Kronecker delta function, s i ∈ {0, 1, . . . , q − 1}, and {h α i } are uncorrelated, quenched random-field variables extracted from a standard normal distribution, i.e., ; ∆ denotes the disorder strength. For q = 2, Eq. (1) maps to the RFIM at coupling J/2 and field strength ∆/ √ 2 [30]. Note that different couplings of random fields to the spins are possible as well as different coupling distributions [31,32,34] but such variations are left for future work. I. Estimates of ∆c, ν, and β/ν according to Eq. (4) as well asγ/ν according to Eq. (9) extracted from scaling collapses of the data for different n as well as the extrapolated data for n → ∞ (Lmin = 24). S1 and S2 are the qualities of the collapses according to (4) and (9), respectively (S ≈ 1 for perfect collapses).
n ∆c 1/ν β/νγ/ν S1 S2 1 1.636(2) 0.837(9) 0.0460 (9)  We perform ground-state calculations for the q = 3 RFPM on simple-cubic lattices of edge length L with periodic boundary conditions. The number of disorder samples ranges from N samp = 50 000 for L = 16 to N samp = 5 000 for L = 96. Approximate GSs are obtained using the algorithm described in Ref. [30] that is based on an embedding of Ising spins into the Potts variables in the spirit of the α-expansion method of Ref. [42]. For each disorder sample, we run our algorithm for n different initial spin configurations and pick the run(s) resulting in the lowest energy as ground-state estimate. The success probability of the resulting approach increases exponentially with n, such that the method becomes exact for n → ∞ [45]. For each sample we determine the order parameter [25] Here, ρ denotes the density of spins in the majority orientation. Also, we measure the bond energy per spin e J (L, ∆, n) = − ij δ si,sj /L 3 [46]. After performing the disorder average, [·] av , we then deduce further quantities such as the Binder cumulant associated to m, av . The employed definitions of the specific heat and magnetic susceptibility will be discussed below. All statistical errors were estimated using the jackknife method [47][48][49]. Estimates for ∆ c and the critical exponents are then extracted from the scaling of observables at sequences of pseudocritical points as well as from scaling collapses [50].
Extrapolation and transition order: To assess the quality of approximation, we first studied the behavior of each quantity as a function of n. We generally find a two-stage behavior, with an initial fast decay followed by a much slower large-n convergence, which is well described by the sum of two power laws [45], where b < e is the asymptotic, slow exponent, e describes the initial fast decay, and O * denotes the limiting value for n → ∞. As was shown elsewhere, this form is quite generic and it holds, in particular, for a certain subset of samples for which exact GSs are known [45]. For such exact samples of size 16 3 we employed our algorithm for up to n = 10 4 runs and found that the residuals with respect to the exact results, i.e., O(n) − O ex for any quantity scale as an −b (1 + cn −e ), with b 0.02 and e 0.5. This behavior is seen to extend to the case where the exact results are not used or known [45]. The value of b is found to be very stable in this regard, such that we fix it for the subsequent fits of our main study reported here, for which n ≤ 100 [51]. We then perform joint fits of the functional form (3) to [m] av , U 4 , and [e J ] av for a common value of the exponent e, yielding extrapolated estimates m * , U * , and e * for any fixed (L, ∆) (see Section S1 of Ref. [52]). We prefix our analysis by a study of the energetic Binder cumulant, whose scaling clearly shows the behavior expected from a continuous transition, as already reported for T > 0 in Ref. [31]; details are provided in Ref.
[52], Section S2. In the following, we perform FSS of all quantities for finite n as well as for n → ∞, in order to determine the transition point and obtain the critical exponents. Magnetization and Binder cumulant: In panel (a) of Fig. 1 we show the extrapolated magnetization m * as a function of the disorder strength ∆ for various lattice sizes L. The expected FSS form [53] implies that a plot of m * (L, ∆)L β/ν against x = (∆ − ∆ c )L 1/ν should yield a collapse of data sets for small |x| for appropriate values of ∆ c , ν, and β/ν. This is shown in Fig. 2(a). The collapse with the least meansquared deviation S from the master curve [50,54] is obtained for ∆ c = 1.606 ± 0.003, 1/ν = 0.723 ± 0.004, and β/ν = 0.0306 ± 0.0023. We also performed FSS of the magnetization for finite n = 1, 5, 10, 50, and 100. As is clear from Table I, there is a weak dependence of the estimates on n with a smooth convergence to the exact limit n → ∞.
Turning to the Binder cumulant, from Fig. 1(b) we  (5), see Table II. For increased clarity, the data for different n are slightly shifted relative to each other. (b) Scaling of the maxima Cs,max(L) as a function of L for different n. The curves correspond to fits of the form (6) to the data (cf. Table II).
see a crossing of U * (L, ∆) = U 4 (L, ∆, n → ∞) in the range ∆ c = [1.59, 1.61], predicting the location of the critical point, where U 4 becomes system-size independent [55,56]. The FSS form of U * is [53] U * (L, ∆) = U (∆ − ∆ c )L 1/ν . As is seen from the rescaled data in Fig. 2(b), a rather clean scaling collapse is achieved for ∆ c = 1.604(2) and 1/ν = 0.720 (6), yielding an alternative, and consistent, set of estimates for ∆ c and 1/ν. Specific heat: Due to the restriction to T = 0 and the uniqueness of the ground state, it is not possible to define the specific heat from a temperature derivative or fluctuation-dissipation relation. Instead, a specific-heat like quantity is given by the derivative C(∆) = ∂[e J (∆)] av /∂∆ [46,57]. Numerically, we compute this using the standard three-point formula at the midpoint [58]. In Fig. 1(d) we show the extrapolated C * (L, ∆). As L is increased, the peaks shift but only weakly vary in height, indicating a small specific-heat exponent α. To determine the latter, we considered additional ∆ values and used parabolic fits C(L, ∆) = a 0 (∆ − ∆ max,C ) 2 + C max to obtain the peak locations ∆ max,C (L) and peak heights C max (L). In a finite system, the singular part of C(L, ∆) scales as [53] C s (L, ∆) = L α/ν C (∆ − ∆ c )L 1/ν . If the maximum of C occurs for argument a 1 , then the positions in ∆ shift as and the maximum value of the singular part of the specific heat C s,max (L) ∼ L α/ν . Figure 3(a) shows our data for ∆ max,C (L) together with fits of the form (5), indicating clear consistency. The estimates of ∆ c and the exponent 1/ν are collected in Table II, where we also indicate the quality Q 1 of these fits [58]. In Fig. 3(b) we present the corresponding peak heights C max (L, n). This plot clearly suggests that α is either positive, but very small, or negative. As a consequence, scaling corrections are relevant and we hence considered the functional form where ω corresponds to the Wegner exponent and C 0 represents a non-singular background term. Since α ≈ 0 would result in a second additive constant in (6), we cannot reliably include all five parameters in the fit. We hence fix C 0 = 0, and the resulting 4-parameter fit yields excellent qualities Q 2 and the estimates of α/ν and ω collected in Table II. We thus conclude that α is very slightly negative or perhaps zero. Susceptibility: We now turn to the connected and disconnected magnetic susceptibilities. Since we operate at T = 0, we cannot use the usual fluctuation-dissipation relation. As outlined in Section S3 of Ref.
[52], we generalize arguments for the RFIM [59] to express the zero-field RFPM susceptibility, for Gaussian random fields as χ µ = [ m µ i h µ i ] av /∆ 2 , where m µ = i δ si,µ /L 3 (see also [13]). We apply a constant external field H to the spin state 1 (i.e., µ = 1) to break the symmetry so that χ displays a peak. As is shown in S3, a minimal field strength ∝ L −3/2 is necessary to break the symmetry, and we choose H(L = 16) = 8 × 10 −2 . Due to the large sample-to-sample fluctuations in χ we do not observe a systematic variation of χ(L, ∆, n) with n, such that instead of an extrapolation we focus on n = 100. Figure 4(a) shows the behavior of χ(L, ∆, n = 100) as a function of ∆ and for various L. To analyze its divergence, we fit a parabola near the peak and obtain the peak positions ∆ max,χ (L) as well as the maxima of the susceptibility χ max (L). FSS predicts that χ(L, ∆) = L γ/ν χ (∆ − ∆ c )L 1/ν , χ max (L) ∼ L γ/ν .   (5) and (6) to the data for the specific heat. Q1 and Q2 refer to the quality-of-fit [58].  Table I. 0.007). The corresponding data are shown in Fig. S5(a) of Ref. [52]. While the value of ∆ c agrees with that obtained from the specific heat (see Table II for n = 100), the estimate for 1/ν is noticeably larger than the values from the magnetization and specific heat, indicative of unresolved scaling corrections. To avoid this problem, we performed additional simulations without external field H at fixed ∆ ∈ [1.602, 1.606], corresponding to the estimate ∆ c = 1.604(2) found above from U * (L, ∆). Here, we find clean scaling, cf. Fig. S5(b), and the fit parameters collected in Table S1, indicating little scaling corrections in χ, and the dependence of γ/ν on ∆ is only weak; we quote γ/ν = 1.51(6) at our best estimate ∆ c = 1.604 for L min = 16.
The additional length scale introduced by the disorder fluctuations results in a different scaling behavior of the static, disconnected susceptibility χ dis = L d [m 2 ] av as compared to the thermal, connected susceptibility χ [60], We use this relation to perform a scaling collapse for the extrapolated χ * dis . As is shown in Fig. 4(b), this leads to an excellent scaling result; the exponentsγ/ν, including those for finite n, are provided in Table I, with γ/ν = 2.9402 (30) for n = ∞. The so-called two-exponent scaling scenario predictsγ = 2γ [61,62]. From our data we find the marginal result (2γ −γ)/ν = 0.08 (6). We can also investigate the validity of the Rushbrooke equality α + 2β + γ = 2 and the modified hyperscaling relation 2 − α = ν(d − θ) [60,63,64], where θ = 2 −η + η = (γ − γ)/ν. It can be inspected from Tables I and II that both relations are well satisfied (within error bars) by the results for infinite n.
Conclusions: We have presented a high-resolution numerical study of the 3D 3-state RFPM via a computationally efficient ground-state method that uses extrapolation in the number n of initial conditions to provide quasi-exact results. With the help of FSS, all critical exponents, including the most elusive specific-heat exponent α, are determined for finite as well infinite n, providing clear evidence for a continuous phase transition in the vicinity of ∆ c = 1.604 (2). Given that we are working at T = 0, it is hence clear that 3 = q < q RF c (3). We are thus able to provide the first comprehensive calculation of critical exponents. As the summary in Table III shows, these are close to but likely distinct from the exponents of the 3D RFIM [12,40]. It will be most interesting to see if this deviation grows stronger as q is increased to 4 and possibly beyond, and if and when the transition turns first order.

ACKNOWLEDGMENTS
We thank N. Fytas for a careful reading of the manuscript. The authors are thankful to Kurt Binder for manifold discussions on topics related to the manuscript, and they dedicate this paper to his memory. MK and MW acknowledge the support of the Royal Society -SERB Newton International fellowship (NIF\R1\180386). All numerical simulations were done on the parallel compute cluster Zeus of Coventry University. The publication of this article was funded by Chemnitz University of Technology and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) -491193532. 1 Institut für Physik, Technische Universität Chemnitz, 09107 Chemnitz, Germany. 2 Department of Physics, Indian Institute of Technology, Hauz Khas, New Delhi -110016, India.

S1. EXTRAPOLATION OF PHYSICAL QUANTITIES
We study the behavior of physical quantities as a function of the number of initial conditions n. For each disorder realization, we run the simulations for at most 100 different initial spin configurations, i.e., n max = 100, and for any n ≤ n max , we determine all observables from the state of minimum energy out of the n approximate ground states.

S2. ORDER OF THE TRANSITION
As was first pointed out by Binder and co-workers [1,2], a useful quantity to study the order of a temperature-driven transition is the energetic cumulant V E defined as Panel (a) of Fig. S4 shows a plot of V E against ∆ for n = 100 with various system sizes in the range 16 ≤ L ≤ 96. It can be seen that this quantity has a minimum at a certain ∆ = ∆ min (L) and the depth of the minima reduces with  (8) . ¡s increasing L. In panel (b), we analyze the minimum of V E , subtracting the trivial limit 2/3, i.e, We see that the behavior of V E,min against L follows a power-law decay with the exponent value 1.958 ± 0.008, i.e., V E,min (L) ∼ L −1.958 (8) . This suggests a second-order phase transition for the three-state RFPM studied here as V E,min → 0 in the thermodynamic limit L → ∞ and V E (∆ = ∆ min ) → 2/3. For a first-order transition, V E (L, ∆) would be expected to saturate at a value less than 2/3 [2].

S3. FLUCTUATION FORMULA FOR THE SUSCEPTIBILITY
We follow the arguments put forward by Schwartz and Soffer [3] for the purpose of establishing a bound on the susceptibility exponent of the RFIM to derive an expression for the zero-field susceptibility of the RFPM. We start from the Hamiltonian of the model in an additional external magnetic field H α , We consider the magnetization, which in general is a vector of q components. We are interested in the disorder average of the derivative ∂ M µ /∂H µ in the limit H µ → 0, where M µ denotes the thermal average of the magnetization. If we writē the disorder average becomes   (14) 0.58 TABLE S1. The exponent γ/ν from fits of the form χ(L, ∆c) ∼ L γ/ν for different estimates ∆c of the critical field strength as a function of the cut-off Lmin (data for n = 100).
of the random fields deep in the ferromagnetic phase is ∆ √ N . Hence a field of strength H ∝ ∆/ √ N is sufficient to break the symmetry.
Here, we stress that the susceptibility χ can also be studied by a standard approach of determining the GS explicitly in the presence of external fields H and fitting [m] av against H to a function y = a 0 + a 1 x + a 2 x 2 . The parameter a 1 then determines χ. While this works relatively well for the RFIM if exact GS algorithms are used [4], for the RFPM the approximate nature of the algorithm can lead to negative values of χ and other artifacts.
Applying the explicit symmetry breaking described above and the parabolic fit procedure described in the main paper, we arrive at estimates for the maxima of the susceptibility which are expected to scale as while the locations of the maxima behave as ∆ max,χ (L) = ∆ c + a 1 L −1/ν . (S10) Our simulation data presented in Fig. S5(a) are consistent with these assumptions, but the somewhat large value for 1/ν hints at the presence of unresolved scaling corrections. For the alternatively performed analysis at fixed random-field strengths ∆ ∈ [1.602, 1.606] results are shown in Fig. S5(b). As the fit parameters in Table S1 imply, the dependence on the fixed ∆ used is quite moderate in this range, and we hence use the central value for ∆ = 1.604, which is γ/ν = 1.51 (6).