# Counting metastable states of Ising spin glasses on arbitrary graphs

## Abstract

Using a field-theoretical representation of the Tanaka-Edwards integral [1] we develop a method to systematically compute the number of 1-spin-stable states (local energy minima) of a glassy Ising system with nearest-neighbor interactions and random Gaussian couplings on an arbitrary graph. In particular, we use this method to determine for -regular random graphs and -dimensional regular lattices for . The method works also for other graphs. Excellent accuracy of the results allows us to observe that the number of local energy minima depends mainly on local properties of the graph on which the spin glass is defined.

## 1Introduction

Glassy systems have non-trivial energy landscapes just as many complex systems observed in nature. The main characteristics of such landscapes is the number of local minima, called also metastable states. Typically, this number grows exponentially with the system size : . The rate of the exponential growth is a fundamental quantity characterizing the complexity of the system. It is however very difficult to calculate, and it has been analytically found in only a few cases: for the Ising model with random Gaussian interactions on a complete graph (SK-model) [1], on the one-dimensional closed chain [3] and for the Ising model with random binary interactions on -regular random graphs [5].

In this paper we describe a method to determine for the Ising model with random Gaussian interactions on an arbitrary graph. In particular we use this method to compute for some -regular graphs. The idea is to express the number of metastable states in terms of the Tanaka-Edwards integral [1] and then to treat this expression as the partition function of a certain statistical field theory. The logarithm of the partition function can be represented as a sum of connected Feynman diagrams which in turn can be generated and summed on a computer, up to a certain order of the perturbative series. Using some general properties of this series we are able to estimate the value already from the first few orders with a very good accuracy. This allows us to observe that the values of are very similar for regular graphs with different topologies, indicating that the number of metastable states depends mainly on local properties of the graph.

## 2Derivation of the statistical field theory

We consider a system of Ising spins, , , residing on nodes of a simple graph described by an adjacency matrix . The graph does not need to be connected. The matrix is an symmetric matrix with if and are connected by an edge or otherwise. The energy of the system is given by:

where the coupling constants ’s are random numbers taken from the standardized Gaussian distribution with zero mean and unit variance.

We are interested in counting the number of local minima of the energy (Equation 1). A local (metastable) minimum is defined as a configuration of spins such that a flip of any single spin increases energy. Such a configuration is also called 1-spin-stable. The number of 1-spin-stable states is given by [1]

where is a step (Heaviside) function. In Ref. [1] it was shown that the averaging over Gaussian couplings leads to the following concise formula:

where

and

Here is called degree of node . The integral (Equation 2) was calculated [1] in the limit of for a complete graph: , yielding the result

with , known also from earlier considerations of the SK model [2].

In this paper we shall propose a systematic method to evaluate this integral also for other graphs. Let us introduce an auxiliary constant to Eq. (Equation 2):

The idea is now to find a systematic way of expanding (Equation 4) in powers of and then to use this series expansion to estimate its value for :

Borrowing some techniques from field theory let us define the following generating function:

which allows for rewriting Eq. (Equation 4) as:

The function has a closed form:

where is the error function and are cumulants of . The coefficients can be easily calculated up to an arbitrary order using a program for symbolic calculations: , , , . Equation (Equation 6) can be graphically represented as a sum of vacuum Feynman diagrams of a field theory with the propagator and -vertices with coupling constants . The Feynman rules to calculate the contribution of a diagram are as follows. To each vertex , at which lines meet, we ascribe a factor . The subscript means that the vertex is decorated by an index which can be thought of as a color taken from a palette of possible colors. A line joining two vertices decorated with colors and contributes a factor . Additionally each diagram has a certain symmetry factor which depends on the shape of the diagram. Finally, one needs to perform the summation over colors.

As usual, the logarithm of Eq. ,

contains only the contribution from connected diagrams. We shall see below that is an extensive quantity in . Thus it is convenient to introduce a density per spin which for large becomes a function of only and can be represented as a power series:

Our goal is to determine , which gives the rate of exponential growth of the number of 1-spin-stable configurations for large .

The coefficients at in the series expansion (Equation 8) come from connected Feynman diagrams with links. In the general case, they must be summed on a computer, because their number grows very fast with . We need a systematic procedure which allows us to sum diagrams in order to calculate the coefficients . In our approach, such a procedure looks as follows: (A) “draw” all possible connected Feynman diagrams with links; (B) calculate their symmetry factors ; (C) decorate each vertex of the diagram with an index ; (D) for each decoration calculate the contribution of the diagram to as:

where is the number of vertices of the diagram, the first product goes over all vertices and the second one over all links of the diagram; (E) sum contributions of all diagrams and decorations.

Clearly, the procedure described above is not efficient if the original graph is sparse since then many propagators are zero. Therefore, many of possible decorations give a zero contribution and one wastes time summing many zeros. In particular, all decorations of a diagram having a self-connecting link give no contribution since . Thus one can omit such diagrams in the sum. One can improve the step (C) of the procedure by concentrating only on decorations which potentially have a chance to contribute. In other words, one should look only for decorations for which propagators do not vanish. This means that and must be neighbors on the graph on which the spins reside. Therefore, instead of the step (C) one should take the step (C’) in which one only checks such decorations which are consistent with the graph structure. This can be done iteratively. First, one assigns a label to one vertex of the diagram. Then one assigns to its neighbors only values such that or equivalently that and are neighbors on the original graph. One repeats this process for neighbors of ’s and so on, and selects only those labels for which the propagator does not vanish. This speeds up the step (E) of the procedure since now the number of decorations is of order , where is the average node degree of the graph.

The procedure (A-C’-E) works for any graph. For a -regular graph one can essentially simplify calculations since in this case , as it stems from Eq. (Equation 3), and thus the elementary contribution (Equation 9) to is

The last product of ’s is either zero or one. Summing over all decorations of a given diagram we get some number of decorations consistent with the graph structure. Each Feynman diagram has now the following contribution:

Actually is the only part of the expression which depends on the graph structure on which spins reside, other quantities can be calculated beforehand. Therefore, we used a C++ program to generate all simple diagrams (without multiple- and self-connections) up to and applied the routine `nauty`

[6] for isomorphism testing to distinguish different diagrams and to determine their symmetry factors . In figure ? we show several examples of diagrams for . Their number grows very fast with : for there are such diagrams, for – and for – . The number of all multi-diagrams, that is diagrams with multiple connections, is much larger but luckily there is no need to generate them.

Any multi-diagram can be obtained from a simple diagram by replacing its links by multiple links with a certain multiplicity . The contribution of all multi-diagrams associated with a simple diagram is:

where and are calculated for . The terms in the sum contribute to the order of the expansion of , and is the total number of links of the multi-diagram. Factors are corrections to the symmetry factor which arise from the fact that one can permute all multiple links joining two vertices without changing the diagram. is the number of links meeting at vertex of the multi-diagram.

## 3A few examples

Let us illustrate how the method works for spins on a graph with vertices. We have , and . As an example we shall calculate the contribution of the linear diagram with links from Fig. ?. Its total contribution is

where the first term accounts for the symmetry factor, the second one is a product of four propagators, the third one is a product of couplings and the fourth one is the combinatorial factor depending on the details of the underlying graph’s structure encoded in the adjacency matrix . There are only two possible assignments of labels and to this diagram, see Fig. ?. All other terms vanish and hence .

This example is somewhat artificial because the number of nodes is small. Now let us suppose that we have a graph where pairs of nodes and are connected, so that the only non-vanishing elements of the adjacency matrix are or . It is a -regular graph. One can easily see that the number of decorations of any Feynman diagram is or . Indeed, vertices of a given diagram can be alternately decorated with two consecutive numbers and . Thus the choice of a single label specifies automatically the whole decoration. Because this label may assume values, we have . But if the diagram has a loop of odd length, one cannot alternately decorate vertices along this loop, so in this case .

Because , one can see that also is proportional to and hence also . The proportionality coefficient can be determined in this case analytically by a straightforward calculation of the integral (Equation 2):

with . We shall use this explicit result to test our method. Using the formula (Equation 11) and setting we find that the first coefficient in the expansion of comes from only one diagram, a line with from Fig. ?, and reads:

The second coefficient is a sum of the previous diagram with doubled line, and the one for :

On the other hand, we can calculate s numerically from the analytic formula (Equation 13):

We see that both and agree perfectly with those obtained before. Higher coefficients can be calculated by performing the sum (Equation 11) on a computer. We checked that all coefficients, up to , agree with those from Eq. (Equation 13).

In the next section we shall make another cross-check by comparing our results with those for the celebrated SK model. We shall see that the method produces correct values of the coefficients as well as of the limiting value .

## 4The SK model

The SK model [7] is the spin glass (Equation 1) on a complete graph, where each node is connected to all other nodes. Before we apply our procedure of diagrams’ summation, let us recall what can be calculated for the SK model using another methods. The integral (Equation 4) can be done in the thermodynamic limit [1], leading to

where is a solution to a saddle-point equation

with . One can find the coefficients by applying Cauchy’s differentiation formula and by integrating Eq. (Equation 15) numerically. This gives:

Let us now calculate ’s using the method described in Secs. II and III. The propagator for a complete graph with vertices is for any pair of . It is a -regular graph with , so as one can see from Eq. (Equation 10), each link introduces a suppression factor . On the other hand, the combinatorial factor contains a power , where is the number of vertices of the diagram. Thus in the thermodynamic limit, a diagram with links and vertices gives a contribution . The exponent is equal to one minus the number of closed loops in the diagram. Therefore, in the limit only tree diagrams give non-vanishing contributions. Our task simplifies therefore to summing only tree graphs. Each tree with links gives the following contribution to :

where is its symmetry factor and is the degree of vertex .

We performed the summation of all tree diagrams up to on a computer and checked that the values ’s obtained in this way agree with those obtained from Eq. (Equation 15). Again, we see that our method gives correct coefficients .

We should, however, remember that our goal is not only to find the coefficients of the expansion but rather which is a sum of infinitely many coefficients. If we naively terminate the series at some : , then e.g. is far away from the true value , because the series is slowly convergent. Therefore, we need to find a method which allows to read off the limiting value from first few coefficients. As we shall see below one can find a very good estimate of the limiting value using some general information about the properties of the series expansion. Let us first observe that the integral (Equation 4) and thus also is convergent only if the matrix

is positive definite. In the Appendix A we show that it is so for and for , and that acquires a zero mode for , so the integral (Equation 4) is divergent for . From this we can conclude that the asymptotic behavior of the coefficients has the following form:

with all , and . From the asymptotic behavior (Equation 18) we can deduce the following formula for :

where is a generalized Riemann Zeta function, and is estimated from the last two coefficients:

With the help of formula (Equation 19) we can predict now the value of for the SK model. To estimate the maximal error we used a method described in Appendix B. Our final result is in an excellent agreement with the analytic result cited above.

## 5Random -regular graphs and Cayley trees

A -regular graph is a graph with all degrees equal to . We say that a regular graph is random if its adjacency matrix is maximally random under the constraints that and for all . If we fix and let , random graphs become sparse and look locally like Cayley trees with degree , because the average density of finite-length loops goes to zero in this limit. Thus for large , instead of computing the coefficients by averaging them over many -regular random graphs one can calculate them for a single Cayley -tree.

The propagator is simply if nodes are connected, and zero otherwise. Unlike for a complete graph, the contribution from diagrams with loops cannot be neglected. Diagrams with multiple connections also do not vanish. In order to compute the contribution of a simple diagram and its multi-linked versions we need to evaluate the sum (Equation 11). The summation can be performed on a computer with only a slight complication as compared to the SK model.

The hardest task in the above is to compute . Recall that for a 1-regular graph, was either or zero. Now is simply equal to the number of different ways a given Feynman diagram can be lied down on a -tree, in such a way that any two neighboring vertices of the diagram are also neighbors on the tree. One expects that the number of such possibilities is of order . The factor comes about since one can lie down one vertex of the diagram anywhere on the tree. But the second has to be located on one of the neighbors of the first one. So each time by finding a position of a vertex we should check neighboring nodes on the tree, which gives the factor . To illustrate this, consider again the linear graph for from Fig. ?. The diagram has vertices. The first vertex can be put anywhere on the tree, but the next one only on one of neighbors of the first vertex. The same for the next one, for which we have again possibilities, etc., so altogether we have . Let us consider now a more complicated example of a square graph from Fig. ? and . Again we have a trivial factor for choosing the position of the first vertex on the -tree. Once it is chosen, we can position remaining vertices only in ways as shown in Fig. ?. So we have which is less than , because of the constraint coming from a closed loop.

Clearly, if a Feynman diagram has no loops then . For a diagram with loops one has to consider constraints on possible decorations. This can be done by enumerating all possible labellings and accepting only those which agree with graph’s structure. It can conveniently be done by a computer program. One point must be clarified here - since all ’s share the same trivial factor coming from possibilities of labeling the initial vertex, in the computer program one can just fix one vertex of the diagram to have some arbitrary label, and consider only decorations consistent with this choice. Next, one multiplies the result by , which then cancels in the definition of so that only numbers independent on remain.

Using this method, we calculated ’s and then the coefficients up to the given order (typically ) for and . The case has been analyzed in Sec. III. The results are summarized in Table ? where we give the values of for with estimated errors and compare them with those obtained by numerical simulations based on enumeration of all metastable states as described in Ref. [8]. To save space, the coefficients are not shown and can be found elsewhere [9].

The agreement with simulations is perfect. For the case there is also a beautiful analytic result [4] to compare with, which gives . As we see from the table, it agrees with our result within the error bars. We observe also that the rate of convergence of the series expansion for grows with . In other words, for smaller one should go to larger order which is however limited by the fast growth of the number of diagrams. This effect is slightly compensated by the fact that the complexity of computing the combinatorial factor is smaller for smaller .

## 6Regular -dimensional lattices

In this section we shall discuss how to calculate for -dimensional regular lattices, namely for (square lattice) and (cubic lattice). The lattices are -regular graphs with , but very special ones. The case , that is a closed chain, gives the same result as a random -regular graph, because the latter always contains at least one long chain, whose contribution to in the limit dominates over the contribution coming from shorter chains. In the general case, the only difference as compared to random regular graphs is that now, while calculating , we shall lie down Feynman diagrams on the -dimensional lattice. Again we see that is proportional to since the first vertex can be put anywhere, but once it is fixed we can put a neighboring one on a neighboring site of the lattice and repeat it iteratively for all remaining vertices of the diagram. The calculated values of are given in Table ?.

## 7Ladders

As a further example we consider another particular type of -regular graphs: graphs which we shall call ladders. A ladder is a graph obtained by stacking above each other infinitely many copies of a -dimensional cube so that corresponding vertices of the copies are aligned on a line and the corresponding vertices of consecutive copies are joined by a link. A ladder is just what one usually would call a ladder except that it is infinitely long. A ladder looks like a bookstand with square-shaped shelves. Such ladders are -regular graphs with . It is interesting to compare for ladders with those for random graphs and regular -dimensional lattices. We adopt the usual scheme. The only thing which changes is again , because now we have to lie down diagrams on the ladders. The final results are presented in Table ?.

Graph | ||||
---|---|---|---|---|

result | simulations | |||

2-reg. graph | 2 | 13 | 0.2414(2) | 0.242(2) |

3-reg. graph | 3 | 12 | 0.22484(4) | 0.226(1) |

2D ladder | 3 | 12 | 0.22568(6) | 0.226(2) |

4-reg. graph | 4 | 11 | 0.21762(2) | 0.219(1) |

2D lattice | 4 | 11 | 0.21808(2) | 0.219(2) |

3D ladder | 4 | 11 | 0.21799(2) | 0.220(4) |

6-reg. graph | 6 | 10 | 0.21101(2) | 0.211(1) |

3D lattice | 6 | 10 | 0.21125(1) | — |

SK model | 11 | 0.199226(5) | 0.199(1) | |

## 8Conclusion

We presented a semi-analytic method to estimate the rate of exponential growth of the number of metastable states of the Ising spin glass on different kinds of graphs. We checked that the method reproduces results known analytically, which are available for a few particular cases. The method is based on a diagrammatic representation of the quantity and on the exact enumeration of all Feynman’s diagrams up to a given order . The accuracy of the method improves with growing but already for close to it yields very precise estimates whose uncertainty varies in the range of order - depending on (see Table ?). The results show that the exponent is determined mainly by the degree . This suggests that the number of local minima depends strongly on local properties of the graph and weakly on its global topology. In other words, an important information about the complexity of the energy landscape of the corresponding spin glass is encoded in the short-range properties of the graph. It would be interesting to test if this also holds for other complex systems.

The method presented here can be applied to any type of graphs. However, the computational complexity of the method and the dependence of the accuracy of on the order has to be tested case by case.

In this paper we calculated for Gaussian ’s. Comparing the results to those for binomial ’s [5] we see that significantly depends on the distribution of ’s. It would be quite interesting to investigate the dependence on the distribution of ’s in a systematic way by calculating for some other continuous distributions of ’s. One can try to do this by applying the Tanaka-Edwards idea to distributions of the type . Another very challenging problem is to calculate higher moments and eventually also .

## Acknowledgments

We thank W. Janke, A. Krzywicki and O. Martin for discussion. We thank EC-RTN Network “ENRAGE” under grant No. MRTN-CT-2004-005616 and the Alexander von Humboldt Foundation for support. Z. B. acknowledges the support from a Marie Curie Actions Transfer of Knowledge project “COCOS”, Grant No. MTKD-CT-2004-517186 and a Polish Ministry of Science and Information Society Technologies Grant 1P03B-04029 (2005-2008).

## Appendix A

In this appendix we shall prove that the matrix is positive definite for . We must show that

is strictly positive for any vector with non-zero length. From Eq. (Equation 17) we have:

Introducing new variables : we can rewrite the right-hand side as:

The formula in the upper line shows that the quadratic form (Equation 20) is non-negative for any and for . It is zero only if and . Therefore, tends to a constant for . On the other hand, the second equation (Equation 22) tells us that for it is positive definite as well. However, for we see that a zero mode appears. Because of the zero mode the Gaussian integration in Eq. (Equation 2) is divergent and consequently for . These two observations indicate that the radius of convergence of given by Eq. (Equation 8) is one, and that has the form:

where and tends to some constant for .

## Appendix B

The estimation of systematic errors like those involved in the calculation of from Eq. (Equation 19) is not an easy task. We believe, that the method described below gives an upper bound on the error of . Let us slightly adjust the notation for the sake of clarity of the discussion. The value of from the left-hand side of Eq. (Equation 19) depends on , for which it has been estimated, so it is convenient to keep trace of in the notation. We will denote the value of the estimate by . Of course it not the same as which is on the right-hand side of Eq. (Equation 19) and which means just a sum of the first coefficients . The method relies on the following observation. We want to compute the deviation but we do not know the limiting value . We can however compute slightly modified differences for all , and say that for being sufficiently large. Now we can plot versus and extrapolate it to to obtain an estimator of which in turns estimates giving the error. Because falls with as a power of or faster, we can estimate from above by fitting a straight line to the points and shifting it so that all points lie below it. By extrapolating it to we get the value and use it to estimate the upper bound for the deviation between and the true value . This procedure is illustrated in Fig. ? for the case of SK model from Sec. III.

How reliable is this method? We checked that intervals obtained from first coefficients always include the value for for all graphs discussed in this paper. We checked also that reducing the error, say, by a factor two would result in many situations for which would lie outside the error bars. This means that the method does not overestimate the error too much.

### References

- F. Tanaka and S. Edwards, J. Phys. F
**10**, 2769 (1980). - A. J. Bray and M. A. Moore, J. Phys. C
**13**, L469 (1980). - T. Li, Phys. Rev. B
**24**, 6579 (1981). - B. Derrida and E. Gardner, J. Physique
**47**, 959 (1986). - D. S. Dean, Eur. Phys. J. B
**15**, 493 (2000). - B. D. McKay,
`http://cs.anu.edu.au/~bdm/nauty`

. - D. Sherrington and S. Kirkpatrick, Phys. Rev. Lett.
**35**, 1792 (1975). - Z. Burda, A. Krzywicki, O.C. Martin, Z. Tabor, Phys. Rev. E
**73**, 036110 (2006). `http://www.physik.uni-leipzig.de/~waclaw/`

`glasses-data.htm`