Disordered Ising model with correlated frustration

TL;DR

This paper investigates a disordered Ising model with correlated frustration, analyzing its critical behavior and universality class through numerical estimates of critical exponents at a specific disorder parameter.

Contribution

It introduces a gauge-invariant disorder distribution with long-range correlations and provides precise estimates of critical exponents for this model.

Findings

Critical exponents estimated: ν = 0.655(15), η_q = 1.05(5).

Frustration exhibits long-range correlations at the studied parameter.

Model belongs to a specific universality class characterized by these exponents.

Abstract

We consider the $\pm J$ Ising model on a cubic lattice with a gauge-invariant disorder distribution. Disorder depends on a parameter $\beta_G$ that plays the role of a chemical potential for the amount of frustration. We study the model at a specific value of the disorder parameter $\beta_G$, where frustration shows long-range correlations. We characterize the universality class, obtaining accurate estimates of the critical exponents: $\nu = 0.655(15)$ and $\eta_q = 1.05(5)$, where $\eta_q$ is the overlap susceptibility exponent.

Tables

Table 1. Table 1: Results (data at β G = β G , c subscript 𝛽 𝐺 subscript 𝛽 𝐺 𝑐 \beta_{G}=\beta_{G,c} ) of the fits to Eq. ( 20 ): the order n 𝑛 n of the expansion is given in the second column. In each fit we only include the data satisfying L ≥ L min 𝐿 subscript 𝐿 min L\geq L_{\rm min} . In the last column we report the sum of the residuals divided by the number of degrees of freedom (DOF) of the fit.

	$n$	$L_{\min }$	${\beta }_c$	$R^{\ast }$	$\nu$	${\chi }^2$/DOF
$R_{\xi }$	5	16	0.225753(12)	0.1726(6)	0.681(9)	1.4
	6	16	0.225755(12)	0.1729(6)	0.676(10)	1.3
	5	20	0.225767(14)	0.1735(8)	0.678(14)	1.5
	6	20	0.225766(14)	0.1736(8)	0.664(15)	1.4
	5	24	0.225756(23)	0.1728(14)	0.678(27)	1.4
	6	24	0.225756(23)	0.1729(14)	0.654(28)	1.3
$U_{q,4}$	7	16	0.225364(12)	4.686(7)	0.653(10)	5.8
	8	16	0.225369(13)	4.687(7)	0.653(10)	5.9
	7	20	0.225446(15)	4.748(10)	0.631(16)	4.2
	8	20	0.225448(15)	4.749(10)	0.634(16)	4.3
	7	24	0.225536(23)	4.812(20)	0.631(33)	1.9
	8	24	0.225537(23)	4.817(20)	0.614(33)	1.9

Table 2. Table 2: Results [data at β G = β G ​ ( L ) subscript 𝛽 𝐺 subscript 𝛽 𝐺 𝐿 \beta_{G}=\beta_{G}(L) ] of the fits to Eq. ( 20 ): the order n 𝑛 n of the expansion is given in the second column. In each fit we only include the data satisfying L ≥ L min 𝐿 subscript 𝐿 min L\geq L_{\rm min} . In the last column we report the sum of the residuals divided by the number of degrees of freedom (DOF) of the fit.

	$n$	$L_{\min }$	${\beta }_c$	$R^{\ast }$	$\nu$	${\chi }^2$/DOF
$R_{\xi }$	3	16	0.225980(12)	0.1157(3)	0.651(14)	1.6
	4	16	0.225977(12)	0.1158(3)	0.649(14)	1.6
	3	20	0.225951(16)	0.1149(4)	0.630(20)	1.8
	4	20	0.225953(17)	0.1150(4)	0.639(22)	1.8
$U_{q,4}$	3	16	0.225629(13)	4.209(10)	0.681(20)	2.0
	4	16	0.225629(13)	4.210(10)	0.670(21)	1.9
	3	20	0.225711(18)	4.288(17)	0.643(30)	1.3
	4	20	0.225708(19)	4.289(18)	0.655(33)	1.2

Table 3. Table 3: Number of samples N s subscript 𝑁 𝑠 N_{s} and number of iterations per sample N m subscript 𝑁 𝑚 N_{m} for the runs at β G = β G , c subscript 𝛽 𝐺 subscript 𝛽 𝐺 𝑐 \beta_{G}=\beta_{G,c} .

$L$	$\beta$	$N_s/{10}^5$	$N_m/{10}^3$	$L$	$\beta$	$N_s/{10}^5$	$N_m/{10}^3$
16	0.2200	3.2	10	24	0.2257	19.2	5.5
	0.2220	3.2	10		0.2258	19.2	5.5
	0.2240	25.6	10		0.2259	19.2	5.5
	0.2245	12.8	10		0.2260	11.5	5.5
	0.2250	12.8	10		0.2262	12.8	5.5
	0.2255	12.8	10		0.2264	12.8	5.5
	0.2260	25.6	10		0.2265	10.2	5.5
	0.2265	12.8	10		0.2267	12.8	5.5
	0.2270	2.56	10		0.2270	3.2	10
	0.2280	3.2	10		0.2280	3.2	10
	0.2300	3.2	10	30	0.2230	3.2	10
20	0.2200	0.64	10		0.2240	3.8	30
	0.2220	0.64	10		0.2245	3.8	30
	0.2240	11.8	10		0.2250	3.8	30
	0.2245	12.8	10		0.2255	19.2	5.5
	0.2250	12.8	10		0.2256	19.2	5.5
	0.2255	19.2	10		0.2257	19.2	5.5
	0.2256	19.2	10		0.2258	19.2	5.5
	0.2257	19.2	5		0.2259	19.2	5.5
	0.2258	19.2	5		0.2260	6.4	5.5
	0.2259	19.2	5		0.2262	6.4	5.5
	0.2260	19.2	5		0.2264	6.4	5.5
	0.2262	12.8	5		0.2265	6.4	5.5
	0.2264	19.2	5		0.2267	6.4	5.5
	0.2265	19.2	5		0.2270	6.4	5.5
	0.2267	19.2	5	34	0.2255	12.8	5.5
	0.2270	16	5		0.2256	12.8	5.5
	0.2280	3.2	10		0.2257	12.8	5.5
	0.2300	3.2	10		0.2258	12.8	5.5
24	0.2200	0.64	10		0.2259	12.8	5.5
	0.2220	0.64	10		0.2260	6.4	5.5
	0.2240	5.12	10		0.2262	6.4	5.5
	0.2245	5.12	10		0.2264	6.4	5.5
	0.2250	5.12	10		0.2265	6.4	5.5
	0.2255	19.2	5.5		0.2267	6.4	5.5
	0.2256	19.2	5.5		0.2270	6.4	5.5

Table 4. Table 4: Number of samples N s subscript 𝑁 𝑠 N_{s} and number of iterations per sample N m subscript 𝑁 𝑚 N_{m} for the runs at β G = β G ​ ( L ) subscript 𝛽 𝐺 subscript 𝛽 𝐺 𝐿 \beta_{G}=\beta_{G}(L) .

$L$	$\beta$	$N_s/{10}^5$	$N_m/{10}^3$	$L$	$\beta$	$N_s/{10}^5$	$N_m/{10}^3$
16	0.2240	6.4	96	20	0.2262	4.48	140
	0.2245	6.4	96		0.2263	4.48	140
	0.2250	6.4	96		0.2265	6.4	140
	0.2255	6.4	96	24	0.2240	2.24	180
	0.2256	6.4	96		0.2245	2.62	180
	0.2257	6.4	96		0.2250	2.82	180
	0.2258	6.4	96		0.2255	6.4	180
	0.2259	6.4	96		0.2256	6.4	180
	0.2260	6.4	96		0.2257	6.4	180
	0.2261	5.12	96		0.2258	6.4	180
	0.2262	5.12	96		0.2259	6.4	180
	0.2263	5.12	96		0.2260	2.62	180
	0.2265	6.4	96		0.2261	2.56	180
20	0.2240	6.4	140		0.2262	2.56	180
	0.2245	6.4	140		0.2263	2.56	180
	0.2250	6.4	140		0.2265	1.92	180
	0.2255	6.4	140	30	0.2255	3.2	30
	0.2256	6.4	140		0.2256	3.2	30
	0.2257	6.4	140		0.2257	4.28	75
	0.2258	6.4	140		0.2258	4.28	75
	0.2259	6.4	140		0.2259	3.2	30
	0.2260	6.4	140		0.2260	3.2	30
	0.2261	4.48	140		0.2261	3.2	30

Equations

$$\mathcal{H}=-\sum _{⟨xy⟩}{J}_{xy}{\sigma }_x{\sigma }_y,$$

$$P(J_{xy})=p\delta (J_{xy}-1)+(1-p)\delta (J_{xy}+1).$$

$${\Pi }_p=J_{{\ell }_1}{J}_{{\ell }_2}{J}_{{\ell }_3}{J}_{{\ell }_4},$$

$$P(\{J_{xy}\})=\frac{1}{Z_G}{e}^{-{\beta }_G{H}_G}{H}_G=-\sum _P{\Pi }_P,$$

$$P(\{J_{xy}\})=\frac{1}{Z_G}\prod _P\left[e^{{\beta }_G}\delta ({\Pi }_P-1)+e^{-{\beta }_G}\delta ({\Pi }_P+1)\right],$$

$$J_{xy}\to J_{xy}^{\prime }={ϵ}_x{J}_{xy}{ϵ}_y,$$

$${\sigma }_x\to {\sigma }_x^{\prime }={ϵ}_x{\sigma }_x.$$

$${\beta }_{G,c}=\frac{1}{2}\ln \coth {\beta }_{I,c}.$$

$$q_x={\sigma }_x^{(1)}{\sigma }_x^{(2)},$$

$$U_{q,4}=\frac{\bigl{[}\langle Q^{4}\rangle\bigr{]}}{\bigl{[}\langle Q^{2}\rangle\bigr{]}^{2}}\qquad Q=\sum_{x}q_{x},$$

$$G_q(x)=\frac{1}{V}\sum _y[⟨q_y{q}_{x+y}⟩]=\frac{1}{V}\sum _y[⟨{\sigma }_y{\sigma }_{x+y}{⟩}^2]{\tilde{G}}_q(k)=\sum _x{e}^{ik\cdot x}{G}_q(x),$$

$${\xi }_q^2=\frac{{\chi }_q-{\chi }_1}{4{\sin }^2(\pi /L){\chi }_1},$$

$$\chi_{1}=\frac{1}{3}\Bigl{[}\tilde{G}_{q}(2\pi/L,0,0)+\tilde{G}_{q}(0,2\pi/L,0)+\tilde{G}_{q}(0,0,2\pi/L)\Bigr{]}.$$

$$R=f_R(t_1{L}^{1/{\nu }_I},t_2{L}^{1/\nu }),$$

$$R=f_R(0,X)X\equiv (\beta -{\beta }_c)L^{1/\nu }.$$

$$C_v=\frac{{\beta }_G^2}{V}\sum _{P,Q}\left[⟨{\Pi }_P{\Pi }_Q⟩-⟨{\Pi }_P⟩⟨{\Pi }_Q⟩\right],$$

$${\beta }_G(L)={\beta }_{G,c}+X_{\max }{L}^{-1/{\nu }_I},$$

$$t_2=\beta -{\beta }_c+a{X}_{\max }{L}^{-1/{\nu }_I},$$

$${\chi }_q(\beta ,L)=L^{2-{\eta }_q}{u}_{\chi }(t_1,t_2)f_{\chi }(t_1{L}^{1/{\nu }_I},t_2{L}^{1/\nu }),$$

$$R(\beta ,L)=R^{\ast }+\sum _{k=1}^n{a}_k{X}^{k}X\equiv (\beta -{\beta }_c)L^{1/\nu },$$

$$R(\beta ,L)=R^{\ast }+\sum _{k=1}^n{a}_k{X}^k+L^{-\omega }\sum _{k=0}^m{b}_k{X}^{k}X\equiv (\beta -{\beta }_c)L^{1/\nu }.$$

$${\beta }_c=0.225754(11)\nu =0.666(8)\omega =1.18(25).$$

$${\beta }_c=0.225823(37)\nu =0.637(14)\omega =0.83(27).$$

$$R(\beta ,L)=R^{\ast }+\sum _{k=1}^n{a}_k{X}^k+L^{{\eta }_q-2}\sum _{k=0}^m{b}_k(\beta -{\beta }_c{)}^{k}X\equiv (\beta -{\beta }_c)L^{1/\nu },$$

$${\beta }_c=0.225742(19)\nu =0.657(11),$$

$${\beta }_c=0.22577(5)\nu =0.655(15)\omega =1.0(3).$$

$$\ln {\chi }_q=(2-{\eta }_q)\ln L+\sum _{k=0}^n{a}_k{X}^k+b(\beta -{\beta }_c).$$

$${\eta }_q=1.05(5).$$

$$\ln {\chi }_q=(2-{\eta }_q)\ln L+\sum _{k=0}^n{a}_k{X}^k+b(\beta -{\beta }_c)+L^{-\omega }\sum _{k=0}^m{c}_k{X}^k,$$

$$\ln {\chi }_q=(2-{\eta }_q)\ln L+\sum _{k=0}^n{a}_k{X}^k+b(\beta -{\beta }_c)+L^{{\eta }_q-2}\sum _{k=0}^m{c}_k(\beta -{\beta }_c{)}^k.$$

Full text

1 Introduction

The $\pm J$ Ising model is a paradigmatic model to study the effects of quenched random disorder and frustration. It is defined by the lattice Hamiltonian

[TABLE]

where $\sigma_{x}=\pm 1$ and the sum is over the nearest-neighbor sites of a lattice. In the following we will focus on the three-dimensional system defined on a simple cubic lattice. Usually, the exchange interactions $J_{xy}$ are uncorrelated quenched random variables, taking values $\pm 1$ with probability distribution

[TABLE]

For $p=1$ we recover the standard ferromagnetic Ising model, while for $p=1/2$ we obtain the bimodal Edwards-Anderson (EA) spin-glass model. The $\pm J$ Ising model is a simplified model [1] for some disordered magnetic materials, such as [2, 3, 4] Fe1-xMnxTiO3 and Eu1-xBaxMnO3 that have a low-temperature glassy phase.

The behavior of the three-dimensional $\pm J$ Ising model is well known. In this model, disorder is a relevant perturbation of the pure Ising fixed point: the introduction of any amount of disorder in the couplings changes the universality class [5, 6]. The $T$-$p$ phase diagram is characterized by a high-temperature paramagnetic phase separated by two different low-temperature phases by a transition line. If we only consider the case $p>1/2$ (the phase diagram is symmetric under $p\to 1-p$), we have a ferromagnetic transition for $p>p^{*}$ and a glassy transition in the opposite case. In the first case, the transition belongs to the randomly-dilute Ising (RDI) universality class [7] characterized by the magnetic critical exponents [8, 9, 10] $\nu_{f}=0.683(2)$ (quite different from the Ising value $\nu_{I}\approx 0.630$ [6]) and $\eta_{f}=0.036(1)$. On the other hand, for any $1-p^{*}<p<p^{*}$, the transition belongs to the same universality class as that of the bimodal Ising spin glass model at $p=1/2$, with [11, 12, 13, 14] $\nu_{SG}=2.56(4)$ and $\eta_{q,SG}=-0.390(4)$. The point $p=p^{*}=0.76820(4)$ is a multicritical point [15, 16, 17, 18]. The renormalization-group (RG) dimensions $y_{1}$ and $y_{2}$ of the relevant operators that control the RG flow close to it are [18] $y_{1}=1.02(5)$, $y_{2}=0.61(2)$. Correspondingly, the temperature and crossover exponents are $\nu=1/y_{2}=1.64(5)$ and $\phi=y_{2}/y_{1}=1.67(10)$, respectively.

In this paper we consider again the $\pm J$ model with Hamiltonian (1), but we use a different probability distribution. We do not directly fix the distribution of the couplings, but rather we specify the amount of frustration present in the system. Moreover, we use a gauge-invariant probability distribution: two sets of bonds that correspond to the same frustration distribution are given the same probability. Frustration is quantified by using the product of the couplings $J_{xy}$ corresponding to the links along an elementary lattice square (plaquette) and the corresponding probability distribution is nothing but the Boltzmann-Gibbs probability for the lattice ${\mathbb{Z}}_{2}$ gauge theory [19, 20] at inverse temperature $\beta_{G}$. In this context the parameter $\beta_{G}$ plays the role of a chemical potential for the frustrated lattice plaquettes. Since the probability distribution depends on the values of the couplings along the different plaquettes, the bond variables on different links are necessarily correlated. If these correlations are very short-ranged (this is the case for $|\beta_{G}|$ small), we expect the behavior to be analogous to that of the usual EA model.111A different correlated distribution was considered in Ref. [21]: the critical behavior at the transition was the same as that in the EA model with coupling distribution (2). On the other hand, as $\beta_{G}$ increases, correlations become strong and we expect a different phase diagram. Here, we focus on a specific value $\beta_{G,c}$ of $\beta_{G}$ that corresponds to a phase transition point of the exchange coupling distribution. At $\beta_{G,c}$ frustration shows long-range power-law correlations. In the temperature-$\beta_{G}$ plane, close to this point, the system shows a multicritical behavior. We study it in detail, identifying the universal features of the new universality class. We also report some preliminary results for the behavior for $\beta_{G}>\beta_{G,c}$. In this range of values of $\beta_{G}$, the model apparently behaves as the pure ferromagnetic Ising model. These results, if confirmed by simulations on larger lattices, imply that the ferromagnetic transition in pure systems is stable with respect to the introduction of a small amount of correlated disorder. This is at variance with what happens for uncorrelated disorder, i.e., in the case one uses distribution (2). In the latter case, the introduction of any amount of antiferromagnetic bonds changes the universality class of the system [5, 6].

The paper is organized as follows. In section 2 we present the model, while in section 3 we define the quantities we measure in the simulation. In section 4 we present our Monte Carlo results at the transition of the gauge model, while some preliminary results in the exchange-coupling low-temperature phase are presented in section 5. We draw our conclusions in section 6. Some computational details are reported in the Appendix.

2 The model

We consider the $\pm J$ Ising model with quenched correlated disorder. The model is defined on a cubic lattice of linear size $L$ in all directions, with Hamiltonian (1). In the usual EA model, the bond couplings on different links are uncorrelated and distributed with probability distribution (2). The presence of negative couplings makes the EA model frustrated. A simple measure of frustration is provided by the product $\Pi_{P}$ of the couplings $J_{xy}$ along an elementary lattice square (plaquette) $P$:

[TABLE]

where $\ell_{1}$, $\ell_{2}$, $\ell_{3}$, and $\ell_{4}$ are the four links belonging to the plaquette $P$. If $\Pi_{P}<0$, the plaquette is frustrated: it is not possible to minimize the local energy on each link of the plaquette, i.e., it is not possible that $J_{xy}\sigma_{x}\sigma_{y}>0$ on all four links $\ell_{1},\ldots,\ell_{4}$.

In this work we consider a different probability distribution for the bond couplings that allows us to tune directly the amount of frustration in the system. We consider

[TABLE]

where the sum is over all lattice plaquettes and $Z_{G}$ is a normalizing factor. The parameter $\beta_{G}$ plays the role of a chemical potential for the frustrated plaquettes. For $\beta_{G}\to\infty$, we obtain an unfrustrated model, while maximal frustration is obtained in the opposite limit $\beta_{G}\to-\infty$. For $\beta_{G}=0$ the model is equivalent to the EA one with $p=1/2$. Note that the probability distribution can also be written as

[TABLE]

that explicitly shows the role of $\beta_{G}$ in controlling frustration.

The model defined by Eq. (4) is the well-known lattice ${\mathbb{Z}}_{2}$ gauge theory [19, 20], which is invariant under the ${\mathbb{Z}}_{2}$ gauge transformations

[TABLE]

for any arbitrary site-dependent $\epsilon_{x}=\pm 1$. Note that also Hamiltonian (1) is gauge invariant, provided the Ising spins are transformed according to

[TABLE]

Using duality [20, 22], for positive values of the coupling $\beta_{G}$, it is possible to map the gauge theory (4) onto the pure Ising model. Therefore, for $\beta_{G}>0$, the gauge model has two different phases separated by a continuous phase transition located at $\beta_{G,c}$, which can be related to the transition inverse temperature $\beta_{I,c}$ in the pure Ising model by

[TABLE]

The pure Ising critical temperature is known with high-precision [23, 24, 25]. Using [24] $\beta_{I,c}=0.221654626(5)$, we obtain $\beta_{G,c}=0.761413292(11)$. It is interesting to note that the distribution (4) is also well defined for negative $\beta_{G}$. In this case the behavior of the gauge model is not related to that of the pure Ising model.

3 Observables

Since the theory is gauge invariant, we should only consider gauge invariant quantities. As usual in studies of glassy systems, we consider observables defined in terms of the overlap variable $q_{x}$

[TABLE]

where the indices in parentheses refer to two replicas of the system with the same bond couplings. Note that the overlap variable $q_{x}$ is gauge invariant and therefore thermal averages of any function of $q_{x}$ are identical for coupling distributions that differ by a gauge transformation. We define the Binder cumulant

[TABLE]

where the angular and square brackets refer to the thermal average and to the quenched average over the bond distribution, respectively. The Binder parameter $U_{q,4}$ plays an important role in our analysis as it is renormalization-group (RG) invariant. We also consider the overlap correlation function

[TABLE]

($V=L^{3}$ is the volume), the corresponding susceptibility $\chi_{q}=\tilde{G}_{q}(0)=[\langle Q^{2}\rangle]/V$ and correlation length

[TABLE]

where

[TABLE]

Note that $R_{\xi}=\xi_{q}/L$ is RG invariant.

4 Multicritical behavior close to the critical point $\beta_{G}$

The behavior of the disordered system depends on both $\beta_{G}$ and $\beta$. Since the gauge model has a transition for $\beta_{G}=\beta_{G,c}$, we expect at least three different universality classes for the behavior of the spin system, one for $\beta_{G}>\beta_{G,c}$, one for $\beta_{G}<\beta_{G,c}$ and one for $\beta_{G}=\beta_{G,c}$. In particular, if the spin system has a transition for $\beta=\beta_{c}$ when $\beta_{G}=\beta_{G,c}$, the point $(\beta_{G,c},\beta_{c})$ is a multicritical point in the parameter space $(\beta_{G},\beta)$. The relevant scaling variables are two linear scaling fields $t_{1}$ and $t_{2}$, that are linear combinations of $\beta_{G}$ and $\beta$. However, the transition in the gauge model is unrelated to the value of $\beta$ and hence one scaling field should simply be $t_{1}=\beta_{G}-\beta_{G,c}$. The second one should be instead $t_{2}=\beta-\beta_{c}+a(\beta_{G}-\beta_{G,c})$, where $a$ is an appropriate constant. In a neighborhood of the multicritical point, any RG invariant quantity $R$ behaves as

[TABLE]

where $\nu_{I}$ is the pure Ising exponent ($\nu_{I}=0.629971(4)$ in the conformal-boostrap approach [25] and $\nu_{I}=0.63002(10)$, $0.62991(9)$ numerically [23, 24]) and $\nu$ is a new critical exponent characterizing the universality class.

To determine the critical behavior of the system, we will study two different cases. First, we perform simulations setting $\beta_{G}=\beta_{G,c}$. In this case, Eq. (14) becomes

[TABLE]

Second, we consider the behavior at the finite-size pseudocritical transition of the gauge model, defined by considering the maximum of the specific heat $C_{v}$,

[TABLE]

where the sum runs over all pairs of plaquettes $P$ and $Q$. In figure 1 we report the specific heat as a function of $\beta_{G}$. For each $L$ in the range $10\leq L\leq 30$, we determine the position $\beta_{G}(L)$ of the maximum, by performing a quadratic fit close to the maximum. The results are then fitted to

[TABLE]

to obtain the constant $X_{\rm max}$. Using [23] $\nu_{I}=0.63002$, we obtain $X_{\rm max}=-0.8080(53)$. The quality of the fit is good: The sum of the residuals ($\chi^{2}$) normalized to the number of degrees of freedom (DOF) is approximately 0.4. Therefore, we perform a second set of simulations varying $\beta_{G}$ with $L$ as in Eq. (17), using $X_{\rm max}=-0.8080$. Note that, for this set of runs, we have

[TABLE]

and, therefore, we can simply set $t_{2}=\beta-\beta_{c}$, neglecting scaling corrections. The comparison of the results of the simulations at fixed $\beta_{G}=\beta_{G,c}$ and at $\beta_{G}(L)$ will provide an important check for the final results. Indeed, we expect quite different $L$ corrections in the two cases, since, for finite $L$, we are sampling a somewhat different set of configurations.

Beside RG invariant quantities, we will also consider the susceptibility $\chi_{q}$, which scales close to the multicritical point as

[TABLE]

where $\eta_{q}$ is a universal exponent and $u_{\chi}$ is a nonlinear scaling field.

4.1 Behavior of the overlap observables

at the gauge critical point

We have performed simulations at $\beta_{G,c}$ for several values of $L$ in the range $16\leq L\leq 34$. Computational details are reported in the Appendix. Results for $R_{\xi}=\xi_{q}/L$ and the Binder parameter $U_{q,4}$ are reported in figure 2. Data show an intersection for $\beta\simeq 0.226$, that gives us a first estimate of the critical temperature. To obtain more precise estimates, we use Eq. (15). Assuming that $X$ is small, we can expand the scaling function in powers of $X$. Then, we fit any RG invariant quantity $R$ to

[TABLE]

where $R^{*}=f_{R}(0,0)$ is a universal quantity. Here $n$ is the order of the polynomial in $X$, which is tuned by looking at the quality of the fit. We have not included scaling corrections. To estimate their role, we repeat the fits for $L_{\rm min}=16,20,24$, each time only including data corresponding to $L\geq L_{\rm min}$.

We first consider $\xi_{q}/L$. For each $L_{\rm min}$ we increase systematically the order $n$ till results and the $\chi^{2}$ per degree of freedom are stable. For all values of $L_{\rm min}$ this is obtained by taking $n\approx 5$-6, see Table 1. The quality of the fits is reasonable and results are stable with respect to $L_{\rm min}$, indicating that scaling corrections are small, compared to the statistical errors. From this analysis we would estimate $\beta_{c}=0.22576(2)$ and $\nu=0.67(3)$.

The estimates obtained analyzing $U_{q,4}$ are worse. Independently of $n$, the quality of the fits is poor, see Table 1. Moreover, results for different values of $L_{\rm\min}$ show a systematic drift, a clear signal of the presence of systematic errors due to the neglected scaling corrections. To obtain more reliable estimates, it is therefore necessary to include scaling corrections. We have thus performed a combined analysis of $R_{\xi}$ and $U_{q,4}$ in which $R_{\xi}$ has been fitted to (20), while $U_{q,4}$ has been fitted to

[TABLE]

Fits are stable only if we take $L_{\rm min}=16$, that is, if we include all data. Using $n=6$ for $R_{\xi}$ and $n=8$, $m=2$ for $U_{q,4}$ we obtain a good-quality fit: $\chi^{2}$/DOF = 1.14. Correspondingly, we obtain the estimates

[TABLE]

Moreover, we have $R_{\xi}^{*}=0.1728(6)$ and $U_{q,4}^{*}=5.20(6)$. The estimates of $\beta_{c}$ and $\nu$ are perfectly consistent with those obtained before from the analysis of $R_{\xi}$ alone. To verify the reliability of the errors we have performed a second fit. Both observables are fitted to (21), with $n=6$, $m=2$ for $R_{\xi}$ and $n=8$, $m=2$ for $U_{q,4}$. We obtain a sligthly better $\chi^{2}$/DOF, 1.10, and

[TABLE]

Moreover, we have $R_{\xi}^{*}=0.182(7)$ and $U_{q,4}^{*}=5.43(20)$. The new estimates are essentially consistent with the results reported in (22), differences being of the order of the sum of the statstical errors.

In fit (21) we have included nonanalytic corrections related to the leading irrelevant operator. Another source of scaling corrections is the analytic background. For $R_{\xi}$ and $U_{q,4}$ such corrections decay as $L^{\eta_{q}-2}$, where $\eta_{q}$ is the susceptibility exponent. Usually, such an exponent is small and analytic corrections are negligible. In the present case, however, as we discuss in Sec. 4.3, $\eta_{q}\approx 1$ and thus they are as important as the nonanalytic ones. We have therefore performed fits to

[TABLE]

where we have only considered the analytic terms in $(\beta-\beta_{c})$ (in principle, the coefficient function in front of $L^{\eta_{q}-2}$ has an expansion in powers of both $(\beta-\beta_{c})$ and $X$). We fix $\eta_{q}=1.05$, see Sec. 4.3. A combined analysis of the two observables, using $n=6$, $m=1$ for $R_{\xi}$ and $n=8$, $m=1$ for $U_{q,4}$ gives $\chi^{2}$/DOF = 1.25 and

[TABLE]

$R_{\xi}^{*}=0.172(3)$ and $U_{q,4}^{*}=5.26(3)$. These results are fully consistent with the previous ones. Collecting all results we end up with the estimates

[TABLE]

It is quite interesting to observe that the critical point position is close to that of the pure Ising model ($\beta_{C,I}\simeq 0.22165$). On the other hand, $R_{\xi}^{*}$ and $U_{q,4}^{*}$ are distinctly different since $R_{\xi}^{*}\approx 0.4$ and $U_{q,4}\approx 2.4$ in the Ising model.

In figure 3 we report the estimates of $R_{\xi}$ and $U_{q,4}$ as a function of $X=(\beta-\beta_{c})L^{1/\nu}$. For $R_{\xi}$ we observe a very good scaling: data are clearly consistent with the predicted scaling behavior. On the other hand, the Binder parameter shows scaling corrections, that appear to be particularly large where the scaling function has a maximum, for $X\approx 0.1$.

4.2 Behavior of the overlap observables

at the finite-size pseudocritical transition

Let us now consider the behavior of the system at the finite-size pseudocritical transition $\beta_{G}(L)$. The data for $R_{\xi}=\xi_{q}/L$ and $U_{q,4}$ are shown in figure 4. Note that in this case we do not have data for $L=34$, we are considering a significantly smaller interval of values of $\beta$, and the the number of samples is significantly smaller than for $\beta_{G}=\beta_{G,c}$. Therefore, we will use the results at $\beta_{G}(L)$ only as a consistency check of those obtained at the critical point, and in particular of the final estimates (26).

As expected, data show an intersection for $\beta\simeq 0.226$, in agreement with the general theory. We analyze $R_{\xi}$ and $U_{q,4}$ as before, performing fits to Eq. (20). Note that the constant $R^{*}$ should take here a value which is different from that it has for $\beta_{G}=\beta_{G,c}$, as here $R^{*}=f_{R}(X_{\rm max},0)$, see Eq. (14). Results are reported in Table 2.

The quality of the fits of $R_{\xi}$ is poor, with no improvement when data with $L=16$ are discarded. While the estimates of $\nu$ are consistent within errors with (26), the estimates of $\beta_{c}$ are significantly larger (by several error bars) than that obtained at the critical point $\beta_{G,c}$. The results obtained from the analysis of $U_{q,4}$ are apparently better. For $L_{\rm min}=20$, the quality of the fits is reasonable and the estimates of $\nu$ and $\beta_{c}$ are both consistent with (26).

We do not have enough data here to perform an analysis that includes scaling corrections with an unconstrained exponent $\omega$, as performed at $\beta_{G,c}$. Stable results are only obtained if we fix $\omega$. We fix it to the central estimate (26), $\omega=1.0$, and we perform a combined analysis of the two observables to (21). Using $n=4$ and $m=2$ for both observables, we obtain $\beta_{c}=0.225857(33)$ and $\nu=0.66(10)$ with $\chi^{2}$/DOF$\approx 1.3$. If we decrease $\omega$ by one error bar, the estimates change by much less than the quoted errors. The estimate of $\beta_{c}$ is essentially consistent within error bars with Eq. (26). This analysis clearly confirms that the discrepancy observed for $\beta_{c}$ in the analysis of $R_{\xi}$ can be ascribed to the neglected scaling corrections.

We show the data as a function of $X=(\beta-\beta_{c})L^{1/\nu}$ in figure 5, using the estimates (26). Scaling corrections are clearly visible, but data are consistent with a universal scaling.

4.3 Overlap susceptibility exponent

It is also interesting to analyze the behavior of the susceptibility $\chi_{q}$ to determine the exponent $\eta_{q}$. For $\beta_{G}=\beta_{G,c}$ we can expand the scaling function defined in (19) in powers of $X$, to obtain

[TABLE]

The last term, proportional to $(\beta-\beta_{c})$ represents the first term in the expansion of the nonlinear scaling field $u_{\chi}$. If we fix $\beta_{c}$ and $\nu$ to the central values (26), we obtain $\eta_{q}=1.055(6)$, 1.059(8), 1.038(15) for $L_{\rm min}=16,20,24$, respectively. If we change $\nu$ by one error bar, results are essentially unchanged, while they vary by less than 0.04 if we change $\beta_{c}$ by the same amount. Thus, we end up with the final estimate

[TABLE]

We have also performed fits including scaling corrections. For the susceptibility there are two sources of corrections. First, there are the nonanalytic corrections that decay as $L^{-\omega}$, as discussed for the RG invariant quantities. Second, there are corrections due to the background that scale as $L^{\eta_{q}-2}$. The exponents $\omega$ and $\eta_{q}-2$ are close and it is not possible to include them both in the fitting function. Therefore, we have performed two different fits:

[TABLE]

and

[TABLE]

Using the central values for $\beta_{c}$, $\nu$, and $\omega$ [see (26)], $L_{\rm min}=16$ and $n=6$, $m=2$, fits (29) and (30) both give $\eta_{q}=1.04(5)$. Results are perfectly consistent with (28).

The value of $\eta_{q}$ is close to the value it assumes in the ferromagnetic Ising model. Indeed, the exponent $\eta_{q}$ for the overlap susceptibility is related to the exponent $\eta_{m}$ for the magnetic susceptibility by $\eta_{q}=1+2\eta_{m}$, so that for the Ising model we have $\eta_{q}=1.0725(2)$ [we use [23] $\eta_{m}=0.03627(10)$].

As a consistency check, we have analyzed the data obtained setting $\beta_{G}=\beta_{G}(L)$. The quality of the fits to (27) is relatively poor, $\chi^{2}$/DOF$\approx 1.9$. We obtain $\eta_{q}=1.226(3)$ and $\eta_{q}=1.217(4)$ for $L_{\rm min}=16,20$, respectively. The results are significantly different from those obtained at $\beta_{G,c}$ at the level of the statistical errors. This clearly indicates that there are significant scaling corrections. We have therefore performed fits including scaling corrections. Using the central values for $\beta_{c}$, $\nu$, and $\omega$ [see (26)], fits (29) and (30) both give $\eta_{q}=1.13(3)$. In both cases $\chi^{2}$/DOF is approximately 1.4. The result is essentially consistent with (28) once we take into account the statistical errors.

4.4 Overlap distribution

To better characterize the nature of the transition points we have determined the overlap distribution defined by

[TABLE]

where $Q=\sum_{x}q_{x}$. The results for $L=20$ and 24 and three slightly different values of $\beta$ at $\beta_{G}=\beta_{G,c}$ are reported in figure 6. The distribution is distinctly different from a Gaussian. There is a significant tail whose importance increases, as $\beta$ increases. This tail is responsible for a critical value $U^{*}_{q,4}$ significantly larger than 3 (the Gaussian value). The function $P(q,\beta)$ should satisfy a general scaling relation. We define the overlap magnetization exponent $\beta_{q}$, which is related to $\eta_{q}$ by the scaling relation

[TABLE]

where $d$ is the space dimension ($d=3$ in our case). Using (28) we find

[TABLE]

In terms of $\beta_{q}$ we have

[TABLE]

where $f_{P}(x,X)$ is a universal (apart from a rescaling of its arguments) function, satisfying

[TABLE]

In figure 6 we report a scaling plot for $\beta=\beta_{c}$ (i.e., for $X=0$). Data for two different values of $L$ fall one on top of the other, confirming the correctness of our estimates of $\eta_{q}$ and $\beta_{q}$.

The overlap distribution has a peculiar shape that seems to indicate the presence of two different behaviors, depending on the disorder distribution. To obtain a better understanding, consider the sample-dependent distribution

[TABLE]

where we only perform the thermal average. The results of figure 6 are consistent with the idea that, for most bond configurations, $P_{0}(q,\{J\})$ is peaked around $q=0$. However, the shape of the tail seems to indicate the presence of rare configurations for which $P_{0}(q,\{J\})$ is quite different. We have looked at the distributions $P_{0}(q,\{J\})$ for several bond distributions and we have indeed identified two different typical shapes, see figure 7. In most of the cases $P_{0}(q,\{J\})$ is clearly peaked around $q=0$, but in a few cases the distribution is broader with a peak for $q\not=0$.

To better single out the presence of two different types of configurations we compute

[TABLE]

As we consider $\langle|Q|\rangle/V$, we expect a sharp peak corresponding to bond configurations such that $P_{0}(q,\{J\})$ is peaked in $q=0$ and a wide tail that gets contributions from configurations with a peak for $q\not=0$. The results, shown in Fig. 8, confirm the general analysis. In figure 8, we also report a scaling plot, that confirms the asymptotic nature of the results.

5 Low-frustration regime

In the low-temperature phase of the gauge model, that is for $\beta_{G}>\beta_{G,c}$, frustration is suppressed. In order to study the intermediate region between $\beta_{G,c}$ and the zero-frustration limit $\beta_{G}=\infty$, we perform simulations for $\beta_{G}=0.9$, considering $16\leq L\leq 30$. For this value of $\beta_{G}$ sample-to-sample fluctuations are small and therefore we obtain significantly more precise results than at $\beta_{G}=\beta_{G,c}$. For $\beta_{G}=0.9$, the fraction of frustrated plaquettes is approximately $0.4\%$, with a weak $L$ dependence.

The results for $\xi_{q}/L$ and $U_{4}^{q}$ have been analyzed as in the previous sections. The $\chi^{2}$ of the fits without scaling corrections is large, so that we have been forced to include them in the analysis. Combined fits of the two observables to (21) (we use $n=8$ and $m=2$ for both observables) give $\beta_{c}=0.222234(2)$, $\nu=0.628(2)$, and $\omega=0.99(15)$. The estimates of $\nu$ and $\omega$ are very close to the pure Ising values [23] $\nu_{I}=0.63002(10)$ and $\omega_{I}=0.832(6)$. These results make us conjecture that the critical behavior of the model for $\beta_{G}=0.9$ is in the same universality class as that of the pure Ising model. If this conjecture is correct, we can obtain a more precise estimate of the critical point, reanalyzing the data fixing $\nu$ and $\omega$ to the Ising values of Ref. [23]. We obtain $\beta_{c}=0.222264(4)$ including all data, and $\beta_{c}=0.222262(9)$ if we only consider data with $L\geq 20$. To provide stronger evidence for the conjecture, we compare in figure 9 the scaling functions in the two models: they should agree, apart from a rescaling of the abscissa. For the Binder cumulant $U_{4,q}$, data fall on top of each other with good accuracy, without the need of any rescaling. On the other hand, the data for $R_{\xi}$ in the two models show a tiny difference that might be due to scaling corrections, that are quite sizable for the range of values we consider, as we discussed.

The results shown here apparently indicate that a small amount of (correlated) frustration does not change the universality class: The system at $\beta_{G}=0.9$ behaves as the pure Ising model. However, the reader should be aware that lattices are quite small and, therefore, conclusive evidence can only be obtained by considering significantly larger lattices. This check should be feasible as the ferromagnetic nature of the transition shoud make cluster algorithms [26] quite efficient.

6 Conclusions

In this paper we have considered a disordered Ising model in which the amount of frustration can be continuously varied in a gauge-invariant way by tuning a single parameter $\beta_{G}$, that plays the role of a chemical potential for the frustrated plaquettes. By changing $\beta_{G}$ we can interpolate between the usual spin-glass Edwards-Anderson model ($\beta_{G}=0$) and the pure Ising model. However, the resulting phase diagram is quite different from that observed in the usual $\pm J$ model with bond distribution (2). In that model, for $p\leq 1/2$, there are only two distinct phases, a ferromagnetic one in the same universality class as that of the randomly dilute Ising (RDI) model and a spin-glass phase separated by a multicritical point. The phase diagram of the model we consider here is apparently richer. For $\beta_{G}=\beta_{G,c}$ we identify a new universality class characterized by the fact that the distribution of frustration is long-range correlated. We compute the correlation-length exponent $\nu$ and the overlap-susceptibility exponent $\eta_{q}$, obtaining

[TABLE]

Given the very small amount of frustration (the fraction of frustrated plaquettes is 2.6% for $\beta_{G}=\beta_{G,c}$), we expect the transition to be ferromagnetic in an appropriate set of variables. Note that, because of gauge invariance, one cannot consider correlations of the spins. For instance, it is trivial to show that $[\langle\sigma_{x}\sigma_{y}\rangle]=\delta_{xy}$. One possibility consists in considering gauge-invariant observables of the form

[TABLE]

where $\mathcal{C}$ is a lattice path that connects sites $x$ and $y$. When $\beta_{G}=\infty$, we recover the usual definition of the magnetic correlation function in an appropriate gauge. Indeed, in this limit $\Pi_{P}=1$ on all plaquettes and therefore by means of a gauge transformation we can set $J_{xy}=1$ on all links. In this gauge the Hamiltonian of the system is that of the pure Ising model and the correlation function becomes $\langle\sigma_{x}\sigma_{y}\rangle$. Observables defined as in (39) were already considered in [27], where it was shown that they had the expected scaling behavior provided the path $\mathcal{C}$ was appropriately chosen.

For $\beta_{G}>\beta_{G,c}$ our preliminary data are consistent with pure Ising behavior, but significant larger systems are needed to draw a definite conclusion. The behavior for $\beta_{G}<\beta_{G,c}$ has not been investigated yet. One possibility is that in the whole interval $0\leq\beta_{G}<\beta_{G,c}$ there is only a spin glass phase. However, at $\beta_{G,c}$ frustration is quite small: the percentage of frustrated plaquettes is only 2.6%, which is significantly smaller than the fraction, 46.2%, at $p^{*}$ for the EA model with disorder distribution (2). Thus, it is possible that in this interval the behavior is analogous to that observed in the $\pm J$ model with uncorrelated disorder: a spin-glass phase for $\beta_{G}<\beta_{G,MC}$ and a ferromagnetic RDI phase for $\beta_{G,MC}<\beta_{G}<\beta_{G,c}$. Finally, it would be interesting to investigate the behavior for $\beta_{G}<0$. For $\beta_{G}=-\infty$ we would obtain the fully frustrated Ising model. For the cubic lattice, all numerical simulations [28, 29, 30, 31, 32] are consistent with the presence of a second-order high-temperature transition (however, so far, the transition has not been explained within the Landau-Ginzburg-Wilson approach, as no appropriate fixed point has been found [33]) and of a first-order low-temperature transition. It is tempting to assume that the two transitions survive if we take $|\beta_{G}|$ large but finite. In this case, the phase diagram would be quite complex as the two transition lines should necessarily merge as $|\beta_{G}|$ decreases.

The simulations presented here have been performed on the Theory cluster at INFN, Sezione di Roma 1.

Appendix A Computational details

The Ising model (1) has been simulated using a standard Metropolis multispin code that deals with 64 systems at the same time, while the bond configurations have been generated using a standard Metropolis algorithm. In practice, we first perform a simulation of the gauge model. Every $N_{\rm it}$ iterations (one iteration consists in proposing an update of every link variable) a bond configuration is saved. Once 64 bond configurations have been generated, they are packed in a single 64-bit word and a simulation of the spin system is performed. Two replicas are simulated together, in order to compute functions of the overlap. The parameters of the different runs are reported in Tables 3 and 4.

The parameter $N_{\rm it}$ has been chosen so that successive bond configurations are approximately independent. For this purpose, we have computed the autocorrelation time [34] $\tau$ for the average energy and for the average of the so-called Polyakov loop $P$. The Polyakov loop is defined as the product of the bond variables associated with a topologically nontrivial line that is parallel to a coordinate direction and connects the opposite boundaries of the lattice. The Polyakov loop is usually considered as the order parameter of the transition in the gauge model. The dynamics of the Polyakov loop is slower than that of the energy and, for $\beta_{G}=\beta_{G,c}$, its correlation times $\tau$ is approximately $26,57,115,171$ for $L=10,14,18,23$ (here $\tau$ is the integrated autocorrelation time as defined in [34]). The results are well fitted by $0.12L^{2.33}$. Note that this is simply a phenomenological fit, which should be correct in the interval of values of $L$ of interest. Indeed, for large $L$ we expect $\tau$ to increase with an exponent very close to 2 [35]. In the calculations at $\beta_{G,c}$ we have taken $N_{\rm it}=\tau$, while in the runs at $\beta_{G}(L)$ we have taken $N_{\rm it}\approx\tau/2$. Since we obtain the bond configurations in a dynamic Monte Carlo simulation, results corresponding to different samples are not uncorrelated. To obtain correct estimates of the statistical errors, they have been determined using a blocking method. Errors have been obtained from the variance of averages over a large (typically of the order of 6400) number of samples. In disordered systems, thermalization is crucial. We discarded in all cases 640 iterations. With this choice there is no evidence of initialization bias. This has been verified in detail in the runs with $\beta_{G}=\beta_{G}(L)$, which are significantly longer than needed.