Counting metastable states of Ising spin glasses on hypercubic lattices

We describe how metastable states of Ising spin glasses can be counted by means of Monte Carlo computer simulations. The method is applied to systems defined on hypercubic lattices in one to six dimensions with up to about 103 spins. It is shown that the number of metastable states obtained for different disorder realizations satisfies a log-normal distribution. We investigate the distribution of energies of metastable states by means of moments and cumulants.


Introduction
One major characteristic of complex glassy systems is the number of metastable states or local energy minima. These are states for which any local modification results in an increase in energy. Often the logarithm of the number of these states is an extensive quantity and not easily accessible. For Ising spin systems with random Gaussian interactions metastable states are usually defined as follows: A spin configuration {s} can be considered a metastable state (MSS) at T = 0 if any flip of a single spin would lead to an increase in energy. Such configurations are also called one-spin-flip stable [1]. Some of them are also two-spin-flip stable, three-spin-flip stable, etc. These higher degrees of stability will not be discussed in this study. An intuitive measure for the complexity of the energy landscape of the model is provided by the number N S of such MSSs. In contrast to systems with geometric frustration this number is in general not known for spin glasses and has been derived analytically only for the one-dimensional spin chain [2] and for a system with infinite dimensions [3], i.e., for the Sherrington-Kirkpatrick mean-field model [4]. However, an expansion of the latter result provides estimates for finite dimensions [3] and further analytical approximations exist for the square and the simple cubic lattice [5]. So far, none of these predictions have been computationally verified.
In this paper we apply Monte Carlo simulations to determine the number of metastable states of the Edwards-Anderson model on lattices of one to six dimensions. The rest of the article is organized as follows. We begin with a short review of previous work in Section 2. In Section 3 we discuss the model and the Contribution to the Topical Issue "Recent Advances in the Theory of Disordered Systems", edited by Ferenc Iglói and Heiko Rieger. a e-mail: Stefan.Schnabel@itp.uni-leipzig.de methods we used. Section 4 is dedicated to the presentation and discussion of our results. Finally in Section 5 we wrap up with a brief summary of our main findings and an outlook to future work.

Previous studies
In one dimension, i.e., for a chain of N spins the problem can easily be solved [2]. An MSS is occupied if for each spin the stronger of its two bonds is satisfied. The weaker one might be broken. Therefore, only bonds that are weaker than the two adjacent bonds can be broken in an MSS, since such bonds are the weaker ones for both involved spins. Hence such bonds act as degrees of freedom in the space of metastable states and since in the limit of large systems one in three bonds is weaker than its two neighbors, it follows that the number of metastable states is N I S = 2 N/3 ≈ e 0.23105N . There is also a solution for the limit of infinite dimensions. Bray and Moore [3] have shown that for the Sherrington-Kirkpatrick (SK) model [4] again in the limit of large systems N ∞ S = e 0.199228N by calculating the distribution of MSSs as function of energy. They also provided an expansion of this result in the inverse coordination number of the lattice z −1 , so that estimates for D-dimensional hypercubic lattices are possible. The identification of the maximum of the approximated distribution function allows one to estimate the number of MSSs and their mean energy. 1 The values thus obtained are listed in Table 1. 1 In reference [3] the distribution of MSSs as function of energy is derived by a first-order expansion in z −1 of the exact result for the SK model. The values in Table 1 are obtained from the numerically determined positions of the maxima of these distributions. Expressions for N S and E /(N √ 2D) in reference [3] (Eqs. (6.8) and (6.9)) are derived from a further first-order expansion and lead to slightly different values.  [5] used an analytic approximation in order to estimate the number of metastable states for several graphs, among which they found that for the square lattice N −1 ln[N II S ] approaches 0.21808 (2) and that for the simple cubic lattice the limit of N −1 ln[N III S ] is 0.21125(1).

Model and method
We consider the D-dimensional Edwards-Anderson (EA) model [6] which is defined by the Hamiltonian Here, ij is an edge of a hypercubic lattice with periodic boundary conditions. Ising spins s i ∈ {−1, 1} on the lattice sites interact via bonds J ij , which are randomly chosen according to a normal distribution: leading to disorder and frustration. Samples with different disorder realizations {J} exhibit different behavior, but these differences are expected to decrease with increasing system size. It is thus justified to rely on disorder averages of observables, where M is the total number of realizations {J} randomly drawn from the disorder distribution (2). The investigation of MSSs of the EA model by means of Monte Carlo simulations is challenging since their number is expected to grow exponentially with system size and they are separated by high energy barriers. For small systems they can easily be found using optimization methods like steepest descent or simulated annealing [7], however, the selection of MSSs so-derived will be biased towards states with a large basin of attraction. Recently, we have introduced an algorithm [8,9] that is able to sample MSSs of the EA model uniformly. Here, we will discuss this method only briefly. For details we refer to our recent publication in reference [8]. Starting from the fact that an MSS within a simulation will be the result of some sort of energy minimization process, the idea is to make this process itself the object of a Monte Carlo simulation and to introduce weights that will remove any bias. To this end we define a combined state {{s}, {ξ}}, where {s} is a spin configuration and {ξ} a set of uniformly distributed random numbers. Both components can easily be modified and the statistical properties of their distribution can be controlled without difficulty. This combined state is then interpreted as a random energy minimization that starts with {s} and uses the numbers ξ to randomly select a spin with positive energy which is then flipped until an MSS is reached. In consequence, the resulting sequence of states is a function of {{s}, {ξ}} and -as it turns out -it is possible to assign weights to each sequence and therefore to each combined state such that all MSSs are produced with equal probability.
Importance sampling of short trajectories in state space has been done in the past, albeit not with the ultimate aim of investigating an ensemble of individual states. Examples can be found in the works of Sun [10] and Hartmann, see, e.g., reference [11].
The problem of counting the total number of states of a system is not often encountered. Usually, this is a known quantity; the number of states -metastable or not -for any Ising spin system comprising N spins, for instance, is 2 N . In most cases only the relative statistical weights of sets of states (e.g., states with one particular energy) are derived from a Monte Carlo simulation. In order to measure the total number of states we need a reference set of which the cardinality is known and to which we relate the whole. In reference [8] we measured the distribution of MSSs as a function of energy which involved 'binning', i.e., the separation of the energy scale into intervals of equal width. If the energy resolution is high enough some intervals will contain exactly one (twice degenerated) MSS, a fact that we could exploit to normalize the entire distribution and count the total number of states. However, this requires a very precise determination of the distribution and was only possible for small systems.
In this study, we propose another approach [12] that relies on the Hamming distance to a specific MSS {s ref } instead of energy. Here δ is the Kronecker delta symbol. Since MSSs of low energy possess larger basins of attraction 2 (BOA) and are, therefore, easier to find, we always used a low-energy reference MSS obtained by a quench in the space of MSSs. We then used a combination of flat-histogram techniques to measure the distribution of MSSs g(d H ) as a function of d H . A variant of the Wang-Landau method [13] was applied in order to get a good estimate where c 1 is an unknown normalization constant. The inverse of g 1 was then used as the weight function for multicanoncial sampling [14,15] that was terminated after completing between 10 and 200 round trips depending on size and dimensionality, i.e., the system was required to travel between 10 and 200 times from d H = 0 to d H = N and back. In this process, a histogram was recorded, with {s} t being the sequence of MSSs generated by the simulation. As the next step, an improved estimate for the distribution was obtained and used for the further analysis. For an individual disorder realization in four dimensions with L = 6 the result is displayed in Figure 1.
and, therefore, we know that g(0) = g(N ) = 1. Since we use Monte Carlo simulations we cannot expect perfect agreement for both values. Instead, one obtains differences that depend on the number of completed round trips. We use the arithmetic average of both values for normalization, which fixes the constant c 2 in (7). The total number of MSSs is obtained by summation, For instance, the system for which the distribution of MSSs is displayed in Figure 1 possesses about N S ≈ 2 × 10 117 MSSs. For all combinations of dimensionality and system size 1000 disorder realizations were generated. During the simulation we also recorded the moments of energy as a function of the Hamming distance, .
The values can be recombined to form the moments of the unbiased distribution of MSSs: for each individual disorder realization. Sampling the space of MSSs one has to be cautious. During test runs for D = 3 we obtained different results (e.g., for average energy) depending on whether straightforward unbiased (simple) sampling or the weighted ensemble described above was used. This means that in combination with at least one of the two ensembles the algorithm is unable to equilibrate. This effect is most likely caused by the different sizes of the BOAs which are much larger if the respective MSS is of low energy. Since for unbiased sampling all individual MSSs are assigned the same weight, the probability to occupy a particular combined state {{s}, {ξ}} within a BOA is inversely proportional to the size of that BOA. Combined states in large BOAs are, therefore, suppressed and such BOAs might become barriers preventing equilibration. Using a quenched reference MSS and a flat distribution of d H as in our weighted ensemble improves the performance in comparison to simple sampling since it increases the weight of some large BOAs.
In reference [8] we measured the distribution of MSSs as a function of energy for D = 3 and where able to directly verify the accuracy of our results for N = 4 3 . We then found that the distributions we obtained depend -when properly rescaled -very little on the size of the system. We see this as evidence that the algorithm is working correctly also for larger sizes. A failure to properly sample the state of MSSs should become more severe with increasing N and an agreement like the one observed would be very unlikely. As we will discuss below values for the average energy obtained in this work agree very well with our older data. Hence, even though we are using a different ensemble, we are confident that the results presented here for D = 3 are valid.
In addition to the tests provided in reference [8] we investigated a spin chain (D = 1) with N = 100. In the spin chain weak bonds, i.e., bonds with an absolute value |J i,i+1 | smaller than those of their neighbors |J i−1,i | and |J i+1,i+2 |, act as the degrees of freedom since only they can be broken within an MSS. While periodic boundaries remove one of these degrees, there is a twofold degeneracy of each bond configuration due to the trivial symmetry {s} ↔ {−s} such that the total number of MSSs is given by two to the power of the number of weak bonds. On the one hand this means that for any given disorder realization log 2 N S is an integer while on the other hand for the disorder average one obtains [log 2 N S ] = N/3. Accordingly, we find that for this system the numerically determined values of the binary logarithm of N S shown in Figure 2 are very close to integer values and [log 2 N I S ]/N = 0.3326(7) or [ln N I S ]/N = 0.2305(5). We draw further confidence in our method for D = 2 by comparing with the analytic approximation by Waclaw and Burda [5] that will be discussed below. Unfortunately, we have no means to independently test our results for D = 4, 5, 6. In the unlikely event that relevant parts of the space of MSSs are not sampled, using the technique described above, one would obtain an estimate for the number N S of MSSs that is too small. That means that our values can at the very least serve as lower bounds.

Results and discussion
We show the results for the distribution of the normalized logarithmic number of states N −1 ln N s in Table 2. As expected we observe that the mean number of MSSs grows exponentially with system size. Only relatively small finite-size effects occur as indicated by the flatness of the curves shown in Figure 3. There, we have also displayed the values of Table 1 that are predicted by the approximation of Bray and Moore [3]. Since the latter is an expansion from the SK model in 1/z = 1/2D, the deviations should decrease with increasing D. It seems, however, that this behavior does not clearly manifest itself for D ≤ 6. In order to make the scaling of the width of the distributions more apparent the standard deviations divided by √ N are shown in Figure 4. The higher-order characteristics skewness and kurtosis are defined as and with the central moments For the distribution of the logarithm of the number of states these values are also provided in Table 2.
The statistical errors of both quantities depend solely on the number of samples M which is identical for all cases (see above). For a Gaussian distribution, the skewness vanishes, S Gauss = 0, and the kurtosis becomes K Gauss = 3. We hence find that our results are consistent with Gaussian distributions for all data sets. The mean values [ln N S ]/N in Table 2 can not directly be compared with the approximations from reference [5]. There, the logarithm of the disorder average was taken while so far we have considered the average of the logarithm. If we presume that ln N S is indeed normally distributed, it follows that and two ways of estimating ln[N S ] are available. The disorder average [N S ] can be taken directly from data or equation (15) can be applied to the values in Table 2.
The respective values are shown in Table 3. Both techniques provide results that agree within statistical errors, which provides additional evidence that the distributions of ln N S and N S are normal and log-normal, respectively. Comparing with the results from reference [5] we find small but systematic deviations. For D = 2 and D = 3 our values for the largest system considered slightly exceed the values presented there: d ln[N II S ]/dN = 0.21808(2) and d ln[N III S ]/dN = 0.21125(1). Since our values tend to increase with system size, it is unlikely that the smaller values from the analytical approximation will be approached in the limit of large N . We point out again that while it is theoretically possible that our results underestimate the exact value, it is hard to see how our estimates could be too large.
As mentioned above the moments of energy of the MSSs E l were calculated for every disorder realization. In this context . . . stands for an unweighted average over all MSSs, i.e., the energies of all MSSs (and only of MSSs) contribute with equal weight. Formally this corresponds to infinite temperature in the ensemble of MSSs. Here, central moments and cumulants characterize the distribution Table 3. Estimates for ln[N S ]/N allowing comparisons with the values 0.21808(2) for D = 2 and 0.21125(1) for D = 3 from reference [5]. of the energies of MSSs for individual disorder realizations and now become subject to disorder averages themselves. It is convenient to use the rescaled energy where z = 2D is the number of nearest neighbors of a spin. Bray and Moore [3] have shown that for the SK model in the limit of large N the average value is = −0.50605. In our recent study [8] we sampled the distribution of MSSs as a function of energy for a simple-cubic lattice directly and found that it is approximated by  Table 4, agree very well. We find that the mean energies converge very well for D = 2, 3, 4 and that the scaling of the standard deviations with increasing system size is also very regular as can be seen in Figure 5. The standard deviations for the three-dimensional system are also in particularly good agreement with the expected value. Using mean values and standard deviations we can attempt to approximate the density of MSSs as a function of energy analogously to equation (17) for D = 2, 4, 5, 6: Comparison of our mean energies [ ] to the values derived from the approximation in reference [3] listed in Table 1 we observe a curious offset: For dimensions D = 3, 4, 5, 6 the former roughly agree with the latter for dimensions D = 2, 3, 4, 5. We see no reason to assume that this is anything but coincidental, though it should be noted that (D + 1) −1 ≈ D −1 with corrections only of order D −2 . The behavior of the higher-order cumulants is not so clear due to the larger errors. For D = 3 assuming (17) to be valid we expect the skewness to vanish like [S ] ∝ N −1/2 and the kurtosis to approach the limit according to [K ] − 3 ∝ N −1 . The scaling behavior of [S ] is displayed in Figure 6 together with fits of the functional form f (N ) = α/ √ N . The systems are still comparatively small and the data are not precise enough to be certain, but we can state that they are consistent with the expected behavior. Finally one might ask the question whether the number of MSSs and their average energies are correlated. This might be expected since both quantities depend on the amount of frustration in the system. However, as can be inspected in Figure 7, we find that there is almost no relation between the two.

Conclusion
We have for the second time employed an algorithm that is designed to perform biased or unbiased random walks in the space of metastable states of the Edwards-Anderson model. Previously, we opted for a more traditional ensemble with the energy as the control parameter [8]. While this allows us to measure the unnormalized distribution of states as function of E, the total number of metastable states is only indirectly accessible and only for small systems. In this study we used a complementary approach. By biasing the ensemble in dependence of the Hamming distance to a reference state, the walker encounters fewer barriers, an easy way of normalization is at hand, and the number of states can reliably be determined not just in two and three, but also in four, five, and six dimensions.
By analyzing cumulants of the distribution we find that for all dimensions and system sizes the logarithmic number of metastable states seems to be normally distributed. Comparing the number of metastable states with the values obtained from the expansion done by Bray and Moore [3] and those derived in an analytical approximation by Waclaw and Burda [5] which are not consistent, we found that our results systematically deviate from both. Except for D = 2 where we find less metastable states than predicted by Bray and Moore [3] our values for [N S ] are above the predictions, although the deviations to [5] are almost insignificant. This is important because even if we assume that the algorithm is unable to find all relevant states the method would merely underestimate the number of metastable states such that the difference between the true values and the predictions becomes even greater. Due to the monotonic increase as function of system size in our data, we also deem it unlikely that this is a finite-size effect and that the analytical approximations will be realized in the thermodynamic limit.
Disorder averages of mean values and standard deviation of the energy of metastable states in three dimensions agree remarkably well with our previous results [8]. We find that disorder averages of skewness and kurtosis approach zero and three, respectively, so that we may conclude that also the energies of the unbiased metastable states are close to normally distributed.
As a future project it would be interesting to modify the bond disorder and study the transition from the spin glass to a ferromagnet.