Numerical results for the Edwards-Anderson spin-glass model at low temperature

We have simulated Edwards-Anderson (EA) as well as Sherrington-Kirkpatrick systems of L^3 spins. After averaging over large sets of EA system samples of 3 =<L =<10, we obtain accurate numbers for distributions p(q) of the overlap parameter q at very low temperature T. We find p(0)/T -->0.233(4) as T -->0. This is in contrast with the droplet scenario of spin glasses. We also study the number of mismatched links --between replica pairs-- that come with large scale excitations. Contributions from small scale excitations are discarded. We thus obtain for the fractal dimension of outer surfaces of q~0 excitations in the EA model d_s -->2.59(3) as T tends to 0. This is in contrast with d_s -->3 as T -->0 that is predicted by mean field theory for the macroscopic limit.


I. INTRODUCTION
Whether spin glasses are complex systems is an important issue. We have discussed this in some detail in Ref. 1, where we gave numerical evidence for fundamental differences between the spin-glass phases of the Edwards-Anderson 2 (EA) and of the Sherrington-Kirkpatrick 3 (SK) models. In short, we studied spikes in probability distributions, p(q), of the overlap parameter, q, that vary widely over different sample systems. The variation of a suitably defined average spike width w over the values of the linear system sizes L we studied was shown to decrease sharply with L in the SK model. Furthermore, rms deviations δp away from p(q) over different system samples (that is, over different realizations of quenched disorder) increase sharply with L. Such behavior is consistent with mean field theory, which predicts the replica symmetry breaking (RSB) scenario in which w → 0 and δp → ∞ in the macroscopic limit. 4,5 Our results 1 for the EA model follow a different trend. Rather, w and δp become, within errors, independent of L in the zero temperature limit. [The statistics of spikes in overlap distributions in different system samples, which have been studied by Yucesoy 6 et al., also point away from an RSB scenario, though this conclusion is criticized in Ref. 7.] Much numerical work on the behavior of the EA model at low temperature stems from the observations of Moore 8 et al., that Monte Carlo simulations had up to then been performed at temperatures that were too close to the critical temperature, and therefore suffered from finite size effects that could be misinterpreted as RSB behavior. Low temperature data for p(q), which was the centerpiece of these considerations, were soon thereafter provided by Katzgraber, Palassini and Young 9 (KPY). A roughly constant value of p(q ∼ 0) over a 3 ≤ L ≤ 8 size range was shown to be consistent with these data. This is as in the RSB, not the droplet scenario 10,11 of spin glasses. We did not report data for p(q) in Ref. [1], because they were essentially the same as KPY's, and our statistical errors did not decisively improve on them. We have since simulated sample sets which are over an order of magnitude larger than KPY's, and cover a range of system sizes which is slightly larger. We are thus able to report here rather accurate data for temperatures as low as 0.16T sg , where T sg is the spin-glass transition temperature.
In the so called trivial-non-trivial (TNT) picture, proposed by Krzakala and Martin 12 and by Palassini and Young 13 (PY), p(q) is size independent in the neighborhood of q = 0, as in the RSB scenario, but the dimensionality d s of outer surfaces of q ∼ 0 excitations is smaller than the dimensionality, d, of the space where spins are embedded. Values of 2.57 d s 2.62 have been calculated by PY, KPY, and by Jörg and Katzgraber. 14 However, Contucci et al. 15 have obtained d s = 3 for the EA model in three dimensions (3D). This would be in accordance with a RSB scenario. These two conflicting results were obtained by different methods. Fractions of mismatched links (FML) between replica pairs are calculated in both methods. All 2.57 d s 2.62 values 9,13,14 were obtained (but see Ref. 16) from the rms deviation of the FML from its mean value (over time and system samples, as well as over all q). On the other hand, d s = 3, was obtained in Ref. 15 from the behavior of the mean FML, f ml (q), for each observed value of q. More specifically, the L → ∞ limit of f ml (q ∼ 0) was studied in Ref. 15 for T 0.5T sg . This limit was argued to be nonzero, which is what one expects of a space fulfilling surface. This conclusion fits with the RSB scenario, and clashes with the ones reached in Refs. 9, 13, and 14.
Most of this paper, is devoted to the fractional number of link mismatches, ∆f ml Q , it costs to create an excitation with a −Q < q < Q value. For a more precise definition of ∆f ml Q , consider first f ml Q , which is the average FML given that q is in a given −Q < q < Q interval. Subtraction from f ml Q of the average FML given that q is not in the −Q < q < Q interval gives ∆f ml Q . Both ∆f ml Q and f ml Q have the same zero temperature limit, arXiv:1304.7641v1 [cond-mat.dis-nn] 29 Apr 2013 but we believe ∆f ml Q is the natural extension to nonzero temperatures of the FML of large-size excitations in the ground state. Whereas f ml Q decreases as T decreases, ∆f ml Q increases. This enables us to bracket very low temperature behavior and confidently make T → 0 extrapolations. These notions stand out clearly in the frustrated box (FB) model which we define below. The plan of this paper is as follows. In Sec. II, we define the models, the spin-overlap and link-overlap parameters, and the simulation procedure. In Sec. III, we report accurate data for p(q) for EA and SK systems at very low temperature. These data show that, in the 3 ≤ L ≤ 10 range, p(q ∼ 0) is independent of L at very low temperatures. The conclusion KPY 9 reached, that the EA model exhibits a clear trend away from the droplet scenario, is thus strengthened. In Sec. IV A we examine the large-scale behavior of the FML in the FB model. This simple model helps to highlight the pitfalls that should be avoided in the interpretation of a nonzero macroscopic limit of FML. In Sec. IV B, we assign a (average) mismatching-link cost, ∆f ml Q , to an excitation with a −Q < q < Q value. Numerical results for ∆f ml Q , which imply a fractal dimension of 2.59(3) for the surface associated to ∆f ml Q , are also given in Sec. IV B. Such a value of d s , smaller than the dimensionality 3D of the space where spins are embedded, is in contradiction with mean field theory predictions, but is as envisioned in the TNT scenario 12,13 of the EA model. We summarize our conclusions in Sec. V.

II. MODELS, DEFINITIONS AND PROCEDURE
In all models we study, an Ising spin sits on each one of the N ≡ L × L × L sites of a simple cubic lattice in three dimensions (3D). We use periodic boundary conditions throughout. In the SK and EA models the interaction energy between a pair of spins at sites i and j is given by J ij σ i σ j . We let J ij = ±1/ √ N randomly, without bias, for all ij site pairs in the SK model. For the EA model, J ij = 0 unless ij are nearest-neighbor pairs, and we draw each nearest-neighbor bond J ij independently from unbiased Gaussian distributions of unit variance.
We let all temperatures be given in units of 1/k B , where k B is Boltzmann's constant. Thus, the transition temperature T sg between the paramagnetic and SG phase of the SK model is T sg = 1. 3,5 For the EA model T sg 0.95. 17 We let σ  for an identical replica, replica (2), of the same system. As usual, we define that is, q is the average (over all sites) spin alignment between the states replicas 1 and 2 are in.  As in Refs. 12 and 13, we define the link-overlap, where N l is the total number of links, and the sum is over all i, j nearest neighbor pairs (of which there are 3N in the nearest neighbor EA model in 3D). The FML between replicas 1 and 2 is given by f ml = (1 − q l )/2. In addition, as in Ref. 15, we define f J ml (q) as the time average of the FML (for a sample with a given set J of bonds) which is observed over all time intervals while the value of the spin-overlap is q. Unfortunately, the dimensionality d s of large scale excitations does not follow straightforwardly from the behavior of f J ml . This is because smaller scale excitations contribute to f J ml (q). More on this can be found in Sec. IV. Let F(q) be any q dependent function, such as p(q) or f ml (q). We define The advantage of working with F Q is that statistical errors for it are smaller than for F(q). Accordingly, most results below are given for p Q and f ml Q rather than for p(q) and f ml (q). How statistical errors on p Q depend on Q is worked out in Appendix A.
In addition, we let A J stand for the value of some observable A on a sample defined by the set J of bonds, and we let A J J stand for the average over samples of J . We make use of the parallel tempered MC method. 18,19 Details on how we apply it to the EA and SK model are as specified in Ref. [1]. However, some details differ. For all sizes of the EA model, temperatures are spaced here by 0.04 (0.08) in the 0.16 ≤ T ≤ 0.48 (0.48 ≤ T ≤ 1.6) range. The rationale for this, as well as checks we perform in order to make sure equilibrium is reached, can be found in Appendix B. Temperatures of all SK systems were evenly spaced by ∆T = 0.04 in the whole 0.12 ≤ T ≤ for all 4 ≤ L ≤ 8 (L = 12). Values for average swap success rates, α, between pairs of EA and SK systems at the lowest two temperatures are given in Table I. Larger values of α are observed for higher temperatures. In the FB model, the smallest value of α, which we give in Table  1, is observed in the critical region. The number, N s , of sample systems we average over, is, as specified in Table I, much larger here than in Ref. [1]. We have tried not to make N s smaller with increasing L. This is because, as we show in Appendix A, statistical errors are independent of L, because of non-self-averaging. (For L = 10, we could only do 20 000 samples. That took some 50 years worth of computer time.)

III. AVERAGE q-DISTRIBUTIONS AT LOW TEMPERATURES
Plots of p(q) vs q are shown in Fig. 1 for EA systems of various sizes at T = 0.16. Plots of p Q /T vs T are shown in Fig. 2 for Q = 1/4 and 1/2. Error sizes are clearly smaller for the larger value of Q. Simulation details, such as sample numbers and running times, are given in Table  I.
For comparison, plots of p Q /T vs T for SK systems of various sizes are shown in Fig. 3 for Q = 1/2 and 1/4. We note that p Q ∼ T if T 0.3T , independently of L, following mean field predictions, 5,20 for the SK model.
For a more accurate picture of how p Q varies with L in the EA model, we show log-log plots of p Q /T vs L, for Q = 1/4 and T = 0.2 in Fig. 4. For comparison, we also show data points from the KPY paper for the same temperature and Q = 0.2. 21 The best fit of p Q ∝ L −θ to the data points shown in Fig. 4 gives θ 0.008. Fits following from letting θ = 0.04 and −0.02 give χ 2 parameters that are over twice as large as the one for θ = 0.008. legend for further details.) For higher temperatures, up to T 0.4, as well as for Q = 1/2 and all T 0.5, all error bars are smaller than the ones shown in Fig. 4, and all best fits of p Q ∝ L −θ to the data give | θ |< 0.01. Thus, future generation of more accurate data that would give θ > 0.04 for the 3 ≤ L ≤ 10 range is rather unlikely.  A. The frustrated-box model We define here a nearest-neighbor Ising model in which most bonds are ferromagnetic. For reasons given below, we term it the frustrated-box (FB) model. Consider plane P 1 , perpendicular to the x axis, at x = 1/2, which cuts all bonds between x = 0, y, z, and 1, y, z sites. Similarly, P 2 at x = L/2 + 1/2, cutting all bonds between x = L/2, y, z and L/2 + 1, y, z. These planes divide the system into two equal portions. In this model, only nearest-neighbor spins interact. All bonds, except the ones that cut across P 1 and P 2 , are of strength 1, that is, ferromagnetic. Half the bonds that cut across both P 1 and P 2 are of strength −1, that is, antiferromagnetic, and the rest are of strength 1. More precisely, all ±1 bonds that cut across both P 1 and P 2 are distributed on a checkerboard pattern. We apply periodic boundary conditions.
In the ground state, all spins within the box (that is, between planes P 1 and P 2 ) are parallel, and so are all spins outside the box. These two spin subsystems can point in the same or opposite directions. Thus, the ground state is (because of invariance under all-spin reversal) four-fold degenerate. The box defined by P 1 and P 2 and the system's boundary is the 3D analog of Toulouse's two-dimensional frustrated plaquettes. 22 Hence, the "frustrated-box" label.
The number of broken bonds in all ground states of the FB model is L 2 , but the number of bond mismatches between two replicas is (in ground states) either 0 or 2L 2 . Thus, f ml = 1/3L but f ml (q = 0) = 2/3L. Plots of f ml Q vs L are shown in Fig. 5 for Q = 1/2 in FB systems at various temperatures. Note f ml Q = 2/3L at T = 1, as expected for T T c 4.5.
Curves for T = 2.4 and 2.8 in Fig. 5  Let us first examine a simple picture of large and small scale excitations. In Fig. 6, the same cross section of an EA sample system of 10 3 spins at T = 0.36 is shown at four different times (consecutive times are at least 10 6 MC sweeps apart) of a single MC simulation. For a 3D picture of the outer surface of a large scale excitation in an EA system at T = 0.16 see Fig. 7.
The general idea is to determine the area of the outer surface of large scale excitations, such as the ones on both right-hand panels of Fig. 6. We intend to do this by subtracting the total surface area of small scale excitations from the total (from small and large excitations) surface area.
For each system sample, we first obtain f J ml (q) for each q by adding 1/2(1 − q l ) whenever q is observed in a given MC run, and we finally divide the result by the number of times q has been observed. Now, the average surface area over all excitations observed in a given sample whenever a value of q is in the −Q < q < Q range is given by, where u J = Q −Q dq p J (q), and by, where v J = |q|>Q dq p J (q), whenever q is outside the −Q < q < Q range. Finally, for the average FML it costs to create an excitation in the −Q < q < Q range, we calculate, We have calculated . . . J in the above equation by each of the following two procedures: (1) giving equal weight to all system samples for which u J > 0.1, and (2) giving each sample a weight proportional to Q −Q dq p J (q). Within statistical errors, we have obtained the same results from these two procedures. from such scaling behavior can be observed, even over this limited range of system sizes, in analogous plots (not shown) for T as small as T 0.3. This effect is more clearly exhibited in Fig. 8(a)

8(b).
In it, plots of ∆f ml Q vs L are shown for the same Q and T as in Fig. 8(a). All data points shown for ∆f ml Q in Fig. 8

V. CONCLUSIONS
We have reported data for p(q) from averages over large sets (numbers are shown in Table I) of EA and SK sys-tems at very low temperature. The data for p Q improves our confidence level in the conclusion that p(q ∼ 0) is nonzero and system-size independent in the EA model. Future generation of more accurate data that would give p(0) ∼ 1/L θ and θ > 0.04 in the 3 ≤ L ≤ 10 range of L values is rather unlikely. Thus, the conclusion KPY had reached, 9 that the EA model exhibits a clear trend away from the droplet scenario, is strengthened. Furthermore, our results are consistent with p(0)/T → 0.233(4) as T → 0 in the EA model (and p(0)/T → 0.51(3) as T → 0 in the SK model).
We have studied the fraction of link mismatches, ∆f ml Q , it costs to create an excitation with a −Q < q < Q value. For a wide range of temperatures in the spin-glass phase, ∆f ml Q seems to vanish, as in the TNT picture, 12 in the macroscopic limit. Data points for ∆f ml Q (for Q = 1/4 and 1/2) are consistent with We show here how statistical errors ∆(N s ) for p Q depend on N s , on system size, on Q and on T for T 1. Consider first the rms deviation δp(q) of the probability density p J (q) from its average over different samples.
It is plotted vs | q | in Fig. 3 of Ref. [1] for EA systems of various sizes at T = 0.1. Because there is no self-averaging, δp p, except near q = 1. In addition, δp(q)/p(q) does not decrease as L increases. This has an unwanted implication, namely, fractional statistical errors in p(q) do not decrease as system size increases if N s remains constant.
To start, let F J (q 1 , q 2 ) ≡ p J (q 1 )p J (q 2 ), and let G(q | Q) be the average of G J (q | Q) over samples. We can then write,

follows. Now, let
define w. We also note, Here, the first term is much smaller than the second one for T 0.5 and all Q 1/2, in both the EA and SK models. This comes from the fact that, whereas p Q vanishes as T → 0, (δp) 2 Q does not. Therefore, ∆ 2 (1) (w/Q)(δp) 2 Q for T 0.5 and all Q 1/2, whence, follows immediately. This is clearly consistent with nonself-averaging. It shows that ∆ is, at least for the values of L we study here, independent of L, for both the EA and SK models. Equation (A5) also shows how much precision is gained by averaging p(q) over −Q < q < Q.
We can substitute into Eq. (A5) the low temperature values of w(δp) 2 Q /p Q from Figs. 11(a) and 11(b) for the EA and SK models, respectively. Further substitutions of p Q from Sec. III give where c 0.2, 0.3 for the EA and SK models, respectively, at low temperatures. This is the desired expression. We show here (i) how we choose the success rate, α, for state swapping (that is, for exchanging spin configurations) between two systems, and (ii) how we checked equilibration was achieved in our simulations.
We first derive an expression for α. In the parallel tempered MC algorithm, 18,19 the probability, P (ss | ∆E), for state swapping to take place between systems 1 and 2, at temperatures T 1 , T 2 , where ∆E = E 2 − E 1 , is given by, P (ss | ∆E) = 1 if ∆E ≤ 0, but P (ss | ∆E) = exp(−∆β∆E) if ∆E > 0. Now, α = d∆E P (ss | ∆E) P (∆E). (B1) In thermal equilibrium, the probability that the energy of systems at T 1 and T 2 differ by ∆E is given by, where, neglecting variations in the specific heat (per spin) c in the T 1 < T < T 2 range, N σ 2 is the mean square energy deviation coming from thermal fluctuations at both T 1 and T 2 , It then follows that, Substitution into Eq. (B1) yields, assuming ∆T [erfc(γ) = (2/ √ π ∞ γ dx exp(−x 2 )]. A choice of α ≈ 0.5 might seem to lead to efficient MC simulations, which, using Eq. (B5), would lead ∆T √ cN /2T ≈ 0.48. Note however that increasing ∆T does make α smaller, but it also implies fewer random steps need be taken by a given state in order to travel from a system at the minimum temperature to one at the maximum temperature. Furthermore, smaller temperature differences imply fewer systems to be simulated, which leads to further computer time saving.
This point is illustrated in Fig. 12(a), where plots of p Q /T vs 10 5 φ/τ s are shown for EA systems of 6 × 6 × 6 spins at T = 0.4. These data points come from tempered MC runs of sets of equally spaced temperatures. Values of the swap success rate, α, between the two systems at the lowest pair of temperatures are given for each ∆T , given by T n+1 − T n , in Fig. 12(a). From the values of φ given in Fig. 12(a), we conclude that values of α as small as 0.2 do not lead significantly slower simulations. Figure Fig. 12(b) is as 12(a) but for L = 8, T min = 0.1, and x axis values are for 10 6 /τ s . Note even an α value as small as 0.04 only slows simulations down by 20%.