Order-parameter fluctuations in Ising spin glasses at low temperatures

We present a numerical study of the order-parameter fluctuations for Ising spin glasses in three and four dimensions at very low temperatures and without an external field. Accurate measurements of two previously introduced parameters, A and G, show that the order parameter is not self-averaging, consistent with a zero-temperature thermal exponent value \theta' \simeq 0, and confirm the validity of the relation G=1/3 in the thermodynamic limit in the whole low-temperature phase, as predicted by stochastic stability arguments.


I. INTRODUCTION
Understanding the low-temperature physics of shortrange spin glasses 1 remains a major unsolved problem. Much of the current debate concentrates on the equilibrium thermodynamics of the Edwards-Anderson model with Ising spins (EAI model), the canonical short-range spin glass. Since analytical approaches pose formidable difficulties, the problem is often studied numerically. However, the existence of large barriers between lowenergy configurations has limited so far numerical calculations to small systems sizes, from which it is hard to draw definite conclusions on the large-volume limit.
Two main issues have been addressed in many numerical studies: the existence and character of the finitetemperature spin-glass transition, and the nature of the low-temperature spin-glass phase. A central quantity of interest in the description of the spin-glass phase is the scaling exponent θ ′ governing the typical energy of the lowest-lying excitations with linear size of order l, which is assumed to scale as E ∼ l θ ′ . In general, θ ′ may be distinct 2,3 from the stiffness exponent θ measured in domain-wall computations 4,5,6,7 , and the stability of the spin-glass phase requires to have θ ′ ≥ 0. In a "many-state" picture 8 such as the replica-symmetrybreaking picture inspired by mean-field theory 1,9,10 , one has θ ′ = 0, hence there are excitations whose energy remains finite (of the order of the coupling strength between two spins) even as their length scale diverges. In a "two-state" picture, such as the droplet model 11,12 , one has θ ′ > 0, hence the energy of large-scale excitations diverges with their size. In this case the identity θ ′ = θ is often assumed.
Both the spin-glass transition and the ordered phase have been usually investigated numerically by computing sample-averaged quantities such as the Binder cumulant 13 or the distribution of the order parameter (OP) and related observables. Recently 14,15 , it was observed that useful information on both issues can be drawn from the sample-to-sample fluctuations of the OP.
In particular, two dimensionless measures of the OP fluctuations were considered 14,15 : A, the normalized fluctuation of the spin-glass susceptibility, and G, a ratio between two cumulants of the OP distribution. These two parameters are related to the Binder cumulant, B, via the relation B = 1 − A/(2G). For a model without time-reversal symmetry (TRS), A, G, and B are given by Eqs. (2), (3), and (4) below. The parameter G serves as a good indicator of the existence of phase transitions 15 in systems lacking TRS (for which B is generally a bad indicator), as recently shown for several systems, including Migdal-Kadanoff spin glasses 16 , RNA folding models 17 , chiral spin systems 18 , and mean-field models such as the SK model (with and without a magnetic field), the infinite-range p-spin model 19 , and the infinite-range Potts model 20 . The parameter A has also been studied before for random diluted models at criticality 21,22 .
In this paper, we investigate the OP fluctuations in the EAI model with Gaussian couplings in three and four dimensions by low-temperature Monte Carlo simulations. We study the case with no external field, which satisfies TRS. In three dimensions (3D), numerical data are available in the literature for A in the high-temperature phase 23 and for G near the critical point 24 , for the "±J" coupling distribution. In four dimensions (4D), G was measured at moderately low temperatures, also for the ±J distribution 14 . Here, we study much lower temperatures than in these studies, in order to reduce crossover effects associated to the critical point 25,26 , which complicate the interpretation of the numerical data at higher temperatures.
A summary of our results is as follows. First, we estimate θ ′ from the system-size dependence of A, finding, for the system sizes we could reach, a small value of θ ′ incompatible with the accepted values of the domain-wall exponent (θ ≃ 0.2 in 3D 4,5,7 and θ ≃ 0.7 in 4D 6 ), and compatible with zero. This agrees with recent determinations of θ ′ from ground-state perturbation methods 2,3,27 and from low-temperature measurements of the OP distribution 10,26,28,29 (which all consider sample-averaged quantities), and supports a picture of the spin-glass phase characterized by two distinct exponents, θ > 0 and θ ′ = 0. The result θ ′ = 0 implies that the OP is not self-averaging in the thermodynamic limit.
Second, we find good evidence that the identity G = 1/3 holds in the whole spin-glass phase in the thermodynamic limit, confirming the validity of sum rules proposed by Guerra 30 and first derived for the SK model, which follow from the property of "replica equivalence" 31,32 .
Third, we find that A and G allow to locate the spinglass transition reasonably well although, as expected due to TRS and as previously numerically observed 24 , B provides a better determination (a much more accurate determination is provided by the correlation length 24 , which we do not investigate here).
We do not study in this paper the surface fractal dimension d s of the excitations, which is the other exponent, besides θ ′ , characterizing the spin-glass phase (in particular, d s = 0 in the standard replica-symmetrybreaking picture 10 , d s > 0 in the droplet model, while the "TNT" picture 2,3 predicts the "mixed" behavior d s > 0, θ ′ = 0).
The rest of the paper is organized as follows. In Sec. II we introduce the different models and observables studied, and discuss the theoretical predictions for these observables. In Sec. III we present and analyze our numerical results for the quantities A, G and B. Finally, in Sec. IV we summarize our conclusions.

II. MODELS, OBSERVABLES, AND THEORETICAL PREDICTIONS
We study the EAI model defined by the Hamiltonian where L D Ising spins S i sit on a (hyper-)cubic lattice in D dimensions with linear size L and periodic boundary conditions in all directions. The couplings J ij are drawn from a Gaussian distribution of zero mean and unit variance. We consider two different models: (i) the case with interactions i, j restricted to nearest neighbors (referred to as NN model) in D = 3 and 4; (ii) the case with interactions restricted to nearest, next-nearest, and nextnext-nearest neighbors (referred to as NNN model) in D = 3, which has a coordination number z = 26.
The NN model has been extensively studied and is known to display a finite-temperature continuous spinglass transition for D ≥ 3. Recent estimates of the critical temperature for Gaussian-distributed couplings are T c = 0.95 ± 0.04 in 3D 33 and T c = 1.80 ± 0.03 in 4D 34 .
The NNN model has been much less studied. In Ref. 23, the 3D case with ±J couplings was considered, but no conclusive evidence of a finite-temperature transition was obtained, the data being compatible with both T c ≈ 3.27 and a zero-temperature singularity. Incidentally, in the NNN case we do not expect a large difference in T c between binary and Gaussian distributions, due to the large coordination number. For the NNN model, we did not consider temperatures as low as for the NN model, but focused on the phase transition region. Another of the results of this paper is a convincing evidence that indeed a finite-temperature transition exists in this model in 3D.
We measure A, G, and B as a function of temperature T and size L using the following definitions: is the spin overlap of two independent systems S a i and S b i with the same random couplings, and ... and (...) stand for thermal and disorder averages, respectively. The OP for a given realization of the disorder is q 2 , therefore A is nothing but the normalized sample-to-sample variance of the OP.
In the paramagnetic phase, T > T c , and for sufficiently large L so that L ≫ ξ (where ξ is the correlation length), the OP follows a Gaussian distribution and all three parameters vanish as 1/L D , as follows from the central limit theorem. Following the terminology of Wiseman and Domany 21 , this means that the OP is strongly selfaveraging.
At T = T c , the correlation length diverges and the central limit theorem cannot be applied. For strongly disordered systems such as spin glasses it is known that the OP is not self-averaging at criticality 21,22 , namely A tends to a finite value in the thermodynamic limit. If A is finite then clearly G must be finite, and standard renormalization-group arguments show that B is also finite at T c . Since B and G are dimensionless and monotonic in T , in plots of these quantities as a function of T , the curves for different values of L must all cross at T = T c , and one can use this to determine T c . From standard finite-size scaling, one can then determine the critical exponent ν. Since much work has been devoted to measuring ν from the standard observables (see for instance Refs. 24,33,35,36,37), and we are primarily interested in the low-temperature phase here, we will not attempt a precise determination.
In the spin-glass phase, T < T c , A is expected 38,39 to vanish linearly with T according to the scaling law where θ ′ is the exponent discussed in the Introduction. This law holds under two hypothesis (both satisfied in the case of continuous couplings studied here): (i) the ground state is unique; (ii) the probability distribution of the energy of the lowest-lying excitations has finite weight at zero energy 38,39 . From the above scaling law we see that if θ ′ > 0, then A vanishes for L → ∞, namely the OP is weakly selfaveraging (where "weakly" indicates 21 that the OP fluctuations vanish more slowly than 1/L d , a consequence of the inequality θ ′ < d). This situation is encountered in the droplet model 11,12 , as discussed in the Introduction, and also in mean-field models with a marginally stable replica-symmetric solution at low temperatures (such as the spherical SK model 40 ). If θ ′ = 0, as in a "manystate" picture, A remains finite in the thermodynamic limit, namely the OP is not self-averaging.
Turning now to G, it is known 31,32 that in the SK model the following relation holds for T < T c : According to Guerra 30 , this relation should hold (for T < T c ) in any model which is "stochastically stable" with respect to a mean-field perturbation and which has a non-self-averaging OP. Under the hypothesis (i) and (ii) above, the more general conjecture has also been made 40 that the above relation holds for T < T c even if the OP is self-averaging. In this case, G would be finite but both the numerator and the denominator in Eq.

III. NUMERICAL RESULTS
We simulated the various models with the parallel tempering technique 41 , which allows to reach significantly lower temperatures than conventional Monte Carlo methods. The parameters of the simulation are given in Table I. Equilibration of the Monte Carlo runs was tested by monitoring all the measured observables on a logarithmic time scale, checking that they had all converged within their statistical errors, and by applying the equilibration test discussed in Ref. 26.

A. Parameter A
In Figures 1, 2, and 3 we show our numerical results for A in the 3D NN, 4D NN, and 3D NNN models, respectively. The vertical lines in Figures 1 and 2 indicate the estimated value of T c . In all cases, the behavior of A resembles that observed in the SK model (see Refs. 14,18,19). At high temperatures, A decreases with L, approximately as 1/L D , showing that the OP is strongly self-averaging in this regime, as expected. Near  T = T c , there is a maximum whose position shifts towards T c as L increases, an effect of finite-size corrections. The shift is modest in the 4D NN model but quite noticeable in the 3D NN model, where even for the largest L the position of the maximum is still significantly larger than T c . In the 3D NNN model, the position of the maximum is also larger than T c (see discussion in Sec. III C on the value of T c in this model), but the shift is less pronounced. The height of the maximum increases with L in all models, indicating that A attains a finite value in the thermodynamic limit (since it is bounded from above), namely that the OP is not self-averaging at T c , as expected. At low temperatures (T < T c ), Figures 1, 2, and 3 show that A is approximately linear in T , in agreement with Eq. (5). Most interestingly, the data for different values of L tend to superimpose to each other. In a scenario with θ ′ > 0, the data should tend to zero for large L in the whole region below T c . In 3D, we see no decrease at all in the data with increasing L, while a modest decrease is observed in 4D.
To analyze in more detail the size dependence of A at low temperatures, in Figure 4 we plot the ratio A/T as function of L at different temperatures for the three models. The straight lines represent the scaling law Eq.(5) assuming θ ′ = θ, and using the estimates of θ from domainwall calculations, θ = 0.2 in 3D 4,5,7 and θ = 0.7 in 4D 6 . No estimates of θ are available for the 3D NNN model, so we use that for the NN model (we expect that θ is a universal exponent equal for both models). Clearly, the data in Figure 4 do not agree with the hypothesis θ ′ = θ for the range of sizes considered, and seem to saturate to a constant value instead. We fitted the data with the form A(T, L)/T = aL −θ ′ , whereθ ′ should be seen as an "effective" exponent which depends on the temperature and which effectively takes into account corrections to the leading scaling behavior, withθ ′ → θ ′ in the limit T L −θ ′ → 0. The fits giveθ ′ varying from 0.03 ± 0.02 (T = 0.7) to 0.00 ± 0.06 (T = 0.2) in the 3D NN model, from 0.30 ± 0.05 (T = 1.0) to 0.003 ± 0.006 (T = 0.32) in the 4D NN model, and from 0.03 ± 0.04 (T = 2.8) to 0.08 ± 0.04 (T = 2.0) in the 3D NNN model. Therefore, in all cases the data is compatible with θ ′ = 0, in agreement with the "many-state" picture, and is statistically incompatible with with θ ′ = θ, in disagreement with the "two-state" picture. As usual, we cannot exclude a crossover 25 to a larger value of θ ′ for larger L. In this case, in the large-volume limit A would be zero at all temperatures, except at T = T c .
A value of θ ′ compatible with zero was also obtained from the OP distribution 10,26,28,29,36 and from direct measurements of the energy of low-lying excitations created by perturbing the ground state 2,3,27 . for a continuous coupling distribution. More importantly, G(T, L) seems to converge to the value 1/3 for large L in the whole low-temperature region, in agreement with Eq. (6). This is particularly evident in the 4D NN model (Figure 6), where G(T, L) has already converged to 1/3 for L = 5 at temperatures below T ≈ T c /2 ≈ 0.9 (the data points above 1/3 are due to incomplete equilibration 19 ). It is quite unlikely that the saturation is a finitesize effect, therefore our results strongly suggest that indeed G = 1/3 in the whole spin-glass phase. As discussed above, if A remains finite as L → ∞ (as indicated by our data), this is an expected consequence of stochastic stability 30 . If A vanishes, instead, our results would support the more general conjecture of Ref. 40.

C. Critical region
In this subsection we comment on the behavior of G and B near the critical temperature, starting from the NN model.
In 3D, the vertical lines in Figure 5 indicates the position of the critical temperature, using the value T c = 0.95±0.04 quoted in Ref. 33, which was obtained from the parameter B measured in a large-scale simulation. One sees that the data for both B and G for different values of L come together as T approaches T c from above, as indicative of a phase transition. Below T c , the data for G separate again in a statistically significant way, while for B one would need a substantially larger statistics (or larger sizes) to see a clear separation, as observed in previous studies 33,37 . For example, at T = 0.82 the separation between the L = 12 and the L = 4 data is 1.4 standard deviations for B and 2.5 standard deviations for G. The small separation below T c is probably due to the vicinity of D = 3 to the lower critical dimension 33,35,37 . We also note that the crossing point of G is at temperatures larger than T c , and close inspection shows that it shifts towards T c from above as L increases. A similar shift was observed for the position of the maximum of A in Figure 1.
In 4D, both B and G display a very clear crossing (see Figure 6), as also observed in previous studies 29,34 . From B we estimate T c = 1.80 ± 0.03 in agreement with the results of Refs. 34,42. This value is indicated by the vertical lines in Figure 6. As in 3D, the crossing point of G is at temperatures larger than T c and shifts towards T c as L increases.
Overall, this confirms that both in 3D and 4D the corrections to scaling are significantly larger for G and A than for B. Since G and A have also much larger statistical errors than B, the latter quantity is to be preferred to G and A for locating T c in models with TRS. As already mentioned, a much more accurate quantity for this purpose is the correlation length, which shows a very clear crossing in 3D 24 , unlike B and G.
Finally, in the 3D NNN model both G and B show a rather clear crossing (see Figure 7). This provides a clear evidence for the existence of a phase transition in 3D Ising spin glasses, confirming recent results for the NN model 24,35,37 that obtained a convincing evidence (especially Ref.24) after the issue had remained unsolved for a long time. The crossing point is at T ≃ 3.3 for B and at somewhat higher temperatures for G, although also here the crossing for G shifts to the left as L increases. From the data for B one might be tempted to conclude that the critical temperature is T c ≃ 3.3. However, if this was the case, the value of B at T c (which is a universal quantity) would be lower in the 3D NNN model than in the 3D NN model, violating universality. This suggests that the actual value of T c is significantly lower than 3.3, despite the clear crossing of B (which would then be strongly affected by scaling corrections), and for this reason we have not indicated the position of T c in Figures 3 and 7. A more detailed analysis 43 clearly shows that indeed T c is significantly lower than 3.3 in this model.  figure) and B (inset) for the 3D NN model, as a function of the temperature and for various system sizes. The vertical lines correspond to Tc = 0.95±0.04, the horizontal line in the main figure corresponds to the limit relation G = 1/3.

IV. CONCLUSIONS
To conclude, we have provided evidence that the order parameter is not self-averaging in the low-temperature phase of the Edwards-Anderson Ising spin glass in 3D and 4D, which implies an exponent θ ′ ≃ 0, in agreement with a "many-state" picture of the spin-glass phase, such as the replica-symmetry-breaking picture or the "TNT" picture. As usual, due to the limited system sizes that are currently reachable in numerical simulations, we cannot exclude that for larger sizes one recovers self-averaging, nevertheless our result is consistent with other studies which used sample-averaged quantities 2,3,26,27 and also found θ ′ ≃ 0. Independently of whether there is selfaveraging or not, we have provided evidence that the identity G = 1/3 holds in the thermodynamic limit in the whole spin-glass phase, a fact that calls for a theoretical explanation in terms of the geometry and energetics of the low-lying excitations. We have confirmed that G and A can be used to locate the spin-glass transition, although in models with time-reversal symmetry the usual sample-averaged parameters provide a better determination. Finally,we have confirmed the existence of a spin-glass phase transition at finite temperature in three dimensions.