# The Ising Spin Glass in dimension four

## Abstract

The critical behaviors of the bimodal and Gaussian Ising spin glass (ISG) models in dimension four are studied through extensive numerical simulations, and from an analysis of high temperature series expansion (HTSE) data of Klein et al. (1991). The simulations include standard finite size scaling measurements, thermodynamic limit regime measurements, and analyses which provide estimates of critical exponents without any consideration of the critical temperature. The higher order HTSE series for the bimodal model provide accurate estimates of the critical temperature and critical exponents. These estimates are independent of and fully consistent with the simulation values. Comparisons between ISG models in dimension four show that the critical exponents and the critical constants for dimensionless observables depend on the form of the interaction distribution of the model.

###### pacs:

75.50.Lk, 05.50.+q, 64.60.Cn, 75.40.Cx## I Introduction

Renormalization Group Theory (RGT) for thermodynamic phase transitions (1) and the Edwards-Anderson model for Ising Spin Glasses (ISGs) (2) were introduced almost simultaneously forty years ago. Ever since it has been tacitly assumed as self-evident that the standard RGT universality rules should apply to ISGs. As far as we know there is no rigorous theoretical proof that this ISG hypothesis holds, though confirmations have been reported a number of times based on numerical data (3); (4); (5); (6). The universality principle states that for all systems within a universality class the critical exponents are strictly identical and do not depend on the microscopic parameters of the model. All ISG models in a given dimension are supposed to be in the same universality class, on the assumption that the form of the interaction distribution is an irrelevant microscopic parameter.

Thus in the family of simple ferromagnets, within a universality class of models having space dimension and spin dimensionality , all models have identical critical properties corresponding to an isolated fixed point in the renormalisation group flow. However, diluted ferromagnets of given have a different (dilution independent) set of critical exponents, with values which correspond to a separate isolated fixed point (7); (8). For a few special cases of spin models in dimension two (discussed for instance in Ref. (9)) the critical behavior is more complicated and corresponds to a line of fixed points rather than an isolated fixed point; critical parameters vary continuously according to motion along the line, produced by a marginal operator. From the empirical ISG data below there appear for the moment to be two possible scenarios : two classes of ISGs (such as models with continuous distributions and those with discrete distributions) or alternatively ISG exponents which vary continuously with a parameter such as the kurtosis of the interaction distribution. In any case it has been stated by authoritative authors that “classical tools of RGT analysis are not suitable for spin glasses” (10); (11); (12) although no explicit theoretical predictions have been made so far concerning the important question of universality.

Here we combine numerical simulation and high temperature series expansion (HTSE) data on the bimodal and Gaussian ISGs in dimension four so as to obtain accurate and reliable values for the critical parameters in this model. We discuss a number of different methods for exploiting numerical data, and show that for each model these are consistent. Comparisons between these and other estimates on ISG models in the same dimension but with different interaction distributions show that the critical exponents and the critical constants for dimensionless parameters depend on the form of the interaction distribution.

## Ii Numerical techniques

The Hamiltonian is as usual

(1) |

with the near neighbor symmetric distributions normalized to . The Ising spins live on simple hyper-cubic lattices with periodic boundary conditions. We have studied the bimodal model with a interaction distribution and the Gaussian interaction distribution model. We will compare with published measurements on these and other d ISGs. We will use the inverse temperature , with the normalization above, or alternatively , to signify the temperature. The spin overlap parameter is defined by

(2) |

where and indicate two copies of the same system and is the number of sites.

The simulations were carried out using the exchange Monte Carlo method for equilibration using so called multi-spin coding, on (up to ) or (for larger ) individual samples at each size. An exchange was attempted after every sweep with a success rate of at least 30%. At least 40 temperatures were used forming a geometric progression reaching down to for the bimodal model and for the Gaussian model. This ensures that our data span the critical temperature region which is essential for the FSS fits. Near the critical temperature the step length was at most . The various systems were deemed to have reached equilibrium when the sample average susceptibility for the lowest temperature showed no trend between runs. For example, for this means about sweep-exchange steps.

After equilibration, at least measurements were made for each sample for all sizes, taking place after every sweep-exchange step. Data were registered for the energy , the correlation length , and for the spin overlap moments. In addition the correlations between the energy and observables were also registered so that thermodynamic derivatives could be evaluated using the relation where is the energy (14). Bootstrap analyses of the errors in the derivatives as well as in the observables themselves were carried out.

## Iii Finite size scaling

ISG simulations are much more demanding numerically than are those on, say, pure ferromagnet transitions with no interaction disorder. The traditional approach to criticality in ISGs has been to study the temperature and size dependence of dimensionless observables, principally the Binder cumulant and the correlation length ratio , in the near-transition region and to estimate the critical temperature and exponents through finite size scaling (FSS) relations after taking means over large numbers of samples. Finite size corrections to scaling must be allowed for explicitly which can be delicate as the range of sizes is generally small. On this FSS approach the estimated values for the critical exponents are very sensitive to the critical inverse temperature estimates. Here we first obtain estimates for and the critical parameters using FSS.

The data for standard dimensionless observables, the Binder cumulant

(3) |

and the correlation length ratio are shown in Figs. 1 and 2 as functions of size at fixed temperatures near . We show data from the present simulations together with data taken from the appropriate figures in Ref. (15). The standard FSS expression for a dimensionless observable , valid in the critical region is

(4) |

where . We plot and against with . The choice of this value, close to the estimate of Ref. (15), will be explained below. By inspection there is good point by point consistency between the present data and those of Ref. (15) where the statistical accuracy was much higher, but where the temperatures studied did not quite span the critical temperature. The consistency between the two data sets implies that full equilibrium has been reached. The data in Figs. 1 and 2 for fixed should tend to straight lines at criticality, curving upwards and downwards respectively at large for higher and lower than ; one can estimate from this criterion. The value is marginally higher than the estimate of Ref. (15). Data we have obtained on other dimensionless parameters lead to very similar estimates for (16).

It is important to consider the value of the irrelevant scaling field thermal correction exponent or equivalently the finite size correction exponent . The first term in the RGT -expansion for the ISG leading irrelevant operator is (17); (18), so , , (see (13) for the analogous site percolation -expansion). These values are only indicative as they are obtained from the leading terms of a series which (like the -expansion series for the other exponents in ISGs) has never been summed. However, the values are qualitatively consistent with the observed and values for , and . FSS estimates in d are and (19) so . HTSE estimates for various ISG models in d are and in d, (21); (20). Simulation estimates in d and d are consistent with these values (15); (16); (22). Potentially there can also be an analytic correction with exponent (36) in d. If an analytic contribution is indeed present in d in addition to the leading conformal correction, the effective exponent of the combined correction will be slightly increased.

The relative strengths of the leading conformal correction and the analytic correction will depend on the observable so the effective exponent for the combination can be expected to vary somewhat from observable to observable within a narrow range of values.

The leading irrelevant scaling field correction is by definition the correction of this type which has the smallest exponent . From the overall agreement between the -expansion and observed values, one can be fully confident that the leading irrelevant operator correction in d is , and hence no hypothetical correction term with a much smaller exponent (which if it existed would modify the behavior at very large ) can exist. Extrapolations of observed FSS data to infinite using a correction exponent are valid, which justifies the natural straight line extrapolations at in Figs. 1 and 2. Similarly valid fits to the ThL data in Figs. 6 and 7 are made with a leading correction term having an exponent . No hypothetical small exponent correction term producing ”reentrant” behavior exists.

For each dimensionless observable , extrapolation to infinite at gives an estimate for a critical parameter characteristic of the universality class. For the bimodal ISG in dimension four these data show and .

With the FSS rule Eq. (4) for a dimensionless observable the critical exponent can be estimated in principle from . To obtain these estimates, the fits are multi-parameter as they involve simultaneous estimates for , , and the parameters , and as well as . Nevertheless in practice reasonably precise values for can be obtained. From analyses of data sets such as those shown in Figs. 1 and 2 on a number of different dimensionless observables, a global estimate for the bimodal model is (16).

## Iv Thermodynamic derivative peak analysis

Near criticality in a ferromagnet for dimensionless observables the heights of the peaks of the thermodynamic derivatives scale for large as (14); (23)

(5) |

The temperature location of the derivative peak similarly scales as

(6) |

As both and vary as at large , a plot of the peak locations against the inverse peak heights tends linearly to at large . These estimates of are independent of the FSS estimates. The observables used for (14) can be for instance the Binder cumulant or the logarithm of the finite size susceptibility .

For an ISG just the same thermodynamic differential peak methodology can be used as in the ferromagnet. As far as we are aware this analysis has not been used previously in the ISG context. In Fig. 3 we show data for peak height locations against the inverse peak height for the derivatives with dimensionless observables , the Binder cumulant , together with and defined by

(7) |

and

(8) |

where the coefficients have been chosen so that the parameters go from 0 at high temperature to 1 at low temperature. In Fig. 3 the linear extrapolations of the three sets of data points lead consistently to , in full agreement with the FSS data.

From the thermodynamic derivatives, ”The critical exponent can be estimated without any consideration of the critical coupling ” (14). As the exponent is close to in d, we choose to plot data for the normalised Binder cumulant derivative against , Fig. 4. Corrections to scaling are visible for small , but for the slope is , which corresponds to , in full agreement with the FSS estimate above.

## V Thermodynamic limit derivatives

Analyses using the standard RGT scaling variable , either in ferromagnets or in spin glasses, are restricted to the critical region and the FSS regime , because tends to diverge at high temperatures so scaling corrections automatically proliferate outside the critical regime. This is unnecessary. A natural scaling variable, which has been used for the scaling of the susceptibility in ferromagnets for more than 50 years, is (24); (25); (26); (27). The Wegner thermodynamic limit (ThL) susceptibility scaling expression (25), which is expressed in terms of and is valid from criticality to infinite temperature,is

(9) |

At criticality becomes identical to and tends to at infinite temperature, so the corrections in the Wegner expression are well behaved for the entire paramagnetic temperature range. The first correction term in the equation is the leading confluent correction and the second an analytic correction. There is a closure condition as . In ISGs as the interaction energy is , not , and the appropriate scaling variable is or . This is the scaling variable which was used from the beginning in ISG HTSE work (28); (29); (20); (21); (30) but curiously not in most analyses of numerical simulations. The Wegner expression for the susceptibility (25) carries over, with the appropriate definition for .

It has been pointed out that in Ising ferromagnets the ThL expression for the second moment correlation length analogous to the Wegner expression for the susceptibility is (30); (31); (27); (22)

(10) |

The factor arises because the generic infinite temperature limit behavior is .

In ISGs the factor becomes , again because replaces (30), so

(11) |

with the ISG . Following a well-established protocol (32); (27) one can define a temperature dependent effective exponent for the susceptibility

(12) |

and for the second moment correlation length

(13) |

which tend to the critical and respectively at criticality, and to and respectively in the high temperature limit, where is the kurtosis of the interaction distribution. ( in the bimodal case and for the Gaussian).

For data in the ThL regime, the Wegner expression with two correction terms Eq. (9) translates exactly into

(14) |

where the second term is an effective subleading correction. There is an equivalent relation for . Fixing , the and plots for all the different sizes are shown in Fig. 5 and Fig. 6. In Fig. 5 the effective evaluated from an explicit sum of the first 15 terms in the HTSE series of Ref. (21) is also shown.

The fit curves are adjusted to those data which are in the ThL regime, and the fit is then extrapolated to criticality to obtain estimates for the critical and . Just as for in the FSS regime, is the lowest correction exponent and so a further hypothetical correction term with a much lower exponent which could modify the form of the fit curve between the ThL data and criticality can be ruled out. The and plots are however rather sensitive to the exact value of .

The fits shown correspond to an ISG temperature dependent susceptibility effectively fitted by

(15) |

and a second moment correlation length

(16) |

so with estimated critical exponents , , , with a leading correction term and further higher order terms.

The ratio

(17) |

does not involve , so can be measured without any knowledge of . Then the as functions of can be extrapolated to for infinite to estimate the critical value of without involving . Alternatively we plot, see Fig. 7, against . At infinite temperature , and towards criticality for the ThL regime data and . Hence the slope of against tends to as criticality is approached. The fit and limiting slope can be estimated without specifying an explicit value of . A satisfactory fit curve with a single effective correction term is with , and adjustable parameters. The values in the fit shown are , , and .

The critical exponents estimated without any knowledge of are therefore from the thermodynamic derivative result Fig. 4 and , and so combining these values. The values are in excellent agreement with the estimates from the extrapolated ThL derivatives of Fig. 5 and Fig. 6 which validates the estimate and the overall methodology.

## Vi High temperature series expansion

In the HTSE approach, series of exact terms are evaluated, which are summed to obtain the temperature dependence of the susceptibility. In ISGs the series for the spin glass susceptibility is

(18) |

where (29); (20), or the same equation with replaced by (21). The terms are exact. The number of terms which can be calculated is limited by practical considerations, and up to now has been an upper limit in ISGs. The series have been analysed by Dlog-Padé, M1 and M2 techniques (29); (20); (21) which are not particularly transparent for a non-specialist. An alternative classical approach is the ratio method (33); (24); (27). For a model where the temperature dependence of the thermodynamic limit (infinite size) susceptibility is (25)

(19) |

with , the corresponding HTSE terms are (27)

(20) |

The ratios between successive terms are

(21) |

In favorable cases the ratio series plotted as a function of can be extrapolated to infinite to estimate the critical temperature from the intercept, , the critical exponent from the initial slope , and the conformal correction exponent from the leading non-linear term. Unfortunately there can be parasitic oscillating terms in the ratios, arising from anti-ferromagnetic poles (27), particularly in simple hypercubic lattices. It turns out that in the ISGs the parasitic terms in the susceptibility series become rather strong when the dimension drops. Explicit lists of terms for in dimension and below were tabulated in Ref. (29).

For the simple hypercubic bimodal ISGs in general dimension Klein et al. (20) tabulated the terms to order 15 not only for the ISG susceptibility series (called in their nomenclature) but also for the higher order susceptibilities and which they define. These susceptibilities have critical exponents and respectively (20); (27). Unfortunately the parasitic terms are so strong in the d susceptibility series that the ratios cannot be readily exploited so as to make accurate estimates of the critical parameters. This explains why the error bars quoted in Ref. (21) for and in d are rather large. However, we have summed the tabulations of (20) for and to obtain the and for these series, and have found that for these higher order susceptibilities the ratios behave very regularly in dimension (and above), as can be seen in Fig. 8. Expression 21 can be applied to these ratios and reliable extrapolations can be made through which one can estimate and . Optimal fits excluding the small values are

(22) |

for , and

(23) |

for . The joint intercept corresponds to , and the initial slopes correspond to and . The leading correction terms have an exponent , and prefactors with negative sign. By inspection any hypothesis of another term (which would then be the true conformal correction to scaling) having a much lower exponent and with a prefactor of the opposite sign can once again be ruled out. There are no statistical errors as all points are exact; systematic errors arising from the extrapolation are small.

The exponents estimated above from the simulations and so quite independently from the HTSE data correspond to and , and so are fully consistent with the HTSE values.

## Vii Gaussian interaction model

We will present equivalent data for the hypercubic ISG model with Gaussian interactions in dimension four. All the discussions above apply equally well to the Gaussian data. Corrections for the Gaussian model are in general much weaker than for the bimodal model. Unfortunately no HTSE higher order susceptibility data are available for the Gaussian model.

The d Gaussian critical inverse temperature from both the FSS data, Figs. 9 and 10, and the thermodynamic derivative peak locations Fig. 11, is estimated to be or . This value is fully consistent with previous FSS estimates from Binder parameter and correlation length ratio measurements (34); (35) and 0.554(3)(6).

The Gaussian critical parameters for the normalised second moment correlation length and the Binder cumulant from Figs. 9 and 10 are and respectively. The value for the Binder cumulant is close to the estimate given in Ref. (6).

The exponent estimated directly form the Binder cumulant thermodynamic derivative peak heights, Fig. 12, is . This value is fully consistent with the estimate given in Ref. (6).

The plot of against with and , Fig. 13, is fitted with a single effective correction :

(24) |

corresponding to a critical exponent . This value is again fully consistent with but more accurate than the estimate given by (6).

The fit to the ThL effective with a fixed , Fig. 14, corresponds to

(25) |

so with a critical exponent and only a very weak high order correction term. The prefactor for the usual leading correction term appears to be accidentally very close to zero. A similar conclusion was reached from the FSS analysis of the d Gaussian model in (6). The fit to the ThL effective with fixed , Fig. 15, corresponds to

(26) |

so a critical exponent and a correction term with a similar exponent to the exponent of the bimodal model leading correction.

The critical constants and exponents for a d diluted bimodal model Ref. (6) were estimated to be , , , and . Each of the diluted bimodal critical values is similar to that for the d Gaussian model, but the diluted bimodal model and the standard bimodal model studied above have quite different critical properties.

## Viii Conclusion

We compare in Table 1 the bimodal and Gaussian interaction distribution model estimates for each of the critical parameters and critical exponents studied above. For all the parameters and exponents the estimates for the two models are quite different, with in each case differences of the order of 10%.

Parameter | Bimodal | Gaussian |
---|---|---|

We conclude that in d ISGs the critical exponents, like the values of the critical ordering parameters, depend on the form of the interaction distribution. ISGs in dimension four with different interaction distributions are clearly not in the same universality class.

## Ix Acknowledgements

The simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at High Performance Computing Center North (HPC2N) and at Chalmers Centre for Computational Science and Engineering (C3SE).

### References

- K.G. Wilson, Rev. Mod. Phys. 47, 773 (1975).
- S. F. Edwards and P. W. Anderson, J. Phys. F: Met. Phys. 5, 965 (1975).
- R. N. Bhatt and A. P. Young, Phys. Rev. B 37, 3707 (1988).
- H. G. Katzgraber, M. Korner, and A. P. Young, Phys. Rev. B 73, 224432 (2006).
- M. Hasenbusch, A. Pelissetto, and E. Vicari, Phys. Rev. B 78, 214205 (2008).
- T. Jörg and H. G. Katzgraber, Phys. Rev. B 77, 214426 (2008).
- H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Muñoz Sudupe, G. Parisi, and J.J. Ruiz-Lorenzo, Phys. Rev. B 58, 2740 (1998).
- P. Calabrese, V. Martín-Mayor, A. Pelissetto, and E. Vicari, Phys. Rev. E 68, 036136 (2003).
- J. L. Cardy, J.Phys. A Math. Gen., 20, L891 (1987).
- G. Parisi, R. Petronzio, and F. Rosati, Eur. Phys. J. B 21, 605 (2001).
- M. Castellana, Eur. Phys. Lett. 95, 47014 (2011).
- M. C. Angelini, G. Parisi, and F. Ricci-Tersenghi, Phys. Rev. B 87, 134201 (2013).
- H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Muñoz Sudupe, Phys. Letts. B 400, 346 (1997).
- A. M. Ferrenberg and D. P. Landau, Phys. Rev. B 41, 5081 (1991).
- R. A. Baños, L. A. Fernandez, V. Martin-Mayor, and A. P. Young, Phys. Rev. B 86, 134416 (2012).
- P. H. Lundow and I. A. Campbell, arXiv:1402.1991.
- C. de Dominicis, private communication (2004).
- A. J. Bray, private communication (2004).
- M. Baity-Jesi et al., Phys. Rev. B 88, 224416 (2013).
- L. Klein, J. Adler, A. Aharony, A. B. Harris, Y. Meir, Phys. Rev. B 43, 11249 (1991).
- D. Daboul, I. Chang and A. Aharony, Eur. Phys. J. B 41, 231 (2004).
- P. H. Lundow and I. A. Campbell, unpublished.
- M. Weigel and W. Janke, Phys. Rev. Lett. 102, 100601 (2009).
- M. E. Fisher and R. J. Burford, Phys. Rev. 156,583 (1967).
- F. Wegner, Phys. Rev. B 5, 4529 (1972).
- S. Gartenhaus and W. S. McCullough, Phys. Rev. B 88, 11688 (1988).
- P. Butera and M. Comi, Phys. Rev. B 65,144431 (2002).
- R. Fisch and A. B. Harris, Phys. Rev. Lett. 38, 785 (1977).
- R. R. P. Singh and S. Chakravarty, Phys. Rev. Lett. 57, 245 (1986).
- I. A. Campbell, K. Hukushima, and H. Takayama, Phys. Rev. Lett. 97, 117202 (2006).
- I. A. Campbell and P. Butera, Phys. Rev. B 78, 024435 (2008).
- J. Kouvel and M. E. Fisher, Phys. Rev. A 136, 1626 (1964).
- C. Domb and M. F. Sykes, Proc. Phys. Soc. A235, 247 (1956).
- G. Parisi, F. Ricci-Tersenghi, and J. J. Ruiz-Lorenzo, J. Phys. A 29, 7943 (1996).
- M. Ney-Nifle, Phys. Rev. B 57, 492 (1998).
- H. G. Ballesteros, L. A. Fernández, V. Martín-Mayor, A. Muñoz Sudupe, G. Parisi, and J. J. Ruiz-Lorenzo, J. Phys. A 32, 1 (1999).
- E. Gardner, J. Phys. 45, 1755 (1984).