# The realizable solutions of the TAP equations

###### Abstract

We show that the only solutions of the TAP equations for the Sherrington-Kirkpatrick model of Ising spin glasses which can be found by iteration are those whose free energy lies on the border between replica symmetric and broken replica symmetric states, when the number of spins is large. Convergence to this same borderline also happens in quenches from a high temperature initial state to a locally stable state where each spin is parallel to its local field; both are examples of self-organized criticality. At this borderline the band of eigenvalues of the Hessian associated with a solution extends to zero, so the states reached have marginal stability. We have also investigated the factors which determine the free energy difference between a stationary solution corresponding to a saddle point and its associated minimum, which is the barrier which has to be surmounted to escape from the vicinity of a TAP minimum or pure state.

## I Introduction

One of the most influential papers in the theory of spin glasses was the paper of Thouless, Anderson and Palmer (TAP) Thouless et al. (1977). They provided a set of coupled equations for the magnetization at site of the Sherrington-Kirkpatrick (SK) Sherrington and Kirkpatrick (1975) model of Ising spin glasses. Since then equations equivalent to those of TAP have been studied for -spin models, which are models for structural glasses, and also for a host of computer science applications Mézard (2017).

The free energy (multiplied by ) associated with a TAP state for the Ising spin SK model is

(1) | |||||

where . The TAP equations themselves are derived from the stationarity equations and take the form

(2) |

The Hessian associated with the stationary points of was studied long ago Bray and Moore (1979); Aspelmeier et al. (2006). It is defined by

where denotes the magnetization at site at a stationary point.

A great deal is already known about the solutions of the TAP equations and their associated Hessians. There are an exponentially large number of solutions for , that is . The complexity of the minima of is defined by

(4) |

where denotes the number of minima of free energy per spin . quantifies the number of solutions when it is exponentially large. It is non-zero over a range of values Bray and Moore (1980, 1981a); Aspelmeier et al. (2004); Cavagna et al. (2004). The solution of lowest free energy per spin is one of the pure states of the Parisi replica symmetry broken solution (RSB) Parisi (1979); Mézard et al. (1987). The pure states are those whose free energies per spin are only above . There is a critical at which the solutions change their nature Bray and Moore (1980, 1981a) . Those solutions at are uncorrelated with each other and their Hessians have a single, nearly null eigenvalue, (whose value vanishes in the limit ), separated by a gap from a band of eigenvalues. As the free energy is reduced towards the gap goes to zero and for all there is no gap between the lowest “null” eigenvalue and the bottom of the band Aspelmeier et al. (2004). The states with have non-trivial RSB overlaps with each other Bray and Moore (1980, 1981a). Those with have trivial (zero) replica symmetric (RS) overlap with each other.

A stationary solution of the TAP equations which corresponds to a minimum will have all the eigenvalues of its associated Hessian non-negative. It turns out that the other stationary points are saddle points with just a single negative eigenvalue Aspelmeier et al. (2004). Every minimum has its associated saddle point, so the complexity of the saddle points is identical to that of the minima. Once over the saddle-point in the direction away from the minimum one goes towards the trivial paramagnetic solution () for which . Thus the free energy landscape of the free energy functional of Eq. (1) is simple; it has a spoke-like arrangement of minimum, associated saddle and at the center of the wheel Aspelmeier et al. (2004); Cavagna et al. (2004). The free energy difference between the free energy at the saddle point and its associated minimum is the barrier which has to be surmounted to escape from the minimum. In this paper we have once more investigated these barriers in order to determine some of the factors which might control their magnitude. The simulations of Billoire et al. Billoire and Marinari (2001); Billoire (2010) indicate that the barriers separating pure states may scale as . Our studies suggest that non-pure state solutions will have much smaller barriers, of , so that it will be possible to escape from their vicinity by thermal fluctuations (in agreement with our earlier study Aspelmeier et al. (2006)).

The chief purpose of this paper is to report a feature of the actual solutions found in numerical work which has not previously been noticed. In numerical work, the free energy minima can be obtained by the iterative map:

(5) |

where is a parameter which controls the approach to the next iterate Aspelmeier et al. (2006). In this work, we set . The value of affects the free energy per spin which is obtained, and instead of obtaining a spread of values of the iteration scheme picks out in the large limit one particular value of , . What we want to point out is that our states at this particular value of have the properties of states at the critical value , which is the borderline between states which are replica symmetric (RS) and those whose replica symmetry is broken (RSB) Bray and Moore (1980, 1981a). is less than . is that associated with all possible minima of the TAP functional, rather than the subset produced by the chosen iteration scheme. In our work the states found lie close to , differing from it by an amount which decreases as becomes large. Even though corresponds to a free energy per spin below , (so it nominally lies in the region where the states would have RSB features), the states produced in the iteration (which are just a subset of all the possible states at ) do not have this feature. Instead the subset of states generated is closer in its properties to those at the borderline itself.

The iteration procedure of Eq. (5) is just one of a large number of ways of solving the TAP equations, but we suspect that any iterative solution will have the same features as those found using Eq. (5). We briefly studied the iterative procedure of Bolthausen Bolthausen (2014) but found it less efficient than that of Eq. (5) in converging to a stationary solution: Starting from some initial state a common feature is just bouncing around without convergence. However, Eq. (5) was more likely to find a solution than the Bolthausen method at large values and we have used it throughout this paper.

What led us to carry out this investigation was the work of Sharma et al. Sharma et al. (2018, 2016). It was found in these papers that a quench from an initial random (i.e. high-temperature) spin configuration by successively putting spins in turn parallel to their local fields until all are so aligned led to a final quenched state which in SK type models lay on the boundary between replica symmetric states and states with RSB. This is what we also find for solutions of the TAP equations; the iterative solution has parallels with the quenching procedure. The number of quenched states at has also been studied as a function of their energy, and there exists a critical energy per spin below which the states have RSB features and above which the states are uncorrelated Bray and Moore (1980, 1981b). A problem with studying the Ising model at , i.e. in the quenched state, is that a Hessian cannot be constructed as the spins take the values , so marginality as indicated by eigenvalues of a Hessian extending down to zero Müller and Wyart (2015), cannot be investigated. A big disadvantage of studying the finite temperature TAP equations is that the values of which can be studied with the TAP equations are much smaller than those which can be handled in a quench. The same self-organized critical features are present in both the solutions of the TAP equations and in the quenched states and presumably the physics behind this is the same Sharma et al. (2018), i.e. somewhat obscure, at least to us. However, for certain aspects of the quenched problem one has some features which are rigorously established; Newman and Stein Newman and Stein (1999) have shown that in the large limit, the quench takes one to a particular energy per spin which is self-averaging, but dependent on the algorithm used. It would be nice if their proof could be extended to the somewhat similar TAP problem, as our work shows that as gets large that there is convergence to a particular free energy .

The details of our numerical work can be found in Sec. II while in Sec. III we present the evidence that in the large limit the TAP solutions which can be found lie at the boundary between solutions whose overlaps are replica symmetric and those whose overlaps are those of broken replica symmetry. Our work on barriers is in Sec. IV. We have fitted the free energy between the minimum and the saddle with a quartic fit and as a consequence can relate the barrier height to the difference in the values of at the minimum and the saddle, and the curvatures at the minimum and the saddle. We then discuss how the barriers between pure states could become of order . Finally in Sec. V we comment upon unresolved issues. In Appendix A we have derived a finite size correction to the position of the Hessian band-edge, which turns out to work well at the rather modest values of which we can study.

## Ii Simulation details

We studied the TAP equations for and spins and 500 bond realizations for each . For each realization we tried to find solutions by iteration according to Eq. (5). As mentionned in Sec. I we used throughout. We chose the temperature as a compromise between having too small a probability of finding any solution at all, as happens for close to , and having , which is the case for close to 0. The latter would lead to complications by causing a very large spread in the eigenvalues of the Hessian, as discussed in Aspelmeier et al. (2006), obfuscating the issues we are focusing on here.

In order to avoid questions of numerical accuracy, which can be very delicate in the complex TAP free energy landscape, we used arbitrary precision arithmetic with 512 binary digits for the final approach to a TAP solution and for subsequent calculations. The final approach is done in terms of the transformed variables by iterating a transformed version of Eq. (5) (see Eq. (12) in Aspelmeier et al. (2006)) and with since the final approach starts off already in the basin of attraction. Use of the transformed variables is necessary because the original may take the values within numerical accuracy upon iterating Eq. (5) directly, in which case the Hessian is ill defined, see Eq. (LABEL:Hessian).

For each solution found we tried to locate the corresponding saddle using the method described in Aspelmeier et al. (2004). We then calculated the eigenvalues and eigenvectors of the Hessian at both the minimum and the saddle (if it was found). Since the Hessian matrices can be very ill-conditioned, it is in this step where the arbitrary precision arithmetic is most useful.

Since the quantities we examine in this work may have strongly non-Gaussian distributions, such as for instance the low-lying eigenvalues of the Hessian, with possibly fat tails, we used the median instead of the mean consistently throughout this work for robustness. Accordingly, all error bars shown are 95% confidence intervals for the median.

## Iii Self-organized criticality

In this section we give the details of why we believe that the solutions of the TAP equations which are found by an iterative process lie at the boundary between replica symmetric and broken replica symmetry solutions. In Fig. 1 we have plotted the free energy per spin of the solutions found at a temperature and with as a function of . The variance decreases strikingly rapidly as increases, suggesting that as , there will be a well-defined limit for the free energy . In Ref. Cavagna et al. (2004) the entire curve was obtained when studying values of up to . Unfortunately as increases this becomes harder and harder to do as the solutions found are approaching and solutions well away from this value are rarely found. Thus the authors of Ref. Cavagna et al. (2004) only succeeded in finding all the solutions by virtue of finite size effects. As grows, the chance of finding solutions well away from rapidly decreases to zero. The rapid decrease of the variance with increasing is very suggestive that the solutions being found do not come from all over the curve, (which would lead to an independent variance) but rather are just those associated with a particular value of the free energy. We do not have any arguments as to why the dependence of and the square root of its variance should vary as ; we only use a as this form arises frequently for finite size effects in the SK model Aspelmeier et al. (2006, 2008). The data for is not even monatonic as a function of which suggests that the values of which we can study are not yet large enough to be in the asymptotic regime for this quantity.

Our contention is not only that the solutions found in an iterative procedure converge to a unique value of the free energy as but the particular free energy converged to is the critical free energy which separates states with vanishing overlaps from those with non-trivial overlaps. We shall refer to this borderline as the RS/RSB critical point. The free energy per spin is the free energy at this borderline when all possible minima of the TAP equations are studied. The subset of these states which we obtain by iteration whose free energies are close to have the features of states at . We have obtained the Hessians associated with the minima obtained by iteration. Fig. 2 shows that the two lowest eigenvalues of the Hessian seem to be both approaching zero as . The smallest eigenvalue is the “null” eigenvalue associated with the broken supersymmetry Aspelmeier et al. (2004); Cavagna et al. (2004). The second eigenvalue lies at the bottom of the band of eigenvalues of the Hessian and for states with should be different from by a finite amount which does not vanish as . Notice that for the smallest value of in the plot, , where one will be sampling states over a wide range of values, one can see that indeed looks quite distinct from .

The bottom of the band can be calculated by considering the matrix defined via , i.e. the Hessian without the projector term. The projector term, being , is only a small perturbation which changes the eigenvalues of only slightly except for the isolated one. We expect that , the second smallest eigenvalue of . If we define , then the band-edge without finite size corrections should be at Bray and Moore (1979); Plefka (2002)

(6) |

This is, however, only an approximation valid for small . In Appendix A we show how to find the exact individual band edge numerically. We have plotted for in Fig. 2 but the agreement with the measued values of is not good, presumably because of finite size effects. In Appendix A we describe how to obtain a finite size correction for the band-edge, which does indeed improve the agreement with .

The results in Fig. 2 show that as increases both and are approaching zero, indicating that the solutions we are finding in this limit are similar to those whose free energy is less than . Below the states are associated with full replica symmetry breaking Bray and Moore (1981a) and would be associated with massless modes so that for all TAP states with one would expect both and to decrease as , just as found in Fig. 2. To see that the convergence is not to a state below but to a state right at the bordeline between RS and RSB states, we have studied the overlaps of the solutions in Fig. 3.

The probability density of overlaps of solutions at a given size was defined as

(7) |

(Note that here is not that of the TAP equation; and denote two distinct solutions). We have also averaged the result over realizations. Fig. 3(a) shows for two different system sizes. Both peak at , a feature which would not be expected when . Right at the expected form of in the large limit. We expect that this peak is broadened by finite size effects so that the data at finite and could be collapsed onto a universal curve by plotting against , but we do not have data at large enough values of to study this. Fig. 3(c) shows the variance of shrinks with for , which is what would be expected if is approaching a delta function at large . For states with the variance of would be expected to approach a non-zero value in the large limit. Fig. 3(b) shows the number of bond realizations for which a given number of solutions was found. For , the most common number of solutions found was zero! For the situation improves, presumably because the larger the value of the more solutions there are to be found. Fig. 3(c) illustrates why the values which we can reach are a long way away from being in the large regime for some quantities, and Fig. 3(b) illustrates how hard it is to get non-trivial solutions of the TAP equations.

The fact that also explains another puzzling feature associated with solving the TAP equations by iteration; the solutions have always been reported from the earliest days as having a Hessian spectrum whose band-edge extended down to zero Bray and Moore (1979), rather than having a finite band gap as expected, for example, at the peak of the or for any . In fact, no finite band gap has ever been clearly seen in numerical studies of the TAP equations.

We conclude that the solutions which are found by iteration are at an RS/RSB border. They are an example of self-organized criticality. The states (solutions) are associated with a Hessian whose eigenvalues extend to zero, and so are also marginal Müller and Wyart (2015) as well as self-organized.

## Iv Barriers

In the SK model, the low-temperature spin glass phase has broken replica symmetry. That means it is associated with pure states, whose free energy per spin differ from each other by terms of Mézard et al. (1987). Escape from a pure state is prevented by large barriers. The simulations in Refs. Billoire and Marinari (2001); Billoire (2010) indicate that the barriers scale with the number of spins as . Unfortunately there seems to be only a little understanding of the origin or form of these barriers Rodgers and Moore (1989). In this paper we shall try to cast some light on them by assuming that TAP solutions whose free energies per spin are within of , the solution of lowest free energy, can be identified as pure states and that the barrier for escaping a pure state can be identified with the free energy difference between the free energy of a TAP minimum and its associated saddle point. Alas, as pointed out in Sec. III the only states which we can find by directly solving the TAP equations are those which are around a free energy on the RS/RSB boundary (i.e. around ) and not those whose free energies lie within of . However, by examining the factors which determine the magnitude of barriers we have been able to understand the features which have to be present for barriers to scale as .

The TAP free energy as a function of and is defined Aspelmeier et al. (2004) via

(8) |

where is the functional of and of Eq. (1) except that now is regarded as an independent variable, unrelated to the ; is a function of the variables , whereas the original TAP free energy depends only on the variables (with defined as ). One can easily show that the stationarity equations for reproduce the TAP equations: . However for these new equations the quantity is in general not equal to the parameter appearing in the equations. However, the additional stationarity equation, forces at stationary points in the full -dimensional space. Therefore at the minimum and the saddle the free-energy functions and have the same and values. By formally eliminating the variables one can obtain the function as a function of and in Fig. 4 we have plots of and as functions of .

To understand how the barrier height, which is the free energy difference between the saddle-point value of the free energy and the minimum value , i.e. , might depend on the values of and the curvatures at the minimum and the saddle, we have used a quartic fit to :

(9) | |||||

where we will relate the coefficients , and to the curvatures at the saddle and the values of at the minimum and at the saddle, . is stationary when

(10) |

is the solution of . The free energy at the saddle is zero and at the minimum

(11) |

The barrier is . The curvature at the saddle is defined as

(12) |

and equals at the saddle , while at the minimum

(13) |

These curvatures at the saddle and at the minimum , where the curvatures and were discussed in Ref. Aspelmeier et al. (2004);

(14) |

evaluated for at values of at the saddle or the minimum.

We can eliminate the coefficients , , and , to get

(15) |

In Fig. 5 we have plotted the observed barrier divided by the right hand side of Eq. (15) to check the accuracy of this equation. It clearly works well for most saddle-minima pairs, but a few are clearly not well-accounted for by the quartic fit of Eq. (9). For these pairs the discrepancy is simply because the neglected higher terms are just not always negligible.

Assuming that the quartic fit provides a good fit to the barrier height we next describe the dependence of the terms in Eq. (15). In Fig. 6 evidence is presented that decreases with as while the curvatures grow like . The variation of means that the saddle becomes very close to the minimum in the large limit. Eq. (15) then implies that the barrier height should be independent. In Ref. Aspelmeier et al. (2006) we showed by varying the iteration parameter that the barriers were -independent, varying as . Hence at the critical free energy between RS/RSB states, the barriers would be expected to be -independent. For pure states , which explains why pure states have barriers of order . Note that at , the free energy associated with our iterative solutions, the barriers are numerically tiny, as can be seen from Fig. 4.

## V Discussion

While the SK model is referred to as a “solvable” model, the finite size corrections to the thermodynamic limit have only been obtained for a few quantities from analytical work. Mostly all that we have are rather unsatisfactory estimates from numerical studies. The same is true of the TAP equations. They become exact in the thermodynamic limit, but finite corrections to them and the dependencies in their solutions have not really emerged from analytical studies. We have suggested that TAP solutions of very low free energies correspond to the pure states of the SK model, in that if one could compute a value of in the pure state it would correspond to that of a TAP solution, but we know of no proof of this possibility.

Our main discovery is that at large values of the solutions of the TAP equations fall at the boundary between states with replica symmetric overlaps and those with overlaps like those of broken replica symmetry. This is like a critical point. These states are associated with massless modes at large and so the solutions found are those of a self-organized marginally stable critical system.

We do not know how this behavior comes about. But as the same behavior arises in quenched states of the SK model, it seems there exists a phenomenon worthy of further study.

## Appendix A Individual band edge and its finite size correction

In this Appendix, we derive the form of the finite size corrections to the individual band-edge which was used in constructing Fig. 2. By individual we mean for a given TAP solution.

The eigenvalue density of , i.e. the Hessian without the projector term, can be obtained from its resolvent

as

(16) |

As explained in Plefka (2002), the resolvent of satisfies the equation

(17) |

in the large limit according to Pastur’s theorem Pastur (1973). The two resolvents are related by , hence satisfies

(18) |

after a change of variables . Using a quadratic approximation to Eq. (17) valid for small , one obtains Plefka (2002)

for small and , where and are defined as in Sec. III. The band edge is thus at ^{1}^{1}1Note that it says in Plefka (2002) for the band edge, which is a typo. Note also that our band edge differs from Plefka’s by a factor of which is hidden in our definition of ..

However, is not always small in our numerical experiments. Hence we refined this approximation by searching numerically for the infimum of real for which Eq. (18) has no appropriate real solution, as this marks the onset of the band of eigenvalues according to Eq. (16). To this end, define and ; Eq. (18) then reads

(19) |

The largest , denoted by , which still allows for a real solution is the maximum of the right hand side for from the interval . This interval follows from the discussion in the appendix of Plefka (2002) about selecting the appropriate solution of Eq. (17). The value of at which is attained is denoted . By numerical optimization both and can easily be found. The individual band edge improves on by going beyond the quadratic approximation but it is still an infinite system result through the use of Pastur’s theorem.

Hence we are now looking for a finite size correction to it. The eigenvalue density , when computed from the full equation (18), still starts off with a square root singularity, i.e.

for close to and some constant . For a system of size the smallest eigenvalue will be roughly determined by the condition

(20) |

such that in our case

so . The second smallest eigenvalue can be calculated in the same way by replacing the right-hand side of Eq. (20) by 2, so .

The constant can be calculated as follows. Tayor expansion to second order of the right hand side of Eq. (19) around the maximum gives

such that

for . On the other hand by definition. Comparison with Eq. (16) shows

This can be evaluated numerically since is already known.

Thus we have calculated the individual band edge and its finite size correction for ,

and for ,

## References

- Thouless et al. (1977) D. J. Thouless, P. W. Anderson, and R. G. Palmer, “Solution of ’Solvable model of a spin glass’,” The Philosophical Magazine: A Journal of Theoretical Experimental and Applied Physics 35, 593–601 (1977).
- Sherrington and Kirkpatrick (1975) David Sherrington and Scott Kirkpatrick, “Solvable Model of a Spin-Glass,” Phys. Rev. Lett. 35, 1792 (1975).
- Mézard (2017) Marc Mézard, “Mean-field message-passing equations in the Hopfield model and its generalizations,” Phys. Rev. E 95, 022117 (2017).
- Bray and Moore (1979) A J Bray and M A Moore, “Evidence for massless modes in the solvable model of a spin glass,” Journal of Physics C: Solid State Physics 12, L441 (1979).
- Aspelmeier et al. (2006) T. Aspelmeier, R. A. Blythe, A. J. Bray, and M. A. Moore, ‘‘Free-energy landscapes, dynamics, and the edge of chaos in mean-field models of spin glasses,” Phys. Rev. B 74, 184411 (2006).
- Bray and Moore (1980) A J Bray and M A Moore, “Metastable states in spin glasses,” Journal of Physics C: Solid State Physics 13, L469–L476 (1980).
- Bray and Moore (1981a) A J Bray and M A Moore, “Metastable states in the solvable spin glass model,” Journal of Physics A: Mathematical and General 14, L377–L383 (1981a).
- Aspelmeier et al. (2004) T. Aspelmeier, A. J. Bray, and M. A. Moore, “Complexity of Ising Spin Glasses,” Phys. Rev. Lett. 92, 087203 (2004).
- Cavagna et al. (2004) Andrea Cavagna, Irene Giardina, and Giorgio Parisi, “Numerical Study of Metastable States in Ising Spin Glasses,” Phys. Rev. Lett. 92, 120603 (2004).
- Parisi (1979) G. Parisi, ‘‘Infinite number of order parameters for spin-glasses,” Phys. Rev. Lett. 43, 1754 (1979).
- Mézard et al. (1987) M. Mézard, G. Parisi, and M. A. Virasoro, Spin Glass Theory and Beyond (World Scientific, Singapore, 1987).
- Billoire and Marinari (2001) Alain Billoire and Enzo Marinari, “Correlation timescales in the Sherrington-Kirkpatrick model,” Journal of Physics A: Mathematical and General 34, L727–L734 (2001).
- Billoire (2010) Alain Billoire, “Distribution of timescales in the Sherrington–Kirkpatrick model,” Journal of Statistical Mechanics: Theory and Experiment 2010, P11034 (2010).
- Bolthausen (2014) Erwin Bolthausen, “An Iterative Construction of Solutions of the TAP Equations for the Sherrington–Kirkpatrick Model,” Communications in Mathematical Physics 325, 333–366 (2014).
- Sharma et al. (2018) Auditya Sharma, Joonhyun Yeo, and M A Moore, “Self-organized critical behavior and marginality in Ising spin glasses,” Journal of Statistical Mechanics: Theory and Experiment 2018, 053302 (2018).
- Sharma et al. (2016) Auditya Sharma, Joonhyun Yeo, and M. A. Moore, “Metastable minima of the Heisenberg spin glass in a random magnetic field,” Phys. Rev. E 94, 052143 (2016).
- Bray and Moore (1981b) A J Bray and M A Moore, “Metastable states in spin glasses with short-ranged interactions,” Journal of Physics C: Solid State Physics 14, 1313–1327 (1981b).
- Müller and Wyart (2015) Markus Müller and Matthieu Wyart, “Marginal Stability in Structural, Spin, and Electron Glasses,” Annual Review of Condensed Matter Physics 6, 177–200 (2015).
- Newman and Stein (1999) C. M. Newman and D. L. Stein, “Metastable states in spin glasses and disordered ferromagnets,” Phys. Rev. E 60, 5244–5260 (1999).
- Aspelmeier et al. (2008) T Aspelmeier, A Billoire, E Marinari, and M A Moore, “Finite-size corrections in the Sherrington Kirkpatrick model,” Journal of Physics A: Mathematical and Theoretical 41, 324008 (2008).
- Plefka (2002) T Plefka, “Modified TAP equations for the SK spin glass,” Europhysics Letters (EPL) 58, 892–898 (2002).
- Rodgers and Moore (1989) G J Rodgers and M A Moore, “Distribution of barrier heights in infinite-range spin glass models,” Journal of Physics A: Mathematical and General 22, 1085–1100 (1989).
- Pastur (1973) L A Pastur, “Spectra of random self adjoint operators,” Russian Mathematical Surveys 28, 1–67 (1973).
- (24) Note that it says in Plefka (2002) for the band edge, which is a typo. Note also that our band edge differs from Plefka’s by a factor of which is hidden in our definition of .