# 1 Introduction

## 1 Introduction

The intuitive picture for spin glasses and other systems with conflicting constraints, for reviews see , is that there exists a large number of degenerate thermodynamic states with the same macroscopic properties but with different microscopic configurations. These states are separated by free-energy barriers in phase space, caused by disorder and frustration. Experimentally this is supported by long characteristic times found in phenomena like the measurement of the remanent magnetization in typical spin glasses such as, e.g., $`(\mathrm{Fe}_{0.15}\mathrm{Ni}_{0.85})_{75}\mathrm{P}_{16}\mathrm{B}_6\mathrm{Al}_3`$ . However, one difficulty of the theory of spin glasses is to give a precise meaning to this classification: No explicit order parameter exists which allows to exhibit the barriers. The way out of this problem appears to use an implicit parametrization, the Parisi overlap parameter $`q`$, which allows to visualize at least some of them . Calculations of the thus encountered barriers in $`q`$ are of major interest. For instance, it is unclear whether the degenerate thermodynamic states are separated by infinite barriers or whether this is just an artifact of mean-field theory.

Before performing numerical calculations of these barriers, one of the questions which ought to be addressed is “What are suitable weight factors for the problem?” The weight factor of canonical Monte Carlo (MC) simulations is $`\mathrm{exp}(\beta E)`$, where $`E`$ is the energy of the configuration to be updated and $`\beta `$ is the inverse temperature in natural units. The Metropolis and other MC methods generate canonical configurations through a Markov process. However, by their very definition the rare-event states associated with the free-energy barriers are suppressed in such an ensemble. Now, it became widely recognized in recent years that MC simulations with a-priori unknown weight factors, like for instance the inverse spectral density $`1/n(E)`$, are also feasible and deserve to be considered, for reviews see . Along such lines progress has been made by exploring \[4–7\] innovative weighting methods for the spin glass problem.

The main idea of the studies \[4–7\] is to avoid getting stuck in metastable low-energy states by using a Markov process which samples the ordered as well as the disordered regions of configuration space in one run. Refreshing the system in the disordered phase clearly benefits the simulations, but the performance has remained below early expectation. One reason appears to be that the direct (i.e. ignoring the dynamics of the system) barrier weights are not affected, such that the simulation slows down due to the tree-like structure of the low-energy spin glass states, see Ref.  for a detailed discussion.

In Ref.  two of the authors introduced a method which focuses directly on enhancing the probability for sampling the barrier regions in the Parisi overlap parameter. Relying on this method, we have meanwhile performed large-scale simulations of the three-dimensional Edwards-Anderson Ising spin glass. First results from this novel investigation are reported in the following. The next section briefly introduces the spin glass model and is followed by an introduction to the method of  in section 3. Technical details of the computer simulation are described in section 4 and some of the obtained physical results are summarized in section 5, followed by conclusions in the final section 6.

## 2 Model

We focus on the three-dimensional ($`3d`$) Edwards-Anderson Ising (EAI) spin glass on a simple cubic lattice of size $`N=L^3`$ with periodic boundary conditions. It is widely considered to be the simplest model to exhibit realistic spin glass behavior and has been the testing ground of Refs. \[4–7\]. The energy is given by

$$E=\underset{ik}{}J_{ik}s_is_k,$$
(1)   
where the sum is over nearest-neighbour sites and the Ising spins $`s_i`$ and $`s_k`$ take values $`\pm 1`$. The exchange coupling constants $`J_{ik}`$ are quenched random variables which are chosen to be $`\pm 1`$ with equal probability. Each fixed assignment of the exchange coupling constants $`J_{ik}`$ defines a realization of the system, and all physical results refer to an average over many such realizations. Early MC simulations of the $`3d`$ EAI model located the freezing temperature at $`\beta _c0.9`$. For a concise review, see Ref. . Recent, very high-statistics canonical simulations  estimate $`\beta _c=0.901\pm 0.034`$, and were interpreted  to improve the evidence in support of a second-order phase transition at $`\beta _c`$.

The Parisi overlap parameter is defined as 

$$q=\frac{1}{N}\underset{i=1}{\overset{N}{}}s_i^1s_i^2,$$
(2)   
where the spins $`s_i^1=\pm 1`$ and $`s_i^2=\pm 1`$ correspond to two independent copies (replica) of the same realization (defined by its couplings $`J_{ik}`$), each with its own time evolution in the MC simulation (realized by different random numbers).

## 3 Multi-overlap algorithm

The method of Ref.  simulates two replica of the same realization in one computer run. In another context this has before been done in . Our basic observation, closely related to multicanonical methods , is that one does still control canonical expectation values at temperature $`\beta ^1`$ when one simulates with a weight function

$$w(q)=\mathrm{exp}\left[\beta \underset{ik}{}J_{ik}(s_i^1s_k^1+s_i^2s_k^2)+S(q)\right].$$
(3)   
Of particular interest is to determine $`S(q)`$ recursively such that the $`q`$-distribution becomes uniform in $`q`$ (“multi-overlap”), and the interpretation of $`S(q)`$ being the microcanonical entropy of the Parisi order parameter. Hence, although an explicit order parameter does not exist, an approach very similar to the multimagnetical  (which is an highly efficient way to sample interface barriers for ferromagnets) exists herewith.

As a measure of the performance of the algorithm we recorded the average number of sweeps that are necessary to perform tunneling events of the form

$$(q=0)(q=\pm 1)\mathrm{and}\mathrm{back}.$$

In the following we refer to this quantity as autocorrelation time.

## 4 Simulation

The multi-overlap algorithm is thus particularly designed for simulations of the interesting region below the freezing temperature where $`P_i(q)`$, the canonical probability density of $`q`$ for realization $`i`$ (additional dependence on lattice size and temperature is implicit), exhibits a rather complicated behavior. The shape of $`P_i(q)`$ ranges from a simple double-peak structure to involved structures of several minima and maxima. A few examples encountered in the simulations of the $`12^3`$ system are shown in Fig. 1. This is the situation which is notoriously difficult to study with standard algorithms. In our EAI study we therefore focussed on simulations at $`\beta =1>\beta _c0.9`$. We investigated lattices of size $`N=L^3`$ with $`L=4`$, 6, 8, and 12, using T3E parallel computers in Grenoble, Berlin, and Jülich. We simulated 4096 different realizations for the smaller systems up to $`L=8`$, and 512 for $`L=12`$, using the Marsaglia random number generator RANMAR  for drawing the $`J_{ik}`$. As a check we also simulated a second set of 4096 realizations for $`L=4`$, 6, and 8, where the $`J_{ik}`$ were generated with a more elaborate version of RANMAR: RANLUX  with luxury level 4. In the simulations themselves we always employed the RANMAR generator due to CPU time considerations.

For each realization the simulation consisted of three steps:

1. Construction of the weight function (3). Here we employed an improved variant of the accumulative iteration scheme discussed in Ref. . Details will be published elsewhere. The iteration was stopped after at least 10 tunnelings (between $`420`$ for $`L=12`$) occurred.
2. Equilibration run. This run was to equilibrate the system for given fixed weight factors.
3. Production run. Each production run of data taking was concluded after at least 20 tunneling events were recorded. To allow for easy (standard) reweighting in temperature we stored besides the Parisi overlap parameter also the energies and magnetizations of the two replica in a time-series file. Since the corresponding autocorrelation times $`\tau `$ vary significantly from realization to realization we used an adaptive data compression routine to make sure that the spacing between measurements was always proportional to $`\tau `$. Using this adaptive scheme we recorded for each realization 65536 measurements.

Let us conclude this section with a few remarks on the performance of the new algorithm. Fitting the estimates of the mean autocorrelation time $`[\tau ]_{\mathrm{av}}`$ to the form $`\mathrm{ln}([\tau ]_{\mathrm{av}})=a+z\mathrm{ln}(N)`$ gives $`z=2.42\pm 0.03`$. As will be discussed below, the implied improvement with respect to barrier calculations is huge. Nevertheless the slowing down is quite off from the theoretical optimum, which is $`z=1`$ for multicanonical simulations . One reason seems to be that we want the $`q`$-distribution to be uniform in the whole admissible interval $`q[1,+1]`$, including the region of $`|q|1`$ that is strongly correlated with ground states and hence difficult to reach by local updates, see for instance . Being content with a smaller region (like the region between the two outmost maxima of $`P_i(q)`$) is expected to give further improvements of the tunneling performance. By also monitoring the encountered minimum, maximum and median tunneling times we observed that the mean values are systematically larger than the median, what means that the tunneling distribution has a rather long tail towards large tunneling times. On the other hand, the effect is not severely hindering our multi-overlap simulations: For the lattice sizes $`L=4`$ to 8 the worst behaved realization took never more than 1% of the entire computer time. Indicating a remarkable increase of complexity, this amount was about 12% for $`L=12`$.

## 5 Results

The data created in this way allows us to calculate a number of physically interesting quantities. Let us first concentrate on the canonical potential barriers in $`q`$ which were in Ref.  defined as

$$B_i=\underset{q=1}{\overset{\mathrm{}q}{}}\mathrm{max}[1,P_i(q)/P_i(q+\mathrm{}q)],$$
(4)   
where $`\mathrm{}q`$ is the step-size in $`q`$. For the double-peak situations of first-order phase transitions  Eq. (4) simplifies to $`B_i=P_i^{\mathrm{max}}/P_i^{\mathrm{min}}`$, where $`P_i^{\mathrm{max}}`$ is the absolute maximum and $`P_i^{\mathrm{min}}`$ is the absolute minimum (for ferromagnets at $`q=0`$) of the probability density $`P_i(q)`$. Our definition generalizes to the situation where several minima and maxima occur due to disorder and frustration. When evaluating (4) from numerical data for $`P_i(q)`$ some care is needed to avoid contributions from statistical fluctuations of $`P_i(q)`$.

Graphically, our values for the $`B_i`$ are presented in Fig. 2. It comes as a surprise that the finite-size dependence of the distributions is very weak. To study this issue further, we have compiled in Table 1 for each lattice size the following informations about our potential barrier results: largest and second largest values $`B_{\mathrm{max}}`$ and $`B_2`$, median values $`B_{\mathrm{med}}`$ and mean values $`[B]_{\mathrm{av}}`$ with their statistical error bars. From this table it becomes obvious, why this investigation could not be performed using canonical methods to which in this context also multicanonical simulations and enlarged ensembles belong, as their weights for those barriers are still canonical. For these methods the slowing down is proportional to the average barrier height $`[B]_{\mathrm{av}}`$, which is already large for $`L=4`$, about 17 thousand, and increases to about 250 thousand for $`L=8`$. The subsequent decrease to about 150 thousand for $`L=12`$ has to be attributed to the smaller number of realizations in this case: Comparison of the extrem values makes only sense when the numbers of realizations match. On $`L=4`$ and $`L=6`$ systems we have performed a number of (very long) canonical simulations to estimate the proportionality constant between barrier height and improvement due to our multi-$`q`$ simulations. Using these result, we estimate that with our computer program a canonical MC calculation of the worst $`L=12`$ barrier alone, the one reported in Table 1, would take about 1000 years on a 500 MHz Alpha processor.

The reader may be puzzled by the very large error bars assigned to the mean barrier values. Their explanation is: The entire mean value is dominated by the largest barrier, which contributes between 41% $`(L=6)`$ and 77% $`(L=8)`$, see $`B_{\mathrm{max}}`$ in the second column of Table 1. Besides $`B_{\mathrm{max}}`$, the second largest value $`B_2`$ is listed in the third column. It may be remarked that most of these worst case barriers exhibit simple double-peak behavior. One lesson from these numbers is that very few of the realizations are responsible for the collapse of canonical simulation methods.

Typical realizations, described by the median results of Table 1, have much smaller tunneling barriers. They turn out to be quite insensitive to the lattice size. This result of an almost constant typical tunneling barrier is consistent with the fact that our tunneling times are rather far apart from their theoretical optimum: Still other reasons than overlap barriers have to be responsible. Therefore, one may question the apparently accepted opinion that these are the typical barriers which are primarily responsible for the severe slowing down of canonical MC simulations.

Even though less detailed than the barrier data also the averaged canonical probability densities $`P(q)=[P_i(q)]_{\mathrm{av}}`$ contain important information on the system. Due to the average over all realizations the $`P(q)`$ exhibit a double-peak shape; see Figs. 3-6. In principle the data stored in the time-series files allow canonical reweighting in a $`\beta `$-range around the simulation point. For disordered systems, however, where many realizations have to be reweighted, great care is needed when estimating the valid reweighting range.

By analyzing the spin glass susceptibility, $`\chi _{\mathrm{SG}}=N[q^2]_{\mathrm{av}}`$, we obtained in Ref.  the best finite-size scaling fit $`\chi _{\mathrm{SG}}L^{\gamma /\nu }`$ at $`\beta =0.88`$ with $`\gamma /\nu =2.37(4)`$ and a goodness-of-fit parameter $`Q=0.25`$. This was corroborated by the curves of the Binder parameter, $`g=(1/2)(3[q^4]_{\mathrm{av}}/[q^2]_{\mathrm{av}}^2)`$, which merge around $`\beta =0.89`$. In the low-temperature phase ($`\beta >\beta _c`$) the curves for different lattice sizes seem to fall on top of each other, but our error bars

were still too large to draw a firm conclusion from this quantity.

Close to the transition temperature $`T_c`$ one expects that the probability densities $`P(q)`$ for different lattice sizes satisfy the finite-size scaling relation

$$P(q)=L^{\beta /\nu }\widehat{P}(L^{\beta /\nu }q,L^{1/\nu }(TT_c)),$$
(5)   
where $`\widehat{P}`$ is a scaling function, and $`\beta `$ and $`\nu `$ are the critical exponents of the order parameter and correlation length, respectively. They are related to $`\gamma /\nu =2\eta `$ by the standard scaling relation $`\beta /\nu =(d\gamma /\nu )/2`$. Using our estimate for $`\gamma /\nu 2.37`$ we thus obtain the value $`\beta /\nu 0.317`$. Figure 3 shows the probability densities $`P(q)`$ of the large-scale simulations, reweighted to the transition temperature $`T_c`$ and rescaled according to (5), using the above exponent estimates. We see that the data for $`L=4`$, 6, and 8 fall almost perfectly onto a common master curve, while the $`L=12`$ data obviously show some deviations, in particular close to the peak. The overall appearance of Fig. 3, therefore, is actually worse than that of the corresponding plot in Ref.  which is based on a much smaller exploratory data set. We strongly suspect that the reason for the deviations of the $`L=12`$ curve can be traced back to problems with the admissible reweighting range for our largest system size. Since for quenched, disordered systems the reweighting range for averaged quantities depends on the common overlap of the reweighting ranges for the individual realizations, it is indeed conceivable that this problem does show up more dramatically for the larger sample of realizations. Given the very good data collapse for the smaller lattices, and with this explanation in mind, we feel that also our new results are compatible with the findings of Ref. . We are currently trying to improve the test of the finite-size scaling prediction (5) by redoing the simulations closer to $`\beta _c`$. Narrowing the $`q`$-range may allow to simulate lattices of size $`L=16`$ and beyond.

The observation that the Binder parameter curves do not splay out at low temperatures suggests that, similar to the 2d XY model, the correlation length $`\xi `$ may be infinite – or at least larger than the simulated lattice sizes. Generalizing the second argument of $`\widehat{P}`$ to $`L/\xi `$ and assuming $`L/\xi =0`$, one thus expects that the $`P(q)`$ should scale also below $`T_c`$. With our simulation set-up, the most reliable test of this conjecture can be done at the simulation point $`T=10.88T_c`$, since then no temperature reweighting is involved. The raw data for $`P(q)`$ at $`T=1`$ are shown in Fig. 5. Notice that $`P(0)`$ is slightly growing with increasing system size. By adjusting the only free parameter, $`\beta /\nu =0.255`$, we obtain the finite-size scaling plot in Fig. 5 which shows a very clear data collapse onto a single master curve for all lattice sizes. Moreover, if we reweight our data to the even lower temperatures $`T=1/1.1`$ (= $`0.80T_c`$) and $`T=1/1.2`$ (= $`0.73T_c`$) we still find

a reasonable data collapse, see Fig. 6, but here again extreme care is needed to control the reliable reweighting range.

Of course, since our lattice sizes are relatively small, we cannot conclude that the correlation length is infinite below $`T_c`$. If $`\xi `$ is large but finite, it is conceivable that we observe an effective scaling behavior as long as $`\xi >L`$. We may conclude, that the correlation length is unusually large ($`\xi >12`$) down to $`0.73T_c`$.

## 6 Conclusions

For the $`3d`$ EAI model at $`\beta =1`$ we have performed a high-statistics calculation for probability distributions depending on the Parisi overlap parameter $`q`$. The results for free-energy barriers in $`q`$ became feasible by using $`q`$-dependent (multi-overlap) weight factors in large-scale MC simulations. Although the tunneling performance is not optimal, the method opens new horizons for spin glass simulations. Using slight modifications of the method (like narrowing the $`q`$-range, including a magnetic field, etc.) will allow us to extend our investigation into various interesting directions, like an improved study of the thermodynamic limit at and below the freezing point, a study of the 4d EAI spin glass model, or $`ϵ`$-physics, where one studies the influence of an interaction term  $`ϵ_{i=1}^Ns_i^1s_i^2=ϵNq`$ in the Hamiltonian (1). In multi-overlap simulations we can obtain expectation values for arbitrary $`ϵ`$-values by reweighting. Physically most interesting is the case where a non-zero magnetic field is combined with a non-zero $`ϵ`$-value.   
W.J. thanks K. Binder for useful discussions and acknowledges support from the Deutsche Forschungsgemeinschaft (DFG) through a Heisenberg Fellowship. The numerical simulations were performed on T3E computers of CEA Grenoble, ZIB Berlin, and HLRZ Jülich. We thank all institutions for their generous support. This research was partially funded by the Department of Energy under Contracts No. DE-FG02-97ER41022 and DE-FG05-85ER2500.