text:
# Metastable configurations on the Bethe lattice

## I Introduction

Despite years of efforts, the nature of the glassy phase of finite dimensional spin glasses is still not clear. It is still debated whether the replica symmetry breaking (RSB) scheme proposed by Parisi to solve the mean field fully connected Sherrington-Kirkpatrick model (SK), and implying the existence of an exponentially large number of pure states with an ultrametric structure, holds in some way for these systems. Part of the difficulty to settle the question is the very poor analytical tractability of finite dimensional systems.

Looking for more realistic – but still tractable – models than SK, much attention was recently paid to spin glasses on random graphs of finite connectivity. These models incorporate the short range nature of interactions, but without the underlying geometry of the finite dimensional models. It has been shown that the replica symmetric (RS) solution of spin glasses on random graphs of finite connectivity is unstable at low temperature mottishaw ; sompo . Unfortunately working out the RSB scheme is far more difficult than for the SK, because the glassy phase can no longer be characterized by the only two spin overlap $`\sigma _a\sigma _b`$ between two distinct replicas $`a`$ and $`b`$, but requires all of the multi spin overlaps. Recently a new approach was suggested in mezpar , based on the cavity method and population algorithms, that allows for a numerical solution at the level of one step RSB (1RSB). The method can be virtually extended to any step of RSB, at the price of increasing computer resources.

In this paper, we use similar ideas to shed a new light onto an old problem: the computation of the number of metastable configurations. A configuration is said to be metastable if its energy cannot be decreased by flipping a single spin. One expects, within the RSB scenario, that a consequence at zero temperature of the exponential number of pure states, is an exponential number of metastable configurations. Beware however that these two concepts are not as obviously linked as it might seem: in systems with finite connectivity the zero temperature limit of the pure states are not the metastable configurations, but are argued to be the configurations stable with respect to an arbitrary finite number of spin flips birolimonasson ; PARISI\_MEZARD\_T0 .

Let us call configurational entropy $`S_C(E)`$ the logarithm, divided by the number of spins, of the number of metastable configurations having an energy density equal to $`E`$. Lots of efforts have been devoted to the study of $`S_C`$ for SK. More recently, attention turned to spin systems on random graphs. On random graphs with fixed finite connectivity, annealed computations were carried out for the binary spin glass dean , and the ferromagnet was addressed in lefevredean . On random graphs of fluctuating finite connectivity, the case of the ferromagnet was solved in berg  – also via a population algorithm, but in a different context as ours –, and the authors give hints how to address the quenched computation of spin glasses.

Our method, which is quite general, enables us to recover all the above results. As new contributions, we carry out quenched computations in the case of the Gaussian or binary spin glass on random graphs of fixed finite connectivity, and we exemplify the implementation of 1RSB.

In the following we stick to the case of random graphs of fixed finite connectivity. For the sake of concision, we display no result about the fluctuating connectivity, but it is easy to generalize our method to this case.

The layout of the paper is the following. In section II we give several definitions and we set the notations. In section III we derive our basic equations in the RS framework. In section IV we analytically solve these equations in the case of a ferromagnet. We compare our results with the microcanonical approach of lefevredean . In section V we turn to the binary spin glass, where the coupling constants are $`\pm 1`$. The population algorithm shows up there. In section VI we address the slightly more involved Gaussian spin glass. Eventually, in section VII we show how to perform 1RSB, and illustrate the algorithm on the case of the spin glass with binary couplings.

## II The system under study

Following mezpar , we call Bethe lattice a random graph with fixed connectivity equal to $`k+1`$, i.e. the number of edges incident to each vertex is exactly $`k+1`$. Locally such a graph has the structure of a Cayley tree, and can be obtained with the following procedure: starting from a root (or ancestor) $`\mathrm{\Psi }`$, one builds a first generation of $`k+1`$ sons, and then successive generations of $`k`$ sons to each vertex, as displayed in Fig.(1). By contrast to the Cayley tree, the Bethe lattice has loops, but small ones are rare: the typical length of a loop is of order $`\mathrm{log}N`$. In practical case one uses the tree-like structure to write down recursive (or ’cavity’ in the language of the physics of disordered systems mpv ) equations, and the existence of loops is enforced by saying that all the spins are equivalent, which is not true for the Cayley tree (in particular the Cayley tree has strong boundary effects).

On each vertex $`A`$ of the Bethe lattice stands a spin, and these spins interact through the Hamiltonian

$$=\underset{\{A,B\}}{}J_{A,B}\sigma _A\sigma _B,$$
(1)
where the sum runs over all the edges. The coupling constants $`J_{A,B}`$ may be fixed, like in the case of the ferromagnet, or they may be quenched random variables in the case of a spin glass.

The local field acting on the spin $`A`$ is $`H_A=_{B𝒱(A)}J_{A,B}\sigma _B`$, where the sum runs over the $`k+1`$ first neighbors of $`A`$. The metastable states are characterized by $`A,H_A\sigma _A0`$. Our aim is to compute the partition function restricted to the metastable configurations:

$$Z=\text{Tr}\left[e^\beta \underset{A}{}\mathrm{\Theta }(H_A\sigma _A)\right],$$
(2)
where $`\mathrm{\Theta }`$ is the Heaviside step-function such that $`\mathrm{\Theta }(u)=1`$ if $`u0`$, 0 otherwise, and Tr stands for the summation over all the possible values of all the spins.

## III The cavity equations

Let us introduce, in close analogy with the construction of the infinite tree in ALDOUS , a simple labeling rule for the vertices of a tree: the sons of the root are labeled $`1,2,\mathrm{}k`$. Then recursively the sons of a vertex labeled $`(i_1\mathrm{}i_p)`$ are labeled $`(i_1\mathrm{}i_p\mathrm{\hspace{0.17em}1}),\mathrm{}(i_1\mathrm{}i_pk)`$. If $`A=(i_1\mathrm{}i_p)`$, we can define the norm of the symbol as $`|A|p`$.

The basic building block of the cavity solution is a branch of the Cayley tree, i.e. a subtree made of a given vertex different from the root and all its descendants. On Fig. (2) is drawn a branch rooted at $`\mathrm{\Phi }`$. Let us call it $`T_\mathrm{\Phi }`$.

The key property of a branch is that it is structurally similar to each of its sub-branches made of a given vertex and all its descendants. This makes a recursive approach possible to study the thermodynamics of a branch. We will see afterwards how to use the results for the Bethe lattice.

### III.1 Computing the metastable states on a branch

As it will be clear later, on a branch it is convenient to require that all the spins are stables, except the root. The Hamiltonian of a branch is

$$_\mathrm{\Phi }=\underset{\{A,B\}T_\mathrm{\Phi }}{}J_{A,B}\sigma _A\sigma _B.$$
(3)
The definition of the partition function is

$$Z_\mathrm{\Phi }=\text{Tr}\left[e^{\beta _\mathrm{\Phi }}\underset{AT_\mathrm{\Phi }\backslash \mathrm{\Phi }}{}\mathrm{\Theta }(H_A\sigma _A)\right].$$
(4)
By contrast to all the other spins in $`T_\mathrm{\Phi }`$, which have $`k+1`$ neighbors, spin $`\mathrm{\Phi }`$ lacks one neighbor. In all the following the local field acting on such a spin with only $`k`$ neighbors will be called a ’cavity’ field and denoted by a small $`h`$ (here $`h_\mathrm{\Phi }`$). We call $`p_\mathrm{\Phi }(h,\sigma )`$ the joint density probability of $`h_\mathrm{\Phi }`$ and the value of spin $`\mathrm{\Phi }`$. One has

$$p_\mathrm{\Phi }(h,\sigma )=\frac{1}{Z_\mathrm{\Phi }}\text{Tr}\left[e^{\beta _\mathrm{\Phi }}\delta (\sigma _\mathrm{\Phi }\sigma )\underset{|A|1}{}\mathrm{\Theta }(H_A\sigma _A)\delta (h_\mathrm{\Phi }h)\right].$$
(5)
We call $`T_1,\mathrm{}T_k`$ the subtrees of $`T_\mathrm{\Phi }`$ engendered by points $`1,\mathrm{},k`$. We can define their Hamiltonians similarly to Eq. (3), and write

$$_\mathrm{\Phi }=\underset{i=1}{\overset{k}{}}_ih_\mathrm{\Phi }\sigma _\mathrm{\Phi }.$$
(6)
Splitting the Tr, equation (5) can be restated as

$$p_\mathrm{\Phi }(h,\sigma )=\frac{1}{Z_\mathrm{\Phi }}e^{\beta h\sigma }\underset{\sigma _1,\mathrm{}\sigma _k}{\text{Tr}}\delta (h_\mathrm{\Phi }h)\underset{|A|2}{\text{Tr}}e^{\beta _{i=1}^k_i}\underset{i=1}{\overset{k}{}}\mathrm{\Theta }(H_i\sigma _i)\underset{|A|2}{}\mathrm{\Theta }(H_A\sigma _A).$$
(7)
Now separate $`H_i`$ into two contributions:

$$H_i=\underset{j=1}{\overset{k}{}}J_{i,ij}\sigma _{ij}+J_{\mathrm{\Phi },i}\sigma _\mathrm{\Phi },$$
(8)
and enforce the fact that the first one is the cavity field acting on spin $`i`$ in the absence of $`\mathrm{\Phi }`$ through the identities

$$1=𝑑h_i\delta \left(\underset{j}{}J_{i,ij}\sigma _{ij}h_i\right).$$
(9)
Plugging (8) and (9) into (7) one gets

$`p_\mathrm{\Phi }(h,\sigma )={\displaystyle \frac{1}{Z_\mathrm{\Phi }}}e^{\beta h\sigma }\underset{\sigma _1,\mathrm{}\sigma _k}{\text{Tr}}\delta (h_\mathrm{\Phi }h){\displaystyle \underset{i=1}{\overset{k}{}}}{\displaystyle 𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Phi },i}\sigma )\sigma _i\right]}`$ (10)
$`\times \underset{|A|2}{\text{Tr}}e^{\beta _{i=1}^k_i}{\displaystyle \underset{|A|2}{}}\mathrm{\Theta }(H_A\sigma _A){\displaystyle \underset{i=1}{\overset{k}{}}}\delta ({\displaystyle \underset{j}{}}J_{i,ij}\sigma _{ij}h_i)`$
Comparing to Eq. (5), it turns out that the second line in the above equation is nothing but the product of $`Z_ip_i(h_i,\sigma _i)`$ for $`i=1,\mathrm{}k`$. So one eventually obtains the following recursion relation:

$$p_\mathrm{\Phi }(h,\sigma )=\frac{Z_1\mathrm{}Z_k}{Z_\mathrm{\Phi }}e^{\beta h\sigma }\underset{\sigma _1,\mathrm{}\sigma _k}{\text{Tr}}\delta (h_\mathrm{\Phi }h)\underset{i=1}{\overset{k}{}}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Phi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i).$$
(11)
By averaging over the coupling constants and the random graphs (operation denoted by $`\overline{\mathrm{}}`$), one defines the probability distributions $`𝒫_\mathrm{\Phi }^{(h,\sigma )}(p)=\overline{\delta (p_\mathrm{\Phi }(h,\sigma )p)}`$. Equation (11) induces a functional relation between $`𝒫_\mathrm{\Phi }^{(h,\sigma )}`$ and the analogous probability distributions $`𝒫_i^{(h,\sigma )}`$ for $`i=1,\mathrm{}k`$. In the Bethe lattice all spins are required to be equivalent so we impose the condition

$$i,𝒫_\mathrm{\Phi }=𝒫_i=𝒫.$$
(12)
This yields a self-consistency equation for $`𝒫`$.

### III.2 Performing measures on the Bethe lattice

Having solved the thermodynamics of a branch, the idea of mezpar  is to describe the Bethe lattice as the result of the merging of several branches. We define two merging procedures:

* Take $`k+1`$ branches, with roots labeled $`i=1,\mathrm{}k`$ and merge them onto a new spin $`\mathrm{\Psi }`$. The resulting tree is the one of Fig. 1.
* Take 2 branches, with roots labeled 1 and 2, and merge them via a new bond of coupling constant $`J`$, without adding any site. See Fig. 3.

Operation (M1) is very similar to the iteration Eq. (11): the differences are that $`k`$ is to be substituted by $`k+1`$, and one must enforce the stability of the spin $`\mathrm{\Psi }`$. After the merging, the joint density probability of the value of spin $`\mathrm{\Psi }`$ and the local field $`H_\mathrm{\Psi }`$ acting on it is $`P_\mathrm{\Psi }(H,\sigma )`$, given by:

$$P_\mathrm{\Psi }(H,\sigma )=\frac{Z_1\mathrm{}Z_{k+1}}{Z}e^{\beta H\sigma }\mathrm{\Theta }(H\sigma )\underset{\sigma _1,\mathrm{}\sigma _{k+1}}{\text{Tr}}\delta (H_\mathrm{\Psi }H)\underset{i=1}{\overset{k+1}{}}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Psi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i).$$
(13)
Summing over all the possible $`(H,\sigma )`$ and taking the logarithm one gets the variation of the free energy of the system during the merging:

$`\mathrm{\Delta }F_1`$ $`=`$ $`{\displaystyle \frac{1}{\beta }}\mathrm{ln}\left[{\displaystyle \underset{\sigma }{}}{\displaystyle 𝑑He^{\beta H\sigma }\mathrm{\Theta }(H\sigma )\underset{\sigma _1,\mathrm{}\sigma _{k+1}}{\text{Tr}}\delta (H_\mathrm{\Psi }H)\underset{i=1}{\overset{k+1}{}}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Psi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i)}\right]`$ (14)
$`=`$ $`{\displaystyle \frac{1}{\beta }}\mathrm{ln}\left[\underset{\sigma ,\sigma _1,\mathrm{}\sigma _{k+1}}{\text{Tr}}e^{\beta H_\mathrm{\Psi }\sigma }\mathrm{\Theta }(H_\mathrm{\Psi }\sigma ){\displaystyle \underset{i=1}{\overset{k+1}{}}}{\displaystyle 𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Psi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i)}\right].`$
As far as operation (M2) is concerned, it is straightforward to derive the variation of the free energy

$$\mathrm{\Delta }F_2=\frac{1}{\beta }\mathrm{ln}\left[\underset{\sigma _1,\sigma _2}{\text{Tr}}e^{\beta J\sigma _1\sigma _2}𝑑h_1𝑑h_2\mathrm{\Theta }\left[(h_1+J\sigma _2)\sigma _1\right]\mathrm{\Theta }\left[(h_2+J\sigma _1)\sigma _2\right]p_1(h_1,\sigma _1)p_2(h_2,\sigma _2)\right].$$
(15)
One can easily deduce the density of free energy $`F`$ of the Bethe lattice from $`\mathrm{\Delta }F_1`$ and $`\mathrm{\Delta }F_2`$. Assuming that each branch has a number of spins equal to $`N`$, the system after the merging (M1) has $`N(k+1)+1`$ spins, and its free energy can by written in two ways:

$$(N(k+1)+1)F=(k+1)F_{branch}+\overline{\mathrm{\Delta }F_1},$$
(16)
where $`F_{branch}`$ is the free energy of one of the branches before the merging. For operation (M2), one has

$$(2N)F=2F_{branch}+\overline{\mathrm{\Delta }F_2}.$$
(17)
Elimination of $`F_{branch}`$ between the above equations yields

$$F=\overline{\mathrm{\Delta }F_1}\frac{k+1}{2}\overline{\mathrm{\Delta }F_2}.$$
(18)
Note that the above derivation holds for any extensive and self-averaging observable. Thus the density of energy $`E`$ of the Bethe lattice is

$$E=\overline{\mathrm{\Delta }E_1}\frac{k+1}{2}\overline{\mathrm{\Delta }E_2},$$
(19)
where $`\mathrm{\Delta }E_1`$, resp. $`\mathrm{\Delta }E_2`$, is the variation of energy under the merging process (M1), resp. (M2):

$`\mathrm{\Delta }E_1`$ $`=`$ $`{\displaystyle \underset{\sigma }{}}{\displaystyle 𝑑HP_\mathrm{\Psi }(H,\sigma )H\sigma }`$ (20)
$`=`$ $`{\displaystyle \frac{\text{Tr}_{\sigma ,\sigma _1,\mathrm{}\sigma _{k+1}}H_\mathrm{\Psi }\sigma e^{\beta H_\mathrm{\Psi }\sigma }\mathrm{\Theta }(H_\mathrm{\Psi }\sigma )_{i=1}^{k+1}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Psi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i)}{\text{Tr}_{\sigma ,\sigma _1,\mathrm{}\sigma _{k+1}}e^{\beta H_\mathrm{\Psi }\sigma }\mathrm{\Theta }(H_\mathrm{\Psi }\sigma )_{i=1}^{k+1}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Psi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i)}},`$
and

$$\mathrm{\Delta }E_2=J\frac{\text{Tr}_{\sigma _1,\sigma _2}\sigma _1\sigma _2e^{\beta J\sigma _1\sigma _2}𝑑h_1𝑑h_2\mathrm{\Theta }\left[(h_1+J\sigma _2)\sigma _1\right]\mathrm{\Theta }\left[(h_2+J\sigma _1)\sigma _2\right]p_1(h_1,\sigma _1)p_2(h_2,\sigma _2)}{\text{Tr}_{\sigma _1,\sigma _2}e^{\beta J\sigma _1\sigma _2}𝑑h_1𝑑h_2\mathrm{\Theta }\left[(h_1+J\sigma _2)\sigma _1\right]\mathrm{\Theta }\left[(h_2+J\sigma _1)\sigma _2\right]p_1(h_1,\sigma _1)p_2(h_2,\sigma _2)}.$$
(21)
These results can be applied to a wide variety of models: up to now we have not specified the coupling constants.

## IV The ferromagnetic case: An exactly solvable model

First, we concentrate on a simple case where all the above equations can be worked out analytically. The previous formalism allows for the computation of the metastable states of the $`k=2`$ Ising ferromagnet. When the couplings all have the same absolute value $`|J|=1`$, the cavity field can take only value in the set $`(2,0,2)`$. Thus integration with respect to the $`h_i`$ in Eq. (11) and other formulas reduces to a finite summation. Moreover, in the case of the Ising ferromagnet (i.e. $`J=1`$) the system is homogeneous, therefore the cavity equations should not depend on the site index and we can define the joint probabilities $`p_0=p(h=2,\sigma =1)`$, $`p_1=p(h=2,\sigma =1)`$, $`p_2=p(h=0,\sigma =1)`$, … , $`p_5=p(h=2,\sigma =1)`$ of Eq. (11) identically for all the lattice sites. A simple enumeration shows that Eq. (11) reduces to the following system of equations:

$`cp_0e^{2\beta }=p_1^2`$ $`cp_1e^{2\beta }=p_1^2+2p_1p_3+p_3^2`$ $`cp_2=2p_1p_2+2p_1p_4`$
$`cp_3=2p_1p_4+2p_3p_4`$ $`cp_4e^{2\beta }=p_2^2+2p_1p_3+p_4^2`$ $`cp_5e^{2\beta }=p_4^2,`$ (22)
together with the normalization condition $`_{i=0}^5p_i=1`$. The symmetries of the system fix some conditions on the values of $`p_i`$. The relations $`p_0=p_5`$, $`p_1=p_4`$, $`p_2=p_3`$ hold in the paramagnetic phase (PA). If the overall $`Z_2`$ symmetry is spontaneously broken the system encounters a ferromagnetic (FM) phase characterized by $`p_2=p_3=p_4=p_5=0`$ (obviously there is also the solution with the opposite magnetization, for which $`p_0=p_1=p_2=p_3=0`$). In the FM case the solution is given by

$$p_0=\frac{1}{1+e^{4\beta }}p_1=\frac{1}{1+e^{4\beta }}c=\frac{1}{e^{2\beta }+e^{6\beta }}.$$
(23)
The solution in the PA case is slightly more involved:

$`p_0`$ $`={\displaystyle \frac{c(x)D(x)^4e^{2\beta }}{16x^2}}p_1={\displaystyle \frac{c(x)D(x)^2}{4x}}p_2={\displaystyle \frac{c(x)D(x)^3}{2x}}`$ (24)
$`c(x)`$ $`={\displaystyle \frac{1}{2}}\left[{\displaystyle \frac{D(x)^4}{64x^5}}+{\displaystyle \frac{D(x)^3}{8x^3}}+{\displaystyle \frac{D(x)^2}{4x}}\right],`$
where $`x=(\frac{e^{2\beta }}{4})^{\frac{1}{3}}`$ and $`D(x)=\sqrt{x^42x}x^2`$. The problem of the iterative stability of both FM and PA solutions can be addressed studying the spectrum of the Jacobian of the system of equations (IV). It turns out that PA solution becomes unstable below a $`T_c^{PA}=\mathrm{log}^1(2\sqrt{\frac{2}{3}})=2.039091`$ and the FM one above $`T_c^{FM}=\frac{2}{\mathrm{log}(2)}=2.885390`$. Once we have calculated the cavity fields we can measure the thermodynamic potential using the merging procedures explained in the former section.

Let us consider first the FM solution.

$`\mathrm{\Delta }E_1=3`$ $`\mathrm{\Delta }F_1=3\left[1{\displaystyle \frac{1}{\beta }}\mathrm{log}(1+e^{4\beta })\right]`$
$`\mathrm{\Delta }E_2=1`$ $`\mathrm{\Delta }F_2={\displaystyle \frac{2}{\beta }}\mathrm{log}(1+e^{4\beta })1.`$ (25)
Inserting Eq. (IV) into Eq. (18) and Eq. (19) we readily obtain $`E=F=\frac{3}{2}`$, i.e. the system is completely frozen in its ground state.

The PA solution is worked analogously.

$`\mathrm{\Delta }E_1={\displaystyle \frac{N_1(x)}{D_1(x)}}`$ $`\mathrm{\Delta }F_1={\displaystyle \frac{1}{\beta }}\mathrm{log}(D_1(x))`$
$`\mathrm{\Delta }E_2={\displaystyle \frac{N_2(x)}{D_2(x)}}`$ $`\mathrm{\Delta }F_2={\displaystyle \frac{1}{\beta }}\mathrm{log}(D_2(x)),`$ (26)
where

$`N_1(x)`$ $`=`$ $`2p_1^3(3e^\beta +e^{3\beta })+6p_1^2p_2(e^{3\beta }+2e^\beta )+6p_1p_2^2(e^\beta +e^{3\beta })`$
$`D_1(x)`$ $`=`$ $`6p_1^3(e^\beta +e^{3\beta })+6p_1^2p_2(3e^{3\beta }+2e^\beta )+6p_1p_2^2(e^\beta +3e^{3\beta })+6p_2^3e^{3\beta }`$
$`N_2(x)`$ $`=`$ $`2p_1^2e^\beta 2(p_1+p_2)^2e^\beta `$
$`D_2(x)`$ $`=`$ $`2p_1^2e^\beta +2(p_1+p_2)^2e^\beta ,`$ (27)
and for the $`p_i`$ we use Eq. (24). In Fig. (4) we display both energy and free energy as a function of the temperature. Starting in the low temperature phase, the system gets trapped in the ferromagnetic solution in which all spins are aligned. This solution is locally stable for the iteration Eq.(IV) up to temperature $`T_c^{FM}`$. We define $`T_c=2.228512`$ as the temperature at which PA and FM solutions have the same energy. Exactly at this temperature we have the coexistence of the two phases of the system. Above $`Tc`$ the PA solution acquires a lower free energy. In Fig. (5) we display the configurational entropy $`S_C(E)`$ as a function of the energy $`E`$. The thin straight line is tangent to the curve exactly at $`S_C(E_c)`$ where $`E_cE(T_c)=1.05613`$, while the dotted curve beneath the tangent is the result of the microcanonical computation of Dean et al. presented in lefevredean . Note that the microcanonical branch can not be obtained in our canonical scheme.

## V The binary spin glass

Let us now turn to a slightly more complicated case, the binary spin glass, where the coupling constants are quenched independent identically distributed (i.i.d.) random variables with the following law:

$$\lambda (J)=\frac{1}{2}(\delta (J1)+\delta (J+1)).$$
(28)
As in the above section, an important simplification occurs: the local field acting on a spin can have only a finite number of values. The cavity fields can have $`k+1`$ values: $`k,k+2,\mathrm{},k`$, whereas the local field acting on a site with $`k+1`$ first neighbors can have $`(k+2)`$ values: $`k1,k+1,\mathrm{}k+1`$. But by contrast to the ferromagnet, the system is no longer homogeneous, and the self consistency equation (12) is a tricky object. To tackle a similar equation, mezpar  suggested a numerical solution via a population algorithm. The basic idea is that a probability distribution can be represented by a large collection of i.i.d. random variables distributed according to it.

In our case, the algorithm works as follows: we use a large population of $`𝒩`$ sites $`i=1,\mathrm{}𝒩`$, each of which is characterized by a finite set of $`2(k+1)`$ numbers, the $`p_i(h,\sigma )`$. After random initialization, one iterates the following sequence:

Algorithm I

* Select $`k+2`$ sites $`i_1,\mathrm{}i_{k+2}`$ at random, and extract some couplings $`J_1,\mathrm{}J_{k+1}`$.
* Perform the merging (M1) of the branches rooted at $`i_1,\mathrm{}i_{k+1}`$ onto the site $`\mathrm{\Psi }=i_{k+2}`$, that is compute $`\mathrm{\Delta }F_1`$ with Eq. (14), $`\mathrm{\Delta }E_1`$ with Eq. (20).
* Perform the merging (M2) of the branches rooted at $`i_1`$ and $`i_2`$, that is compute $`\mathrm{\Delta }F_2`$ with Eq. (15), $`\mathrm{\Delta }E_2`$ with Eq. (21), using for instance the coupling $`J_1`$.
* Update the population: perform the merging of the branches rooted at $`i_1,\mathrm{}i_k`$ onto the site $`\mathrm{\Phi }=i_{k+2}`$, that is substitute $`p_{k+2}(h,\sigma )`$ by the result of Eq. (11).

The algorithm converges in a stochastic sense: after a sufficient number of iterations, the $`p_i(h,\sigma )`$ are distributed according to $`𝒫`$. So one can compute the average of the thermodynamic quantities with respect to the couplings by ’time’ averaging of the measures performed in (ii) and (iii).

In the practical implementation of the algorithm we have scanned values of $`𝒩`$ ranging from $`200`$ to $`4000`$. A careful finite size scaling analysis shows that results are really mildly dependent on $`𝒩`$ as soon as $`𝒩>1000`$, and the asymptotic extrapolation is always compatible, within statistical errors, with the biggest size we have simulated. The errors are calculated with standard binning procedure, discarding the first half of the simulation. In this way we have a complete control on the “thermalization” of the algorithm. Most of the simulations have been performed also for the same system removing the stability condition on the local fields. We will always refer to these data as AC (all configurations) to distinguish them from the metastable (MS) set of data. For most of the simulations we have performed from $`10^3\times 𝒩`$ up to $`10^4\times 𝒩`$ iterations of Algorithm I, which we have verified to be enough both for thermalization and numerical accuracy.

The configurational entropy $`S_C`$ is the MS entropy.

One should also note that in the case of the Bethe lattice the choice of the distribution Eq. (28) is not crucial and the same results hold in the case of a distribution $`\lambda (J)=p\delta (J1)+(1p)\delta (J+1)`$ regardless how small $`p`$ is and even in the limiting case of a purely anti-ferromagnetic system. This observation allows for the analytic computation of the critical temperature and the critical energy of this model. Following the route specified in Sec. IV we calculated the equivalent of the system of equations (IV) for purely anti-ferromagnetic couplings. The stability analysis of the paramagnetic solution gives the same temperature of the FM case, i.e. $`T_c=\mathrm{log}^1(2\sqrt{\frac{2}{3}})`$ and $`E_cE(T_c)=15/141.071429`$. Below this temperature (energy) the replica symmetry is spontaneously broken.

In figure (6) we display the logarithm of the number of metastable states as a function of the energy. In the inset we zoom the results for the lowest energy. It is interesting to note that in the region $`E<1.22`$ both AC and MS entropy seem to behave linearly on $`E`$, and a linear fit works really well. The two curves meet at $`E=1.273(5)`$ and $`S=0.0172(5)`$. This energy is clearly compatible with the replica symmetric value of the ground state $`E_0^{RS}=23/18=1.277777`$ PARISI\_MEZARD\_86 . We have an excess of entropy at $`E_0^{RS}`$ which is proportional to the density of spin fous in the ground state, i.e. the number of spins with zero local field – such a spin can be flipped without changing the energy.

In figure (7) we compare our data with the annealed approximation presented in lefevredean , while in inset we zoom the difference $`S_{ann}(E)S(E)`$ in the region $`E<E_c`$. Above this value the annealed approximation is believed to be exact, and in fact our data fall on the analytic curve within the error bars, while below $`E_c`$ the two curves split. It is interesting to note that the splitting, barely visible in the main panel for $`E>1.2`$, seems to be exponential, at least near $`E_c`$ (note the logarithmic scale on the $`y`$ axis of inset).

Let us stress here that it is straightforward to adapt the algorithm to the case of random graphs of fluctuating connectivity of finite mean $`c`$. In such a graph, given two points, there is an edge connecting them with probability $`c/N`$, and no edge with probability $`1c/N`$. Thus the number of first neighbors of a given site is a random variable distributed according to a Poisson law of mean $`c`$. To implement this, one must extract the number $`k+1`$ according to this law in step (i). The rest is unchanged, except the important fact that step (iv) is no longer to be performed with $`k`$ branches, but with $`k+1`$ branches, like (ii) and that the $`k+1`$ term in Eqs. (18) and (19) should be replaced by $`c`$ MEZ\_ZEC .

## VI The Gaussian spin glass

Our method can tackle even more complicated cases, such as the Gaussian spin glass, in which the coupling constants are quenched i.i.d random variables whose law is a Gaussian with unit variance and zero mean.

The local fields are now continuous variables. We choose to represent the probability distribution $`p_i(h,\sigma )`$ on a site by a large population of $``$ couples $`(h,\sigma )`$ distributed according to $`p_i`$.

So on each site $`i=1,\mathrm{}𝒩`$, we have a population $`(h_i^\nu ,\sigma _i^\nu ),\nu =1,\mathrm{}`$. The algorithm’s layout is similar to algorithm I:

Algorithm II

* Select $`k+2`$ sites $`i_1,\mathrm{}i_{k+2}`$ at random, and extract some couplings $`J_1,\mathrm{}J_{k+1}`$.
* Perform the merging of the branches rooted at $`i_1,\mathrm{}i_{k+1}`$ onto the site $`\mathrm{\Psi }=i_{k+2}`$. In the context of this algorithm, the Eq.(14) and (20) become

$$\mathrm{\Delta }F_1=\frac{1}{\beta }\mathrm{ln}\left[\underset{\nu =1}{\overset{}{}}\underset{\sigma =\pm 1}{}e^{\beta H_\mathrm{\Psi }^\nu \sigma }\mathrm{\Theta }(H_\mathrm{\Psi }^\nu \sigma )\underset{j=1}{\overset{k+1}{}}\mathrm{\Theta }\left[(h_{i_j}^\nu +J_j\sigma )\sigma _{i_j}^\nu \right]\right]$$
(29)
$$\mathrm{\Delta }E_1=\frac{_{\nu =1}^{}_{\sigma =\pm 1}H_\mathrm{\Psi }\sigma e^{\beta H_\mathrm{\Psi }^\nu \sigma }\mathrm{\Theta }(H_\mathrm{\Psi }^\nu \sigma )_{j=1}^{k+1}\mathrm{\Theta }\left[(h_{i_j}^\nu +J_j\sigma )\sigma _{i_j}^\nu \right]}{_{\nu =1}^{}_{\sigma =\pm 1}e^{\beta H_\mathrm{\Psi }^\nu \sigma }\mathrm{\Theta }(H_\mathrm{\Psi }^\nu \sigma )_{j=1}^{k+1}\mathrm{\Theta }\left[(h_{i_j}^\nu +J_j\sigma )\sigma _{i_j}^\nu \right]},$$
(30)
where $`H_\mathrm{\Psi }^\nu =_{j=1}^{k+1}J_j\sigma _{i_j}`$.
* Perform the merging of the branches rooted at $`i_1`$ and $`i_2`$, using for instance the coupling $`J_1`$. The formulas (15) and (21) read:

$$\mathrm{\Delta }F_2=\frac{1}{\beta }\mathrm{ln}\left[\underset{\nu =1}{\overset{}{}}e^{\beta J_1\sigma _{i_1}\sigma _{i_2}}\mathrm{\Theta }\left[(h_{i_1}^\nu +J_1\sigma _{i_2}^\nu )\sigma _{i_1}^\nu \right]\mathrm{\Theta }\left[(h_{i_2}^\nu +J_1\sigma _{i_1}^\nu )\sigma _{i_2}^\nu \right]\right]$$
(31)
$$\mathrm{\Delta }E_2=J_1\frac{_{\nu =1}^{}\sigma _{i_1}\sigma _{i_2}e^{\beta J_1\sigma _{i_1}\sigma _{i_2}}\mathrm{\Theta }\left[(h_{i_1}^\nu +J_1\sigma _{i_2}^\nu )\sigma _{i_1}^\nu \right]\mathrm{\Theta }\left[(h_{i_2}^\nu +J_1\sigma _{i_1}^\nu )\sigma _{i_2}^\nu \right]}{_{\nu =1}^{}e^{\beta J_1\sigma _{i_1}\sigma _{i_2}}\mathrm{\Theta }\left[(h_{i_1}^\nu +J_1\sigma _{i_2}^\nu )\sigma _{i_1}^\nu \right]\mathrm{\Theta }\left[(h_{i_2}^\nu +J_1\sigma _{i_1}^\nu )\sigma _{i_2}^\nu \right]}.$$
(32)
* Update the population, performing the merging of the branches rooted at $`i_1,\mathrm{}i_k`$ onto the site $`\mathrm{\Phi }=i_{k+2}`$: one wants to substitute the $`(h_\mathrm{\Phi }^\nu ,\sigma _\mathrm{\Phi }^\nu ),\nu =1,\mathrm{}`$ with a new population distributed according to $`p_\mathrm{\Phi }(h,\sigma )`$ of Eq. (11). This is performed in two steps

+ First one builds a list $`(\stackrel{~}{h}^\nu ,\stackrel{~}{\sigma }^\nu ),\nu =1,\mathrm{}`$ as follows. Start with $`\nu =1`$. Enter a loop in which you extract randomly an $`\nu ^{}\{1,\mathrm{}\}`$ and a spin value $`\sigma =\pm 1`$, until you have the stability condition $`j\{1,\mathrm{}k\},(h_{i_j}^\nu ^{}+J_j\sigma )\sigma _{i_j}^\nu ^{}0`$. Then set $`\stackrel{~}{h}^\nu =_{j=1,\mathrm{}k}J_j\sigma _{i_j}^\nu ^{}`$ and $`\stackrel{~}{\sigma }^\nu =\sigma `$. Increment $`\nu `$, and enter the loop again.
    + The list $`(\stackrel{~}{h}^\nu ,\stackrel{~}{\sigma }^\nu )`$ is not the one to overwrite the $`(h_\mathrm{\Phi }^\nu ,\sigma _\mathrm{\Phi }^\nu )`$, because one must enforce the presence of the factor $`e^{\beta h\sigma }`$ in Eq. (11). Thus the elements of the list are to be reweighted: a suitable new population $`(h_\mathrm{\Phi }^\nu ,\sigma _\mathrm{\Phi }^\nu ),\nu =1,\mathrm{}`$ is obtained by repeating $``$ times the process of picking an element in the list $`(\stackrel{~}{h}^\nu ,\stackrel{~}{\sigma }^\nu ),\nu =1,\mathrm{}`$ with a probability proportional to $`e^{\beta \stackrel{~}{h}^\nu \stackrel{~}{\sigma }^\nu }`$.

We implemented the above algorithm in the case $`k=2`$, using the values $`𝒩=1000`$ and $`=1000`$. The results are presented in Fig. (8). By contrast with the binary spin glass, the ground state is not degenerated. So both entropies AC and MS go to 0 at the energy density of the ground state $`E_{GS}`$. Using the algorithm of mezpar  at low temperature and extrapolating the results to $`T=0`$, one finds $`E_{GS}=1.12(1)`$, which is compatible with our data.

## VII Breaking the replica symmetry

Let us now see how our method can be generalized to perform RSB.

### VII.1 General considerations

So far we have assumed there is a single pure state, stable under the iteration process. It is well known that this is actually not true for the spin glass on the Bethe lattice mezpar . One should allow for an infinity of pure states as described by the hierarchical continuous RSB scheme proposed by Parisi mpv . As a first approximation, one can implement the 1RSB scheme. Briefly speaking, one assumes that there exists an infinity of pure states labeled by $`\alpha =1,\mathrm{}+\mathrm{}`$. The free energies of the states on one branch are independent identically distributed random variables, with an exponential density

$$\rho (F)=\mathrm{exp}(\beta x(FF^R))$$
(33)
where $`F^R`$ is a reference free energy, and $`x`$ is Parisi’s parameter mpv . In this approach one computes the average values of the observables as a function of $`x`$. The physical value of $`x`$ is the one which maximizes the free energy.

If one considers a branch rooted at site $`\mathrm{\Phi }`$, in each pure state $`\alpha `$ the probability distribution $`p_\mathrm{\Phi }^\alpha (h,\sigma )`$ is different, and the iteration Eq. (11) is valid only inside a given pure state. The self-consistency equation to be solved is consequently more complicated than in the replica symmetric case. Let us write $`𝐩_i(h,\sigma )=(p_i^\alpha (h,\sigma ))_{\alpha =1,\mathrm{}}`$. One is interested in the set of probability distributions $`𝒬_\mathrm{\Phi }^{(h,\sigma )}(𝐩)=\overline{\delta (𝐩_\mathrm{\Phi }(h,\sigma )𝐩)}`$, which are functionals of the analogous quantities $`𝒬_i^{(h,\sigma )}`$, and one asks for all of them to be equal. A comprehensive theoretical description of the extension of the cavity method to 1RSB can be found in mezpar .

Let us turn to the description of the algorithm. We use a population of $`𝒩`$ sites $`i=1,\mathrm{}𝒩`$, $``$ states $`\alpha =1,\mathrm{}`$. On each site, each state is characterized by its own distribution $`p_i^\alpha (h,\sigma )`$. Depending on the distribution of the coupling constants, $`p_i^\alpha (h,\sigma )`$ is represented as a set of $`2(k+1)`$ numbers (case of binary couplings), or again as a population of $``$ couples $`(h,\sigma )`$ (case of Gaussian couplings).

The RS algorithm I is nested into the present algorithm: it is used to implement the iteration (11) and compute the expectation values of the observables inside each state. Then a meta-algorithm takes into account all the states, with appropriate weights, to update the population and compute the global expectation values. Two new observables are required with respect to the RS case. First the variation of the free energy during the iteration process described by Eq. (11):

$$\mathrm{\Delta }F_{iter}=\frac{1}{\beta }\mathrm{ln}\left[\underset{\sigma ,\sigma _1,\mathrm{}\sigma _k}{\text{Tr}}e^{\beta h_\mathrm{\Phi }\sigma }\underset{i=1}{\overset{k}{}}𝑑h_i\mathrm{\Theta }\left[(h_i+J_{\mathrm{\Phi },i}\sigma )\sigma _i\right]p_i(h_i,\sigma _i)\right].$$
(34)
Second the derivative of the free energy $`F`$ with respect to $`x`$. By derivation of Eq. (14) and (15), one gets

$$\frac{dF}{dx}=\frac{F}{x}+d^1\frac{k+1}{2}d^2,$$
(35)
where

$$d^1=\frac{1}{x}\frac{_\alpha \mathrm{\Delta }F_1^\alpha e^{\beta x\mathrm{\Delta }F_1^\alpha }}{_\alpha e^{\beta x\mathrm{\Delta }F_1^\alpha }},$$
(36)
and similarly for $`d^2`$.

An iteration of the algorithm goes as follows:

Algorithm III

* Perform step (i) of algorithm I.
* For each state $`\alpha `$, perform steps (ii), (iii) and (iv) (computing en passant $`\mathrm{\Delta }F_{iter}^\alpha `$), of algorithm I. One gets quantities bearing $`\alpha `$ as a superscript: $`\mathrm{\Delta }F_1^\alpha ,\mathrm{}`$.
* Reweight the states: this is a crucial step, motivated in mezpar . The states with low $`\mathrm{\Delta }F_{iter}^\alpha `$ have to be favored. Thus one picks up $``$ times an element in the list of distributions $`p_\mathrm{\Phi }^\alpha ,\alpha =1,\mathrm{}`$ with a probability proportional to $`\mathrm{exp}(\beta x\mathrm{\Delta }F_{iter}^\alpha )`$. The resulting list overwrites the one at site $`\mathrm{\Phi }`$. (Note that the elements one picks up are composite objects, i.e. either a set of $`2(k+1)`$ numbers or a population of $``$ couples).
* Compute the global average values of the observables $`\mathrm{\Delta }F_1\mathrm{}`$. The formulas can be found in mezpar :

$$\mathrm{\Delta }F_1=\frac{1}{\beta x}\mathrm{ln}\left[\frac{1}{}\underset{\alpha }{}e^{\beta x\mathrm{\Delta }F_1^\alpha }\right]$$
(37)
$$\mathrm{\Delta }E_1=\frac{_\alpha \mathrm{\Delta }E_1^\alpha e^{\beta x\mathrm{\Delta }F_1^\alpha }}{_\alpha e^{\beta x\mathrm{\Delta }F_1^\alpha }},$$
(38)
and their obvious analogues for $`\mathrm{\Delta }F_2`$ and $`\mathrm{\Delta }E_2`$.

Note that it is essential that the sites and couplings be the same for all the $`\alpha `$ in step (ii) above.

In practice, we found that the implementation of this algorithm in the case of the Gaussian couplings requires too many computer resources to reach a satisfying accuracy. So we limited ourselves to the case of binary couplings.

### VII.2 Application to the binary spin glass

We considered the binary spin glass of section V, at temperature $`T=0.5`$. Our implementation used the values $`𝒩=1000`$ and $`=1000`$. We computed $`dF/dx`$ for several values of $`x`$. The result is presented is Fig. (9). To determine the value $`x^{}`$ of $`x`$ where $`dF/dx`$ is zero, we fitted the curve by a polynomial of degree 4, and found its roots. We got $`x^{}=0.19(1)`$. Then we measured the thermodynamic quantities for $`x=x^{}`$, taking into account the uncertainty on $`x^{}`$:

$$F=1.280\pm 0.001,E=1.265\pm 0.001,S=0.029\pm 0.001.$$
(39)
This is to be compared to the output of the replica symmetric algorithm of section V:

$$F=1.2816\pm 0.0008,E=1.2685\pm 0.0006,S=0.0274\pm 0.0008.$$
(40)
The overall improvement is small, particularly on $`F`$. It is more obvious on $`E`$ and $`S`$ but still only of the order of $`10^3`$.

## VIII Conclusion

In this paper we have described a general procedure to compute the number of 1-spin flip stable configurations on the Bethe lattice. Given some integer $`k2`$, the method – at least conceptually – can be easily generalized to compute the number of $`k`$-spin flips stable configurations, i.e. whose energies can not be decreased by flipping a number of spins ranging from 1 to $`k`$. The practical difficulty is that the recursion relations can no longer involve only quantities related to the root of a branch, but must also take into account the $`k1`$ successive generations of spins. In the case $`k=2`$ it remains quite straightforward to write down the equations, but we have not worked out their solution. It might be interesting to do so in order to clarify the nature of the zero temperature limit of the pure states: recently Mézard and Parisi PARISI\_MEZARD\_T0  presented a computation at 1RSB level of the number of locally ground states (LGS), i.e. configurations stable with respect to $`k`$-spin flips with the number $`k`$ going to infinity with the size of the system in some unprecised way, and they found some surprising features. There is a need for more precise definitions, and one of the points at issue is how these LGS are related to $`k`$-spin flips stable configurations when $`k+\mathrm{}`$.

ACKNOWLEDGMENTS

We acknowledge very useful discussions with F. Ricci-Tersenghi.