arXiv:0704.0427v1  [cond-mat.dis-nn]  3 Apr 2007
The 3D ±J Ising model at the ferromagnetic transition line
Martin Hasenbusch,1 Francesco Parisen Toldin,2 Andrea Pelissetto,3 and Ettore Vicari1
1 Dipartimento di Fisica dell’Università di Pisa and INFN, Pisa, Italy.
2 Scuola Normale Superiore and INFN, Pisa, Italy.
3Dipartimento di Fisica dell’Università di Roma “La Sapienza” and INFN, Roma, Italy.
(Dated: November 1, 2018)
Abstract
We study the critical behavior of the three-dimensional ±J Ising model [with a random-exchange
probability P (Jxy) = pδ(Jxy−J)+(1−p)δ(Jxy+J)] at the transition line between the paramagnetic
and ferromagnetic phase, which extends from p = 1 to a multicritical (Nishimori) point at p =
pN ≈ 0.767. By a finite-size scaling analysis of Monte Carlo simulations at various values of p
in the region pN < p < 1, we provide strong numerical evidence that the critical behavior along
the ferromagnetic transition line belongs to the same universality class as the three-dimensional
randomly-dilute Ising model. We obtain the results ν = 0.682(3) and η = 0.036(2) for the critical
exponents, which are consistent with the estimates ν = 0.683(2) and η = 0.036(1) at the transition
of randomly-dilute Ising models.
PACS numbers: 75.10.Nr, 75.40.Cx, 75.40.Mg, 64.60.Fr
http://arxiv.org/abs/0704.0427v1
I. INTRODUCTION
The ±J Ising model has played an important role in the study of the effects of quenched
random disorder and frustration on Ising systems. It is defined by the lattice Hamiltonian
H±J = −
Jxyσxσy, (1)
where σx = ±1, the sum is over the nearest-neighbor sites of a simple cubic lattice, and the
exchange interactions Jxy are uncorrelated quenched random variables, taking values ±J
with probability distribution
P (Jxy) = pδ(Jxy − J) + (1− p)δ(Jxy + J). (2)
For p = 1 we recover the standard Ising model, while for p = 1/2 we obtain the usual
bimodal Ising spin-glass model.
The phase diagram of the three-dimensional (3D) ±J Ising model is sketched in Fig. 1.
The high-temperature phase is paramagnetic for any p. The low-temperature phase depends
on the value of p: it is ferromagnetic for small values of 1 − p, while it is a spin-glass
phase with vanishing magnetization for sufficiently large values of 1 − p. The different
phases are separated by transition lines, which meet at a multicritical point N located
along the so-called Nishimori line.1,2,3 The spin-glass transition has been mostly studied at
the symmetric point p = 1/2, see, e.g., Refs. 3,4 and references therein. The spin-glass
transition line extends up to the Nishimori multicritical point,2 located at5,6,7,8 pN ≈ 0.767.
For larger values of p, the transition is ferromagnetic, up to p = 1 where one recovers the pure
Ising model, and therefore a transition in the Ising universality class. At the ferromagnetic
transition line, for pN < p < 1, the critical behavior is expected to belong to a different
universality class.
An interesting hypothesis, which has already been put forward in Refs. 3,9, is that the
ferromagnetic transition of the ±J Ising model belongs to the 3D randomly-dilute Ising
(RDIs) universality class (see, e.g., Refs. 10,11 for reviews on randomly-dilute spin mod-
els). A representative of the RDIs universality class is the randomly site-dilute Ising model
(RSIM) defined by the lattice Hamiltonian
Hd = −J
ρx ρy σxσy, (3)
glass
ferro
1 − p
FIG. 1: Sketch of the phase diagram of the 3D ±J Ising model in the T -p plane.
where ρx are uncorrelated quenched random variables, which are equal to 0, 1 with proba-
bility
P (ρx) = pδ(ρx − 1) + (1− p)δ(ρx). (4)
For p < 1 and above the percolation threshold of the spins (pperc ≈ 0.3116081(13) on a cubic
lattice12), the RSIM undergoes a continuous phase transition between a disordered and a
ferromagnetic phase, whose nature is independent of p. This transition is definitely different
from the usual Ising transition: for instance, the correlation-length critical exponent13,14,15,16
ν = 0.683(2) differs from the Ising value17,18 ν = 0.63012(16). The RDIs universality class is
expected to describe the ferromagnetic transition in generic diluted ferromagnetic systems.
For instance, it has been verified that also the randomly bond-diluted Ising model (RBIM)
belongs to the RDIs universality class.13,19 These results do not necessarily imply that also
the ±J Ising model has an RDIs ferromagnetic transition line. Indeed, while the RSIM (3)
has only ferromagnetic exchange interactions, the ±J Ising model is frustrated for any value
of p < 1. Therefore, the ferromagnetic transition in the ±J Ising model belongs to the RDIs
universality class only if frustration is irrelevant, a fact that is not obvious and should be
carefully investigated.
Reference 9 investigated the issue by means of a Monte Carlo (MC) renormalization-
group (RG) study, claiming that the ±J Ising model belongs to the same RDIs universality
class as the RSIM and the RBIM. It should be noted however that the quoted estimate for
the correlation-length exponent at the ferromagnetic transition, ν = 0.658(9), is close to
but not fully consistent with the RDIs value ν = 0.683(2).13 Another numerical MC work20
investigated the nonequilibrium relaxation dynamics of the ±J Ising model and showed an
apparent nonuniversal dynamical critical behavior along the ferromagnetic transition line.
These results are not conclusive and further investigation is called for to clarify this issue.
In this paper we focus on the transition line of the 3D ±J Ising model between the
paramagnetic and the ferromagnetic phase. We investigate the critical behavior by means
of MC simulations at various values of p in the region pN < p < 1. Our finite-size scaling
(FSS) analysis provides a strong evidence that the critical behavior of the 3D ±J Ising
along the ferromagnetic line belongs to the 3D RDIs universality class. For example, we
obtain ν = 0.682(3) and η = 0.036(2), which are in good agreement with the presently most
accurate estimates13 ν = 0.683(2) and η = 0.036(1) for the 3D RDIs universality class.
The paper is organized as follows. In Sec. II we summarize some FSS results which are
needed for the analysis of the MC data, and describe our strategy to check whether the
transition belongs to the RDIs universality class. In Sec. III we describe the MC simula-
tions. In Sec. IV we report the results of the FSS analysis. Finally, in Sec. V we draw our
conclusions. In App. A we report the definitions of the quantities we compute.
II. STRATEGY OF THE FINITE-SIZE SCALING ANALYSIS
In this work we check whether the ferromagnetic transition line in the 3D ±J Ising models
belongs to the RDIs universality class. For this purpose, we present a FSS analysis of MC
data for various values of p in the region 1 > p > pN ≈ 0.767. We follow closely Ref. 13,
which studied the ferromagnetic transition line in the 3D RSIM and RBIM and provided
strong numerical evidence that these transitions belong to the same RDIs universality class.
We refer to Ref. 13 for notations (a short summary is reported in App. A) and a detailed
discussion of FSS in these disordered systems.
According to the RG, in the case of periodic boundary conditions and for L → ∞, where
L is the lattice size, a generic RG invariant quantity R at the critical temperature 1/βc
behaves as
R(L, β = βc) = R
1 + c11L
−ω + c12L
−2ω + · · ·+ c21L
−ω2 + · · ·
, (5)
where R∗ is the universal infinite-volume limit and ω and ω2 are the leading and next-
to-leading correction-to-scaling exponents. In RDIs systems scaling corrections play an
important role,16,21 since ω is quite small. Indeed we have ω = 0.33(3) and ω2 = 0.82(8)
in the 3D RDIs universality class.13 These slowly-decaying scaling corrections make the
accurate determination of the universal asymptotic behavior quite difficult.
Instead of computing the various quantities at fixed Hamiltonian parameters, we keep
a RG invariant quantity R fixed at a given value Rf .
22 This means that, for each L, we
determine the pseudocritical inverse temperature βf(L) such that
R(β = βf(L), L) = Rf . (6)
All interesting thermodynamic quantities are then computed at β = βf (L). The pseudocrit-
ical inverse temperature βf(L) converges to βc as L → ∞. The value Rf can be specified
at will, as long as Rf is taken between the high- and low-temperature fixed-point values
of R. The choice Rf = R
∗ (where R∗ is the critical-point value) improves the conver-
gence of βf to βc for L → ∞; indeed βf − βc = O(L
−1/ν) for generic values of Rf , while
βf −βc = O(L
−1/ν−ω) for Rf = R
∗. This FSS method has already been applied to the study
of the critical behavior of N -vector spin models,22,23 and of randomly-dilute Ising models.13
As in Ref. 13, we perform a FSS analysis at fixed Rξ ≡ ξ/L = 0.5943, which is very close to
the fixed-point value R∗ξ = 0.5944(7) of Rξ at βc. Given any RG invariant quantity R, such as
the quartic cumulants U4 and U22, we consider its value at fixed Rξ, i.e., R̄(L) = R(L, βf(L)).
For L → ∞, R̄(L) behaves as R(L, βc):
R̄(L) = R̄∗
1 + b11L
−ω + b12L
−2ω + · · ·+ b21L
−ω2 + · · ·
, (7)
where the coefficients bij depend on the Hamiltonian. The derivative R̄
′ with respect to β
of a generic RG invariant quantity R behaves as
R̄′(L) = aL1/ν
1 + a11L
−ω + a12L
−2ω + · · ·+ a21L
−ω2 + · · ·
. (8)
Finally, the FSS of the magnetic susceptibility χ is given by13
χ̄(L) ≡ χ(L, β = βf (L)) = eL
1 + e11L
−ω + e12L
−2ω + · · ·+ e21L
−ω2 + · · ·
+ eb (9)
where eb represents the background contribution.
A standard RG analysis, see, e.g., Ref. 13, shows that the amplitudes of the O(L−kω)
scaling corrections are proportional to uk3 (with a universal coefficient), where u3 is the
leading irrelevant scaling field with RG dimension y3 = −ω. Hamiltonians such that u3 = 0—
we call them improved Hamiltonians—have a faster approach to the universal asymptotic
behavior, because the O(L−kω) scaling corrections vanish: b1k = a1k = e1k = 0 in Eqs. (7),
(8), and (9). In this case the leading scaling corrections are proportional to u4L
−ω2 , where u4
is the next-to-leading irrelevant scaling field and y4 = −ω2 is its RG dimension. In Ref. 13
it was shown that the RSIM for p = p∗ = 0.800(5) and the RBIM for p = p∗ = 0.56(2)
are improved. Since scaling fields are analytic functions of the Hamiltonian parameters, u3
must be proportional to p − p∗ close to p = p∗, i.e. u3 ≈ c3(p − p
∗). Therefore, since the
coefficients b1k, a1k, and e1k that appear in Eqs. (7), (8), and (9) are proportional to u
3, we
b1k, a1k, e1k ∼ (p− p
∗)k. (10)
Beside the quantities defined in App. A, we also consider observables—in analogy with
the previous terminology, we call them improved quantities—characterized by the fact that
the leading scaling correction proportional to L−ω (approximately) vanishes in any model
belonging to the RDIs universality class.13 We consider the combination of quartic cumulants
Ūim = Ū4 + 1.3Ū22, (11)
and improved estimators of the critical exponent ν defined as
R′ξ,im ≡ R̄
d , U
4,im ≡ Ū
d (12)
(Ūd is defined in App. A). In Ref. 13 we showed that, if the transition belongs to the RDIs
universality class, the leading scaling correction proportional to L−ω of these improved ob-
servables is suppressed. More precisely, we showed that the universal ratio of the amplitudes
of the leading scaling correction in Ūim and Ū4 satisfies
|b11,Ūim/b11,Ū4 | .
, (13)
while the one for the quantities R′ξ,im and R̄
ξ is bounded by
|a11,R′
/a11,R̄′
. (14)
The remaining scaling corrections are of order L−2ω and L−ω2. These improved observables
are particular useful to check whether the transition in a given system belongs to the 3D
RDIs universality class.
10 100
FIG. 2: Exponential autocorrelation time τ of the magnetic susceptibility for a mixture of
Metropolis and cluster updates as discussed in the text, at p = 0.87. The dotted line shows
the result of a fit to τ = cLz: this fit gives z ≈ 1.6.
To summarize: in order to check whether the ferromagnetic transition of the 3D ±J
Ising model belongs to the RDIs universality class, we perform a FSS analysis at fixed
Rξ = 0.5943, and check if the results for the critical exponents and other universal quantities
are consistent with those obtained for the RDIs universality class, which is characterized
by13 critical exponents ν = 0.683(2) and η = 0.036(1), by the leading and next-to-leading
scaling-correction exponents ω = 0.33(3) and ω2 = 0.82(8) and by the universal infinite-
volume values of the quartic cumulants Ū∗22 = 0.148(1), Ū
im = 1.840(4), and Ū
d = 1.500(1).
Notice that the fact that we fix Rξ = 0.5943 does not introduce any bias in our FSS analysis.
III. MONTE CARLO SIMULATIONS
We performed MC simulations of Hamiltonian (1) with J = 1 for p =
0.94, 0.90, 0.883, 0.87, 0.83, 0.80, close to the critical temperature on cubic lattices of size
L3 with periodic boundary conditions, for a large range of lattice sizes: from L = 8 to
L = 80 for p = 0.883, 0.87, to L = 64 for p = 0.94, 0.90, 0.83, and to L = 48 for p = 0.80. We
chose values of p not too close to p = 1: indeed, as p → 1 we expect crossover effects due to
the presence of the Ising transition for p = 1 and, therefore, that the asymptotic behavior
sets in only for large values of L. We return to this point later.
We used a Metropolis algorithm and multispin coding.24 In the simulation nbit systems
evolve in parallel, where nbit = 32 or nbit = 64 depending on the computer that is used.
For each of these nbit systems we use a different set of couplings Jxy. This allows us to
perform 64 parallel simulations on a 64-bit machine, and therefore to gain a large factor
in the efficiency of the MC simulations. We used high-quality random-number generators,
such as the RANLUX25 or the twister26 generators.27 Using the twister random-number
generator, we need about 1.2 × 10−9 seconds for one Metropolis update of a single spin on
an Opteron processor running at 2 GHz. Our simulations took approximately 3 CPU years
on an Opteron (2 GHz) processor.
It is worth mentioning that cluster algorithms, such as the Swendsen-Wang cluster28
and the Wolff single-cluster29 algorithm, show significant slowing down in the ±J Ising
model. At the earlier stage of this work we performed some simulations of the ±J Ising
model at p = 0.87 using the algorithm used in Ref. 13 to simulate the RSIM and the
RBIM. There we used a combination of Metropolis, Swendsen-Wang cluster,28 and Wolff
single-cluster29 updates. More precisely, each updating step consisted of 1 Swendsen-Wang
update, 1 Metropolis update, and L single-cluster updates. In all cases the exponential
autocorrelation times τ of the magnetic susceptibility was small: τ . 1 in units of the
above updating step, even for the largest lattice sizes considered, i.e. L = 192. In the ±J
Ising model at p = 0.87 autocorrelation times are much larger. In Fig. 2 we plot estimates
of τ as obtained from the magnetic susceptibility. They show a clear evidence of critical
slowing down: τ ∼ Lz with z ≈ 1.6. Such a value of z should be compared with the
dynamic exponent of Swendsen-Wang and Wolff cluster algorithms in the RSIM, which is
much smaller:30 z . 0.5. These results show that cluster algorithms behave differently in
the ±J Ising model, likely due to frustration. They suggest that frustration is relevant for
the cluster dynamics.
Taking also into account the computer time required by the cluster algorithms, we then
turned to a multispin Metropolis algorithm. This turns out to be much more effective at the
lattice sizes considered, although it has a larger dynamic exponent z & 2, see, e.g., Ref. 30
and references therein. We also mention that the autocorrelation time significantly increases
with decreasing p (keeping L fixed). For example, for L = 48 it increases by approximately
a factor of 10 from p = 0.90 to p = 0.80. This represents a major limitation to perform
simulations for large lattices close to the multicritical point.
For each lattice size we considered Ns disorder samples, with Ns decreasing with in-
creasing L, from Ns & 10
6 for L = 8 to Ns & 2 × 10
4 for the largest lattices. For each
disorder sample, we collected a few hundred independent measurements at equilibrium. The
averages over disorder are affected by a bias due to the finite number of measures at fixed
disorder.13,31 A bias correction is required whenever one considers the disorder average of
combinations of thermal averages. We used the formulas reported in App. B of Ref. 13.
Errors were computed from the sample-to-sample fluctuations and were determined by using
the jackknife method.27
Our FSS analysis is performed at fixed Rξ ≡ ξ/L. In order to determine expectation
values at fixed Rξ, one needs the values of the observables as a function of β in some
neighborhood of the inverse temperature βrun used in the simulation. In Ref. 13 we used the
reweighting method for this purpose. This requires that the observables and, in particular,
the values of the energy are stored at each measurement. For the huge statistics like those we
have for the smaller values of L, this becomes unpractical. Therefore, we used here a second-
order Taylor expansion, determining O(β, L) from O(βrun, L)+aO(β−βrun)+ bO(β−βrun)
The coefficients aO and bO are obtained from appropriate expectation values as in Ref. 23.
Since their computation involves disorder averages of products of thermal averages, we have
implemented in all cases an exact bias correction, using the formulas of Ref. 13. Derivatives
with respect to β are then obtained as O′(β, L) = aO + 2bO(β − βrun). Of course, this
method requires |βrun − βf | to be sufficiently small. We have carefully checked the results
by performing, for each L and p, runs at different values of β.
The MC estimates of the quantities introduced in Sec. II and in App. A at fixed Rξ ≡
ξ/L = 0.5943 are available on request.
IV. FINITE-SIZE SCALING ANALYSIS
In this section we present the results of our FSS analysis of the MC data at fixed Rξ =
0.5943.
0.0 0.1 0.2 0.3 0.4 0.5
-0.33
p=0.90
p=0.883
p=0.87
p=0.83
p=0.80
FIG. 3: MC estimates of Ū22 versus L
−ω with ω = 0.33 for different values of p. The dotted lines
show results of fits to c0 + c1L
−ε1 + c2L
−ε2 , fixing c0 = 0.148, ε1 = 0.33, and ε2 = 0.82. In the
RDIs universality class Ū∗22 = 0.148(1).
A. Renormalization-group invariant quantities
In Fig. 3 we show the MC estimates of Ū22 versus L
−ω with ω = 0.33(3), which is the
leading scaling exponent of the RDIs universality class. The data vary significantly with p
and L. This p and L dependence is always consistent with the existence of the expected
next-to-leading scaling corrections, i.e. with a behavior of the form
Ū22 = Ū
22 + c1L
−ε1 + c2L
−ε2 , (15)
where Ū∗22, ε1 and ε2 are fixed to the RDIs values:
13 Ū∗22 = 0.148, ε1 = 0.33 and ε2 =
0.66, 0.82. The fits corresponding to ε2 = 0.82 are shown in Fig. 3. Note that in most of
the cases it is crucial to include a next-to-leading correction. Only for p = 0.90 the data are
well fitted by taking only the leading scaling correction.
An unbiased estimate of ω can be obtained from the difference of data at different values
of p, i.e. by considering
Ū22(p1;L)− Ū22(p2;L) ≈ cL
−ω. (16)
Linear fits of the logarithm of these differences give results in reasonable agreement with
0.00 0.05 0.10 0.15 0.20
-0.82
p=0.90
p=0.883
p=0.87
p=0.83
p=0.80
FIG. 4: MC estimates of Ūim versus L
−0.82. The filled square on the vertical axis corresponds to
the RDIs estimate13 Ū∗im = 1.840(4).
the RDIs estimate ω = 0.33(3), especially when only data corresponding to L ≥ Lmin = 24
are used. For Lmin = 24 [Lmin = 32], we obtain ω = 0.27(2) [ω = 0.27(3)] from the data at
p1 = 0.83 and p2 = 0.90, ω = 0.19(5) [ω = 0.31(9)] from those at p1 = 0.883 and p2 = 0.90,
and ω = 0.29(3) [ω = 0.25(4)] from the results at p1 = 0.83 and p2 = 0.883. We also fitted
the difference Ū22− 0.148 at p = 0.90 to cL
−ε (for this value of p next-to-leading corrections
are apparently very small, see Fig. 3). We obtain ω = 0.35(3) [ω = 0.39(6)] for Lmin = 24
[Lmin = 32].
The results of the above-reported fits of Ū22 show that the leading scaling corrections
proportional to L−ω vanish for p ≈ 0.883. Note that, close to p∗, the relevant next-to-
leading scaling corrections should be those proportional to L−ω2 with ω2 ≈ 0.82. Indeed,
according to Eq. (10), the coefficient of those proportional to L−2ω is of order (p − p∗)2,
i.e. b12 ≈ b̄12(p − p
∗). Therefore, b12 is small if b̄12 = O(1) (we checked this numerically).
This applies to the FSS at p = 0.87 and 0.90, where the L−2ω corrections can be neglected,
although in these two cases we cannot neglect the leading L−ω correction whose coefficient is
proportional to p−p∗. An analysis of the leading scaling corrections at p = 0.87, 0.883, 0.90,
assuming the RDIs values Ū∗22 = 0.148(1) and ω = 0.33(3) (we perform combined fits to (15)
10 20 30 40 50
p=0.900, ε=0.82
p=0.900, ε=0.66
p=0.883, ε=0.82
p=0.870, ε=0.82
p=0.870, ε=0.66
p=0.830, ε=0.82
p=0.830, ε=0.66
p=0.800, ε=0.66
FIG. 5: Estimates of Ū∗im as obtained by fits to Ū
im+ cL
−ε, versus the minimum lattice size Lmin
allowed in the fits. Some data are slightly shifted along the x-axis to make them visible. The
dotted lines correspond to the RDIs estimate13 Ū∗im = 1.840(4).
with ε1 = ω) gives the estimate
p∗ = 0.883(3), (17)
which is approximately in the middle of the ferromagnetic line, i.e. 1 − p∗ ≈ (1 − pN)/2.
We performed a similar analysis for Ūd, obtaining a consistent estimate of p
∗. Thus, the ±J
Ising model for p = 0.883 is approximately improved. Therefore, at p = 0.883, fits of the
data assuming O(L−ω2) leading scaling corrections should provide reliable results.
As discussed in Sec. II, a useful quantity to perform stringent checks of universality within
the RDIs universality class is the combination Ūim of quartic cumulants reported in Eq. (11).
For this quantity the scaling corrections proportional to L−ω are small, cf. Eq. (13), and
thus the dominant corrections should behave as L−2ω, with 2ω ≈ 0.66. As already discussed,
for values of p close to p∗, such as p = 0.87, 0.883, 0.90, also the L−2ω term is expected to be
small and thus the dominant corrections should scale as L−ω2 with ω2 ≈ 0.82. In Fig. 4 we
show the MC results for Ūim for various values of p. Fig. 5 shows results of fits to
Ū∗ + cL−ε, (18)
0 20 40 60 80
 L, L
ε=0.82
ε=0.74
ε=0.90
FIG. 6: Estimates of Ū∗d as obtained by fits of Ūd at p = 0.883, to Ū
d + cL
−ε. Some data are
slightly shifted along the x-axis to make them visible. The dotted lines correspond to the RDIs
estimate13 Ū∗d = 1.500(1).
with ε = 0.66, 0.82. We obtain Ū∗im = 1.840(3)[3], 1.842(2)[1], 1.845(2)[3] respectively for
p = 0.90, 0.883, 0.87, fixing ε = ω2 = 0.82(8) (the error in brackets is related to the uncer-
tainty of ω2) and using data with L ≥ 32; moreover we obtain Ū
im = 1.847(3)[2], 1.840(7)[1]
respectively for p = 0.83, 0.80, fixing ε = 0.66(6) and using data with L ≥ 24. For all values
of p the results are in good agreement with RDIs estimate13 Ū∗im = 1.840(4). They provide
strong support to a RDIs critical behavior along the ferromagnetic line.
A further stringent check of universality comes from the analysis of the data Ūd at
p = 0.883, because the data of Ūd are very precise due to a cancellation of the statisti-
cal fluctuations.13 Since the model is improved, the L−kω scaling corrections are negligible
and the large-L behavior is approached with corrections of order L−ω2 , ω2 = 0.82(8). We
thus fit the data to Eq. (18) with ε = 0.74, 0.82, 0.90. In Fig. 6 we show the results. We
obtain Ū∗d = 1.5001(1)[15], 1.5004(2)[9], 1.5003(10)[6] (the error in brackets is related to the
uncertainty of ω2 = 0.82(8)) for Lmin = 12, 24, 48 respectively. Moreover, by fitting the data
to Ū∗d +c1L
−ε1+c2L
−ε2 with ε1 = 0.33 and ε2 = 0.82, we obtain Ū
d = 1.5006(7), 1.500(3) for
Lmin = 12, 24 respectively. These results are in perfect agreement with the RDIs estimate
10 15 20 25 30
0.665
0.670
0.675
0.680
0.685
0.690
R’ξ,im , ε=0.82
R’ξ,im , ε=0.66
 , ε=0.82
R’ξ , ε=0.82
R’ξ , ε=0.33
 , ε=0.82
p=0.883
FIG. 7: Estimates of the critical exponent ν, as obtained by fits of R̄′ξ, Ū
ξ,im, and U
4,im at
p = 0.883. Lmin is the minimum lattice size allowed in the fits. Some data are slightly shifted along
the x-axis to make them visible. The dotted lines correspond to the estimate ν = 0.682(3).
Ū∗d = 1.500(1). Such an agreement is also confirmed by the analysis of the data of Ū22, for
example a fit to Eq. (18) with ε = ω2 = 0.82(8) gives Ū
22 = 0.1486(8)[3] for Lmin = 32, to
be compared with the RDIs estimate13 Ū∗22 = 0.148(1).
We have not shown results for values of p too close to 1, for p > 0.90 say, because they
are affected by crossover effects due to presence of the Ising transition for p = 1, as it also
occurs in randomly dilute Ising models.16,21,32 For instance, for p = 0.94 the data are not
compatible with a behavior of the form (15) with Ū∗22 fixed to the RDIs value. Our data
that correspond to lattice sizes L ≤ 64 apparently converge to a smaller value, consistently
with the expected crossover from pure to random behavior (in pure systems Ū∗22 = 0). The
same quantitative differences are observed in the RSIM and in the RBIM close to the Ising
transition. This suggests that in FSS analyses up to L ≈ 100 the asymptotic RDIs behavior
can only be observed for p . 0.94.
B. Critical exponents
The correlation-length exponent ν can be estimated by fitting the derivative of Rξ and U4
to the expression (8). Accurate estimates are only obtained for improved Hamiltonians. For
generic models, as shown in Ref. 13, good estimates are only obtained by using improved
estimators, such as those reported in Eq. (12).
We analyze the data at p = 0.883, which is a very good approximation of the improved
value p∗ = 0.883(3). In Fig. 7 we report several results for the critical exponent ν, obtained
by analyzing R̄′ξ, Ū
4, and their improved versions R
ξ,im and U
4,im. We show results of fits of
their logarithms to
lnL+ a+ bL−ε, (19)
fixing ε to several values. Since the Hamiltonian is approximately improved, scaling cor-
rections are expected to decrease as L−ω2 with ω2 = 0.82(8). Since p = 0.883 is only
approximately equal to p∗, one may be worried of the residual leading scaling corrections
that are small but do not vanish exactly. Improved estimators should provide the most
reliable results since the leading scaling corrections are additionally suppressed.
As can be seen in Fig. 7, the results obtained by using R̄′ξ and Ū
4 and ε = 0.82 are
perfectly consistent with those obtained from their improved versions. This confirms that
the Hamiltonian is improved. Fits of R̄′ξ to (19) with ε = 0.33 do not provide stable results.
The results approach the values obtained in the other fits only when increasing the minimum
size Lmin allowed in the fit. This is expected, since the L
−ω corrections should be negligible
with respect to the L−ω2 ones. In conclusion, our final estimate of the correlation-length
exponent is
ν = 0.682(3), (20)
which includes all results (with their errors) of the fits of R̄′ξ, Ū
ξ,im, U
4,im to Eq. (19)
with ε = 0.82(8) and Lmin = 16, 24. Estimate (20) is in perfect agreement with the most
precise RDIs estimate ν = 0.683(2).
Estimate (20) is also confirmed by the analysis of the data at the other values of p. Fig. 8
shows results obtained by fitting the logarithm of R′ξ,im to the function (19) for other values
of p. They are definitely consistent with the result obtained at p = 0.883. Results for
p = 0.80 are not shown because the available data are not sufficient to get reliable results.
10 15 20 25 30
p=0.900, R’ξ,im, ε=0.82
p=0.900, R’ξ,im, ε=0.66
p=0.900, U’
, ε=0.82
p=0.870, R’ξ,im, ε=0.82
p=0.870, R’ξ,im, ε=0.66
p=0.870, U’
, ε=0.82
p=0.830, R’ξ,im, ε=0.66
p=0.830, U’
, ε=0.66
FIG. 8: Results from ν obtained by fitting R′ξ,im to aL
1/ν(1+bL−ε). Some data are slightly shifted
along the x-axis to make them visible. The dotted lines correspond to the estimate obtained at
p = 0.883, i.e. ν = 0.682(3).
In order to estimate the critical exponent η, we analyze the FSS of the magnetic suscep-
tibility χ̄, cf. Eq. (9). We fit it to aL2−η + b (where b represents a constant background
term), to aL2−η(1 + cL−ε), and to aL2−η(1 + cL−ε) + b (more precisely, we fit ln χ̄ to the
logarithm of the previous expressions). The results at p = 0.883 are shown in Fig. 9, versus
the minimum size Lmin allowed in the fits. We obtain the estimate
η = 0.036(2), (21)
which includes all results obtained for Lmin & 16. This estimate agrees with the most precise
RDIs estimate η = 0.036(1). Fig. 10 shows results for the other values of p. Again, they are
in good agreement.
C. The critical temperature
The critical temperature can be estimated by extrapolating the estimates of βf at Rξ =
0.5943, cf. Eq. (6). Since we have chosen Rξ = 0.5943 ≈ R
ξ = 0.5944(7),
13 we expect
10 15 20 25 30
+cost
ε=0.82
ε free
+cost & ε=0.82
p=0.883
FIG. 9: Estimates of the critical exponent η, obtained by fitting ln χ̄ at p = 0.883, to a + (2 −
η) lnL+bLη−2 (denoted by +cost), a+(2−η) ln L+a1L
−ε, and to a+(2−η) ln L+a1L
−ε+bLη−2.
Some data are slightly shifted along the x-axis to make them visible. The dotted lines correspond
to the final estimate η = 0.036(2).
in general that βf − βc = O(L
−1/ν−ω). For p = 0.883, since the model is approximately
improved, the leading scaling corrections are related to the next-to-leading exponent ω2.
Thus, in this case βf −βc = O(L
−1/ν−ω2). This behavior is nicely observed in Fig. 11, which
shows βf(L) at p = 0.883 vs L
−1/ν−ω2 with 1/ν + ω2 ≈ 2.28. A fit to βc + aL
−1/ν−ω2 gives
βc = 0.300611(1).
For the other values of p we expect βf − βc = O(L
−1/ν−ω) with 1/ν + ω ≈ 1.79. Linear
fits of βf (L) (for L ≥ Lmin with Lmin sufficiently large to give an acceptable χ
2) give the
estimates βc = 0.25544(2) for p = 0.94, βc = 0.285285(5) for p = 0.90, βc = 0.313748(1) for
p = 0.87, βc = 0.365459(5) for p = 0.83, βc = 0.42501(3) for p = 0.80. We finally recall
that18 βc = 0.22165452(8) for p = 1 (the standard Ising model), and that
7 βc = 0.5967(11)
at the multicritical Nishimori point at pN = 0.7673(3). In Fig. 12 we plot the available
estimates of the critical temperature Tc ≡ 1/βc in the region 1 ≥ p ≥ pN .
The estimates of Tc shown in Fig. 12 hint at a smooth linear behavior for small values of
w ≡ 1− p, close to the Ising point at w = 0. This can be explained by some considerations
10 15 20 25 30 35
p=0.900, +cost
p=0.900, ε free
p=0.870, +cost
p=0.870, ε free
p=0.830, +cost
p=0.830, ε free
FIG. 10: Estimates of the critical exponent η, obtained by fitting χ̄ for various values of p. See the
caption of Fig. 9 for an explanation of the fits. Some data are slightly shifted along the x-axis to
make them visible. The dotted lines correspond to the result η = 0.036(2) obtained at p = 0.883.
on the multicritical behavior around the Ising point at w = 0. The Ising critical behavior
at w = 0 is unstable against the RG perturbation induced by quenched disorder at w > 0,33
which leads to the RDIs critical behavior. Indeed such a perturbation has a positive RG
dimension yw at the Ising fixed point:
10,34 yw = αIs/νIs = 2/νIs − 3 where αIs and νIs are the
Ising specific-heat and correlation-length critical exponents, and therefore17 yw = 0.1740(8).
Thus, in the absence of an external magnetic field, beside the scaling field ut related to the
temperature, there is another relevant scaling field uw associated with the quenched disorder
parameter w ≡ 1− p. General RG scaling arguments21,35 show that the singular part of the
free energy for w → 0 behaves as
Fsing ∼ u
2−αIs
t F (X), X = uwu
t , (22)
where φ = ywνIs = αIs = 0.1096(5) is the crossover exponent, and F (X) is a crossover
scaling function which is universal (apart from normalizations). The scaling fields ut and
uw depend on the parameters of the model. In general, we expect
ut = t+ kw, (23)
0.000 0.002 0.004 0.006
-1/ν−ω2
0.3006
0.3007
FIG. 11: Estimates of βf (L) at p = 0.883 versus L
−(1/ν+ω2) for 1/ν + ω2 ≈ 2.28. The dashed line
corresponds to a linear fit of the data for L ≥ 12.
where t ≡ T/TIs − 1, TIs is the critical temperature of the Ising model, and k is a constant.
No such mixing between t and w occurs in uw, since uw vanishes for w = 0. Hence, we can
take uw = w. The system has a critical transition for w > 0 at Tc(w). Since the singular
part of the free energy close to a critical point behaves as (T −Tc)
2−α (α = −0.049(6) is the
specific-heat exponent of the RDIs universality class), we must have F (Xc) = 0, where Xc
is the value of X obtained by setting T = Tc(w) (see, e.g., Ref. 36 and references therein).
Hence, we obtain
w [Tc(w)/TIs − 1 + kw]
= Xc, (24)
and therefore
Tc(w)/TIs − 1 = (w/Xc)
1/φ − kw + · · · , (25)
where the dots indicate higher-order terms. This expression provides the w dependence of
the critical temperature for w small. Note that the nonanalytic term in Eq. (25) is suppressed
with respect to the analytic ones, because 1/φ ≈ 9.1. Thus, Tc(w) ≈ TIs(1 − kw + O(w
Since Tc(w) < TIs, we can also infer that k > 0. From the results for Tc(w) we estimate
k ≈ 2.2 for the ±J Ising model.
0.0 0.1 0.2
Ising
multicritical point
FIG. 12: The critical temperature Tc ≡ 1/βc vs 1− p.
V. CONCLUSIONS
In this paper we have studied the critical behavior of the 3D ±J Ising model at the
transition line between the paramagnetic and the ferromagnetic phase, which extends from
p = 1 to a multicritical (Nishimori) point at p = pN ≈ 0.767. We presented a FSS analysis of
MC simulations at various values of p in the region pN < p < 1. The results for the critical
exponents and other universal quantities are consistent with those of the RDIs universality
class. For example, we obtained ν = 0.682(3) and η = 0.036(2), which are in good agreement
with the presently most accurate estimates13 ν = 0.683(2) and η = 0.036(1) for the 3D RDIs
universality class. Therefore, our FSS analysis provides a strong evidence that the critical
behavior of the 3D ±J Ising along the ferromagnetic line belongs to the 3D RDIs universality
class.
We also note that the random-exchange interaction in the ±J Ising model gives rise
to frustration, while the RDIs universality class describes transitions in generic diluted
Ising systems with ferromagnetic exchange interactions. This implies that frustration is
irrelevant at the ferromagnetic transition line of the 3D ±J Ising model. Moreover, the
observed scaling corrections are consistent with the RDIs leading and next-to-leading scaling
correction exponents ω = 0.33(3) and ω2 = 0.82(8). This indicates that frustration does not
introduce new irrelevant perturbations at the RDIs fixed point with RG dimension yf & −1.
APPENDIX A: NOTATIONS
We define the two-point correlation function
G(x) ≡ 〈σ0 σx〉, (A1)
where the overline indicates the quenched average over the Jxy probability distribution.
Then, we define the corresponding susceptibility χ ≡
xG(x) and the correlation length ξ
G̃(0)− G̃(qmin)
q̂2minG̃(qmin)
, (A2)
where qmin ≡ (2π/L, 0, 0), q̂ ≡ 2 sin q/2, and G̃(q) is the Fourier transform of G(x). We also
consider quantities that are invariant under RG transformations in the critical limit. Beside
the ratio
Rξ ≡ ξ/L, (A3)
we consider the quartic cumulants U4, U22 and Ud defined by
, (A4)
U22 ≡
µ22 − µ2
Ud ≡ U4 − U22,
where
µk ≡ 〈 (
k〉 . (A5)
We also define corresponding quantities Ū4, Ū22, and Ūd at fixed Rξ = 0.5943. Finally, we
consider the derivative R′ξ of Rξ, and U
4 of U4, with respect to β ≡ 1/T , which allow one to
determine the critical exponent ν.
1 H. Nishimori, Prog. Theor. Phys. 66, 1169 (1981).
2 P. Le Doussal and A.B. Harris, Phys. Rev. Lett. 61, 625 (1988).
3 N. Kawashima and H. Rieger, in Frustrated Spin Systems, edited by H.T. Diep (World Scientific,
Singapore, 2004); cond-mat/0312432.
4 H. Katzgraber, M. Körner, and A.P. Young, Phys. Rev. B 73, 224432 (2006).
5 Y. Ozeki and H. Nishimori, J. Phys. Soc. Japan 56, 3265 (1987).
6 R.R.P. Singh, Phys. Rev. Lett. 67, 899 (1991).
7 Y. Ozeki and N. Ito, J. Phys. A 31, 5451 (1998).
8 Refs. 5,6,7 report the estimates pN = 0.767(2), pN = 0.7656(20), and pN = 0.7673(3) respec-
tively.
9 K. Hukushima, J. Phys. Soc. Japan 69, 631 (2000).
10 A. Pelissetto and E. Vicari, Phys. Rept. 368, 549 (2002).
11 R. Folk, Yu. Holovatch, and T. Yavors’kii, Uspekhi Fiz. Nauk 173, 175 (2003) [Phys. Usp. 46,
175 (2003)].
12 H. G. Ballesteros, L. A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and
J. J. Ruiz-Lorenzo, J. Phys. A 32, 1 (1999).
13 M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp.
P02016 (2007).
14 P. Calabrese, V. Mart́ın-Mayor, A. Pelissetto, and E. Vicari, Phys. Rev. E 68, 036136 (2003).
15 A. Pelissetto and E. Vicari, Phys. Rev. B 62, 6393 (2000).
16 H.G. Ballesteros, L.A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and J.J.
Ruiz-Lorenzo, Phys. Rev. B 58, 2740 (1998).
17 M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E 65, 066127 (2002).
18 Y. Deng and H.W.J. Blöte, Phys. Rev. E 68, 036125 (2003).
19 W. Janke, in Proceedings of the XXIII International Symposium on Lattice Field Theory,
Dublin, July 2005, POS(LAT2005)018
20 N. Ito, Y. Ozeki, and H. Kitatani, J. Phys. Soc. Jpn. 68, 803 (1999).
21 P. Calabrese, P. Parruccini, A. Pelissetto, and E. Vicari, Phys. Rev. E 69, 036120 (2004).
22 M. Hasenbusch, J. Phys. A 32, 4851 (1999).
23 M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B 74, 144506 (2006);
Phys. Rev. B 63, 214503 (2001).
24 See, e.g., S. Wansleben, J.B. Zabolitzky, and C. Kalle, J. Stat. Phys. 37, 271 (1984); G. Bhanot,
D. Duke, and R. Salvador, Phys. Rev. B 33, 7841 (1986).
25 M. Lüscher, Comput. Phys. Commun. 79, 100 (1994).
26 The SIMD-oriented fast Marsenne twister random number generator has been introduced by
M. Matsumoto and M. Saito. Details can be found in M. Saito, Master Thesis (2007) and at
http://www.math.sci.hiroshima-u.ac.jp/∼m-mat/MT/emt.html.
27 In order to make the use of these expensive (in terms of CPU-time) generators affordable,
we employed the same sequence of random numbers for the update of all nbit systems (for
the initialization of the configurations at the beginning of the simulation we used independent
random numbers for each of the systems). This may give rise to a statistical correlation among
the nbit systems. This effect is probably small and we have not detected it. Anyway, in order to
ensure a correct estimate of the statistical error, all nbit systems that use the same sequence of
random numbers have been put in the same bin in our jackknife analysis.
28 R.H. Swendsen and J-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
29 U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
30 D. Ivaneyko, J. Ilnytskyi, B. Berche, and Yu. Holovatch, Physica A 370, 163 (2006).
31 H.G. Ballesteros, L.A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and J.J.
Ruiz-Lorenzo, Nucl. Phys. B 512, 681 (1998).
32 The crossover exponent from pure Ising to RDIs critical behavior is the Ising specific-heat
exponent34 αIs, see also Sec. IVC. This implies that the crossover scaling variable in the FSS at
Tc is given by the combination X = cwL
αIs/νIs , where w = 1 − p, αIs/νIs = 0.1740(8), and c is
a normalization constant. When w → 0, strong crossover effects are expected for X . 1, which
corresponds to L . (cw)−5.75. The RDIs asymptotic critical behavior is observed for X ≫ 1.
33 A. B. Harris, J. Phys. C 7, 1671 (1974).
34 A. Aharony, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. Lebowitz
(Academic Press, New York, 1976), Vol. 6, p. 357.
35 I.D. Lawrie and S. Sarbach, in Phase Transitions and Critical Phenomena, Vol. 9, edited by C.
Domb and J. Lebowitz (Academic, London, 1984).
36 A. Pelissetto and E. Vicari, cond-mat/0702273.

  We study the critical behavior of the three-dimensional $\pm J$ Ising model
[with a random-exchange probability $P(J_{xy}) = p \delta(J_{xy} - J) + (1-p)
\delta(J_{xy} + J)$] at the transition line between the paramagnetic and
ferromagnetic phase, which extends from $p=1$ to a multicritical (Nishimori)
point at $p=p_N\approx 0.767$. By a finite-size scaling analysis of Monte Carlo
simulations at various values of $p$ in the region $p_N<p<1$, we provide strong
numerical evidence that the critical behavior along the ferromagnetic
transition line belongs to the same universality class as the three-dimensional
randomly-dilute Ising model. We obtain the results $\nu=0.682(3)$ and
$\eta=0.036(2)$ for the critical exponents, which are consistent with the
estimates $\nu=0.683(2)$ and $\eta=0.036(1)$ at the transition of
randomly-dilute Ising models.

arXiv:0704.0427v1  [cond-mat.dis-nn]  3 Apr 2007
The 3D ±J Ising model at the ferromagnetic transition line
Martin Hasenbusch,1 Francesco Parisen Toldin,2 Andrea Pelissetto,3 and Ettore Vicari1
1 Dipartimento di Fisica dell’Università di Pisa and INFN, Pisa, Italy.
2 Scuola Normale Superiore and INFN, Pisa, Italy.
3Dipartimento di Fisica dell’Università di Roma “La Sapienza” and INFN, Roma, Italy.
(Dated: November 1, 2018)
Abstract
We study the critical behavior of the three-dimensional ±J Ising model [with a random-exchange
probability P (Jxy) = pδ(Jxy−J)+(1−p)δ(Jxy+J)] at the transition line between the paramagnetic
and ferromagnetic phase, which extends from p = 1 to a multicritical (Nishimori) point at p =
pN ≈ 0.767. By a finite-size scaling analysis of Monte Carlo simulations at various values of p
in the region pN < p < 1, we provide strong numerical evidence that the critical behavior along
the ferromagnetic transition line belongs to the same universality class as the three-dimensional
randomly-dilute Ising model. We obtain the results ν = 0.682(3) and η = 0.036(2) for the critical
exponents, which are consistent with the estimates ν = 0.683(2) and η = 0.036(1) at the transition
of randomly-dilute Ising models.
PACS numbers: 75.10.Nr, 75.40.Cx, 75.40.Mg, 64.60.Fr
http://arxiv.org/abs/0704.0427v1
I. INTRODUCTION
The ±J Ising model has played an important role in the study of the effects of quenched
random disorder and frustration on Ising systems. It is defined by the lattice Hamiltonian
H±J = −
Jxyσxσy, (1)
where σx = ±1, the sum is over the nearest-neighbor sites of a simple cubic lattice, and the
exchange interactions Jxy are uncorrelated quenched random variables, taking values ±J
with probability distribution
P (Jxy) = pδ(Jxy − J) + (1− p)δ(Jxy + J). (2)
For p = 1 we recover the standard Ising model, while for p = 1/2 we obtain the usual
bimodal Ising spin-glass model.
The phase diagram of the three-dimensional (3D) ±J Ising model is sketched in Fig. 1.
The high-temperature phase is paramagnetic for any p. The low-temperature phase depends
on the value of p: it is ferromagnetic for small values of 1 − p, while it is a spin-glass
phase with vanishing magnetization for sufficiently large values of 1 − p. The different
phases are separated by transition lines, which meet at a multicritical point N located
along the so-called Nishimori line.1,2,3 The spin-glass transition has been mostly studied at
the symmetric point p = 1/2, see, e.g., Refs. 3,4 and references therein. The spin-glass
transition line extends up to the Nishimori multicritical point,2 located at5,6,7,8 pN ≈ 0.767.
For larger values of p, the transition is ferromagnetic, up to p = 1 where one recovers the pure
Ising model, and therefore a transition in the Ising universality class. At the ferromagnetic
transition line, for pN < p < 1, the critical behavior is expected to belong to a different
universality class.
An interesting hypothesis, which has already been put forward in Refs. 3,9, is that the
ferromagnetic transition of the ±J Ising model belongs to the 3D randomly-dilute Ising
(RDIs) universality class (see, e.g., Refs. 10,11 for reviews on randomly-dilute spin mod-
els). A representative of the RDIs universality class is the randomly site-dilute Ising model
(RSIM) defined by the lattice Hamiltonian
Hd = −J
ρx ρy σxσy, (3)
glass
ferro
1 − p
FIG. 1: Sketch of the phase diagram of the 3D ±J Ising model in the T -p plane.
where ρx are uncorrelated quenched random variables, which are equal to 0, 1 with proba-
bility
P (ρx) = pδ(ρx − 1) + (1− p)δ(ρx). (4)
For p < 1 and above the percolation threshold of the spins (pperc ≈ 0.3116081(13) on a cubic
lattice12), the RSIM undergoes a continuous phase transition between a disordered and a
ferromagnetic phase, whose nature is independent of p. This transition is definitely different
from the usual Ising transition: for instance, the correlation-length critical exponent13,14,15,16
ν = 0.683(2) differs from the Ising value17,18 ν = 0.63012(16). The RDIs universality class is
expected to describe the ferromagnetic transition in generic diluted ferromagnetic systems.
For instance, it has been verified that also the randomly bond-diluted Ising model (RBIM)
belongs to the RDIs universality class.13,19 These results do not necessarily imply that also
the ±J Ising model has an RDIs ferromagnetic transition line. Indeed, while the RSIM (3)
has only ferromagnetic exchange interactions, the ±J Ising model is frustrated for any value
of p < 1. Therefore, the ferromagnetic transition in the ±J Ising model belongs to the RDIs
universality class only if frustration is irrelevant, a fact that is not obvious and should be
carefully investigated.
Reference 9 investigated the issue by means of a Monte Carlo (MC) renormalization-
group (RG) study, claiming that the ±J Ising model belongs to the same RDIs universality
class as the RSIM and the RBIM. It should be noted however that the quoted estimate for
the correlation-length exponent at the ferromagnetic transition, ν = 0.658(9), is close to
but not fully consistent with the RDIs value ν = 0.683(2).13 Another numerical MC work20
investigated the nonequilibrium relaxation dynamics of the ±J Ising model and showed an
apparent nonuniversal dynamical critical behavior along the ferromagnetic transition line.
These results are not conclusive and further investigation is called for to clarify this issue.
In this paper we focus on the transition line of the 3D ±J Ising model between the
paramagnetic and the ferromagnetic phase. We investigate the critical behavior by means
of MC simulations at various values of p in the region pN < p < 1. Our finite-size scaling
(FSS) analysis provides a strong evidence that the critical behavior of the 3D ±J Ising
along the ferromagnetic line belongs to the 3D RDIs universality class. For example, we
obtain ν = 0.682(3) and η = 0.036(2), which are in good agreement with the presently most
accurate estimates13 ν = 0.683(2) and η = 0.036(1) for the 3D RDIs universality class.
The paper is organized as follows. In Sec. II we summarize some FSS results which are
needed for the analysis of the MC data, and describe our strategy to check whether the
transition belongs to the RDIs universality class. In Sec. III we describe the MC simula-
tions. In Sec. IV we report the results of the FSS analysis. Finally, in Sec. V we draw our
conclusions. In App. A we report the definitions of the quantities we compute.
II. STRATEGY OF THE FINITE-SIZE SCALING ANALYSIS
In this work we check whether the ferromagnetic transition line in the 3D ±J Ising models
belongs to the RDIs universality class. For this purpose, we present a FSS analysis of MC
data for various values of p in the region 1 > p > pN ≈ 0.767. We follow closely Ref. 13,
which studied the ferromagnetic transition line in the 3D RSIM and RBIM and provided
strong numerical evidence that these transitions belong to the same RDIs universality class.
We refer to Ref. 13 for notations (a short summary is reported in App. A) and a detailed
discussion of FSS in these disordered systems.
According to the RG, in the case of periodic boundary conditions and for L → ∞, where
L is the lattice size, a generic RG invariant quantity R at the critical temperature 1/βc
behaves as
R(L, β = βc) = R
1 + c11L
−ω + c12L
−2ω + · · ·+ c21L
−ω2 + · · ·
, (5)
where R∗ is the universal infinite-volume limit and ω and ω2 are the leading and next-
to-leading correction-to-scaling exponents. In RDIs systems scaling corrections play an
important role,16,21 since ω is quite small. Indeed we have ω = 0.33(3) and ω2 = 0.82(8)
in the 3D RDIs universality class.13 These slowly-decaying scaling corrections make the
accurate determination of the universal asymptotic behavior quite difficult.
Instead of computing the various quantities at fixed Hamiltonian parameters, we keep
a RG invariant quantity R fixed at a given value Rf .
22 This means that, for each L, we
determine the pseudocritical inverse temperature βf(L) such that
R(β = βf(L), L) = Rf . (6)
All interesting thermodynamic quantities are then computed at β = βf (L). The pseudocrit-
ical inverse temperature βf(L) converges to βc as L → ∞. The value Rf can be specified
at will, as long as Rf is taken between the high- and low-temperature fixed-point values
of R. The choice Rf = R
∗ (where R∗ is the critical-point value) improves the conver-
gence of βf to βc for L → ∞; indeed βf − βc = O(L
−1/ν) for generic values of Rf , while
βf −βc = O(L
−1/ν−ω) for Rf = R
∗. This FSS method has already been applied to the study
of the critical behavior of N -vector spin models,22,23 and of randomly-dilute Ising models.13
As in Ref. 13, we perform a FSS analysis at fixed Rξ ≡ ξ/L = 0.5943, which is very close to
the fixed-point value R∗ξ = 0.5944(7) of Rξ at βc. Given any RG invariant quantity R, such as
the quartic cumulants U4 and U22, we consider its value at fixed Rξ, i.e., R̄(L) = R(L, βf(L)).
For L → ∞, R̄(L) behaves as R(L, βc):
R̄(L) = R̄∗
1 + b11L
−ω + b12L
−2ω + · · ·+ b21L
−ω2 + · · ·
, (7)
where the coefficients bij depend on the Hamiltonian. The derivative R̄
′ with respect to β
of a generic RG invariant quantity R behaves as
R̄′(L) = aL1/ν
1 + a11L
−ω + a12L
−2ω + · · ·+ a21L
−ω2 + · · ·
. (8)
Finally, the FSS of the magnetic susceptibility χ is given by13
χ̄(L) ≡ χ(L, β = βf (L)) = eL
1 + e11L
−ω + e12L
−2ω + · · ·+ e21L
−ω2 + · · ·
+ eb (9)
where eb represents the background contribution.
A standard RG analysis, see, e.g., Ref. 13, shows that the amplitudes of the O(L−kω)
scaling corrections are proportional to uk3 (with a universal coefficient), where u3 is the
leading irrelevant scaling field with RG dimension y3 = −ω. Hamiltonians such that u3 = 0—
we call them improved Hamiltonians—have a faster approach to the universal asymptotic
behavior, because the O(L−kω) scaling corrections vanish: b1k = a1k = e1k = 0 in Eqs. (7),
(8), and (9). In this case the leading scaling corrections are proportional to u4L
−ω2 , where u4
is the next-to-leading irrelevant scaling field and y4 = −ω2 is its RG dimension. In Ref. 13
it was shown that the RSIM for p = p∗ = 0.800(5) and the RBIM for p = p∗ = 0.56(2)
are improved. Since scaling fields are analytic functions of the Hamiltonian parameters, u3
must be proportional to p − p∗ close to p = p∗, i.e. u3 ≈ c3(p − p
∗). Therefore, since the
coefficients b1k, a1k, and e1k that appear in Eqs. (7), (8), and (9) are proportional to u
3, we
b1k, a1k, e1k ∼ (p− p
∗)k. (10)
Beside the quantities defined in App. A, we also consider observables—in analogy with
the previous terminology, we call them improved quantities—characterized by the fact that
the leading scaling correction proportional to L−ω (approximately) vanishes in any model
belonging to the RDIs universality class.13 We consider the combination of quartic cumulants
Ūim = Ū4 + 1.3Ū22, (11)
and improved estimators of the critical exponent ν defined as
R′ξ,im ≡ R̄
d , U
4,im ≡ Ū
d (12)
(Ūd is defined in App. A). In Ref. 13 we showed that, if the transition belongs to the RDIs
universality class, the leading scaling correction proportional to L−ω of these improved ob-
servables is suppressed. More precisely, we showed that the universal ratio of the amplitudes
of the leading scaling correction in Ūim and Ū4 satisfies
|b11,Ūim/b11,Ū4 | .
, (13)
while the one for the quantities R′ξ,im and R̄
ξ is bounded by
|a11,R′
/a11,R̄′
. (14)
The remaining scaling corrections are of order L−2ω and L−ω2. These improved observables
are particular useful to check whether the transition in a given system belongs to the 3D
RDIs universality class.
10 100
FIG. 2: Exponential autocorrelation time τ of the magnetic susceptibility for a mixture of
Metropolis and cluster updates as discussed in the text, at p = 0.87. The dotted line shows
the result of a fit to τ = cLz: this fit gives z ≈ 1.6.
To summarize: in order to check whether the ferromagnetic transition of the 3D ±J
Ising model belongs to the RDIs universality class, we perform a FSS analysis at fixed
Rξ = 0.5943, and check if the results for the critical exponents and other universal quantities
are consistent with those obtained for the RDIs universality class, which is characterized
by13 critical exponents ν = 0.683(2) and η = 0.036(1), by the leading and next-to-leading
scaling-correction exponents ω = 0.33(3) and ω2 = 0.82(8) and by the universal infinite-
volume values of the quartic cumulants Ū∗22 = 0.148(1), Ū
im = 1.840(4), and Ū
d = 1.500(1).
Notice that the fact that we fix Rξ = 0.5943 does not introduce any bias in our FSS analysis.
III. MONTE CARLO SIMULATIONS
We performed MC simulations of Hamiltonian (1) with J = 1 for p =
0.94, 0.90, 0.883, 0.87, 0.83, 0.80, close to the critical temperature on cubic lattices of size
L3 with periodic boundary conditions, for a large range of lattice sizes: from L = 8 to
L = 80 for p = 0.883, 0.87, to L = 64 for p = 0.94, 0.90, 0.83, and to L = 48 for p = 0.80. We
chose values of p not too close to p = 1: indeed, as p → 1 we expect crossover effects due to
the presence of the Ising transition for p = 1 and, therefore, that the asymptotic behavior
sets in only for large values of L. We return to this point later.
We used a Metropolis algorithm and multispin coding.24 In the simulation nbit systems
evolve in parallel, where nbit = 32 or nbit = 64 depending on the computer that is used.
For each of these nbit systems we use a different set of couplings Jxy. This allows us to
perform 64 parallel simulations on a 64-bit machine, and therefore to gain a large factor
in the efficiency of the MC simulations. We used high-quality random-number generators,
such as the RANLUX25 or the twister26 generators.27 Using the twister random-number
generator, we need about 1.2 × 10−9 seconds for one Metropolis update of a single spin on
an Opteron processor running at 2 GHz. Our simulations took approximately 3 CPU years
on an Opteron (2 GHz) processor.
It is worth mentioning that cluster algorithms, such as the Swendsen-Wang cluster28
and the Wolff single-cluster29 algorithm, show significant slowing down in the ±J Ising
model. At the earlier stage of this work we performed some simulations of the ±J Ising
model at p = 0.87 using the algorithm used in Ref. 13 to simulate the RSIM and the
RBIM. There we used a combination of Metropolis, Swendsen-Wang cluster,28 and Wolff
single-cluster29 updates. More precisely, each updating step consisted of 1 Swendsen-Wang
update, 1 Metropolis update, and L single-cluster updates. In all cases the exponential
autocorrelation times τ of the magnetic susceptibility was small: τ . 1 in units of the
above updating step, even for the largest lattice sizes considered, i.e. L = 192. In the ±J
Ising model at p = 0.87 autocorrelation times are much larger. In Fig. 2 we plot estimates
of τ as obtained from the magnetic susceptibility. They show a clear evidence of critical
slowing down: τ ∼ Lz with z ≈ 1.6. Such a value of z should be compared with the
dynamic exponent of Swendsen-Wang and Wolff cluster algorithms in the RSIM, which is
much smaller:30 z . 0.5. These results show that cluster algorithms behave differently in
the ±J Ising model, likely due to frustration. They suggest that frustration is relevant for
the cluster dynamics.
Taking also into account the computer time required by the cluster algorithms, we then
turned to a multispin Metropolis algorithm. This turns out to be much more effective at the
lattice sizes considered, although it has a larger dynamic exponent z & 2, see, e.g., Ref. 30
and references therein. We also mention that the autocorrelation time significantly increases
with decreasing p (keeping L fixed). For example, for L = 48 it increases by approximately
a factor of 10 from p = 0.90 to p = 0.80. This represents a major limitation to perform
simulations for large lattices close to the multicritical point.
For each lattice size we considered Ns disorder samples, with Ns decreasing with in-
creasing L, from Ns & 10
6 for L = 8 to Ns & 2 × 10
4 for the largest lattices. For each
disorder sample, we collected a few hundred independent measurements at equilibrium. The
averages over disorder are affected by a bias due to the finite number of measures at fixed
disorder.13,31 A bias correction is required whenever one considers the disorder average of
combinations of thermal averages. We used the formulas reported in App. B of Ref. 13.
Errors were computed from the sample-to-sample fluctuations and were determined by using
the jackknife method.27
Our FSS analysis is performed at fixed Rξ ≡ ξ/L. In order to determine expectation
values at fixed Rξ, one needs the values of the observables as a function of β in some
neighborhood of the inverse temperature βrun used in the simulation. In Ref. 13 we used the
reweighting method for this purpose. This requires that the observables and, in particular,
the values of the energy are stored at each measurement. For the huge statistics like those we
have for the smaller values of L, this becomes unpractical. Therefore, we used here a second-
order Taylor expansion, determining O(β, L) from O(βrun, L)+aO(β−βrun)+ bO(β−βrun)
The coefficients aO and bO are obtained from appropriate expectation values as in Ref. 23.
Since their computation involves disorder averages of products of thermal averages, we have
implemented in all cases an exact bias correction, using the formulas of Ref. 13. Derivatives
with respect to β are then obtained as O′(β, L) = aO + 2bO(β − βrun). Of course, this
method requires |βrun − βf | to be sufficiently small. We have carefully checked the results
by performing, for each L and p, runs at different values of β.
The MC estimates of the quantities introduced in Sec. II and in App. A at fixed Rξ ≡
ξ/L = 0.5943 are available on request.
IV. FINITE-SIZE SCALING ANALYSIS
In this section we present the results of our FSS analysis of the MC data at fixed Rξ =
0.5943.
0.0 0.1 0.2 0.3 0.4 0.5
-0.33
p=0.90
p=0.883
p=0.87
p=0.83
p=0.80
FIG. 3: MC estimates of Ū22 versus L
−ω with ω = 0.33 for different values of p. The dotted lines
show results of fits to c0 + c1L
−ε1 + c2L
−ε2 , fixing c0 = 0.148, ε1 = 0.33, and ε2 = 0.82. In the
RDIs universality class Ū∗22 = 0.148(1).
A. Renormalization-group invariant quantities
In Fig. 3 we show the MC estimates of Ū22 versus L
−ω with ω = 0.33(3), which is the
leading scaling exponent of the RDIs universality class. The data vary significantly with p
and L. This p and L dependence is always consistent with the existence of the expected
next-to-leading scaling corrections, i.e. with a behavior of the form
Ū22 = Ū
22 + c1L
−ε1 + c2L
−ε2 , (15)
where Ū∗22, ε1 and ε2 are fixed to the RDIs values:
13 Ū∗22 = 0.148, ε1 = 0.33 and ε2 =
0.66, 0.82. The fits corresponding to ε2 = 0.82 are shown in Fig. 3. Note that in most of
the cases it is crucial to include a next-to-leading correction. Only for p = 0.90 the data are
well fitted by taking only the leading scaling correction.
An unbiased estimate of ω can be obtained from the difference of data at different values
of p, i.e. by considering
Ū22(p1;L)− Ū22(p2;L) ≈ cL
−ω. (16)
Linear fits of the logarithm of these differences give results in reasonable agreement with
0.00 0.05 0.10 0.15 0.20
-0.82
p=0.90
p=0.883
p=0.87
p=0.83
p=0.80
FIG. 4: MC estimates of Ūim versus L
−0.82. The filled square on the vertical axis corresponds to
the RDIs estimate13 Ū∗im = 1.840(4).
the RDIs estimate ω = 0.33(3), especially when only data corresponding to L ≥ Lmin = 24
are used. For Lmin = 24 [Lmin = 32], we obtain ω = 0.27(2) [ω = 0.27(3)] from the data at
p1 = 0.83 and p2 = 0.90, ω = 0.19(5) [ω = 0.31(9)] from those at p1 = 0.883 and p2 = 0.90,
and ω = 0.29(3) [ω = 0.25(4)] from the results at p1 = 0.83 and p2 = 0.883. We also fitted
the difference Ū22− 0.148 at p = 0.90 to cL
−ε (for this value of p next-to-leading corrections
are apparently very small, see Fig. 3). We obtain ω = 0.35(3) [ω = 0.39(6)] for Lmin = 24
[Lmin = 32].
The results of the above-reported fits of Ū22 show that the leading scaling corrections
proportional to L−ω vanish for p ≈ 0.883. Note that, close to p∗, the relevant next-to-
leading scaling corrections should be those proportional to L−ω2 with ω2 ≈ 0.82. Indeed,
according to Eq. (10), the coefficient of those proportional to L−2ω is of order (p − p∗)2,
i.e. b12 ≈ b̄12(p − p
∗). Therefore, b12 is small if b̄12 = O(1) (we checked this numerically).
This applies to the FSS at p = 0.87 and 0.90, where the L−2ω corrections can be neglected,
although in these two cases we cannot neglect the leading L−ω correction whose coefficient is
proportional to p−p∗. An analysis of the leading scaling corrections at p = 0.87, 0.883, 0.90,
assuming the RDIs values Ū∗22 = 0.148(1) and ω = 0.33(3) (we perform combined fits to (15)
10 20 30 40 50
p=0.900, ε=0.82
p=0.900, ε=0.66
p=0.883, ε=0.82
p=0.870, ε=0.82
p=0.870, ε=0.66
p=0.830, ε=0.82
p=0.830, ε=0.66
p=0.800, ε=0.66
FIG. 5: Estimates of Ū∗im as obtained by fits to Ū
im+ cL
−ε, versus the minimum lattice size Lmin
allowed in the fits. Some data are slightly shifted along the x-axis to make them visible. The
dotted lines correspond to the RDIs estimate13 Ū∗im = 1.840(4).
with ε1 = ω) gives the estimate
p∗ = 0.883(3), (17)
which is approximately in the middle of the ferromagnetic line, i.e. 1 − p∗ ≈ (1 − pN)/2.
We performed a similar analysis for Ūd, obtaining a consistent estimate of p
∗. Thus, the ±J
Ising model for p = 0.883 is approximately improved. Therefore, at p = 0.883, fits of the
data assuming O(L−ω2) leading scaling corrections should provide reliable results.
As discussed in Sec. II, a useful quantity to perform stringent checks of universality within
the RDIs universality class is the combination Ūim of quartic cumulants reported in Eq. (11).
For this quantity the scaling corrections proportional to L−ω are small, cf. Eq. (13), and
thus the dominant corrections should behave as L−2ω, with 2ω ≈ 0.66. As already discussed,
for values of p close to p∗, such as p = 0.87, 0.883, 0.90, also the L−2ω term is expected to be
small and thus the dominant corrections should scale as L−ω2 with ω2 ≈ 0.82. In Fig. 4 we
show the MC results for Ūim for various values of p. Fig. 5 shows results of fits to
Ū∗ + cL−ε, (18)
0 20 40 60 80
 L, L
ε=0.82
ε=0.74
ε=0.90
FIG. 6: Estimates of Ū∗d as obtained by fits of Ūd at p = 0.883, to Ū
d + cL
−ε. Some data are
slightly shifted along the x-axis to make them visible. The dotted lines correspond to the RDIs
estimate13 Ū∗d = 1.500(1).
with ε = 0.66, 0.82. We obtain Ū∗im = 1.840(3)[3], 1.842(2)[1], 1.845(2)[3] respectively for
p = 0.90, 0.883, 0.87, fixing ε = ω2 = 0.82(8) (the error in brackets is related to the uncer-
tainty of ω2) and using data with L ≥ 32; moreover we obtain Ū
im = 1.847(3)[2], 1.840(7)[1]
respectively for p = 0.83, 0.80, fixing ε = 0.66(6) and using data with L ≥ 24. For all values
of p the results are in good agreement with RDIs estimate13 Ū∗im = 1.840(4). They provide
strong support to a RDIs critical behavior along the ferromagnetic line.
A further stringent check of universality comes from the analysis of the data Ūd at
p = 0.883, because the data of Ūd are very precise due to a cancellation of the statisti-
cal fluctuations.13 Since the model is improved, the L−kω scaling corrections are negligible
and the large-L behavior is approached with corrections of order L−ω2 , ω2 = 0.82(8). We
thus fit the data to Eq. (18) with ε = 0.74, 0.82, 0.90. In Fig. 6 we show the results. We
obtain Ū∗d = 1.5001(1)[15], 1.5004(2)[9], 1.5003(10)[6] (the error in brackets is related to the
uncertainty of ω2 = 0.82(8)) for Lmin = 12, 24, 48 respectively. Moreover, by fitting the data
to Ū∗d +c1L
−ε1+c2L
−ε2 with ε1 = 0.33 and ε2 = 0.82, we obtain Ū
d = 1.5006(7), 1.500(3) for
Lmin = 12, 24 respectively. These results are in perfect agreement with the RDIs estimate
10 15 20 25 30
0.665
0.670
0.675
0.680
0.685
0.690
R’ξ,im , ε=0.82
R’ξ,im , ε=0.66
 , ε=0.82
R’ξ , ε=0.82
R’ξ , ε=0.33
 , ε=0.82
p=0.883
FIG. 7: Estimates of the critical exponent ν, as obtained by fits of R̄′ξ, Ū
ξ,im, and U
4,im at
p = 0.883. Lmin is the minimum lattice size allowed in the fits. Some data are slightly shifted along
the x-axis to make them visible. The dotted lines correspond to the estimate ν = 0.682(3).
Ū∗d = 1.500(1). Such an agreement is also confirmed by the analysis of the data of Ū22, for
example a fit to Eq. (18) with ε = ω2 = 0.82(8) gives Ū
22 = 0.1486(8)[3] for Lmin = 32, to
be compared with the RDIs estimate13 Ū∗22 = 0.148(1).
We have not shown results for values of p too close to 1, for p > 0.90 say, because they
are affected by crossover effects due to presence of the Ising transition for p = 1, as it also
occurs in randomly dilute Ising models.16,21,32 For instance, for p = 0.94 the data are not
compatible with a behavior of the form (15) with Ū∗22 fixed to the RDIs value. Our data
that correspond to lattice sizes L ≤ 64 apparently converge to a smaller value, consistently
with the expected crossover from pure to random behavior (in pure systems Ū∗22 = 0). The
same quantitative differences are observed in the RSIM and in the RBIM close to the Ising
transition. This suggests that in FSS analyses up to L ≈ 100 the asymptotic RDIs behavior
can only be observed for p . 0.94.
B. Critical exponents
The correlation-length exponent ν can be estimated by fitting the derivative of Rξ and U4
to the expression (8). Accurate estimates are only obtained for improved Hamiltonians. For
generic models, as shown in Ref. 13, good estimates are only obtained by using improved
estimators, such as those reported in Eq. (12).
We analyze the data at p = 0.883, which is a very good approximation of the improved
value p∗ = 0.883(3). In Fig. 7 we report several results for the critical exponent ν, obtained
by analyzing R̄′ξ, Ū
4, and their improved versions R
ξ,im and U
4,im. We show results of fits of
their logarithms to
lnL+ a+ bL−ε, (19)
fixing ε to several values. Since the Hamiltonian is approximately improved, scaling cor-
rections are expected to decrease as L−ω2 with ω2 = 0.82(8). Since p = 0.883 is only
approximately equal to p∗, one may be worried of the residual leading scaling corrections
that are small but do not vanish exactly. Improved estimators should provide the most
reliable results since the leading scaling corrections are additionally suppressed.
As can be seen in Fig. 7, the results obtained by using R̄′ξ and Ū
4 and ε = 0.82 are
perfectly consistent with those obtained from their improved versions. This confirms that
the Hamiltonian is improved. Fits of R̄′ξ to (19) with ε = 0.33 do not provide stable results.
The results approach the values obtained in the other fits only when increasing the minimum
size Lmin allowed in the fit. This is expected, since the L
−ω corrections should be negligible
with respect to the L−ω2 ones. In conclusion, our final estimate of the correlation-length
exponent is
ν = 0.682(3), (20)
which includes all results (with their errors) of the fits of R̄′ξ, Ū
ξ,im, U
4,im to Eq. (19)
with ε = 0.82(8) and Lmin = 16, 24. Estimate (20) is in perfect agreement with the most
precise RDIs estimate ν = 0.683(2).
Estimate (20) is also confirmed by the analysis of the data at the other values of p. Fig. 8
shows results obtained by fitting the logarithm of R′ξ,im to the function (19) for other values
of p. They are definitely consistent with the result obtained at p = 0.883. Results for
p = 0.80 are not shown because the available data are not sufficient to get reliable results.
10 15 20 25 30
p=0.900, R’ξ,im, ε=0.82
p=0.900, R’ξ,im, ε=0.66
p=0.900, U’
, ε=0.82
p=0.870, R’ξ,im, ε=0.82
p=0.870, R’ξ,im, ε=0.66
p=0.870, U’
, ε=0.82
p=0.830, R’ξ,im, ε=0.66
p=0.830, U’
, ε=0.66
FIG. 8: Results from ν obtained by fitting R′ξ,im to aL
1/ν(1+bL−ε). Some data are slightly shifted
along the x-axis to make them visible. The dotted lines correspond to the estimate obtained at
p = 0.883, i.e. ν = 0.682(3).
In order to estimate the critical exponent η, we analyze the FSS of the magnetic suscep-
tibility χ̄, cf. Eq. (9). We fit it to aL2−η + b (where b represents a constant background
term), to aL2−η(1 + cL−ε), and to aL2−η(1 + cL−ε) + b (more precisely, we fit ln χ̄ to the
logarithm of the previous expressions). The results at p = 0.883 are shown in Fig. 9, versus
the minimum size Lmin allowed in the fits. We obtain the estimate
η = 0.036(2), (21)
which includes all results obtained for Lmin & 16. This estimate agrees with the most precise
RDIs estimate η = 0.036(1). Fig. 10 shows results for the other values of p. Again, they are
in good agreement.
C. The critical temperature
The critical temperature can be estimated by extrapolating the estimates of βf at Rξ =
0.5943, cf. Eq. (6). Since we have chosen Rξ = 0.5943 ≈ R
ξ = 0.5944(7),
13 we expect
10 15 20 25 30
+cost
ε=0.82
ε free
+cost & ε=0.82
p=0.883
FIG. 9: Estimates of the critical exponent η, obtained by fitting ln χ̄ at p = 0.883, to a + (2 −
η) lnL+bLη−2 (denoted by +cost), a+(2−η) ln L+a1L
−ε, and to a+(2−η) ln L+a1L
−ε+bLη−2.
Some data are slightly shifted along the x-axis to make them visible. The dotted lines correspond
to the final estimate η = 0.036(2).
in general that βf − βc = O(L
−1/ν−ω). For p = 0.883, since the model is approximately
improved, the leading scaling corrections are related to the next-to-leading exponent ω2.
Thus, in this case βf −βc = O(L
−1/ν−ω2). This behavior is nicely observed in Fig. 11, which
shows βf(L) at p = 0.883 vs L
−1/ν−ω2 with 1/ν + ω2 ≈ 2.28. A fit to βc + aL
−1/ν−ω2 gives
βc = 0.300611(1).
For the other values of p we expect βf − βc = O(L
−1/ν−ω) with 1/ν + ω ≈ 1.79. Linear
fits of βf (L) (for L ≥ Lmin with Lmin sufficiently large to give an acceptable χ
2) give the
estimates βc = 0.25544(2) for p = 0.94, βc = 0.285285(5) for p = 0.90, βc = 0.313748(1) for
p = 0.87, βc = 0.365459(5) for p = 0.83, βc = 0.42501(3) for p = 0.80. We finally recall
that18 βc = 0.22165452(8) for p = 1 (the standard Ising model), and that
7 βc = 0.5967(11)
at the multicritical Nishimori point at pN = 0.7673(3). In Fig. 12 we plot the available
estimates of the critical temperature Tc ≡ 1/βc in the region 1 ≥ p ≥ pN .
The estimates of Tc shown in Fig. 12 hint at a smooth linear behavior for small values of
w ≡ 1− p, close to the Ising point at w = 0. This can be explained by some considerations
10 15 20 25 30 35
p=0.900, +cost
p=0.900, ε free
p=0.870, +cost
p=0.870, ε free
p=0.830, +cost
p=0.830, ε free
FIG. 10: Estimates of the critical exponent η, obtained by fitting χ̄ for various values of p. See the
caption of Fig. 9 for an explanation of the fits. Some data are slightly shifted along the x-axis to
make them visible. The dotted lines correspond to the result η = 0.036(2) obtained at p = 0.883.
on the multicritical behavior around the Ising point at w = 0. The Ising critical behavior
at w = 0 is unstable against the RG perturbation induced by quenched disorder at w > 0,33
which leads to the RDIs critical behavior. Indeed such a perturbation has a positive RG
dimension yw at the Ising fixed point:
10,34 yw = αIs/νIs = 2/νIs − 3 where αIs and νIs are the
Ising specific-heat and correlation-length critical exponents, and therefore17 yw = 0.1740(8).
Thus, in the absence of an external magnetic field, beside the scaling field ut related to the
temperature, there is another relevant scaling field uw associated with the quenched disorder
parameter w ≡ 1− p. General RG scaling arguments21,35 show that the singular part of the
free energy for w → 0 behaves as
Fsing ∼ u
2−αIs
t F (X), X = uwu
t , (22)
where φ = ywνIs = αIs = 0.1096(5) is the crossover exponent, and F (X) is a crossover
scaling function which is universal (apart from normalizations). The scaling fields ut and
uw depend on the parameters of the model. In general, we expect
ut = t+ kw, (23)
0.000 0.002 0.004 0.006
-1/ν−ω2
0.3006
0.3007
FIG. 11: Estimates of βf (L) at p = 0.883 versus L
−(1/ν+ω2) for 1/ν + ω2 ≈ 2.28. The dashed line
corresponds to a linear fit of the data for L ≥ 12.
where t ≡ T/TIs − 1, TIs is the critical temperature of the Ising model, and k is a constant.
No such mixing between t and w occurs in uw, since uw vanishes for w = 0. Hence, we can
take uw = w. The system has a critical transition for w > 0 at Tc(w). Since the singular
part of the free energy close to a critical point behaves as (T −Tc)
2−α (α = −0.049(6) is the
specific-heat exponent of the RDIs universality class), we must have F (Xc) = 0, where Xc
is the value of X obtained by setting T = Tc(w) (see, e.g., Ref. 36 and references therein).
Hence, we obtain
w [Tc(w)/TIs − 1 + kw]
= Xc, (24)
and therefore
Tc(w)/TIs − 1 = (w/Xc)
1/φ − kw + · · · , (25)
where the dots indicate higher-order terms. This expression provides the w dependence of
the critical temperature for w small. Note that the nonanalytic term in Eq. (25) is suppressed
with respect to the analytic ones, because 1/φ ≈ 9.1. Thus, Tc(w) ≈ TIs(1 − kw + O(w
Since Tc(w) < TIs, we can also infer that k > 0. From the results for Tc(w) we estimate
k ≈ 2.2 for the ±J Ising model.
0.0 0.1 0.2
Ising
multicritical point
FIG. 12: The critical temperature Tc ≡ 1/βc vs 1− p.
V. CONCLUSIONS
In this paper we have studied the critical behavior of the 3D ±J Ising model at the
transition line between the paramagnetic and the ferromagnetic phase, which extends from
p = 1 to a multicritical (Nishimori) point at p = pN ≈ 0.767. We presented a FSS analysis of
MC simulations at various values of p in the region pN < p < 1. The results for the critical
exponents and other universal quantities are consistent with those of the RDIs universality
class. For example, we obtained ν = 0.682(3) and η = 0.036(2), which are in good agreement
with the presently most accurate estimates13 ν = 0.683(2) and η = 0.036(1) for the 3D RDIs
universality class. Therefore, our FSS analysis provides a strong evidence that the critical
behavior of the 3D ±J Ising along the ferromagnetic line belongs to the 3D RDIs universality
class.
We also note that the random-exchange interaction in the ±J Ising model gives rise
to frustration, while the RDIs universality class describes transitions in generic diluted
Ising systems with ferromagnetic exchange interactions. This implies that frustration is
irrelevant at the ferromagnetic transition line of the 3D ±J Ising model. Moreover, the
observed scaling corrections are consistent with the RDIs leading and next-to-leading scaling
correction exponents ω = 0.33(3) and ω2 = 0.82(8). This indicates that frustration does not
introduce new irrelevant perturbations at the RDIs fixed point with RG dimension yf & −1.
APPENDIX A: NOTATIONS
We define the two-point correlation function
G(x) ≡ 〈σ0 σx〉, (A1)
where the overline indicates the quenched average over the Jxy probability distribution.
Then, we define the corresponding susceptibility χ ≡
xG(x) and the correlation length ξ
G̃(0)− G̃(qmin)
q̂2minG̃(qmin)
, (A2)
where qmin ≡ (2π/L, 0, 0), q̂ ≡ 2 sin q/2, and G̃(q) is the Fourier transform of G(x). We also
consider quantities that are invariant under RG transformations in the critical limit. Beside
the ratio
Rξ ≡ ξ/L, (A3)
we consider the quartic cumulants U4, U22 and Ud defined by
, (A4)
U22 ≡
µ22 − µ2
Ud ≡ U4 − U22,
where
µk ≡ 〈 (
k〉 . (A5)
We also define corresponding quantities Ū4, Ū22, and Ūd at fixed Rξ = 0.5943. Finally, we
consider the derivative R′ξ of Rξ, and U
4 of U4, with respect to β ≡ 1/T , which allow one to
determine the critical exponent ν.
1 H. Nishimori, Prog. Theor. Phys. 66, 1169 (1981).
2 P. Le Doussal and A.B. Harris, Phys. Rev. Lett. 61, 625 (1988).
3 N. Kawashima and H. Rieger, in Frustrated Spin Systems, edited by H.T. Diep (World Scientific,
Singapore, 2004); cond-mat/0312432.
4 H. Katzgraber, M. Körner, and A.P. Young, Phys. Rev. B 73, 224432 (2006).
5 Y. Ozeki and H. Nishimori, J. Phys. Soc. Japan 56, 3265 (1987).
6 R.R.P. Singh, Phys. Rev. Lett. 67, 899 (1991).
7 Y. Ozeki and N. Ito, J. Phys. A 31, 5451 (1998).
8 Refs. 5,6,7 report the estimates pN = 0.767(2), pN = 0.7656(20), and pN = 0.7673(3) respec-
tively.
9 K. Hukushima, J. Phys. Soc. Japan 69, 631 (2000).
10 A. Pelissetto and E. Vicari, Phys. Rept. 368, 549 (2002).
11 R. Folk, Yu. Holovatch, and T. Yavors’kii, Uspekhi Fiz. Nauk 173, 175 (2003) [Phys. Usp. 46,
175 (2003)].
12 H. G. Ballesteros, L. A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and
J. J. Ruiz-Lorenzo, J. Phys. A 32, 1 (1999).
13 M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp.
P02016 (2007).
14 P. Calabrese, V. Mart́ın-Mayor, A. Pelissetto, and E. Vicari, Phys. Rev. E 68, 036136 (2003).
15 A. Pelissetto and E. Vicari, Phys. Rev. B 62, 6393 (2000).
16 H.G. Ballesteros, L.A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and J.J.
Ruiz-Lorenzo, Phys. Rev. B 58, 2740 (1998).
17 M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E 65, 066127 (2002).
18 Y. Deng and H.W.J. Blöte, Phys. Rev. E 68, 036125 (2003).
19 W. Janke, in Proceedings of the XXIII International Symposium on Lattice Field Theory,
Dublin, July 2005, POS(LAT2005)018
20 N. Ito, Y. Ozeki, and H. Kitatani, J. Phys. Soc. Jpn. 68, 803 (1999).
21 P. Calabrese, P. Parruccini, A. Pelissetto, and E. Vicari, Phys. Rev. E 69, 036120 (2004).
22 M. Hasenbusch, J. Phys. A 32, 4851 (1999).
23 M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B 74, 144506 (2006);
Phys. Rev. B 63, 214503 (2001).
24 See, e.g., S. Wansleben, J.B. Zabolitzky, and C. Kalle, J. Stat. Phys. 37, 271 (1984); G. Bhanot,
D. Duke, and R. Salvador, Phys. Rev. B 33, 7841 (1986).
25 M. Lüscher, Comput. Phys. Commun. 79, 100 (1994).
26 The SIMD-oriented fast Marsenne twister random number generator has been introduced by
M. Matsumoto and M. Saito. Details can be found in M. Saito, Master Thesis (2007) and at
http://www.math.sci.hiroshima-u.ac.jp/∼m-mat/MT/emt.html.
27 In order to make the use of these expensive (in terms of CPU-time) generators affordable,
we employed the same sequence of random numbers for the update of all nbit systems (for
the initialization of the configurations at the beginning of the simulation we used independent
random numbers for each of the systems). This may give rise to a statistical correlation among
the nbit systems. This effect is probably small and we have not detected it. Anyway, in order to
ensure a correct estimate of the statistical error, all nbit systems that use the same sequence of
random numbers have been put in the same bin in our jackknife analysis.
28 R.H. Swendsen and J-S. Wang, Phys. Rev. Lett. 58, 86 (1987).
29 U. Wolff, Phys. Rev. Lett. 62, 361 (1989).
30 D. Ivaneyko, J. Ilnytskyi, B. Berche, and Yu. Holovatch, Physica A 370, 163 (2006).
31 H.G. Ballesteros, L.A. Fernández, V. Mart́ın-Mayor, A. Muñoz Sudupe, G. Parisi, and J.J.
Ruiz-Lorenzo, Nucl. Phys. B 512, 681 (1998).
32 The crossover exponent from pure Ising to RDIs critical behavior is the Ising specific-heat
exponent34 αIs, see also Sec. IVC. This implies that the crossover scaling variable in the FSS at
Tc is given by the combination X = cwL
αIs/νIs , where w = 1 − p, αIs/νIs = 0.1740(8), and c is
a normalization constant. When w → 0, strong crossover effects are expected for X . 1, which
corresponds to L . (cw)−5.75. The RDIs asymptotic critical behavior is observed for X ≫ 1.
33 A. B. Harris, J. Phys. C 7, 1671 (1974).
34 A. Aharony, in Phase Transitions and Critical Phenomena, edited by C. Domb and J. Lebowitz
(Academic Press, New York, 1976), Vol. 6, p. 357.
35 I.D. Lawrie and S. Sarbach, in Phase Transitions and Critical Phenomena, Vol. 9, edited by C.
Domb and J. Lebowitz (Academic, London, 1984).
36 A. Pelissetto and E. Vicari, cond-mat/0702273.
