The critical behavior of three-dimensional Ising spin glass models
aa r X i v : . [ c ond - m a t . d i s - nn ] S e p The critical behavior of three-dimensional Ising spin glass models
Martin Hasenbusch, Andrea Pelissetto, and Ettore Vicari Institut f¨ur Theoretische Physik, Universit¨at Leipzig,Postfach 100 920, D-04009 Leipzig, Germany. Dipartimento di Fisica dell’Universit`a di Roma “La Sapienza” and INFN, Roma, Italy. Dipartimento di Fisica dell’Universit`a di Pisa and INFN, Pisa, Italy. (Dated: October 24, 2018)
Abstract
We perform high-statistics Monte Carlo simulations of three-dimensional Ising spin-glass modelson cubic lattices of size L : the ± J (Edwards-Anderson) Ising model for two values of the disorderparameter p , p = 0 . p = 0 . L = 28 and L = 20, respectively), and the bond-dilutedbimodal model for bond-occupation probability p b = 0 .
45 (up to L = 16). The finite-size behaviorof the quartic cumulants at the critical point allows us to check very accurately that these modelsbelong to the same universality class. Moreover, it allows us to estimate the scaling-correctionexponent ω related to the leading irrelevant operator: ω = 1 . p b = 0 . p b = 0 .
35 (up to L = 10) and of the Ising spin-glass model with Gaussian bond distribution (up to L = 8) also support the existence of a uniqueIsing spin-glass universality class. A careful finite-size analysis of the Monte Carlo data which takesinto account the analytic and the nonanalytic corrections to scaling allows us to obtain precise andreliable estimates of the critical exponents ν and η : we obtain ν = 2 . η = − . PACS numbers: 75.10.Nr, 64.60.Fr, 75.40.Cx, 05.10.Ln . INTRODUCTION The Ising model with random ferromagnetic and antiferromagnetic couplings is a sim-plified model for disordered uniaxial magnetic materials which show glassy behavior insome region of their phase diagram, such as Fe − x Mn x TiO and Eu − x Ba x MnO ; see, e.g.,Refs. 2,3,4. The random nature of the short-ranged interactions is mimicked by nearest-neighbor random bonds. This model is also interesting per se , since it provides a laboratoryin which the combined effect of quenched disorder and frustration can be studied.It is now well established that three-dimensional Ising spin-glass models present a high-temperature paramagnetic phase and, for some values of the parameters, a low-temperatureglassy phase (if the frustration is small, the low-temperature phase is ferromagnetic). Thetwo phases are separated by a continuous phase transition, which is expected to haveuniversal features, i.e. to belong to a universality class which is independent of the de-tails of the model and, in particular, of the disorder distribution. Several numericalworks have addressed these issues,considering various Ising spin-glass models, characterized by different disorder distribu-tions, with or without dilution. Over the years many estimates of the critical exponentshave been obtained. We mention the most recent ones for the correlation-length expo-nent ν : ν = 2 . ν = 2 . ν = 1 . ν = 1 . ν = 2 . ν = 1 . obtained from simulations of the symmetric model with bimodal distribution; ν = 2 . p b = 0 . ν = 2 . and ν = 2 . for the symmetric model with Gaussian disorder distribution; ν = 2 . which is ex-pected to be in the same Ising spin-glass universality class. Moreover, the analysis ofdifferent quantities has often given different estimates of the same critical exponent, evenin the same model. For instance, recent Monte Carlo (MC) studies find significant dis-crepancies among the estimates of the exponent ν obtained from the finite-size scaling (FSS)at T c of the temperature derivative of ξ/L , of the Binder cumulant, and of the overlap sus-ceptibility. For the bimodal Ising model Ref. 30 quotes ν = 2 . ν = 2 . ν = 1 . which allows us to predict the corrections to the asymptotic criticalbehavior for the different quantities. In particular, we show that the analytic dependenceof the relevant scaling fields on the model parameters, such as the temperature, may giverise to scaling corrections that decay as powers of L − /ν , where L is the linear size of thelattice. Since ν ≈ .
45 in Ising spin-glass systems, they decay quite slowly and may give riseto systematic deviations which are difficult to detect, given the small interval of values of L which can be considered in MC simulations. Their presence explains some inconsistencies inthe standard analyses of MC data reported in the literature. Thus, it is crucial to take scalingcorrections into account for an accurate study of the critical behavior, for a robust check ofuniversality among different models, and for reliable estimates of universal quantities suchas the critical exponents.In this paper we report a high-statistics MC study of different Ising spin-glass models.2e consider the ± J Ising model for two values of the disorder parameter, the bond-dilutedsymmetric bimodal model with various values of the dilution, and also the model withGaussian disorder distribution. We determine the FSS behavior of several RG invariantquantities, such as the ratio R ξ ≡ ξ/L ( ξ is the second-moment correlation length) andthe quartic cumulants defined from the overlap variables. We verify with good precisiontheir independence on the model and disorder distribution, providing an accurate evidenceof universality. Then, we obtain an estimate of the leading correction-to-scaling exponent ω : ω = 1 . ω . We obtain ν = 2 . ,η = − . . (1)Then, using scaling and hyperscaling relations we obtain β = ν (1 + η ) / . ,γ = (2 − η ) ν = 5 . ,α = 2 − ν = − . . (2)In this work we extend the results presented in Ref. 32. First, we have significantly increasedthe statistics of the large- L data for the bimodal symmetric Ising model and we have addedsome data for other diluted models and for the bimodal model with Gaussian distributedcouplings. Second, we present a much more detailed analysis of the critical-point data anda new analysis of the high-temperature data obtained in the parallel-tempering simulations.This allows us to check the universality of the critical behavior of the correlation lengthand of the susceptibility in the high-temperature phase. Finally, we discuss the extended-scaling scheme, which is inspired by the high-temperature expansion. As already notedin Ref. 29, this scheme shows an apparent improvement of the scaling behavior with respectto the naive approach in which scaling corrections are simply neglected, at least for somequantities, e.g., the overlap susceptibility. However, as we shall show, such an improvementis only marginal for the purpose of obtaining accurate results. Indeed, this requires to takeinto account the analytic and nonanalytic scaling corrections predicted by the RG theory.The paper is organized as follows. In Sec. II we define the models we investigate andthe quantities that are computed in the MC simulation. In Sec. III we derive the FSSpredictions of the RG theory which are the basis of our FSS analyses. Some details arereported in App. A. In Sec. IV and in App. B we give some details on the MC simulations.In Sec. V A we discuss universality, verifying that the infinite-volume limit of the quarticcumulants and of R ξ ≡ ξ/L is independent of the model. In Sec. V B we compute the leadingcorrection-to-scaling exponent ω . In Sec. VI we compute the critical exponents ν and η andthe critical temperature for the different models by using the data close to the critical point.In Sec. VII we present a global analysis of all available high-temperature data obtained inour parallel-tempering MC simulations. We again determine the critical exponents and showthat the FSS behavior of several quantities is universal. Moreover, we discuss the extended-scaling scheme of Ref. 29. In Sec. VIII we compute the high-temperature zero-momentumquartic couplings. Finally, in Sec. IX we present our conclusions.3 IsT para
MGP glassyN line
RDI ferro1 − p ? ISG
FIG. 1: Phase diagram of the three-dimensional ± J (Edwards-Anderson) Ising model in the T - p plane for p ≥ /
2, i.e., 1 − p ≤ /
2. The phase diagram is symmetric under p → − p . II. ISING SPIN GLASS SYSTEMSA. The ± J Edwards-Anderson Ising model and its phase diagram
We consider the ± J Edwards-Anderson Ising model on a simple cubic lattice of linearsize L with periodic boundary conditions. The corresponding Hamiltonian is H = − X h xy i J xy σ x σ y , (3)where σ x = ±
1, the sum is over the nearest-neighbor lattice sites, and the exchange interac-tions J xy are uncorrelated quenched random variables with probability distribution P ( J xy ) = pδ ( J xy −
1) + (1 − p ) δ ( J xy + 1) . (4)The usual bimodal Ising spin glass model, for which [ J xy ] = 0 (brackets indicate the averageover the disorder distribution), corresponds to p = 1 /
2. For p = 1 / J xy ] = 2 p − = 0 , (5)and ferromagnetic (or antiferromagnetic) configurations are energetically favored. Note thatthe free energy and also the correlations of the overlap variables that we shall define beloware invariant under p → − p and thus we can always assume p ≥ / T - p phase diagram of the three-dimensional ± J Ising model is sketched in Fig. 1.The high-temperature phase is paramagnetic for any p . The nature of the low-temperaturephase depends on the value of p : it is ferromagnetic for small values of 1 − p , while it isglassy with vanishing magnetization for sufficiently large values of 1 − p . The three phasesare separated by transition lines, which meet at a magnetic-glassy multicritical point (MGP),4sually called Nishimori point, which is located along the so-called Nishimori line defined by the relation ( p ≥ / T = T N ( p ) , T N ( p ) = 2ln p − p . (6)On the Nishimori line the magnetic and the overlap two-point correlation functions are equal.The paramagnetic-ferromagnetic (PF) transition line starts at the Ising transition at p = 1 and extends up to the MGP at p = p ∗ . For p = 1 there is the standard Ising transitionat T Is = 4 . p = 1 is a multicritical point withcrossover exponent φ = α Is , where α Is = 0 . > p > p ∗ belongs to the randomly-dilute Ising (RDI) universalityclass, characterized by the correlation-length critical exponent ν = 0 . η = 0 . T - p plane are T ∗ = 1 . p ∗ = 0 . thethermal exponent ν = 1 . φ = 1 . η = − . p =1 / p = 1 − p ∗ = 0 . p → − p of the phasediagram). A reasonable hypothesis is that the PG critical behavior is independent of p , asfound in mean-field models. Hence, for any p ∗ > p > − p ∗ the PG transition is expectedto belong to the same universality class (named ISG in Fig. 1) as that of the bimodal Isingspin-glass model at p = 1 /
2. The critical behavior along this transition line is the mainsubject of this paper. As we shall see, the universality hypothesis is fully confirmed by ourFSS analyses at p = 0 . p = 0 . p the following inequality holds: | [ h σ x σ y i T ] | ≤ [ |h σ x σ y i T N ( p ) | ] , (7)where the subscripts indicate the temperature of the thermal average, and T N ( p ) is thetemperature along the Nishimori line, defined in Eq. (6). This relation shows that ferromag-netism can only exist in the region p > p ∗ and that the system is maximally magnetizedalong the Nishimori line. Ref. 50 (see also Refs. 41,51) also argues that the FG transitionline coincides with the line p = p ∗ , from T = T ∗ to T = 0. Recent numerical investiga-tions of the two-dimensional ± J model have shown that this conjecture is not exact,though quantitative deviations are small: at T = 0 the critical value p c where ferromag-netism disappears is definitely larger than p ∗ , indicating a reentrant transition line. In threedimensions Ref. 55 quotes p c = 0 . p ∗ = 0 . ν = 1 . in which ferromagnetism and glassy order coexist, is found in mean-field models such asthe infinite-range Sherrington-Kirkpatrick model. Its presence has been confirmed in the ± J Ising model on a Bethe lattice. However, there is no evidence of a mixed phase in the ± J Ising model on a cubic lattice and in related models. In particular, the numericalresults of Ref. 55 show that the onset of the glassy behavior at T = 0 occurs close to thepoint where the ferromagnetic phase disappears, and are consistent with a single transition5ithin numerical precision. Nevertheless, the existence of such a mixed phase is still an openproblem, as discussed in Ref. 58. B. Other Ising spin-glass models
We also consider the bond-diluted bimodal Ising model (BDBIM) defined by Hamiltonian(3) with bond probability distribution P ( J xy ) = p b (cid:20) δ ( J xy − J ) + 12 δ ( J xy + J ) (cid:21) + (1 − p b ) δ ( J xy ) . (8)A PG transition occurs for sufficiently small values of 1 − p b , i.e. for p b > p SG . Whileinvestigations at T = 0 indicate that p SG should be identified with the bond-percolationpoint ( p perc = 0 . ), recent investigations of thecritical behavior close to the percolation point suggest that p SG is larger than p perc . The model can be extended by considering the distribution P ( J xy ) = p b [ pδ ( J xy − J ) + (1 − p ) δ ( J xy + J )] + (1 − p b ) δ ( J xy ) . (9)In this case, for p b > p SG ( p SG may depend on p ) we expect a T - p phase diagram similarto the one sketched in Fig. 1 for the ± J Ising model without dilution, with a PF and a PGtransition line meeting at a multicritical point.A PG transition is also observed in the random-bond Ising spin-glass model with Gaussianbond distribution: P ( J xy ) = 1 √ π e − J xy / . (10)This transition is expected to be in the same universality class as that of the bimodal Isingmodel. C. Overlap thermodynamic quantities
In this work we focus on the critical behavior of the overlap parameter q x ≡ σ (1) x σ (2) x , (11)where the spins σ ( i ) x belong to two independent replicas with the same disorder realization { J xy } . The corresponding correlation function is G ( x ) ≡ [ h q q x i ] = [ h σ σ x i ] , (12)where the angular and square brackets indicate the thermal average and the quenched averageover { J xy } , respectively. We define the susceptibility χ and the second-moment correlationlength ξ as χ ≡ X x G ( x ) , (13) ξ ≡
14 sin ( p min / e G (0) − e G ( p ) e G ( p ) , (14)6here p = ( p min , , p min ≡ π/L , and e G ( q ) is the Fourier transform of G ( x ).We also define the RG invariant quantities R ξ ≡ ξ/L, (15) U ≡ [ µ ][ µ ] , (16) U ≡ [ µ ] − [ µ ] [ µ ] , (17)where µ k ≡ h ( X x q x ) k i . (18)We call them phenomenological couplings and denote them by R in the following.In the high-temperature paramagnetic phase, we also consider the zero-momentum quarticcouplings G ≡ − χ ξ χ , (19) G ≡ − χ ξ χ , (20)where χ ≡ L (cid:0) [ µ ] − µ ] (cid:1) , (21) χ ≡ L (cid:0) [ µ ] − [ µ ] (cid:1) . (22)The critical limit T → T + c of the zero-momentum quartic couplings G and G is universal. III. FINITE-SIZE SCALING
In this section we summarize some basic results concerning FSS, which allow us to un-derstand the role of the analytic and nonanalytic scaling corrections. We consider two Isingspin-glass systems coupled by an interaction h X x q x = h X x σ (1) x σ (2) x , (23)where h is a constant external field. The model is defined on a cubic lattice of linear size L with periodic boundary conditions.By applying standard RG arguments we expect the disorder-averaged free-energy densityto be the sum of a regular part and a singular part: F ( β, h, L ) = F reg ( β, h, L ) + F sing ( β, h, L ) , (24)where β ≡ /T . The regular part is expected to depend on L only through exponentiallysmall terms, while the singular part encodes the critical behavior. The starting point of FSS7s the scaling behavior of the singular part of the free-energy density (see, e.g., Refs. 36,47,64,65): F sing ( β, h, L ) = L − d F ( u h L y h , u t L y t , { v i L y i } ) , (25)where d is the space dimension, u h and u t are the scaling fields associated with h andthe reduced temperature t ∼ − β/β c (their RG dimensions are y h = ( d + 2 − η ) / y t = 1 /ν , respectively), and v i are irrelevant scaling fields with y i <
0. At the critical pointwe have u t ( t = 0 , h = 0) = 0 and u h ( t = 0 , h = 0) = 0, while, generically, we expect v i ( t = 0 , h = 0) = 0. Since y i <
0, for large L the free energy can be expanded in powers of { v i L y i } . Therefore, we can write F sing ( β, h, L ) = L − d f ( u h L y h , u t L y t ) + v ω L − d − ω f ω ( u h L y h , u t L y t ) + . . . (26)where the leading nonanalytic correction-to-scaling exponent ω is related to the RG dimen-sion y ω of the leading irrelevant scaling field v ω ≡ v , ω = − y ω . The scaling fields areanalytic functions of the system parameters—in particular, of h and t —and are expectednot to depend on L . Note also that the size L is expected to be an exact scaling field forperiodic boundary conditions. For a general discussion of these issues, see Ref. 65, Sec. IIIof Ref. 66, and references therein. In general, u t and u h can be expanded as u h = h ¯ u h ( t ) + O ( h ) , ¯ u h ( t ) = a h + a t + O ( t ) , (27) u t = c t t + c t + c h + c h t + O ( t , h , h t ) , (28)where we used the fact that the free energy is symmetric under h → − h . In the expansion of u h,t around the critical point h, t = 0, the terms beyond the leading ones give rise to analytic scaling corrections.The scaling behavior of zero-momentum thermodynamic quantities can be obtained byperforming appropriate derivatives of F with respect to h . For instance, for the overlapsusceptibility at h = 0 we obtain χ ( β, L ) = ∂ F ∂h (cid:12)(cid:12)(cid:12)(cid:12) h =0 = L − η ¯ u h ( t ) g ( u t L y t ) (cid:2) v ω L − ω g ω ( u t L y t ) + . . . (cid:3) + g reg ( β ) . (29)The function g reg ( β ) represents the contribution of the regular part F reg ( β, h, L ) of the free-energy density and is L independent (apart from exponentially small terms). It gives rise toa correction proportional to L η − . Analogous formulae hold for the 2 n -point susceptibilities.The FSS of the phenomenological couplings is given by R ( β, L ) = r ( u t L y t ) + v ω r ω ( u t L y t ) L − ω + . . . = R ∗ + r ′ (0) c t tL y t + . . . + c ω L − ω + . . . , (30)where R ∗ ≡ r (0), c ω = v ω r ω (0), and the second line holds only very close to the criticalpoint, for tL y t ≪
1. A proof of Eq. (30) for the phenomenological couplings U , U , and R ξ is presented in App. A.The thermal RG exponent y t = 1 /ν is usually computed from the FSS of the derivative R ′ of a phenomenological coupling R with respect to β at β c . Using Eq. (30) one obtains R ′ ≡ ∂R∂β = L y t ( ∂ β u t ) (cid:2) r ′ ( u t L y t ) + v ω L − ω r ′ ω ( u t L y t ) + · · · (cid:3) . (31)8ne may also consider the derivative χ ′ ≡ dχ/dβ of the susceptibility χ . From Eq. (29) weobtain χ ′ = L − η + y t ¯ u h ∂ β u t (cid:8) g ′ ( u t L y t ) + v ω L − ω [ g ′ ( u t L y t ) g ω ( u t L y t ) + g ( u t L y t ) g ′ ω ( u t L y t )] (cid:9) +2 L − η ¯ u h ∂ β ¯ u h g ( u t L y t ) + · · · + g ′ reg ( β ) . (32)Note that the second term in the right hand side gives rise to scaling corrections proportionalto L − y t = L − /ν , while the background term g ′ reg ( β ) leads to corrections proportional to L − y t − η .At T = T c , setting t = 0 in the above-reported equations, we obtain: R = R ∗ + c ω L − ω + . . . , (33) χ = cL − η (1 + c ω L − ω + . . . ) , (34) R ′ = cL /ν (1 + c ω L − ω + . . . ) , (35) χ ′ = cL − η +1 /ν (1 + c ω L − ω + . . . + c a L − /ν + . . . ) . (36)Note that, unlike the temperature derivative R ′ of a RG-invariant quantity, χ ′ also presentsan L − /ν scaling correction, due to the analytic dependence on t of the scaling field u h (forthis reason we call it analytic correction). Since, as we shall see, in the Ising spin-glass case1 /ν ≈ . ω ≈ .
0, the scaling corrections in χ ′ decay significantly more slowly thanthose occurring in R ′ . This makes the ratio χ ′ χ ∼ L /ν (37)unsuitable for a precise determination of ν and explains the significant discrepancies observedin Ref. 30.Instead of computing the various quantities at fixed Hamiltonian parameters, one canalso consider FSS keeping a phenomenological coupling R fixed at a given value R f . Thismeans that, for each L , one determines β f ( L, R f ), such that R ( β = β f ( L, R f ) , L ) = R f , andthen considers any quantity at β = β f ( L, R f ). The value R f can be specified at will, as longas R f is taken between the high- and low-temperature fixed-point values of R . For R f = R ∗ ,where R ∗ is defined in Eq. (33), β f converges to β c as β f − β c ∼ L − /ν , (38)while for R f = R ∗ we have β f − β c ∼ L − /ν − ω . (39)Indeed, if u t,f ( L, R f ) is the value of u t for β = β f ( L, R f ), we obtain from Eq. (30) u t,f L y t = B ( R f ) + v ω B ω ( R f ) L − ω + · · · (40)where, using Eq. (30), r ( B ( x )) = x , (41) B ω ( x ) = − r ω ( B ( x )) r ′ ( B ( x )) . (42)9ow, if R f = R ∗ , we have B ( R f ) = 0, which implies u t,f ∼ L − y t − ω , hence Eq. (39). Other-wise, B ( R f ) is different from zero and we obtain the behavior (38).If we now substitute relation (40) into Eqs. (29), (30), (31), and (32), we obtain theexpansion of the different quantities at fixed R f , which we denote by adding a bar: given O ( β, L ), we define ¯ O ( L, R f ) ≡ O [ β f ( L, R f ) , L ]. For R f = R ∗ , since u t,f ∼ L − y t − ω wereobtain Eqs. (33), (34), (35), and (36), with different coefficients, of course. If R f = R ∗ ,we must be more careful. The behavior of another phenomenological coupling R α does notchange qualitatively and we still have¯ R α ( L, R f ) ≈ ¯ R ∗ α + c α L − ω + . . . , (43)where ¯ R ∗ α = r α [ B ( R f )] is universal. It depends on R f and satisfies ¯ R ∗ α = R ∗ α for R f = R ∗ .Also ¯ χ ′ behaves as it does at fixed T = T c , i.e. it follows Eq. (36). On the other hand, ¯ χ and¯ R ′ present additional analytic corrections. Indeed, since [see Eq. (27)]¯ u h,f = a h + a c t B ( R f ) L − y t + · · · (44)(nonanalytic O ( L − ω ) corrections have been neglected), Eq. (29) gives¯ χ α ( L, R f ) = L − η (cid:20) a h + 2 a a h c t B ( R f ) L − y t + O ( L − y t ) (cid:21) g ( B ( R f ))[1 + O ( L − ω )] . (45)If R f = R ∗ , B ( R f ) = 0, and thus analytic corrections occur. Note that, if R f is close to R ∗ ,since B ( R ∗ ) = 0, we have B ( R f ) ∼ R f − R ∗ . (46)Hence, in this case the analytic corrections are small, of order ( R f − R ∗ ) L − /ν . In general,corrections of order L − k/ν have amplitudes proportional to ( R f − R ∗ ) k . IV. MONTE CARLO SIMULATIONS
In the MC simulations we employed the Metropolis algorithm, the random-exchangemethod (often called parallel-tempering or multiple Markov-chain method), and multispincoding. See App. B for details on their implementation.We simulated the ± J Ising model at p = 0 . L = 3-14,16,20,24,28, at p = 0 . L = 3-12,14,16,20, and the BDBIM at p b = 0 .
45 for L = 4-12,14,16. We averagedover a large number N s of disorder samples: N s ≈ . · up to L = 12, N s / ≈ , , , , ,
18, respectively for L = 13 , , , , ,
28 in the case of the ± J Ising model at p = 0 .
5. Similar statistics were collected at p = 0 . L = 20 wherethe statistics were approximately 1/3 of those for p = 0 . p b = 0 . p b = 0 .
35 and the Ising modelwith Gaussian distributed couplings, but only performed simulations for small values of L : L = 4 , , ,
10 for the BDBIM and L = 4 , , , L and model we performed parallel-tempering runs. This allowed us to estimatethe different quantities in a large interval [ β min , β max ]. To fix β max we used the results ofRef. 30, which provided the best estimates of R ∗ ξ at the time we started our simulations:10 ABLE I: Quartic cumulant ¯ U at fixed R ξ = 0 .
63 for the ± J Ising model at p = 0 . p = 0 . p b = 0 . , . , .
35, and for the Ising spin-glass model with Gaussian bonddistribution. L ± J p =0 . ± J p =0 . BDBIM p b =0 . BDBIM p b =0 . BDBIM p b =0 . Gaussian4 1.48231(6) 1.46813(5) 1.49036(8) 1.48480(6) 1.49164(9) 1.49145(13)5 1.48985(6) 1.47597(6) 1.49853(8) 1.4996(2)6 1.49446(6) 1.48193(6) 1.50300(9) 1.49618(9) 1.50788(9) 1.5033(2)7 1.49753(6) 1.48642(6) 1.50544(9)8 1.49987(6) 1.48984(6) 1.50714(9) 1.50082(9) 1.51320(13) 1.5063(5)9 1.50136(6) 1.49260(6) 1.50815(9)10 1.50273(6) 1.49478(6) 1.50889(9) 1.50382(11) 1.5146(2)11 1.50383(6) 1.49665(6) 1.50946(9)12 1.50469(6) 1.49781(7) 1.50984(13)13 1.50541(11)14 1.50618(10) 1.50030(12) 1.5103(3)16 1.50702(13) 1.50220(13) 1.5113(3)20 1.5081(3) 1.5048(5)24 1.5089(4)28 1.5108(13) β max was chosen so that R ξ ( β max , L ) ≈ .
63. Onlyin the most recent simulations (two runs with L = 20 and 24) of the ± J Ising model at p = 0 . β max such that R ξ ( β max , L ) ≈ .
66. We checked thermalization by usingthe recipe outlined in Ref. 30, see App. B for some details.As already emphasized in Refs. 46,69, in high-precision MC studies of random systemsone should be careful when computing disorder averages of products of thermal expectations,for instance the quartic cumulant U defined in Eq. (17). Indeed, naive estimators have abias that may become significantly larger than the statistical error if N s is large. We use(essentially) bias-free estimators, defined as discussed in Ref. 46; some details are given inApp. B.In total, the MC simulations took approximately 40 years of CPU time on a single coreof a 2.4 GHz AMD Opteron processor. V. UNIVERSALITY AND CORRECTION-TO-SCALING EXPONENT ω A. Analysis of the quartic cumulants at fixed R ξ In this section we investigate the universality of the critical behavior of Ising spin-glassmodels by comparing the limit for L → ∞ of U and U computed at fixed R ξ , denoted by¯ U and ¯ U , respectively. As discussed in Sec. III, for sufficiently large L , ¯ U and ¯ U areexpected to behave as ¯ U = ¯ U ∗ + c L − ω , (47)11 ABLE II: Quartic cumulant ¯ U at fixed R ξ = 0 .
63 for the ± J Ising model at p = 0 . p = 0 . p b = 0 . , . , .
35, and for the Ising spin-glass model with Gaussian bonddistribution. L ± J p =0 . ± J p =0 . BDBIM p b =0 . BDBIM p b =0 . BDBIM p b =0 . Gaussian4 0.13714(6) 0.13955(6) 0.13581(9) 0.13602(7) 0.14635(10) 0.13522(15)5 0.14088(7) 0.14087(6) 0.13935(10) 0.14014(24)6 0.14277(7) 0.14181(7) 0.14166(10) 0.14193(10) 0.14501(10) 0.14227(24)7 0.14392(6) 0.14262(7) 0.14308(10)8 0.14478(7) 0.14318(7) 0.14415(10) 0.14414(10) 0.14597(15) 0.1434(5)9 0.14522(7) 0.14380(7) 0.14470(10)10 0.14561(7) 0.14418(7) 0.14518(10) 0.14536(12) 0.14586(24)11 0.14595(7) 0.14453(7) 0.14566(10)12 0.14618(7) 0.14465(7) 0.14605(14)13 0.14650(11)14 0.14671(10) 0.14531(13) 0.1462(3)16 0.14675(14) 0.14553(14) 0.1469(4)20 0.1469(4) 0.1458(5)24 0.1477(5)28 0.1496(14) where the constants ¯ U ∗ are universal (model independent) but depend on the fixed value of R ξ ; ω is the leading scaling-correction exponent.In Tables I and II we report the estimates of ¯ U and ¯ U at fixed R ξ = 0 .
63 for differentmodels. Without performing any analysis, one can immediately observe that the resultsobtained for the different models are very close, and appear to approach the same large- L limit as L increases. For instance, the estimates of ¯ U for the largest lattices differ at mostby 0.5%, while those of ¯ U vary by a few percent. This already provides strong support touniversality.For a more detailed analysis, let us first consider the three models for which we have mostdata, the ± J model at p = 0 . p = 0 .
7, and the BDBIM at p b = 0 .
45. The MC estimatesof ¯ U ( L ) are shown in Fig. 2 versus 1 /L . The results for the ± J Ising model at p = 0 . p = 0 . L → ∞ . Theysupport the universality of the critical behavior and show the presence of scaling correctionswith an exponent ω ≈ .
0. In the case of the BDBIM with p b = 0 .
45, the data apparentlyshow a faster approach to the same infinite-volume limit. A fit of ¯ U ( L ) to ¯ U ∗ + cL − ε usingall data with L > ε ≈
2. However, for sufficiently large L , L &
10 say, scaling corrections are consistent with a linear 1 /L behavior, as in the othermodels. The corresponding amplitude | c | is quite small, at least a factor of two smaller thanin the undiluted bimodal model. This means that for L .
10 next-to-leading corrections toscaling dominate and determine the apparent approach of the data to the infinite-volumelimit. Since the ratios of the amplitudes of the leading scaling corrections are universal, asmall | c | implies that the leading nonanalytic scaling correction is small for any quantity.Thus, in this model the approach of any thermodynamic quantity to the critical limit should12 .00 0.05 0.10 0.15 0.20 0.25 L -1 U p=0.5p=0.7p b =0.45 FIG. 2: (Color online) Estimates of ¯ U ( L ) vs L − , for the ± J model at p = 0 . p = 0 .
7, andthe BDBIM at p b = 0 .
45. The dotted lines are drawn to guide the eye. L - U p=0.5p=0.7p b =0.45 FIG. 3: (Color online) Estimates of ¯ U ( L ) vs L − , for the ± J model at p = 0 . p = 0 .
7, andthe BDBIM at p b = 0 .
45. The dotted lines are drawn to guide the eye. be faster than in the ± J models without dilution, as already noted in Ref. 28, at least forsufficiently large lattice sizes, where next-to-leading corrections can be neglected. Similarresults are obtained for ¯ U , see Fig. 3. However, in this case next-to-leading correctionsappear to be more important, since the 1 /L behavior is observed for somewhat larger valuesof L .In order to estimate the universal infinite-volume values ¯ U ∗ and ¯ U ∗ , we perform fits of¯ U to ¯ U ∗ + c L − ε . The exponent ε is either a free parameter, or is set equal to 1, whichcorresponds, as we discuss below, to our best estimate of the leading scaling-correction13 ABLE III: Results of the fits of ¯ U and ¯ U at fixed R ξ = 0 .
63 to ¯ U ∗ + cL − ε for several valuesof L min , the minimum lattice size allowed in the fits. DOF is the number of degrees of freedom ofthe fit. The fits labelled “ p = 0 . p = 0 .
7” are fits in which the data for p = 0 . p = 0 . U ∗ and ¯ U ∗ .Model L min ¯ U ∗ χ / DOF ¯ U ∗ χ / DOF ± J , p = 0 . ε ε = 1 8 1.5144(2) 0.8 0.1490(2) 0.612 1.5139(4) 0.6 0.1488(4) 0.9 ± J , p = 0 . ε ε = 1 8 1.5143(2) 2.3 0.1478(2) 0.912 1.5153(5) 0.1 0.1482(5) 0.8BDBIM, p b = 0 .
45 free ε ε = 1 8 1.5153(3) 0.8 0.1496(3) 0.412 1.5148(12) 1.2 0.1489(13) 1.0 ± J , p = 0 . p = 0 . ε = 1 8 1.5143(1) 1.2 0.1485(1) 2.112 1.5145(3) 0.9 0.1486(3) 0.816 1.5133(9) 0.4 0.1489(10) 1.1 ε = 1 . exponent ω . In Table III we report the results. They are all consistent with the estimates¯ U ∗ = 1 . , ¯ U ∗ = 0 . . (48)The accuracy and stability of the fits can also be appreciated in Fig. 4, where the fit resultsare plotted versus the minimum lattice size L min allowed in the fits. We can thus concludethat the estimates of the quartic cumulants for the ± J Ising model and the BDBIM at p b = 0 .
45 are fully consistent with universality. The relative differences are approximatelyone per mille in the case of ¯ U and one per cent for ¯ U .We also considered the BDBIM for the other dilution values, i.e., for p b = 0 . p b =0 .
35, though in this case we have data only up to L = 10. The estimates of ¯ U are shownversus 1 /L in Fig. 5. For p b = 0 . /L and are largerthan at p b = 0 .
45. For the scaling-correction amplitude c defined in Eq. (47) we obtainthe estimate c ≈ − .
10, to be compared with the result c ≈ − .
11 for the ± J model at p = 0 . c ≈ − .
05 for the diluted model at p b = 0 .
45. For p b = 0 .
35 we do not haveenough data to estimate c ; however, the data reported in Fig. 5 apparently approach theinfinite-volume limit faster than at p b = 0 .
45. As indicated by the MC data for p b . .
27 ofRef. 63, the scaling corrections should increase as p b further decreases. This suggests that14
10 15 20 L min U p=0.5, free ε p=0.5, ε=1.0 p=0.7, free ε p=0.7, ε=1.0 p b =0.45, free ε p=0.5 & p=0.7, ε=1.0 * L min U p=1/2, free ε p=1/2, ε=1.0 p=0.7, free ε p=0.7, ε=1.0 p b =0.45, free ε p=0.5 & p=0.7, ε=1.0 * FIG. 4: Estimates of ¯ U ∗ and ¯ U ∗ , obtained by fitting the quartic cumulants ¯ U at fixed R ξ = 0 . U ∗ + cL − ε . Here L min is the minimum lattice size allowed in the fits. The fits labelled “ p = 0 . p = 0 .
7” are fits in which the data for p = 0 . p = 0 . U ∗ = 1 . U ∗ = 0 . the optimal model—the one with the smallest leading scaling corrections—corresponds to adilution in the range 0 . . p ∗ b . .
4: the model with p b = 0 .
35 should be close to the optimalone.The estimates of ¯ U and ¯ U for the Ising spin-glass model with Gaussian bond distributionare shown in Figs. 5 and 6: they are very close, and clearly approach, the estimates (48).The results for L = 8 differ by 0.5% ( ¯ U ) and 3% ( ¯ U ) from the asymptotic value.In conclusion, the results for the ± J model, for the BDBIM, and for the model withGaussian distributed couplings provide a very strong evidence of universality for the criticalbehavior of Ising spin-glass models along the PG transition line.15 .00 0.05 0.10 0.15 0.20 0.25 L -1 U p b =0.45p b =0.7p b =0.35 Gaussian
FIG. 5: Estimates of ¯ U ( L ) for the BDBIM for several values of p b and for the Ising spin-glass modelwith Gaussian bond distribution. The dotted lines correspond to the estimate (48), ¯ U ∗ = 1 . p b = 0 . p b = 0 .
45 show the linearbehavior of the data for their largest lattices. L -1 U U __ FIG. 6: ¯ U (below) and ¯ U (above) for the Ising spin-glass model with Gaussian bond distribution.The dotted lines correspond to the estimates (48), ¯ U ∗ = 1 . U ∗ = 0 . B. The RG exponent of the leading irrelevant operator
The analyses of ¯ U and ¯ U also give estimates of ω . The most precise ones are obtainedfrom fits of ¯ U . In Fig. 7 we show the estimates of ω obtained from fits of ¯ U to ¯ U ∗ + c p L − ω and from fits of the difference ¯ U ( p = 0 . L ) − ¯ U ( p = 0 . L ) to bL − ω . We estimate ω = 1 . . (49)Consistent but less precise results are obtained from fits of ¯ U . The result (49) is consistentwith the less precise estimates ω = 0 . +0 . − . and ω = 1 . L min ω U (p=0.5;L)U (p=0.7;L)U (p=0.5;L) & U (p=0.7;L)U (p=0.5;L)-U (p=0.7;L) FIG. 7: (Color online) Estimates of the leading correction-to-scaling exponent ω versus the min-imum lattice size L min used in the fit. The fits labelled “ p = 0 . p = 0 .
7” are fits in whichthe data for p = 0 . p = 0 . ω = 1 . respectively.In order to verify that scaling corrections are controlled by a single correction-to-scalingterm with exponent ω = 1 and we have not determined an effective exponent arising fromseveral almost degenerate exponents, we check that the ratio c /c is universal, where c is the scaling-correction amplitude appearing in Eq. (47), as predicted by standard RGarguments. In order to determine this ratio, we fit together the available estimates of ¯ U ( L )for the ± J Ising models at p = 0 . p = 0 . U ( L ) = ¯ U ∗ + c p L − ε , fixing ε = 1 andtaking ¯ U ∗ and the two amplitudes c p =0 . and c p =0 . as free parameters. Fits of the data for L ≥ L min = 10 give¯ U ∗ = 1 . , c ( p = 0 .
5) = − . , c ( p = 0 .
7) = − . , (50)( χ / DOF ≈ .
3, where DOF is the number of degrees of freedom of the fit), and¯ U ∗ = 0 . , c ( p = 0 .
5) = − . , c ( p = 0 .
7) = − . , (51)( χ / DOF ≈ . c /c = 0 . p = 0 . , (52) c /c = 0 . p = 0 . , (53)which are in good agreement. These results are quite stable with increasing L min . For L min = 12 we obtain c /c = 0 . c /c = 0 . p = 0 . p = 0 .
7, while for L min = 14 we find c /c = 0 . c /c = 0 . χ / DOF .
1. Moreover, the variation of the ratio c /c with respect to a change ofthe exponent ε by ± . ω , ω = 1 . O ( L − . )corrections as the corrections arising from the leading irrelevant scaling field.17 .00 0.02 0.04 0.06 0.08 L -1.4 β cross R ξ U FIG. 8: (Color online) Crossing point β cross ( L ) for the ± J Ising model at p = 0 .
5, from several pairs L, L , as obtained from R ξ and U , versus L − /ν − ω ∼ L − . . The dotted lines show the result of afit to β c + cL − . , which takes into account only the data satisfying L ≥ L min = 8—it correspondsto L − . ≈ . y -axis corresponds to the final estimate β c = 0 . L min β c FIG. 9: (Color online) Estimates of β c as obtained from fits of R ξ and U to Eq. (56) with ω = 1,for the ± J Ising model at p = 0 .
5. All fits have χ / DOF < .
5. The lines correspond to the finalestimate β c = 0 . .0 0.1 0.2 0.3 0.4 0.5 T c PF IsingPF RDIMNPPG ISGN line
FIG. 10: (Color online) The values of the critical temperature T c versus 1 − p for the varioustransitions of the ± J Ising model. The estimates at the paramagnetic-ferromagnetic transitionsare taken from Ref. 43 for the Ising transition, from Ref. 44 for the RDI transition, and from Ref. 42for that along the N line. VI. CRITICAL TEMPERATURE AND EXPONENTSA. Estimates of β c In order to determine the critical temperature, we perform a standard FSS analysis of R ξ and the quartic cumulants. Figure 8 shows the crossing points β cross ( L ) determined bysolving the equation R ( β cross , L ) = R ( β cross , L ) , (54)for the ± J Ising model at p = 0 .
5, using R ξ and U . In the large- L limit β cross ( L ) is expectedto converge to β c as β cross ( L ) − β c ∼ L − /ν − ω ∼ L − . , (55)since ν ≈ .
45 and ω ≈ .
0. We perform a combined fit of R ξ and U to β c + c L − ε with ε = − .
4, taking β c and the two amplitudes c R ξ and c U as free parameters. Using onlythe results with L ≥ L min = 8, we obtain β c = 0 . χ / DOF ≈ . R ξ and U are correlated the error is only indicative). The corresponding linesare drawn in Fig. 8. We also consider larger values of L min . We find that the estimates of β c vary significantly (much more than the statistical error) when L min varies, indicating thatthe systematic error due to the additional scaling corrections is significantly larger than thestatistical one. Equivalently, one can perform fits to R ( L, β c ) = R ∗ + bL − ω , (56)taking R = R ξ or U . The results are shown in Fig. 9. Taking into account the dependenceof the results on L min and their variation as ω varies by one error bar we obtain the finalestimates β c = 0 . , T c = 1 . . (57)19 L min ν wsc, R ξ =0.654 wsc, R ξ =0.63ε=1.0, R ξ =0.654ε=1.1, R ξ =0.654ε=0.9, R ξ =0.654ε=1.0, R ξ =0.63 FIG. 11: (Color online) Estimates of the exponent ν from the analysis of ¯ R ′ ξ , at R ξ = 0 . R ξ = 0 . ± J Ising model at p = 0 .
5, obtained by fitting the data without scalingcorrections (wsc), Eq. (59), and with scaling corrections, Eq. (60). We only show the results of thefits which satisfy χ / DOF . .
5, up to values of L min where the errors blow up. The dotted linescorrespond to the final estimate ν = 2 . We regard the error as conservative. We also consider the pseudocritical β f ( L ) defined inSec. III at fixed R ξ , which converges to β c as L → ∞ , cf. Eqs. (38) and (39). The results areconsistent with the estimate (57). The analysis of the crossing points and the fits to Eq. (56)also provide estimates of the limit L → ∞ of R ξ and U at the critical point. We obtain R ∗ ξ = 0 . , U ∗ = 1 . . (58)Concerning the other models, similar analyses of R ξ and U give β c = 0 . T c =1 . ± J model at p = 0 .
7, and β c = 1 . T c = 0 . p b = 0 .
45. The corresponding values of R ∗ ξ and U ∗ are consistent with the estimates (58).Finally, in Fig. 10 we collect all results for the critical temperature T c of the ± J Isingmodel at its various PF and PG transitions as a function of the disorder parameter p . B. The critical exponent ν We estimate the critical exponent ν from the finite-size behavior of R ′ ξ ≡ dR ξ /dβ and U ′ ≡ dU /dβ at a fixed value R ξ,f of R ξ . As we discussed in Sec. III, the best choicefor R ξ,f is R ∗ ξ , otherwise ¯ R ′ ξ and ¯ U ′ present analytic corrections. In the present case, theestimate (58) of R ∗ ξ is not very precise, hence the corrections of order L − /ν ∼ L − . cannotbe completely neglected. However, since their amplitude is proportional to ( R ξ − R ∗ ξ ), we canassume that they are small for R ξ,f close to R ∗ ξ . Thus, in order to estimate the systematicerror induced by them, we proceed as follows. We repeat the analysis of ¯ R ′ at fixed R ξ using two different values of R ξ and neglecting in both cases the L − /ν corrections: we use R ξ = 0 .
630 and R ξ = 0 . R ∗ ξ = 0 . L min ν wsc, R ξ =0.654 wsc, R ξ =0.63ε=1.0, R ξ =0.654ε=1.0, R ξ =0.63 FIG. 12: (Color online) Estimates of the exponent ν from the analysis of ¯ U ′ , at R ξ = 0 .
63 and R ξ = 0 . ± J Ising model at p = 0 .
5, obtained by fits without scaling corrections (wsc),Eq. (59), and fits with scaling corrections, Eq. (60). We only show the results of the fits whichsatisfy χ / DOF . .
5, up to values of L min where the errors blow up. The dotted lines correspondto ν = 2 . both values of R ξ we determine an estimate of ν . Our final result is obtained by linearinterpolation to R ξ = 0 . R ξ = 0 .
630 and R ξ = 0 . ν is obtained from fits of ¯ R ′ ξ and ¯ U ′ toln ¯ R ′ = a + 1 ν ln L, (59)ln ¯ R ′ = a + 1 ν ln L + bL − ǫ , (60)with ǫ = ω = 1 . R ′ ξ for the ± J Ising model with p = 0 . L min allowed in the fits. They arequite stable with respect to L min , depend very little on the fixed value R ξ , and change onlyslightly as ω varies by ± .
1, corresponding to one error bar. In particular, the fit of ¯ R ′ ξ at R ξ = 0 .
654 to Eq. (60) gives ν = 2 . ν = 2 . L ≥ L min = 8 ( L min = 10);here the error in brackets takes into account the uncertainty on ω . Analogously, the fit ofthe data at R ξ = 0 .
630 gives ν = 2 . ν = 2 . L min . Inboth cases χ / DOF ≈ . χ / DOF ≈ . ν = 2 . . R ξ = 0 .
645 for L min = 8 ,
10. The systematic error due to the analytic corrections issmall, 0.01 and 0.04 for the two values of L min . The results from fits of ¯ U ′ , which are shownin Fig. 12, favor a smaller value of ν , although in substantial agreement. Indeed, fits with L min = 10 give ν = 2 . R ξ = 0 .
654 and ν = 2 . R ξ = 0 . ν = 2 . R ′ ξ , approximately 0.07. As our final result we take ν = 2 . , (61)21 L min ν wsc, R ξ ’wsc, U ’ ε=1.0, R ξ ’ ε=1.0, U ’ FIG. 13: (Color online) Estimates of the exponent ν from the analysis of ¯ R ′ ξ and ¯ U ′ , at R ξ = 0 . ± J Ising model at p = 0 .
7, obtained by fitting the data without scaling corrections (wsc)and with scaling corrections, Eq. (60). We only show the results of the fits with χ / DOF . .
5, upto values of L min where the errors blow up. The dotted lines correspond to ν = 2 . which is consistent with the results obtained for ¯ U ′ and ¯ R ′ ξ .Results at R ξ = 0 .
63 for the ± J model at p = 0 . p b = 0 .
45, shown inFigs. 13 and 14, respectively, are fully consistent with the estimate (61). For the ± J model at p = 0 .
7, fits with ε = 1 . L min = 8 of ¯ R ′ ξ and ¯ U ′ at fixed R ξ = 0 .
63 give ν = 2 . ν = 2 . p b = 0 .
45 and again for R ξ = 0 . R ′ ξ to (60) assuming ε = 2 . L min = 7 they give ν = 2 . U ′ with ε = 2 . ν = 2 . L min = 7. This slight discrepancy is probably due to the fact that scaling correctionswith exponent ω = 1 are small in the BDBIM at p b = 0 .
45 but not completely negligible:thus, the fits assume that corrections vanish faster than they really do, obtaining a slightlyincorrect asymptotic estimate.
C. The critical exponent η We estimate the exponent η by analyzing the susceptibility χ at fixed R ξ . Its criticalbehavior is discussed in Sec. III: at fixed R ξ , it behaves as¯ χ ( L, R ξ ) = aL − η (cid:2) a a ( R ξ − R ∗ ξ ) L − /ν + · · · + a ω L − ω + . . . + a b L − η + · · · (cid:3) , (62)where the O ( L − η ) term is the background contribution, cf. Eq. (29). Since R ξ is not exactlyequal to R ∗ ξ , we must take into account the O ( L − /ν ) corrections. As for ν , the systematicerror they induce is estimated by comparing the results at R ξ = 0 .
630 and R ξ = 0 . R ∗ ξ = 0 . L min ν wsc, R ξ ’wsc, U ’ ε=2.0, R ξ ’ ε=2.0, U ’ FIG. 14: (Color online) Estimates of the exponent ν from the analysis of ¯ R ′ ξ and ¯ U ′ , at R ξ = 0 . p b = 0 .
45, obtained by fitting the data without scaling corrections (wsc) and withscaling corrections, Eq. (60). We only show the results of the fits with reasonable values of χ / DOF,up to values of L min where the errors blow up. The dotted lines correspond to ν = 2 . L min -0.42-0.41-0.40-0.39-0.38-0.37-0.36 η R ξ =0.654, wsc R ξ =0.654, ε=1.0 R ξ =0.654, ε=2.4 R ξ =0.63, wsc R ξ =0.63, ε=1.0 R ξ =0.63, ε=2.4 FIG. 15: (Color online) Estimates of the exponent η from the analysis of ¯ χ ( L, R ξ ) for R ξ = 0 . R ξ = 0 . ± J Ising model at p = 0 .
5. Only results of fits with χ / DOF < . η = − . corrections, i.e. to ln χ = c + (2 − η )ln L, (63)ln χ = c + (2 − η )ln L + c ε L − ε , (64)23ith ε = 1 (the leading scaling correction arising from irrelevant scaling fields) and ε =2 . ≈ − η (to check the relevance of the background term, which might be large for smalllattices) at R ξ = 0 .
63 and R ξ = 0 . η for the ± J model at p = 0 .
5. As already mentioned, the systematic error due to the neglected L − /ν correctionsis estimated from the difference of the estimates at R ξ = 0 .
63 and R ξ = 0 . ν , we obtain η = − . { } , where the first erroris the statistical one while the second error takes into account the systematic error due tothe residual O ( L − /ν ) corrections and corresponds to the uncertainty on our estimate of R ∗ ξ . Slightly smaller but perfectly consistent results are obtained in the analyses of the datafor the other models. For example, in the case of the ± J Ising model at p = 0 .
7, a fit ofthe data at R ξ = 0 .
63 to Eq. (64) with ε = 1 . η = − . L ≥ L min = 10( χ / DOF ≈ . p b = 0 .
45, a fit of the data at R ξ = 0 .
63 toEq. (64) with ε = 2 . η = − . L ≥ L min = 9 ( χ / DOF ≈ . η = − . . (65)We finally discuss the behavior of the derivative χ ′ ≡ dχ/dβ of the susceptibility, which insome numerical works, see, e.g. Ref. 30, has led to contradictory results for the exponent ν .We show in the following that this discrepancy is essentially due to the analytic O ( L − /ν )corrections discussed in Sec. III, cf. Eq. (36). At fixed R ξ , ¯ χ ′ is expected to behave as¯ χ ′ ( L, R ξ ) = bL σ (cid:2) b a L − /ν + · · · + b ω L − ω + . . . + b b L − σ + · · · (cid:3) , (66)where σ ≡ /ν + 2 − η . Using our best estimates of ν and η , we obtain σ = 2 . χ and the derivative of the phenomenological quantities, the amplitude b a of the O ( L − /ν ) corrections does not vanish at T c and thus at R ξ = R ∗ ξ ; therefore, it is not expectedto be small when R ξ ≈ R ∗ ξ . The O ( L − σ ) term in Eq. (66) is the background contribution,cf. Eq. (32). In Fig. 16 we show the estimates of σ obtained in fits of ln χ ′ toln χ ′ = c + σ ln L + c ε L − ε , ln χ ′ = c + σ ln L + c L − ε + c L − ε , (67)for several values of ε, ε , ε . The results are consistent with the expected value σ ≈ . O ( L − /ν ) corrections are taken into account. VII. HIGH-TEMPERATURE ESTIMATES
In the parallel-tempering simulations we have collected a large amount of data for severalvalues of β in the high-temperature phase. They have not been used in the analyses wehave presented in the previous Sections. Here, we shall present analyses that consider allhigh-temperature results. They fully confirm the critical-point estimates discussed aboveand provide additional evidence that all models belong to the same universality class.24 L min σ ε=0.4ε=1.0ε =0.4, ε =1.0ε =0.4, ε =2.8 FIG. 16: (Color online) Estimates of the critical exponent σ ≡ /ν +2 − η , as obtained by analyzing χ ′ at R ξ = 0 . σ = 2 . (β c - β) L ξ/ L L=40L=32L=28L=24L=20L=18L=16L=14L=13L=12L=11L=10L=9L=8
FIG. 17: (Color online) Plot of ξ/L vs ( β c − β ) L /ν with β c = 0 .
902 and ν = 2 .
45. Data for the ± J model with p = 0 . A. Finite-size scaling behavior of ξ and estimates of ν We determine ν from the FSS behavior of the correlation length. The starting point isthe FSS equation ξ ( β, L ) L = f ( u t L /ν ) + v ω ( β ) L ω f ω ( u t L /ν ) + · · · , (68)where f ( x ) and f ω ( x ) are universal (once the normalization of the scaling fields has beenfixed) and have a regular expansion in powers of x . In order to ensure a finite infinite-volume25imit, for x → ∞ they behave as f ( x ) ∼ x − ν , f ω ( x ) ∼ x ν ( ω − . (69)The scaling field u t is an analytic function of β that vanishes at the critical point. Hence, ithas an expansion u t = β c − β + b ( β c − β ) + O [( β c − β ) ] . (70)In Fig. 17 we report ξ ( β, L ) /L versus ( β c − β ) L /ν using the data for the ± J model at p = 0 . β c and ν obtained in the previous Sections: β c = 0 .
902 and ν = 2 . u t on β .In order to obtain quantitative estimates of ν , we follow Ref. 71 and perform three differentfits of ξ ( β, L ) /L . In the first fit we neglect the nonanalytic scaling corrections, set u t = β c − β ,and fit the data to (fit a) ξ ( β, L ) L = P n ( x ) − ν/n , x ≡ ( β c − β ) L /ν , (71)where P n ( x ) is a polynomial of degree n . We have chosen polynomials for computationalconvenience, but any other complete set of functions can be used. Note that the specificchoice (71) of fitting function satisfies the large- x behavior (69). In this fit, the free parame-ters are the ( n + 1) coefficients of the polynomial P n ( x ), the critical inverse temperature β c ,and ν . The critical-point value R ∗ ξ is equal to P n ( x = 0) − ν/n .In order to take into account the analytic corrections (fit b), we slightly modify the Ansatz(71), setting x ≡ β c − β + b ( β c − β ) and taking b as an additional free parameter. Finally,we consider the nonanalytic scaling corrections. We use the Ansatz (fit c) ξ ( β, L ) L = (cid:20) P n ( x ) + 1 L ω (1 + ax ) ων Q n ( x ) (cid:21) − ν/n , x ≡ ( β c − β ) L /ν , (72)where P n ( x ) and Q n ( x ) are both polynomials of degree n , and a is a free parameter. We cancheck that Ansatz (72) has the correct infinite-volume limit. For L → ∞ at fixed β , i.e., for x → ∞ , we have P n ( x ) ≈ p n x n , Q n ( x ) ≈ q n x n , and ξ ( β, L ) ≈ L (cid:20) p n x n + 1 L ω ( ax ) ων q n x n (cid:21) − ν/n ≈ Lx − ν [ p n + a ων q n ( β c − β ) ων ] − ν/n ∼ ( β c − β ) − ν (cid:2) k ( β c − β ) ∆ (cid:3) (73)where ∆ ≡ ων . In these fits n must be chosen so that P n provides an accurate parametriza-tion of the scaling function. We have tried several values of n , until the χ of the fit doesnot change significantly as n increases by 1. In practice, we have taken n between 6 and 10.In fit (72) corrections to scaling are parametrized by a polynomial that has the same degreeas that parametrizing the leading behavior. Our data are not so precise and corrections arenot so large to require such a large number of parameters. To reduce the number of freeparameters we have taken n = 6 , Q n ( x ) = q + q x + · · · + q n x n . (74)26 .55 0.60 0.65 0.70 0.75 β min ν L min = min = min = min = FIG. 18: (Color online) Estimates of ν for the ± J model at p = 0 .
5. We report estimates obtainedby fitting ξ/L to the Ansatz (71) for several values of L min and β min . The dotted lines correspondto the estimate ν = 2 . This choice, which is completely arbitrary, significantly reduces the number of free parame-ters, but still allows us to parametrize accurately (at the level of the statistical errors) thescaling corrections.In order to take into account the additional scaling corrections which are not taken intoaccount by our fit Ans¨atze, we have repeated each fit several times; each time we only includedata such that β ≥ β min and L ≥ L min . In Fig. 18 we report the estimates of ν obtained infits of ξ/L to Eq. (71) for the ± J model at p = 0 .
5. Corrections to scaling are clearly visiblefor small values of β min and L min . Indeed, at fixed L min the estimates decrease, becomingapproximately independent of β min for β min & .
65. The dependence on L min is instead tiny,and all results with L min ≥
10 are consistent within errors.The results of fits that also take into account the analytic and the nonanalytic corrections(with ω = 1) are reported in Fig. 19. As far as ν is concerned, no significant differencesare observed and in all cases the results become independent of β min for β min & .
65. Allresults (with their statistical errors) lie in the interval 2 . ≤ ν ≤ .
5, and are therefore inperfect agreement with the estimate ν = 2 . β c . Indeed, while analyses without nonanalytic scalingcorrections give estimates of β c that cluster around 0.895, those that take the correctionsinto account give values which are somewhat larger. In any case, all results are consistentwith the estimate β c = 0 . R ∗ ξ . Analyses without nonanalytic scaling corrections give R ∗ ξ = 0 . R ∗ ξ = 0 . R ∗ ξ = 0 . b that appears in theexpansion (70) of u t : they vary somewhat with β min and give approximately 0 . b . . ν and determine the FSS curves also for the27 .55 0.60 0.65 0.70 0.75 β min ν L min = fit b L min = fit b L min = fit c L min = fit d L min = fit a β min β c L min = fit b L min = fit b L min = fit c L min = fit d L min = fit a FIG. 19: (Color online) Estimates of ν (top) and β c (bottom) for the ± J model at p = 0 .
5. Wereport estimates obtained by fitting ξ/L to Eq. (71) (fit a), to Eq. (71) with x ≡ [( β c − β ) + b ( β c − β ) ] L /ν (fit b), to Eq. (72) with ω = 1 (fit c), and to the extended-scaling Ansatz (86) (fit d). Thedotted lines correspond to the estimates ν = 2 . β c = 0 . ± J model at p = 0 . p b = 0 .
45. We report here only estimates of ν obtained by fitting the data to Eq. (71), since we do not have data precise enough to allowfor a detailed study of the scaling corrections. In any case, the results for p = 0 . ν . For the ± J modelat p = 0 . ν = 2 . ν = 2 . β min = 0 .
59 and L min = 7 [ L min = 9]; ν = 2 . ν = 2 . β min = 0 .
68 and the same values of L min . For the BDBIMat p b = 0 .
45 we obtain: ν = 2 . ν = 2 . β min = 0 .
82 and again L min = 7[ L min = 9]; ν = 2 . ν = 2 . β min = 1 .
02. These results are in very goodagreement with the estimate ν = 2 . .0 0.2 0.4 0.6 0.8 1.0 (β c - β) L ξ/ L L=20L=16L=14L=11L=10L=9L=8L=7L=6L=5L=4 (β c - β) L ξ/ L L=16L=14L=11L=10L=9L=8L=7L=6L=5L=4
FIG. 20: (Color online) Plot of ξ/L vs ( β c − β ) L /ν for the ± J model at p = 0 . p b = 0 .
45 (bottom). We use ν = 2 . β c = 0 .
87 (top), and β c = 1 .
54 (bottom). Thedashed line is the curve R ξ ( cx ), where R ξ ( x ) is the function (76), which is estimated in fits of ξ/L for the ± J model at p = 0 .
5, and c is a model-dependent constant: c ≈ .
026 for the ± J model at p = 0 . c ≈ . p b = 0 . In order to verify the universality of the FSS curves we first fitted the data for the ± J model at p = 0 . ν = 2 . β c = 0 . L ≥ β ≥ .
62, we obtain ξL = R ξ ( x ) x ≡ ( β c − β ) L /ν , (75)with (this expression is valid in the interval of values of x for which we have data, i.e., for0 ≤ x . . R ξ ( x ) = (cid:0) . . x − . x + 1926 . x .0 0.2 0.4 0.6 ξ/ L U p=0.5, L=24p=0.7, L=20p=0.7, L=16p=0.5, p b =0.45, L=16p=0.5, p b =0.45, L=14 FIG. 21: (Color online) Plot of U vs ξ/L for the ± J model at p = 0 . p b = 0 . − . x + 88711 . x − . x + 446776 . x − . x + 243001 . x − . x (cid:1) − . . (76)The function R ξ ( x ) is universal apart from a rescaling of its argument. Thus, if we plot ξ/L vs x ≡ ( β c − β ) L /ν in any model, the data should fall on the curve R ξ ( cx ), where c is amodel-dependent constant. In Fig. 20 we report the data for the ± J model at p = 0 . p b = 0 .
45. The results show very good scaling and fall on top of the curvecomputed from the data of the bimodal model at p = 0 .
5, confirming universality.We finally consider the cumulant U . In the critical limit U should be a universal functionof ξ/L , independent of the model. Corrections scale as L − ω h ( ξ/L ), where h ( x ) is alsouniversal, apart from a multiplicative constant. The numerical estimates of U and ξ/L arereported in Fig. 21. We consider the ± J model at p = 0 . p = 0 . p d = 0 .
45. We only consider the largest lattices, so that nonanalytic scaling corrections arenot visible on the scale of the figure (a detailed study of the L − ω corrections for ξ/L = 0 . B. Finite-size scaling of the susceptibility and estimates of η We now turn to the determination of the critical exponent η . The starting point isEq. (29), which we rewrite as χ ( β, L ) = L − η ¯ u h ( β ) g ( u t L /ν ) (cid:20) v ω ( β ) L ω g ω ( u t L /ν ) + . . . (cid:21) . (77)This equation is not very convenient since it involves u t , hence the critical temperature, andthe exponent ν . To reduce the number of unknown parameters, we note that Eq. (68) can30 .0 0.1 0.2 0.3 0.4 0.5 0.6 ξ/ L χ ξ − + η L=40L=32L=28L=24L=20L=18L=16L=14L=13L=12L=11L=10L=9L=8 ξ/ L u h − χ ξ − + η L=40L=32L=28L=24L=20L=18L=16L=14L=13L=12L=11L=10L=9L=8 | FIG. 22: (Color online) Plot of χξ η − (upper panel) and of χξ η − ˜ u − h (lower panel) vs ξ/L with η = − . C ( ξ/L ), as estimated in fitsof χ to Eq. (80). Data for the ± J model at p = 0 . be inverted to give u t L /ν = F ( ξ/L ) + v ω ( β ) L ω F ω ( ξ/L ) + . . . (78)Inserting in Eq. (77) we obtain the scaling form χ ( β, L ) = ξ − η ¯ u h ( β ) C ( ξ/L ) (cid:2) v ω ( β ) ξ − ω C ω ( ξ/L ) + . . . (cid:3) . (79)In Eq. (79) we have singled out ξ − η and ξ − ω instead of L − η and L − ω . With this choice C ( x ) and C ω ( x ) are regular for x →
0. Note the presence of the function ¯ u h ( β ). In theFSS limit ξ → ∞ , L → ∞ , at fixed ξ/L , we always have β → β c , so that asymptotically itshould be possible to replace ¯ u h ( β ) with the constant ¯ u h ( β c ). Therefore, this function givesrise to scaling corrections, that we have named analytic corrections in the previous sections.31 .55 0.60 0.65 0.70 0.75 β min -0.44-0.42-0.40-0.38-0.36 η L min = fit a L min = fit a L min = fit a L min = fit b L min = fit b FIG. 23: (Color online) Estimates of η for the ± J model at p = 0 .
5. We report estimates obtainedby fitting χ to the Ansatz (80) (fit a) and (81) (fit b) for several values of L min and β min . Thedotted lines correspond to the estimate η = − . In order to understand their relevance for our data, in Fig. 22 (upper panel) we plot ξ η − χ versus ξ/L for the ± J model at p = 0 .
5. It is evident that the data do not fall onto asingle curve: the scaling-field term ¯ u h ( β ) varies significantly with β and therefore cannot beneglected.The previous discussion indicates that, in order to estimate accurately the exponent η , itis essential to include the analytic corrections in the fitting function. We perform two fits.In the first one (fit a) we neglect the nonanalytic scaling corrections and considerln χξ = − η ln ξ + P n ( ξ/L ) + Q m ( β ) , (80)where P n ( x ) and Q m ( x ) are polynomials of degree n and m , respectively. Moreover, werequire Q m (0) = 0 in order to avoid the presence of two constant terms. As before, n and m are varied till the quality of the fit does not change significantly by varying the parametersby 1. In practice, we take 6 ≤ n, m ≤
10. To include the scaling corrections, we also considerln χξ = − η ln ξ + P n ( ξ/L ) + Q m ( β ) + ξ − ω S p ( ξ/L ) , (81)where S p ( x ) is a polynomial of degree p : we take p ≤
3. Note that here, as we already didin the analysis of ξ/L , we neglect the β dependence of the scaling field v ω .The results of the fits for the ± J model at p = 0 . β min and L min .For instance, for β min = 0 .
55 we obtain η = − . L min = 8), η = − . L min = 14), while for β min = 0 .
75 we have η = − . L min = 8), η = − . L min = 14). Fits with nonanalytic scaling corrections are less stable. We observe significantfluctuations that indicate that the data are not precise enough to be sensitive to this type of32 .55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 β u h β min = min =8 β min = min =14 β min = min =8 β min = min =14( β/β c ) η/2 _ FIG. 24: (Color online) Scaling-field function ˜ u h ( β ) for the ± J model at p = 0 .
5, as determinedin the fits of χ to the Ansatz (80), for different values of β min and L min . It is normalized suchthat ˜ u h ( β c ) = 1. We also report the approximation ˜ u h ( β ) ≈ ( β/β c ) − η/ , which is used in theextended-scaling scheme (see the discussion in Sec. VII C). scaling corrections. This is consistent with what is observed at the critical point: While theresults for η depend strongly on the chosen value for R ξ,f , indicating that the analytic scalingcorrections are important, essentially no dependence is observed on the nonanalytic ones;see, e.g., Fig. 15. In any case, all results are consistent with the estimate η = − . u h ( β ) that appears in Eqs. (77) and (79). InFig. 24 we plot ˜ u h ( β ) = ¯ u h ( β ) / ¯ u h ( β c ) [˜ u h ( β ) is normalized so that ˜ u h ( β c ) = 1] as obtainedin the different fits. The results corresponding to different values of L min and β min agreenicely, supporting the scaling Ansatz (79). A simple expression which reproduces the resultsreported in Fig. 24 is˜ u h ( β ) = 1 + 0 . − β/ . . − β/ . , (82)which is valid for 0 . ≤ β ≤ . u h ( β ) has been determined, we can compute thescaling function ˜ C ( ξ/L ) = ¯ u h ( β c ) C ( ξ/L ) by considering L η − χ ˜ u h ( β ) − . Such a combina-tion is shown in Fig. 22 (lower panel): all points fall on top of each other, confirming thevalidity of the FSS Ansatz. Moreover, as expected, we find that ˜ C (0) is finite and ˜ C ( ξ/L ) isapproximately constant for ξ/L . .
15, two properties which are not obvious from the upperpanel of Fig. 22. These conclusions are consistent with the FSS results for χ (2 L, β ) /χ ( L, β )(in this quantity the analytic function ¯ u h ( β ) cancels out) reported in Refs. 20,28, whichshow that this ratio has a tiny dependence on ξ/L up to ξ/L ≈ .
15 (for ξ/L = 0 .
15 wehave χ (2 L, β ) /χ ( L, β ) ≈ . L → ∞ . Since the rate of convergence is very slow (at fixed ξ/L data converge as L − /ν ), it is clear that such an asymptotic behavior can only be observedon enormously large lattices. Thus, in order to estimate η , it is crucial to take the function33 u h ( β ) into account.The function ˜ C ( x ) is universal, apart from a model-dependent multiplicative constant.We write it as ˜ C ( x ) = b Γ( x ) Γ(0) = 1 , (83)where Γ( x ) is universal. A fit of the data reported in Fig. 22 givesΓ( x ) = 1 + 5 . y − . y + 1516 . y − . y + 50105 . y − . y ,y ≡ exp( − /x ) , (84)and b ≈ . x ≡ ξ/L . R ∗ ξ ≈ . ± J model at p = 0 . p b = 0 .
45. Here we only present results corresponding to fits to the Ansatz (80). Our dataare not precise enough to allow us to perform fits which include the nonanalytic scalingcorrections. The results are consistent with the estimate η = − . ± J modelat p = 0 .
7, we obtain η = − . β min = 0 .
59) and η = − . β min = 0 .
68) for L min = 7, and η = − . β min = 0 .
59) and η = − . β min = 0 .
68) for L min = 9.For the BDBIM, we obtain η = − . − . L min = 7 and 9 and any β min in the range [0 . , . χξ η − ˜ u − h . In Fig. 25 we report this quantity for the ± J model at p = 0 . p b = 0 .
45 (middle panel), using in both cases η = − .
375 and our bestestimate of ˜ u h . As expected all points fall onto a single curve. If universality holds, thesecurves should be parametrized as b Γ( x ), where Γ( x ) is given in Eq. (84) and b is a model-dependent constant. In the case of the ± J model (upper panel of Fig. 25) we observe a verygood agreement, while for the BDBIM at p b = 0 .
45 (middle panel) some discrepancies occurfor ξ/L & .
4. There are two reasons for them. First, the function ˜ u h is not precisely knownfor β ≈ β c : in this range of values of β it varies slightly (10%) with β min and L min . Second,the plot depends on η . If we use η = − . u h , we obtain the lower panel of Fig. 25. Discrepancies are nowsignificantly reduced.It is worth noting that the functions ¯ u h are approximately the same in the three modelswe study, if one considers them as a function of the reduced temperature t ≡ − β/β c .For 0 . t . . t which is probed by our simulations—the ratio¯ u h, model1 ( t ) / ¯ u h, model2 ( t ) is constant, within our precision, for any pair of models. This resultis somewhat unexpected within RG theory, because these functions are not universal.Finally, let us comment on the FSS approach of Ref. 73, applied to spin-glass systems inRefs. 20,28. In this approach one considers the ratio χ (2 L, β ) /χ ( L, β ). This choice has asignificant advantage: the scaling-field function ¯ u h ( β ) cancels out, so that the leading scalingcorrections are the nonanalytic ones. As we have shown here, they are quite small, so thatvery good scaling is observed and reliable infinite-volume estimates are obtained. Analyticscaling corrections come in again when considering the critical limit of the infinite-volumeresults χ ∞ ( β ). Indeed, since ∆ ≡ ων ≈ .
45, for β → β c the analytic corrections dominate: χ ∞ ( β ) = ( β c − β ) − γ [ b + b ( β c − β ) + b ( β c − β ) + b ∆ ( β c − β ) ∆ + · · · ] . (85)34 .0 0.1 0.2 0.3 0.4 0.5 0.6 ξ/ L u h − χ ξ − + η L=20L=16L=14L=11L=10L=9L=8L=7 | ξ/ L u h − χ ξ − + η L=16L=14L=12L=11L=10L=9L=8L=7 | ξ/ L u h − χ ξ − + η L=16L=14L=12L=11L=10L=9L=8L=7 | FIG. 25: (Color online) Plot of χξ η − ˜ u − h vs ξ/L for the ± J model at p = 0 . η = − . p b = 0 .
45 with η = − .
375 (middle panel) and η = − . b Γ( ξ/L ), where Γ( ξ/L ) is defined in Eq. (84) and b is a model-dependent constant. We use b = 1 . ξ/L ≤ . .0 0.1 0.2 0.3 0.4 0.5 0.6 ξ/ L ( β / β c ) − η χ ξ − + η L=40L=32L=28L=24L=20L=18L=16L=14L=13L=12L=11L=10L=9L=8 β min -0.48-0.46-0.44-0.42-0.40-0.38-0.36 η L min = fit c L min = fit c L min = fit d L min = fit d L min = fit d L min = fit e FIG. 26: (Color online) Extended-scaling results for the ± J model at p = 0 .
5. In the upper panelwe report β − η χξ η − vs ξ/L . The dashed line is the universal curve ˜ C ( ξ/L ), as estimated in fits of χ to Eq. (80). In the lower panel we report estimates of η obtained in three different fits for severalvalues of L min and β min . Fit c uses the Ansatz (89), fit d and e are defined in text, below Eq. (89).The lines correspond to the estimate η = − . C. Extended-scaling scheme
In this Section we consider the extended-scaling scheme introduced in Ref. 29. It consistsin a particular choice of scaling variables, which, according to Ref. 29, should somehowdecrease scaling corrections and thus allow a faster convergence to the critical limit. Let usconsider first ξ/L . In this scheme the appropriate fit Ansatz is ξ ( β, L ) L = P n ( x ) − ν/n , x ≡ ( β c − β )( L/β ) /ν , (86)36here P n ( x ) is a polynomial of degree n . The results for the ± J model at p = 0 . β c and ν , wehave ( β c − β ) β − /ν = 1 . β c − β )[1 − . β c − β ) + . . . ] . (87)Thus, fit (86) is essentially equivalent to a fit with analytic corrections (fit b) with b = − . b is small—hence, this change of the scaling variable does not have muchinfluence on the final results—and is close to what we obtain numerically, though not fullyconsistent (we predict 0 . b . . x ≡ ( β c − β ) L /ν , whichshould be the natural variable in spin-glass systems, given the symmetry under β → − β . These fits are significantly worse than the previous ones for β . .
70. For larger values, nosignificant differences are observed. These results can be understood by noting ( β c − β ) =1 . β c − β )[1 − . β c − β ) + . . . ]. Thus, this choice of scaling variable corresponds toassuming b = − .
55 in Eq. (70), which is significantly larger than what we find numerically.Therefore, if we use ( β c − β ) as approximate thermal scaling field, the analytic corrections—in this case they are proportional to ( β c − β ) —are more important than in the case in which u t is simply approximated with β c − β .The extended-scaling scheme can also be applied to the analysis of the susceptibility. Itamounts to consider the scaling Ansatz χ ( β, L ) = β η − ξ − η C ( ξ/L ) . (88)In Fig. 26 we show β − η χξ η − versus ξ/L for the ± J model at p = 0 .
5. Scaling is betterthan that observed in the upper panel of Fig 22: the scatter of the data points is signifi-cantly reduced, indicating that β η − approximates the scaling-field term ¯ u h better than aconstant. However, the rescaled data are still far from the asymptotic curve (the dashedline) determined numerically above, indicating that the residual analytic scaling correctionsare also in this case not negligible. This is better understood, by comparing the function ˜ u h ,as determined in the fits, with the approximation ( β/β c ) η/ − , which follows from Eq. (88).As can be seen from Fig. 24, the approximate expression proposed in Ref. 29 has the correctqualitative shape, but differs significantly from the quantitative point of view. For thesereasons, we do not expect the scaling Ansatz (88) to be particularly useful to estimate η from our data.To understand the role of the residual analytic scaling corrections on the determinationsof η , we fit the data for the ± J model at p = 0 . χβ ξ = − η ln ξβ + P n ( ξ/L ) , (89)where P n ( x ) is a polynomial of degree n . The results are reported in Fig. 26. They varystrongly with β min and L min , indicating that scaling corrections are sizable and not negligible.As a test we have also repeated the fits to Eqs. (80) and (81) replacing ln χ/ξ with ln χβ /ξ and ln ξ with ln ξ/β (we call fit d and fit e the fits corresponding to Eqs. (80) and (81),respectively). As expected, the results are identical to those obtained in fit b and fit c,respectively. Indeed, the fits only differ in the parametrization of the analytic function¯ u h . Note that the extended-scaling approximation for ¯ u h is not analytic in β , since β =0 is a branching point. This is, however, irrelevant in practice, since we are looking forapproximations of ¯ u h in the interval 0 . . β . .
9, which is quite far from β = 0.37 III. THE ZERO-MOMENTUM QUARTIC COUPLINGS G AND G IN THEHIGH-TEMPERATURE PHASE
In this section we consider the zero-momentum quartic couplings G and G defined inEq. (19) and (20) and estimate their infinite-volume critical value defined by G ∗ = lim β → β + c lim L →∞ G ( L, β ) . (90)We consider the ± J Ising model at p = 0 . β = 0 .
55 and L = 12, 14, 16, 18, and 32, and for β = 0 .
625 and L = 16, 18, 28, 32, and 40.We combine results obtained in the random-exchange runs that we performed for our FSSstudy around β c and results obtained in standard MC simulations for these two values of β . First, we investigate the infinite-volume limit. The correlation length converges rapidly(see Fig. 27). For instance, for β = 0 .
55 we obtain ξ = 1 . , . , . L = 18 , ,
32, while for β = 0 .
625 we have ξ = 3 . , . L = 32 and L =40: for L/ξ &
10 the results vary by less than 0.1%, indicating that the difference from theirthermodynamic limit is within 0.1%. Thus, we can take as infinite-volume estimates thoseobtained on the largest lattices. The quartic couplings show larger finite-size corrections. Asshown in Figs. 27 the infinite-volume limit is approximately reached for
L/ξ &
12, withinour statistical precision. Indeed at β = 0 .
625 we find G = 90 . , . , . G = − . , − . , − . L = 18 , ,
32, respectively. We can thustake the estimates on the largest lattices as infinite-volume estimates.In the critical limit we expect G ( L = ∞ , β ) = G ∗ + c ξ − ω , (91)with ω = 1 . β = 0 .
55 and β = 0 .
625 leads us to theestimates G ∗ = 90 . , G ∗ = − , (92)where the error takes also into account the effects of the O ( ξ − ω ) scaling corrections, whichis roughly estimated from the difference of the infinite-volume results at the two values of β . These results can be compared with those obtained in Ref. 27 for the random-anisotropyHeisenberg model in the limit of infinite anisotropy: G ∗ = 88(8) and G ∗ = − IX. CONCLUSIONS
In this paper we discuss the critical behavior of three-dimensional Ising spin-glass systemswith the purpose of verifying universality, clarifying the role of scaling corrections, anddetermining the critical exponents. More precisely, our results can be summarized as follows.i) By using the RG we derive the behavior for L → ∞ and β → β c of several quantitieswhich are routinely measured in MC simulations. In particular, we show that theanalytic dependence of the scaling fields on the model parameters may give rise tocorrections which behave as L − /ν ∼ L − . . If they are neglected, FSS analyses give38
10 15 20 L/ ξ ξ β=0.55β=0.625 L/ ξ G β=0.55β=0.625 L/ ξ -21-20-19-18-17-16-15 G β=0.55β=0.625 FIG. 27: (Color online) Estimates of ξ ( L, β ) (top), G ( L, β ) (middle), and G ( L, β ) (bottom)versus
L/ξ for β = 0 .
55 and β = 0 . ξ ∞ ≈ .
79 and ξ ∞ ≈ .
27. Results for the ± J model at p = 0 .
5. The dotted lines correspond tothe estimates G ∗ = 90 . G ∗ = − ν is typically large.ii) We determine the leading nonanalytic correction-to-scaling exponent ω . We obtain ω = 1 . ω ≈ .
8, see Ref. 47. The exponent ω is alsosignificantly larger than that at the PF transition, which occurs for small frustration: ω = 0 . iii) We accurately verify universality. A careful analysis of the quartic cumulants U and U at fixed R ξ = 0 .
63 shows that their limit for L → ∞ is independent of the modeland of the disorder distribution. The results obtained in the different models differat most by approximately one per mille in the case of ¯ U and one per cent for ¯ U .They support the existence of a unique Ising spin-glass universality class. Universalityis also supported by the FSS analyses of ξ and χ in the high-temperature phase. Weverify that the FSS curves for these two quantities are independent of the model.iv) We determine the critical exponents. For this purpose we perform analyses at thecritical point and analyses which take into account all high-temperature data. Resultsare consistent, once the analytic and the nonanalytic corrections are taken into account.Moreover, they do not depend on the model and disorder distribution. Again, thissupports the universality of the paramagnetic-glassy transition. We obtain ν = 2 . , η = − . . (93)Using scaling and hyperscaling relations, we find β = ν (1 + η ) / . γ =(2 − η ) ν = 5 . α = 2 − ν = − . ν are reported in the introduction and in Table I of Ref. 30. Some of themost recent ones are close to our final estimate. For the exponent η , we quote here the mostrecent results: η = − . − . η = − . and η = − . Theyare all in substantial agreement with our result, which is however significantly more precise.We also mention the estimate β = 0 . ν = 2 . η = − . L = 20 ,
24 for the ± J model at p = 0 .
5. The new resultshave slightly shifted the estimates of β c , R ∗ ξ , and of the critical exponents. Second, we havebeen more conservative. With our present error bars, the estimates are fully consistent withthe results of all analyses for the three models we considered.We also analyzed our data by using the extended-scaling scheme proposed in Ref. 29. Thisapproach might partly take into account the scaling corrections arising from the analytic de-pendence of the scaling fields on the reduced temperature but it neglects the nonanalyticcorrections arising from the irrelevant operators. In some cases, for instance for the over-lap susceptibility, this scheme shows an apparent improvement of the scaling behavior withrespect to the naive approach in which the analytic corrections are simply neglected. How-ever, the approximate expressions which follow from the extended-scaling scheme are not40ufficiently precise for a high-precision study of the critical behavior: if one aims at accurateestimates, it is necessary to determine the corrections directly from the data. APPENDIX A: FINITE-SIZE BEHAVIOR OF THE PHENOMENOLOGICALCOUPLINGS
We now provide a detailed proof of Eq. (30) for the phenomenological couplings U , U ,and R ξ ≡ ξ/L . First, we present a simple physical argument that clarifies the origin ofthe function ¯ u h ( t ) introduced in Eq. (27), and then, a more formal argument based on theWegner expansion. A short and physical proof goes as follows. First, note that χ is not RG invariant and, inparticular, it is not invariant under field redefinitions. Given a theory with fields φ , define χ φ = Z d x [ h φ ( x ) φ (0) i ] . (A1)If we change variables φ = Zψ , we obtain χ φ = Z Z d x [ h ψ ( x ) ψ (0) i ] = Z χ ψ . (A2)Thus, under a field redefinition χ → Z χ . The function ¯ u h ( t ) is exactly the field redefinitionthat relates the bare lattice fields to the renormalized fields. Hence, we expect at criticality χ = ¯ u h L − η . (A3)A t dependent prefactor is absent in the phenomenological couplings, because these quantitiesare RG invariant and in particular, they do not depend on the normalization of the fields.This is quite obvious from their definitions. In the continuum limit they can be written as U ≡ R dx dx dx dx [ h φ ( x ) φ ( x ) φ ( x ) φ ( x ) i ] (cid:2)R dx dx h φ ( x ) φ ( x ) i (cid:3) , (A4) U ≡ R dx dx dx dx [ h φ ( x ) φ ( x ) ih φ ( x ) φ ( x ) i ] (cid:2)R dx dx h φ ( x ) φ ( x ) i (cid:3) − , (A5) ξ ≡ R dx dx ( x − x ) [ h φ ( x ) φ ( x ) i ] R dx dx h φ ( x ) φ ( x ) i . (A6)The formula we report for ξ correspond to Eq. (14) for L → ∞ , disregarding corrections oforder L − . Since the number of fields in the denominator is equal to the number of fields inthe numerator in all expressions, each quantity is invariant under field redefinitions.This discussion clarifies the origin of the t -dependent prefactor in Eq. (A3) and is at thebasis of the statement, which is well accepted in the literature, that U ( T c ), U ( T c ), and ξ ( T c ) /L approach universal constants as L → ∞ . The presence of an analytic prefactorwould violate universality. Indeed, imagine that a phenomenological coupling R behaves as R ( β, L ) = ¯ u h ( t ) some power g ( u t L /ν ) + scaling corrections (A7)41lose to the critical point. Since ¯ u h ( t ) is model dependent, at the critical point R ( β c , L )would not be universal for L → ∞ . Hence, on purely physical grounds, there must be no t -dependent prefactor in any phenomenological coupling, hence no analytic corrections.This short discussion shows that R ξ and χ/L − η are conceptually two very different ob-jects. The first quantity is RG invariant and shows a universal FSS behavior. The secondquantity is not RG invariant and, for instance, χ ( T c ) /L − η converges to a model-dependent constant as L → ∞ .Eq. (30) can also be derived from the usual Wegner’s scaling expression for the free energy.Let us first consider U . Using Eq. (26), we find[ µ ] = L d ∂ F ∂h (cid:12)(cid:12)(cid:12)(cid:12) h =0 = L y h ¯ u h f (4) (0 , u t L y t ) + · · · (A8)[ µ ] = L d ∂ F ∂h (cid:12)(cid:12)(cid:12)(cid:12) h =0 = L y h ¯ u h f (2) (0 , u t L y t ) + · · · (A9)[the dots correspond to nonanalytic scaling corrections and bulk contributions, and thederivatives refer to the first variable appearing in the scaling function f ( x, y )] so that U = f (4) (0 , u t L y t ) f (2) (0 , u t L y t ) + · · · , (A10)which proves Eq. (30).To discuss U one should generalize Wegner’s scaling expression (see Sec. 3.1 of Ref. 46for a detailed discussion). Define Z ( β, h, L ) as the partition function of two systems atinverse temperature β defined on a lattice of size L coupled by an interaction h X x σ x σ x . (A11)Then, consider F ( β, h , h , L ) = L − d [ln Z ( β, h , L ) ln Z ( β, h , L )] . (A12)A scaling Ansatz like Eq. (25) allows one to obtain an expression analogous to that obtainedfor U and to prove Eq. (30) for U .In order to determine the scaling behavior of R ξ we consider a momentum-dependentmagnetic field. The argument goes as follows: Define Z ( β, h, L, p ) as the partition functionof two systems at inverse temperature β defined on a lattice of size L coupled by an interac-tion h P x σ x σ x cos( p · x ). Then, consider the corresponding disorder-averaged free-energydensity F ( β, h, L, p ) = L − d [ln Z ( β, h, L, p )] . (A13)Under RG transformations L → λL , momenta scale as p → p/λ , so that the singular partof the free-energy density scales as F sing ( β, h, L, p ) = L − d f ( pL, u h ( h, t, p ) L y h , u t ( h, t, p ) L y t ) , (A14)where we have neglected the nonanalytic scaling corrections and now the scaling fields dependalso on p . Taking derivatives with respect to h and then setting h = 0, we obtain for thetwo-point function [of course u h ( h, t, − p ) = u h ( h, t, p )] e G ( p ) = ¯ u h ( t, p ) L − η f (2) ( pL, , u t (0 , t, p ) L y t ) , (A15)42here we write as before u h ( h, t, p ) = h ¯ u h ( t, p ) + O ( h ) , (A16)and we have neglected subleading terms. For p →
0, because of the cubic symmetry of thelattice, we have ¯ u h ( t, p ) = ¯ u h ( t ) + O ( p ) , (A17) u t (0 , t, p ) = u t (0 , t ) + O ( p ) , (A18)where ¯ u h ( t ) and u t (0 , t ) are the usual (zero-momentum) scaling fields. Hence, for p → p , we can express e G ( p ) in terms of the scaling fields thatappear for p = 0: e G ( p ) = ¯ u h ( t ) L − η f (2) ( pL, , u t (0 , t ) L y t ) + O ( p ) + · · · (A19)In the definition (14) of the correlation length ξ we should consider p = q ∼ /L . Thus,disregarding terms of order L − we have e G (0) − e G ( p ) e G ( p ) = f (2) (0 , , u t (0 , t ) L y t ) f (2) ((2 π, , , , u t (0 , t ) L y t ) − u t (0 , t ) L y t ) . (A20)Neglecting again corrections of order 1 /L , we have ξ L = 14 π Φ( u t (0 , t ) L y t ) , (A21)which proves Eq. (30). If we consider the corrections to scaling, this derivation shows that R ξ behaves essentially as U and U . The only difference is the presence of correctionsdue to the momentum-dependence of the scaling fields and to the specific definition of thecorrelation length: they scale as L − , L − , . . . , L − ω − , . . . Since ω ≈
1, in Eq. (30) theyrepresent additional subleading corrections and can thus be neglected. This allows us toconsider R ξ , U , and U on the same footing. APPENDIX B: SOME TECHNICAL DETAILS ON THE MC SIMULATIONS
In our MC simulations we implement the standard Metropolis algorithm with a sequentialupdate of the spins. We use a multispin implementation, in which n bit = 64 systems aresimulated in parallel. For each of them we use a different set of bonds { J xy } .For the random numbers we use the SIMD-oriented Fast Mersenne Twister (SFMT) generator. In particular, we use the genrand res53() function that produces double-precisionoutput. Independent random numbers are employed to generate the starting configurationsfor each disorder realization and in the parallel-tempering updates. In the latter case veryfew random numbers are used and thus, it takes virtually no extra time to use individualrandom numbers for each of the n bit systems which are simulated in parallel. On the otherhand, in order to save CPU time, we use the same sequence of random numbers for thelocal Metropolis update of any of the n bit systems. Though this choice does not lead towrong estimates of the expectation values, it might create a statistical correlation amongthe n bit systems. However, since each of the n bit systems correspond to a different set of43ond couplings J xy we expect this effect to be negligible. Nevertheless, in order to ensure acorrect estimate of the statistical error, in our jackknife analysis we put all n bit systems thatuse the same sequence of random numbers in the same bin. In order to compute overlapobservables, we performed runs for two systems with the same set { J xy } in parallel. In ourMC simulations a single Metropolis update of a single spin takes about 1 . × − secondson an Opteron CPU running at 2 GHz (this should be compared with the speed of thededicated computer Janus, the fastest computer simulating discrete spin models, whichtakes 2 × − seconds to update an Ising spin).To reduce autocorrelations we used the random-exchange or parallel-tempering method. To this end, we divided the interval [ β min , β max ] into N β − β . The pa-rameter β max was chosen such that ξ ( β max ) /L ≈ .
63 in most cases; in the latest runs weconsidered larger values, such that ξ ( β max ) /L ≈ .
66. The parameter β min was chosen suchthat ξ ( β min ) ≪ L . We computed the observables in the neighborhood of β max by using theirsecond-order Taylor expansion around β max . The coefficients of the expansion were estimatedby measuring appropriate correlators with the energy.An elementary update unit consists in n met Metropolis sweeps followed by a replica ex-change of all pairs of systems at nearby temperatures. The different systems were sequentiallyvisited, starting from those at β min and β min + ∆ β . As a candidate for the exchange, weconsidered one of two replicas with equal probability. The acceptance probability for theexchange is min[1 , exp( − ∆ β ∆ H )]. Since the measurement of the energy in our implemen-tation costs more CPU time than a Metropolis sweep, we chose n met ≫ β . The computation of disorder averages of products of thermal expectations requires par-ticular care. Indeed, naive estimators show a bias which may become larger than statisticalerrors. To avoid the problem we consider essentially bias-free estimators, defined followingRef. 46. For this purpose we divide the measurement phase of the run into 12 intervals.Between each pair of subsequent interval there is a decorrelation phase. In total, the runconsists of the following phases: Eq , D M , D M , ..., D M . After some tests, we fixedthe number of update steps for each of them. The equilibration phase Eq corresponds to20 n temp elementary update units; the measurement phases M i correspond to n temp updateunits, while the length of D i is n temp for i = 7 and 5 n temp for i = 7. Recall that eachelementary update unit corresponds to n met Metropolis sweeps of all systems and to one fulltempering sweep.The presence of different measurement phases allows us to define bias-free quantities. Todefine [ h A ih B i ] (for instance, this is relevant for the computation of U ) we average overthe samples the quantity 12 × × X i =1 12 X j =7 [ µ ( A ) i µ ( B ) j + µ ( A ) j µ ( B ) i ] , (B1)where µ ( A ) i is the average of estimates of A obtained in the measurement phase M i . Analo-gously, to compute [ h A ih B ih C i ] (these correlators are necessary to compute the coefficientsof the Taylor expansions around a given value of β ), we average over the samples the quantity13! × X i =1 8 X j =5 12 X k =9 [ µ ( A ) i µ ( B ) j µ ( C ) k + 5 permutations] . (B2)44 ABLE IV: Summary of the parameters for the runs at p = 0 . N β is the number of β values usedin the parallel-tempering simulation. The rest of the notation is explained in the text. The CPUtime refers to a single core of a dual core Opteron CPU running at 2.4 GHz. L samples/64 n met n temp N β β min β max CPU time(days)4 100000 5 40 5 0.58 0.925 100000 5 50 5 0.58 0.9086 100000 5 50 5 0.58 0.9018 17 119103 5 80 5 0.58 0.899 38 100000 5 80 5 0.58 0.8975 49 110850 10 100 8 0.55 0.8962 2710 100681 10 150 8 0.55 0.896 5011 109779 10 300 10 0.54 0.8955 18312 106812 10 400 10 0.54 0.8955 30813 38282 10 600 10 0.54 0.8955 21014 31600 50 200 10 0.62 0.8955 36116 24331 10 1000 20 0.52 0.895 83020 1542 20 2000 32 0.5125 0.895 65820 2291 50 1500 20 0.625 0.91 114624 717 25 2500 32 0.5125 0.895 82624 1627 50 2000 20 0.625 0.91 187428 285 60 2500 20 0.6575 0.895 782
In order to check equilibration, and decorrelation for the bias correction, we followed thesuggestion of Ref. 30. We doubled the length of the run until the estimates of all observableswere consistent within error bars. We performed this check only for the observables at β max ,because these are expected to be the most difficult ones for equilibration and decorrelation.Starting from disordered configurations, we determined the number of update steps n half that are needed to reach (averaged over samples) half of the equilibrium value of the overlapsusceptibility. In total, the equilibration consisted of at least 100 n half update steps. Usingthese methods to check equilibration, we came up with the choices summarized in Table IV.The parameters are not highly tuned, since we had the CPU time available on short notice.The runs that were done later have typically a larger β min than those done earlier. The runfor L = 28 is a bit at the edge of the criterion given above for equilibration. However, giventhe rather small number of samples ( N s = 18240), we are quite confident that the estimatesare correct within the quoted error bars. S. F. Edwards and P. W. Anderson, J. Phys. F , 965 (1975). A. Ito, H. Aruga, E. Torikai, M. Kikuchi, Y. Syono, H. Takei, Phys. Rev. Lett. , 483 (1986). K. Gunnarsson, P. Svedlindh, P. Nordblad, L. Lundgren, H. Aruga, and A. Ito, Phys. Rev. B , 8199 (1991). S. Nair and A. K. Nigam, Phys. Rev. B , 214415 (2007). A. T. Ogielski and I. Morgenstern, Phys. Rev. Lett. , 928 (1985). A. T. Ogielski, Phys. Rev. B , 7384 (1985). W. L. McMillan, Phys. Rev. B , 340 (1985). A. J. Bray and M. A. Moore, Phys. Rev. B , 631 (1985). R. N. Bhatt and A. P. Young, Phys. Rev. Lett. , 924 (1985). R. R. P. Singh and S. Chakravarty, Phys. Rev. Lett. , 245 (1986). J. D. Reger and A. Zippelius, Phys. Rev. Lett. , 3225 (1986). R. N. Bhatt and A. P. Young, Phys. Rev. B , 5606 (1988). E. Marinari, G. Parisi, and F. Ritort, J. Phys. A , 2687 (1994). N. Kawashima and A. P. Young, Phys. Rev. B , R484 (1996). L. W. Bernardi, S. Prakash, and I. A. Campbell, Phys. Rev. Lett. , 2798 (1996). D. I˜niguez, G. Parisi, and J. J. Ruiz-Lorenzo, J. Phys. A , 4337 (1996). B. A. Berg and W. Janke, Phys. Rev. Lett. , 4771 (1998). E. Marinari, G. Parisi, and J. J. Ruiz-Lorenzo, Phys. Rev. B , 14852 (1998). P. O. Mari and I. A. Campbell, Phys. Rev. E , 2653 (1999). M. Palassini and S. Caracciolo, Phys. Rev. Lett. , 5128 (1999). H. G. Ballesteros, A. Cruz, L. A. Fern´andez, V. Mart´ın-Mayor, J. Pech, J. J. Ruiz-Lorenzo,A. Taranc´on, P. T´ellez, C. L. Ullod, and C. Ungil, Phys. Rev. B , 14237 (2000). P. O. Mari and I. A. Campbell, Phys. Rev. B , 184409 (2002). T. Nakamura, S.-i. Endoh, and T. Yamamoto, J. Phys. A , 10895 (2003). D. Daboul, I. Chang, and A. Aharony, Eur. Phys. J. B , 231 (2004). M. Pleimling and I. A. Campbell, Phys. Rev. B , 184429 (2005). S. Perez Gaviro, J. J. Ruiz-Lorenzo, and A. Taranc´on, J. Phys. A , 8567 (2006). F. Parisen Toldin, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp. P06002 (2006). T. J¨org, Phys. Rev. B , 224431 (2006). I. A. Campbell, K. Hukushima, and H. Takayama, Phys. Rev. Lett. , 117202 (2006). H. G. Katzgraber, M. K¨orner, and A. P. Young, Phys. Rev. B , 224432 (2006). J. Machta, C. M. Newman, and D. L. Stein, J. Stat. Phys. , 113 (2008); arXiv:0805.0794. M. Hasenbusch, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp. L02001 (2008). F. Belletti, M. Cotallo, A. Cruz, L. A. Fern´andez, A. Gordillo-Guerrero, M. Guidetti, A. Maio-rano, F. Mantovani, E. Marinari, V. Mart´ın-Mayor, A. Mu˜noz Sudupe, D. Navarro, G. Parisi,S. Perez-Gaviro, J. J. Ruiz-Lorenzo, S. F. Schifano, D. Sciretti, A. Taranc´on, R. Tripiccione, J.L. Velasco, D. Yllanes (the Janus Collaboration), Comp. Phys. Comm. , 208 (2008). K. H. Chen and T. C. Lubensky, Phys. Rev. B , 2106 (1977). F. Liers, J. Lukic, E. Marinari, A. Pelissetto, and E. Vicari, Phys. Rev. B , 174423 (2007). F. J. Wegner, in
Phase Transitions and Critical Phenomena , edited by C. Domb and M. S. Green(Academic Press, New York, 1976), Vol. 6. I. A. Campbell, K. Hukushima, and H. Takayama, Phys. Rev. B , 134421 (2007). H. Nishimori, Prog. Theor. Phys. , 1169 (1981). A. Georges, D. Hansel, P. Le Doussal, and J. Bouchaud, J. Phys. (Paris) , 1827 (1985). P. Le Doussal and A. B. Harris, Phys. Rev. Lett. , 625 (1988); Phys. Rev. B , 9249 (1989). H. Nishimori,
Statistical Physics of Spin Glasses and Information Processing: An Introduction (Oxford University Press, Oxford, 2001). M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, Phys. Rev. B , 184202 (2007). Y. Deng and H. W. J. Bl¨ote, Phys. Rev. E , 036125 (2003). M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, Phys. Rev. B , 094402 (2007). M. Campostrini, A. Pelissetto, P. Rossi, and E. Vicari, Phys. Rev. E , 066127 (2002). M. Hasenbusch, F. Parisen Toldin, A. Pelissetto, and E. Vicari, J. Stat. Mech.: Theory Exp.P02016 (2007). A. Pelissetto and E. Vicari, Phys. Rep. , 549 (2002). G. Toulouse, J. Physique Lettres , 447 (1980). N. Kawashima and H. Rieger, in
Frustrated Spin Systems , edited by H. T. Diep (World Scientific,Singapore, 2004); cond-mat/0312432. H. Nishimori, J. Phys. Soc. Japan , 3305 (1986). H. Kitatani, J. Phys. Soc. Japan , 4049 (1992). C. Wang, J. Harrington, and J. Preskill, Ann. Phys. , 31 (2003). C. Amoruso and A. K. Hartmann, Phys. Rev. B , 134425 (2004). M. Picco, A. Honecker, and P. Pujol, J. Stat. Mech.: Theory Expt. P09006 (2006). A. K. Hartmann, Phys. Rev. B , 3617 (1999). It can be shown rigorously that the Nishimori line never intersects the spin-glass phase, H.Kitatani, J. Phys. Soc. Japan , 2070 (1994). Since we must also have p F G ≤ p ∗ , the mixedphase, if it exists, should be confined to the region below the Nishimori line and on the left ofthe line p = p ∗ (see Fig. 1). D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett. , 1792 (1975). T. Castellani, F. Krzakala, and F. Ricci Tersenghi, Eur. Phys. J. B , 99 (2005). F. Krzakala and O. C. Martin, Phys. Rev. Lett. , 267202 (2002). S. Boettcher and E. Marchetti, Phys. Rev. B , 100405(R) (2008). S. Boettcher and J. Davidheiser, Phys. Rev. B , 214432 (2008). C. D. Lorenz and R. M. Ziff, Phys. Rev. E , 230 (1998). T. J¨org and F. Ricci-Tersenghi, Phys. Rev. Lett. , 177203 (2008). V. Privman, in
Finite Scaling and Numerical Simulations of Statistical Systems , edited byV. Privman (World Scientific, Singapore, 1990). J. Salas and A. D. Sokal, J. Stat. Phys. , 551 (2000). M. Campostrini, M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B , 144506 (2006). M. Hasenbusch, J. Phys. A , 4851 (1999). C. J. Geyer in