{'fraction_non_alphanumeric': 0.03938447025529155, 'fraction_numerical': 0.009413753119902522, 'mean_word_length': 4.108835341365462, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The cavity method at zero temperature

In this note we explain the use of the cavity method directly at zero temperature, in the case of the spin glass on a Bethe lattice. The computation is done explicitly in the formalism equivalent to 'one step replica symmetry breaking'; we compute the energy of the global ground state, as well as the complexity of equilibrium states at a given energy. Full results are presented for a Bethe lattice with connectivity equal to three.


I. INTRODUCTION
The cavity method, initially invented to deal with the Sherrington Kirkpatrick model of spin glasses [1,2], is a powerful method to compute the properties of ground states in many condensed matter and optimization problems. It was originally developed at finite temperature, but if one is mainly interested in the optimization problem concerning the structure of ground states, it turns out to be easier to apply it directly at zero temperature, where many of the concepts can be explained in a more straightforward way. The aim of this note is to discuss the use of the cavity method at zero temperature. In order to be definite, we will study the concrete case of a spin glass on the Bethe lattice (the precise definition of the Bethe lattice which we use for spin glasses is given in sect. II), although the scope of application of the method is much wider; indeed the best applications of the formalism we present here are those dealing with optimization problems [3], like for instance the K-satisfiability problem [4] or some of its variants [5].
The cavity method is in principle equivalent to the replica method, but it turns out to have a much clearer and more direct interpretation, that allows in practice to find solutions to some problems which remain rather difficult to understand in the replica formalism: the replica approach is very elegant and compact, but it is more difficult to get an intuitive feeling of what is going on. Also, the cavity approach deals with usual probabilistic objects, and can lend itself to rigorous studies [7,8]. We shall present it here at two successive levels of approximation. The first one, corresponding in replica language to the replica symmetric (RS) solution, is an easy one and has already been studied a lot. The main aim of the paper is to explain in some details how one can solve the problems at the level of a 'one step replica symmetry breaking' (1RSB). For years, this was only possible for systems with infinite range interactions. We have recently found a general procedure allowing to get this 1RSB solution for problems defined on the Bethe lattice [9], and here we explain in some details its use directly at zero temperature. We shall give an explicit example of this 1RSB solution for the Bethe lattice spin glass with connectivity equal to three.
The paper is organised as follows: in sect. II we introduce the Bethe lattice spin glass, the basic ideas of iterative techniques, and the main questions which one wants to answer. Sect. III presents the cavity method when one assumes that there is only one local ground state (the RS approximation). The 1RSB solution is described in sect. IV. The paper ends with some remarks (sect. V), and is complemented by two appendices: Appendix A shows how the present zero temperature work is connected to the limit T → 0 of finite temperature analyses, and how this can be used to estimate the validity of a given level of RSB. Appendix B explains how the cavity equations can be interpreted from a variational point of view. Appendix C contains some miscellaneous comments on the definition of the local ground states.

II. GENERALITIES
We consider an Ising spin glass model on a special class of random lattices which we call Bethe lattice. This is defined here as as a random lattice with fixed connectivity k + 1. Such graphs locally look like a portion of a Cayley tree, but they display loops at large distances (of length of order log(N )/ log(k), where N is the number of points in the lattice). In unfrustrated problems, the Bethe lattice is defined as the interior of a large Cayley tree, but this definition is not appropriate for frustrated systems like spin glasses, where the presence of loops is essential to insure the existence of frustration. It has been argued in [9] that the natural generalization of the Bethe lattice to frustrated systems is provided by the random lattices with fixed connectivity which we study here. Notice that in these systems, the frustration and the disorder are due to the presence of loops, and thus occur only on large scales. Locally, the structure around almost all point is tree-like.
Spins are located at the vertices of the graph, and interact with neighbouring spins with exchange couplings. In the spin glass case the Hamiltonian is defined as: (1) The sum is over all links of the lattice (there are k + 1 links incoming onto each site i). For each link < ij > the coupling J ij is an independent random variable chosen with the same probability distribution P (J). One of our aims is to compute, in the infinite N limit, the value of the energy density of the global ground state (GGS), which is the configuration of Ising spins σ i = ±1 which minimizes the Hamiltonian. More precisely the ground state energy of a N spin system, averaged over the distribution of samples (i.e. both over the choices of random graphs, and the values of the couplings) will be denoted by E N . We want to compute We are also interested in computing, in the same limit, the number of local ground states (LGS) with a given energy, where a LGS is defined as a configuration whose energy cannot be decreased by flipping a finite number of spins (when N → ∞).
The basic locally tree-like structure is best exploited by an iterative method that has been called, in the context of spin glasses, the cavity method. Let us introduce an intermediate object which is a spin glass model with N spins, on a slightly different random lattice, where q randomly chosen 'cavity' spins have only k neighbours, while the other N − q spins all have k + 1 neighbours (see fig.1). We call such a graph a G N,q 'cavity graph'. The cavity spins are fixed, their values are σ 1 , ..., σ q . The GGS energy of the corresponding spin glass model obviously depends on the values of the cavity spins. 1. An example, for the case k = 2, of a GN,6 cavity graph where q = 6 randomly chosen cavity spins have 2 neighbours only. All the other N − 6 spins outside the cavity are connected through a random graph such that every spin has k + 1 = 3 neighbours.
While our primary interest is in the ground state configurations on G N,0 graphs, the intermediate construction of G N,q is helpful. The basic operations which one can perform on cavity graphs are the following: • Iteration: By adding a new spin σ 0 of fixed value into the cavity, connecting it to k of the cavity spins say σ 1 , ..., σ k , and optimizing the values of these k spins, one changes a G N,q into a G N +1,q−k+1 graph: • Link addition: By adding a new link between two randomly chosen cavity spins σ 1 , σ 2 , and optimizing the values of these 2 spins one changes a G N,q into a G N,q−2 graph: • Site addition: By adding a new spin σ 0 into the cavity, connecting it to k + 1 of the cavity spins say σ 1 , ..., σ k+1 , and optimizing the values of the k + 2 spins σ 1 , ..., σ k+2 , one changes a G N,q into a G N +1,q−k−1 graph: In particular, if one starts from a G N,2(k+1) cavity graph and perform k + 1 link additions, one gets a G N,0 graph, i.e. our original spin glass problem with N spins. Starting from the same G N,2(k+1) cavity graph and performing 2 site additions, one gets a G N +2,0 graph, i.e. our original spin glass problem with N + 2 spins. Therefore the variation in the GGS energy when going from N to N + 2 sites (E N +2 − E N ) is related to the average energy shifts ∆E (1) for a site addition, and ∆E (2) for a link addition, through: Using the fact that the total energy is asymptotically linear in N , the energy density of the ground state is finally An intuitive interpretation of this result (for even k + 1) is that in order to go from N to N + 1 one should remove (k + 1)/2 links (the energy for removing a link is minus the energy for adding a link) and then add a site. There are several alternative derivations and expressions of this energy density [10], but the previous one is the simplest one for our purpose. When q/N ≪ 1, generically, the distance on the lattice between two generic cavity spins is large (it diverges logarithmically in the large N limit). It is thus reasonable to assume that the various cavity spins become uncorrelated. This is the basic assumption of the RS solution, although it is in general not explicitly formulated. Several papers have worked out this RS solution in details at all temperatures [11][12][13][14][15][16][17][18][19]. Here we are presenting for completeness the zero temperature version of this RS approach, the modifications induced in the 1RSB case will be explained below. This RS assumption amounts to saying that the GGS energy of a G N,q spin glass can be written as an additive function of the values of the q cavity spins: The quantities h i , that we call the local cavity fields, depend on the sample. When considering the ensemble of random cavity graphs, they are independent identically distributed (iid) random variables, and their distribution is denoted P(h). The computation of P(h) will be a crucial step in our approach. The reader should notice that in general there is no simple expression for the local cavity fields of the type: h i ≡ k J ik σ k : These fields are related to the difference in energy of two GGS with flipped cavity spins, and these two GGS may in principle differ in an arbitrarily large number of spins. Indeed the quantity E({σ}) is computed by minimizing the energy as function of the other N − q spins for fixed values of the q cavity spins. The key hypothesis of the RS treatment is that the GGS of the spin glass before and after any of the previous graph operations (e.g. iteration) are related. Equivalently one should assume that the perturbation corresponding to the variation of one of the cavity spins remains localized and it does not propagate to the whole lattice. Under this hypothesis (whose consistency can be checked and, as we shall see later, is not always correct) a self consistent equation for this distribution is easily found by considering the iteration procedure (see fig.3). We add a new spin in the cavity, fix its value to σ 0 , and connect it to the spins σ 1 , ..., σ k through a new set of coupling constants J 1 , ..., J k . Each of the spins σ i , i ∈ {1, ..., k}, sees a local magnetic field h i + J i σ 0 , and has to point into the direction of this field in order to minimize its energy (if the local field is zero this spin is free to point in any direction). The energy of the link between σ 0 and σ i is thus: where an elementary computation gives: This shows in particular that the new local field on site 0, i.e. the coefficient of σ 0 in the expression for the GGS energy, is given by [20]: This implies a recursion relation for the distribution P(h): (Throughout this paper, we use a generic notation E J [·] to indicate the average over the distribution of all coupling constants and over random graphs). Let us suppose that the previous equation is sufficient to determine the distribution P(h). If we know this distribution we can compute the energy. The average energy shift due to a site addition is and the average energy shift due to a link addition is These expressions give the ground state energy density through (7).

B. Application to the ±J model
In order to obtain the RS solution, one must first solve the functional equation (12) for P(h). This is particularly easy in the case where the couplings are taken from the distribution P (J) = (1/2)(δ(J − 1) + δ(J + 1)). The definition of the local fields implies that they are integers, and the functions u and a take simpler forms: where S(h) = 0 if h = 0, S(h) = Sign(h) otherwise. In other words the function S(h) is −1 or +1 depending on the sign of h and it is zero for h = 0. We can write the most general P(h) in the form Furthermore the symmetry of the J distribution, together with (12), implies that P(h) is symmetric (p r = p −r ). The property of the function S(h) implies that what matters is only the sign of h. The probabilities p + = p − for having a positive or negative h are obviously given by (1 − p 0 )/2, where in our notation (16) p 0 is probability of having a zero h. The self consistency condition (12) now becomes a closed equation for p 0 , (where C q k = k!/(q!(q − k)!)) and the other coefficients are given by If we specialize for simplicity to a Bethe lattice with k = 2, the solution is p 0 = 1/3, p 1 = 2/9, p 2 = 1/9, and the ground state energy is E = −23/18 ≈ −1.278. All these results were obtained in [17].
We would like to comment here on the fact that the fields are integers. It was noticed already long ago that, even in the simple RS case, the iteration equations (12) admit some solutions which are not distributed only on integers, but have a continuous part; in fact a lot of efforts have been devoted to find 'the best' among these solutions. Moreover, if a computation is done at finite temperature, where no ambiguity is present, and the results are extrapolated at zero temperature, a continuous part is found in this RS case. A detailed computation shows that if, working at zero temperature, one starts from fields that are integer plus a small correction, this correction is amplified under the iteration procedure until the final distribution is no more concentrated near the integers. In our opinion it is quite likely that this instability is a signal that the RS solution is incorrect. A similar effect can be found on the 1RSB solution but it is weaker (see appendix A). We believe that whenever one reaches the correct solution (which in the present case should be full RSB), this artefact will disappear. Roughly speaking the rational for this belief is the following: the fact that a small variation of the field in one point propagates and leads to a large effect throughout the whole lattice is precisely the so called "replicon" instability (that would make the whole approach inconsistent). The RSB solution was invented just to cure this problem and we know that it is very successful in doing this [2].

A. Basic hypotheses of the cavity method
It is well known that the above result is wrong [21][22][23][24][25]9], because the hypothesis of continuity of the global ground state when we add new spins or links to the graph is incorrect.
A mechanism leading to the breaking of this continuity hypothesis, which is not taken into account in the RS approximation, is the existence of several local ground states (LGS, not to be confused with GGS). By a LGS here we mean a state whose energy cannot be lowered by flipping a finite number of spins. Of course such a definition applies only to the N = ∞ limit. In principle one should define the LGS for a finite system as being stable with respect to flipping a number of spins less than f (N ), where f (N ) is a well chosen increasing function of N , diverging at large N . Determining the most appropriate form of f is a difficult problem. In its simplest version, the cavity method assumes the existence of a function f (N ) such that there are many ground states, but does not try to prove it. Controlling this point would allow one to turn the cavity into a full mathematical proof, and has been done so far in only a few non trivial cases [7,8]). Here we shall not try to provide this kind of rigorous treatment, we just want to show how the cavity method works, within its usual assumptions, and check it through the validity of its results. Some further comments on the definition and counting of LGS are displayed in appendix C. Although the precise definition of LGS is a very subtle issue, we want to underline that there exist at least one model, rather closely related to our Bethe lattice spin glass problem, where the present one step RSB cavity method can be applied [5]. In this system the LGS can be characterized and enumerated fully and the results can be checked versus exact computations done using completely different methods [5,6].
In presence of several ground states, the assumption is that there is a one-to-one correspondence among the LGS before and after the addition of spins or links (at least for the LGS with low energies). Equivalently we assume that the perturbation due to the change of the value of a cavity spin propagates (in the limit N going to infinity) only to an infinitesimal fraction of the lattice. Therefore it is possible to write an iteration procedure for the whole population of LGS with given energy. However it may well be that the order of the LGS energies change during the graph operations, and the GGS after iteration is not the same LGS as the one before. The problem is to take into account these 'level crossings', which is not done in the RS solution of the previous section, and is done in the RSB solutions. Below we shall describe one such solution, which is equivalent to what is called in replica language the 1RSB solution. One is forced to follow a large population of the LGS of lowest energy, large enough so that one can be sure to obtain the GGS when iterating.
Let us state here the postulates of the cavity method, at the level equivalent to the 1RSB solution in the replica language.
• The cavity spins are uncorrelated within one given LGS. Labelling by α a LGS of a G N,q spin glass, its energy E α is given by: The previous equation would be wrong for the GGS because the GGS may correspond to different LGS depending on the value of the cavity spins σ. This extreme sensitivity of the GGS to the variation of a single variable is typical of hard optimization problems.
• The energies E α 0 of the LGS of low energy (near to that of the GGS) are assumed to be iid variables with a distribution given by a Poisson process of density where E ref is a reference energy, which is near to the GGS energy, and µ is a parameter whose physical interpretation will be explained in sect. IV D. This hypothesis is compatible with the idea that LGS are extremes of the energies, and with the Gumbel universality class for extremes [26,27].
• On a given site i, the local cavity fields in the various states, h α i , are iid variables taken from the same distribution P i (h). However the distribution P i (h) fluctuates from site to site, so that the correct order parameter is a functional Q[P (h)] giving the probability, when one picks up a site at random to find on this site a cavity field distribution P i (h) = P (h) [28]. Moreover the cavity fields and the LGS energies are not correlated.
These hypotheses are by no means evident and are likely to be wrong in many models where subtle correlations exist among the different cavity field. This more complex situation can appear, for instance, in cases where higher order breaking of the replica symmetry is present and it will not be discussed here. One should note that, while the 1RSB which we discuss here has a nice interpretation in terms of LGS having independent energies (and thus a Gumbel distribution), the higher order RSB solutions involve correlated variables and thus describe new universality classes of extremes. In the spin glass case which we study here, we expect that an infinite order of RSB will be needed to solve the problem (as is well known for the large k case [29,2], which is equivalent to the Sherrington Kirkpatrick (SK) model [30]). The 1RSB which we explore is thus an approximation to the problem, which usually produces better results than the RS approach. In the limit where k goes to infinity it is know the the GGS energy per spin (divided by a normalization factor k −1/2 ) is -.798 in the RS case, -.766 in the 1RSB , to be compared with the conjectured correct value -.763.
It is interesting to notice that the RS case of only one LGS can be obtained in this framework by fixing the local distribution to be P i (h) = δ(h − a i ), so that all the LGS are automatically equal. On a given site i, this distribution is fixed by the single number a i . The various a i 's are iid, taken from a distribution P(a) which satisfies exactly the RS recursion relation (12), and the result is µ independent, as it should. Alternatively in the limit where µ → 0, the number of different LGS does no more increase with the energy, the gap between the two lowest LGS diverges, and thus level crossings can be neglected. We shall see that in this limit we again recover the RS value of the average energy.

B. Self consistency under iteration
Let us study the effect of the iteration, starting from these hypotheses. When iterating, the new local cavity field, in one given LGS labeled by α, is given by the same expression as (11): and the energy of the LGS α is globally shifted from E α The functions u and a are given in (10). In this case the energy shifts and the local fields on the new spin σ 0 are correlated. Given a site (and thus fixing the couplings J i ), the pairs (h α 0 , ∆E α 0 ) are iid, taken from a distribution P 0 (h 0 , ∆E 0 ) which can be read-off from the iteration equations above: LGS energy after adding the new spin σ 0 . We first note that E α 0 ′ are iid, described again by a Poisson process with an exponential density: LGS with lowest energies, we need to keep only the states with the lowest new energy). The joint distribution R 0 (h 0 , E ′ 0 ) of the local field and the new LGS energy is given by: where the constant C being fixed in such a way that P ′ 0 (h 0 ) is a normalized probability distribution. The crucial point is the following: for each given LGS the recursion for the probability distribution of the h's is the naive one, i.e. the previous one with µ = 0. However if we look at the recursion for the probability distribution of the h's, conditioned to a given value of the LGS energy we find that it depends crucially on µ. The probability distribution of the cavity field in a generic LGS is different from the probability distribution of the cavity field for a given value of the LGS energy. The two quantities may be different because we are looking at the LGS of low energy which are for large N an infinitesimal fraction of the total number of local ground states (indeed we shall see later that the total number of local ground states increases exponentially with N ).
This distribution P ′ 0 of the cavity field on the new site depends on the added site through the choice of the various P i 's and of the coupling constants J i . When one averages over these choices, the distribution of P ′ 0 (h) must be equal to the original Q[P (h)]: when this self-consistency equation is satisfied, the hypotheses of the cavity method on the N spin system are found to be true again in the N + 1 spin system obtained after iteration.

C. Computation of the energy
In order to compute the GGS energy, we must see the effect of adding a site or adding a link within this 1RSB scenario.
Site addition: the energy shift during a site addition is computed through the same method as (24), the extra term coming from the fact that one optimizes the value of the added spin: Link addition: Note that one gets back the RS expressions when µ → 0. Once the distribution Q[P (h)] of the cavity field distributions is known, one can compute the energy shifts ∆E 1 , ∆E 2 which depend on µ, and deduce from them the energy function Φ(µ) defined as in (7): One finally finds: where dQ[P ] ≡ Q[P ]dP denotes the probability distribution of the probability P (the integral is done over all the probability distributions P with weight Q[P ]). Explicitly, one has: Notice that all this program can be carried out for any given µ. In order to fix the parameter µ which gives the increase rate of the number of LGS, one must compute Φ as a function of µ and maximize it. One justification for that procedure can be to go back to the replica formulation: µ turns out to be the zero temperature limit of m/T , where m is the breakpoint in Parisi's order parameter function at the 1RSB level [2]. The necessity of maximising over m is a well known feature of the replica method (in the limiting case k → ∞ it can be rigorously proven [32] that the value of Φ(µ) is a lower bound to the correct result so that it is natural that the preferred value of µ is obtained by maximizing Φ(µ)), but this does not provide a clear physical reason. The next section will explain the physical meaning of µ, and explain why µ must be equal to µ * , chosen such that Φ(µ) is maximum, in order to compute the GGS energy.
It is also interesting to notice that the whole self-consistency procedure obtained by the iteration equation (26) is a variational procedure which can be deduced from the stationarity condition of the functional Φ(µ) given in (31) with respect to changes of Q[P ]. This is shown in appendix B. The existence of this variational formulation is crucial in simplifying the computation of all the derivatives with respect to the various parameters, because only explicit derivatives must be considered.

D. Computing the complexity
In presence of many LGS, one can be interested in knowing their number. We shall make the basic assumption that for large N , the typical number N N (E) of LGS with a given energy E behaves as The function Σ(ǫ) is a positive function called the complexity of the problem. In the range of energy densities ǫ where Σ(ǫ) > 0, the log of the number of LGS is supposed to be a selfaveraging function, so that for almost all sample (in the large N limit), (1/N ) log(N N (E)) is given by Σ. Technically, given a definition of the LGS (see the discussion in app. C), one can thus define the complexity as Σ(ǫ) = lim N →∞ (1/N )E J log(1 + N N (N ǫ)). In general, for ground states α with a fixed energy density ǫ = E α /N , the local cavity fields h α i are iid variables taken from the distribution P (ǫ) i (h), which now depends on ǫ, and fluctuates from site to site. The corresponding order parameter is a functional Q (ǫ) [P ] giving the probability, when one picks up a site at random to find on this site a cavity field distribution (for states with energy density ǫ): P Our goal here is to check the self consistency of these hypotheses under iteration, and to determine the 'complexity' function Σ. One expects that Σ(ǫ) will vanish at some value ǫ 0 , which gives the GGS energy density, and the analysis of the previous section applies for LGS having energy close to N ǫ 0 (plus terms of order one). In this context the value of U can be found by looking for the solution of the equation Let us first study the effect of the iteration. The number of LGS at energy E after iteration is given by We need to compute the distribution of local cavity fields P (ǫ) 0 (h) at a fixed GS energy density E/N = ǫ. Expanding to first order in the small shift ∆E/N in (34), we obtain: where Notice that we get back exactly the same expression as in the previous section (see (26)), but now µ has a well defined meaning: it is the derivative of the complexity with respect to the energy density ǫ. When one varies ǫ, the value of µ changes. The whole formalism of the previous section can be used for any µ. Its result will give the distribution of cavity fields for the states α which have energy density ǫ, which is related to µ through (36). The functional probability distribution found in this way is Q[P (ǫ) (h)]. In order to compute the complexity function Σ(ǫ), we must see the effect of adding 2 sites or adding k + 1 links to a G N,2k+2 graph within this 1RSB scenario. Let us call Σ 2k+2 (ǫ) the complexity function of the G N,2k+2 graph from which we start.
Site addition: By adding one site, we go to a G N +1,k+1 graph with complexity Σ k+1 (ǫ). To compute this, we notice that when one site is added and the corresponding spin is optimised, the energy of the LGS is shifted by a value ∆E, and the corresponding probability distribution P site (∆E) is: After the site addition, the new complexity Σ k+1 (ǫ) is: Link addition: By adding one link, we go to a G N,2k graph with complexity Σ 2k (ǫ). The probability distribution for the corresponding energy shift is: After the link addition, the new complexity Σ 2k (ǫ) is: Going from a N site graph to a N + 2 site graph: As before, we need to add two sites and take away k + 1 links. When performing two times the site addition, we go to a graph G N +2,0 where the number of GS is while performing k + 1 times the link addition leads to a graph G N,0 where the number of LGS is exp (N Σ(E/N )). Adding up the effect of the link additions and the site additions, one obtains: where the average E J also implies an averaging over the choices of P i (h i ) appearing in P site and P link (the functions P i have to be taken from the functional distribution Q[P (h)]). Comparing to the expression for the energy density derived previously (31), one gets simply: So the knowledge of the function Φ(µ) allows to reconstruct the complexity function Σ(ǫ) through a Legendre transform [31]. The usual relations for such a transform, show that the condition of maximisation of Φ(µ) is equivalent to having Σ = 0. This proves that the GGS energy density is found by maximising Φ(µ) with respect to µ. It also gives the practical way of deducing the complexity from the knowledge of Φ(µ). Because of the structure of (43), we call Φ(µ) the zero temperature free energy of the problem. In fact the previous equations could also be written as where α labels the LGS and E α is the total energy of the LGS. In other words Φ(µ) is the free energy (at a temperature 1/µ) of a system where the sum is done only on the LGS and not over all the spin configuration (as in the case of the usual free energy). One could also apply the same techniques used at finite temperature [9] to rederive the previous expression for Φ(µ) and Σ(µ) in a slightly different way. Let us give here for completeness the full expressions used in the computation of the complexity. Once one has solved the iteration equations at a given µ, and obtained the functional probability distribution Q[P ], one can deduce Φ(µ) from (31). The energy density ǫ can then be obtained as: where we have defined the function D µ (x) = xe µx . The previous equations can be rewritten with an evident meaning of the symbols as: Finally the complexity is where Φ is given in (31) and ǫ is given in (46). The smart reader may cast some doubts on the consistency of the whole approach: our first hypothesis was the existence of a one-to-one correspondence between the ground states when increasing N by two units (and therefore the number of ground states should no depend on N ) and from this hypothesis we arrive to the conclusion that the number of ground states increase with N as exp(N Σ(ǫ)). We believe that the hypothesis of one-to-one correspondence of the LGS is valid only in a low energy region, where the computation that we have presented is correct, and it should be modified in the region where Σ(ǫ) is near to its maximum.

E. Factorized case
A particularly simple case is when we assume that the distribution P i (h) are i independent (i.e. the Q[P (h)] is a functional δ function). This case, which was first studied in [33], and developed for the Bethe lattice spin glass in [34], is named the 'factorized case' because of the special pattern of RSB to which it leads. It is simple because the order parameter is a single function P (h) (fixed from the self-consistency relation (26) by imposing P ′ 0 (h) = P 1 (h) = P 2 (h) = ... = P k (h)), and the RS equation are only slightly modified. However one should note that, in general, one may expect a P i (h) which fluctuates: this is obviously the case whenever the connectivity fluctuates, but this may also happen in the case of the fixed connectivity random graphs which we study here. Some special models where the factorized Ansatz gives an exact solution have been studied recently [35,36].
F. Application to the ±J model

Factorized case
Our first task is to compute the distribution Q[P (h)]. We will do it here first in the simple factorized case where P i (h) does not depend on the site (so that Q[P (h)] is a functional δ). The distribution P i (h) ≡ P (h) still has the form of a sum of delta peaks on the integers P (h) = r p r δ(h − r). This factorized solution can be stable only if P (h) is symmetric (p r = p −r ). The equations for the weights of the δ peaks p 0 , p 1 , ..., p k , taking into account the free energy shift, are: (the function S has already been defined as S(h) = 0 if h = 0, S(h) = Sign(h) otherwise). The quantity A is a normalisation constant fixed by the condition One can notice the similarity of (49) with the RS solution (16). With respect to this RS solution, the self consistency equations are only modified by the reweighting term in exp(µ...). However the physical interpretation of P (h) studied here and the P(h) of the RS solution are totally different. The present P (h) gives, for one site i, the probability that in a LGS α chosen at random, with energy density ǫ fixed by the value of µ, the field h α i is equal to h. In the RS case, there is only one LGS, and P(h) governs the fluctuation of h i from site to site.

The full 1RSB solution
Let us now go beyond the factorized approximation. On a given site j the probability distribution of the cavity field is It is parametrized by a vector of weights, p j = (p j −k , ..., p j k ) which can fluctuate from one site j to the next. Because we work at zero temperature and the fields take integer values, the probability depends on a finite number (2k) of parameters, and the full order parameter is not a functional, but a function R( p) of the vector of weights p = (p −k , ..., p k ), which is given in the limit of large N by: Let us perform one iteration, by adding a spin σ 0 in the cavity, connecting it to the k cavity sites j 1 , ..., j k , and optimizing the values of σ j1 , ..., σ j k . The distribution of the cavity field on site 0 again takes the form (55), with: where A 0 is a site-dependent normalisation constant, fixed by r p 0 r = 1. The energy shifts are given by (27) and (28), and are also site dependent. The iteration (57) offers a way to obtain the order parameter R( p) by a method of population dynamics [9]. The method amounts to following a population of N vectors p 1 , ..., p N . Each iteration of the dynamics corresponds precisely to the iteration of the cavity method: one picks up k vectors p j1 , ..., p j k at random in the population, one computes the new vector p 0 according to (57), and substitutes the new factor p 0 in the population, in the place of a randomly chosen vector p i . This stochastic process yields, after some transient regime, a population whose distribution is stationary, and the probability distribution of p inside this population is nothing but the R( p) order parameter. The energy shifts ∆E (1) and ∆E (2) obtained when adding one site or one link are computed as before, and they are then averaged over many iterations. This allows to compute, for any µ, the zero T free energy Φ(µ) using (29), which is shown in fig. 4. It has a maximum at a value µ ∼ .7, and predicts a GGS energy density U ≃ −1.2717.
One should notice that in this problem the quantitative effect of RSB on the GGS energy density are rather small. The full 1RSB result U ≃ −1.2717 differs from the RS result U ≃ −1.2777 by 5.10 −3 , and the factorized solution gives a very good approximation U ≃ −1.2723, precise at 5.10 −4 . Therefore one can expect that the quantitative effects of higher order RSB on the GGS energy will be tiny. However, the qualitative consequences of RSB are clearly visible on the fact that the probability distribution of cavity fields is site dependent. This is exemplified by fig. 5, which gives the site to site fluctuations P (p 0 ) = (1/N ) i δ p i 0 − p 0 of the probability of having a zero cavity field. We have also checked that, while the individual cavity field distributions P j ( p) are not symmetric under field reversal (i.e. P j (p r ) = P j (p −r )), the full order parameter is statistically symmetric (i.e. the site to site fluctuations of p r are identical to those of p −r ).
Using the previous results one can compute the complexity function Σ(ǫ). In fig.6 we show the results in the region where the complexity is positive. In order to obtain these results we only need the function Φ(µ) in the region to the left of its maximum, where dΦ/dµ ≥ 0. The complexity is rather similar in the factorized method and in the true solution, apart from some small variation of the ground state energy. However it shows a somewhat unusual form which we wish to comment. One sees two branches. The right branch is concave and goes from ǫ = U , the ground state energy where the complexity vanishes, to ǫ = ǫ M , the maximal energy beyond which one does not find any local ground state, which corresponds to a value µ = µ M . It is obtained for µ ∈ [µ M , µ * ], where µ * is the point where Φ(µ) reaches its maximum, corresponding to the GGS energy U . The second branch is convex, and interpolates between the RS solution (obtained at µ = 0) and the maximal complexity point (obtained at µ = µ M ). This second branch does not seem to have a direct physical interpretation and in this context can be simply ignored. On the other hand it must be present insofar as the µ → 0 limit of our RSB solution gives back the RS solution. Clearly a better understanding of this second branch would be welcome.
It is curious that this type of complexity curve has not been discussed in the literature so far, to our knowledge. In the rest of this section we want to present a short comment on some apparent discrepancy between our results and some widely used results on the SK model [30]. In the SK model which corresponds to the k → ∞ limit of our problem, the states can be defined at any temperature as solutions of TAP equations [37,2]. Bray and Moore [38] have computed the number of such solutions, and they have deduced that in the zero temperature limit, the corresponding complexity becomes equal to the complexity of one spin flip stable (1SF) states, which has been computed in [38,39]. In fig.7 we show the complexity of 1SF states, and also the complexity at the 1RSB level computed through the use of the Legendre transform of the free energy as in (43). This 1RSB result is indeed very different, qualitatively and quantitatively, from the complexity of 1SF states, and also different from the general shape found by Bray and Moore at finite temperatures.
It seems to us that this difference can be explained by two arguments: -The 1SF configurations are generically not stable to several spin flips (in contradiction to what is claimed in [40]) and in particular a change on the cavity fields may propagate through the whole system -The free energy functions of [38] admit also other new saddle points which have not yet been considered. Clearly a full discussion of these points goes much beyond the scope of the present paper; much more work is needed in order to understand the relationship between TAP states, the q-spin flip stable states, and the Legendre transform construction in spin glasses [41].

V. CONCLUSION
With respect to our previous work [9] on the cavity method for the Bethe lattice spin glass, the fact of working directly at zero temperature allows to simplify the discussion of the cavity method in several aspects 1 . On the one hand, one can discuss the physics directly in terms of ground states. The physics of RSB comes into play when there exist many LGS and their energies may cross when one adds a new spin. At a more technical level, the cavity fields are integers, and it is thus much simpler to parametrize their probability distributions. This allows in some cases to get analytic solutions, but also when one has to resort to the numerical solution of the cavity equations using the population dynamics, this aspects simplifies a lot the procedure: one can follow a population of the local distributions P i (h), because each of them is easily parametrized, while in the finite T approach one must represent each P i (h) by a population of M fields. Computationally, this zero T approach is thus much more efficient. This is very interesting for discussing optimization problems, and in fact it has already allowed for a full solution of the Ksat problem [4]. We expect that it can have many other applications in this context. If we study the thermodynamics of the model at finite temperature, we can use a very similar approach where the free energy of valleys is used instead of the energy of the configuration. A similar cavity method can be developed at all temperatures, as was done in [9]. In general one would expect that the finite free energy density computed in this way be continuous at T = 0, which means that one should get back the results of the present paper by using the T → 0 limit of the finite T computations.
Strictly speaking this does not happens for the model with J = ±1. The reason is the following. Strictly at T = 0 if the local fields have integer values when we start the iteration, they are then always integers. On the contrary, starting from the same initial condition at T = 0, after a few iterations the local fields are integer plus corrections of order T . However it is possible (and this is what happens in the RS and in the 1RSB cases) that the coefficient of the corrections proportional to T increases exponentially with the iteration number so that, after a number of iterations which is proportional to | ln T |, one obtains a distribution which is no more concentrated near the integers.
In our opinion it is quite likely that this instability is an artefact of not using the full RSB scheme and it should disappear when more and more precise computations, with higher levels of RSB, will be done. In order to substantiate this point we have computed the distribution of the cavity fields in the J = ±1 case for k = 2 at very small T (i.e. T = 0.01) both in the RS case and in the 1RSB case. The results are shown in fig. 8. It is evident that in the 1RSB case the distribution is much more concentrated on the integers.

APPENDIX B: VARIATIONAL FORMULATION
It is interesting to notice that all the equations of stability under iteration of the cavity procedure, giving the probability distributions of local fields, can be obtained from a variational principle.
In the replica symmetric case, let us define, for a generic probability distribution R(h), the functional where Φ (1) and Φ (2) are given by: and where dP (h) = P (h)dh. One can show, following a method similar to the one given in [9], that the stationarity equation of this functional coincides with the equation (12) giving the field distribution P (h). Furthermore, on this solution, one can use the stationarity equation in order to derive: where ∆E (1) , ∆E (1) are given in eq. (13,14). This shows that, on the solution, The variational formulation also exists in the 1RSB case. In this case, we need to introduce a functional Φ[Q[P (h)]] (a somewhat complicated object). We introduce first the two expressions for Φ (1) and Φ (2) and then construct as before: We must now compute the functional derivative with respect to Q. We first notice that Q must be a normalized distribution so that This constraint will be enforced using a Lagrange multiplier λ. Let us consider the functional derivative of Φ rsb [Q] with respect to Q[P ], evaluated on an arbitrary function P (h) = P * (h). We find for the first piece: and for the second one: Finally, the equation: should be true, with an appropriate choice of the Lagrange multiplier λ, for any P * .
In the first piece we can use: Let us then define the probability distribution P 0 (h 0 ) through: where C[{P }] is a functional of the P i with i ∈ {1, ..., k}. As Q[P ] satisfies the 1RSB self-consistency equation, we know that the distribution P 0 (h 0 ) defined in (72) is generated with a probability Q[P 0 ]. Therefore we get, on the 1RSB self-consistent-under-iteration solution for Q[P ]: Clearly this is equal to −(k + 1)/2 times the derivative of Φ (2) rsb , for any P * , up to some constant that compensates with the Lagrange multiplier λ. The curious reader may find interesting to know that Notice that the expression (46) for the internal energy can be directly derived using the relation d(µΦ(µ))/dµ = ǫ. Indeed where the derivative d dµ * act only on the explicit dependence on µ. The second term in the previous equation is zero when Q satisfies the equilibrium condition. If we compute the derivative of eq. (42) we get the result contained in equation (46). It is interesting to notice that all the equations of stability under iteration of the cavity procedure, giving the probability distributions of local fields, can be obtained from a variational principle. Moreover the fact that the free energy we use is a variational one is crucial in simplifying the computation of all the derivatives with respect to the various parameters, because only explicit derivatives must be considered.

APPENDIX C: MISCELLANEOUS COMMENTS ON THE DEFINITION OF THE STATES
In this appendix we give some further comments on the definition of the LGS and of the complexity. The subject is rather delicate from the mathematical point of view and if one is not precise it is easy to make wrong statements and to find contradictions. Our aim here is just to point out some of the subtleties which one should take into account.
In the definition of LGS we have asked that the configuration should be stable with respect to k spin flips with the number k going to infinity with N in some unprecised way. We could start by defining a LGS of order k as a configuration that is stable with respect to k spin flips. In this way one could define the complexities Σ k (N ) as (1/N ) log(Number of k-spin-flips-stable configurations). The complexity of LGS could be defined as An alternative natural definition is to define LGS as stable to changing an infinitesimal fraction of all spins: The question of the equivalence among various possible definitions of Σ, and in particular whether Σ a = Σ b , is a well posed mathematical problem and may have a non trivial answer.
One case in which one encounters problems when studying the large k, large N limit is the SK model. Although it is quite possible that these problems are an artefact of the infinite number of neighbours in the SK model and they may be not present on the Bethe lattice, let us mention them here briefly. In the SK model one can compute [38,39] the complexity of one spin stable configurations, lim N →∞ Σ 1 (N ). A detailed computation shows that (at least at not too low energies) the probability distribution of the local fields h in these states P (h) goes to a constant when h goes to zero. It is easy to deduce that with probability one these states are not 2-spins stable and therefore they do not correspond to LGS, with any reasonable definition. This argument might seem to imply that lim N →∞ Σ 2 (N ) < lim N →∞ Σ 1 (N ). However this is not true: we can count the limit for large N of the complexity,Σ(δ), for the "δ-stable configurations" such that Irrespectively of the precise value of δ, a simple computation shows thatΣ(δ) is continuous at δ = 0, whileΣ(0) = lim N →∞ Σ 1 (N ). On the other hand δ-stable configurations are also stable for k spin flips (at least when k ≪ δ √ N ). This argument shows that for the SK model lim N →∞ Σ k (N ) is k independent. These configurations are likely not to be LGS if we ask for configurations that are stable with respect to k spin flips with k ∝ N 1/2 but they should be k-spin-flips-stable if k ∝ N γ with γ < 1/2. It appears rather likely that there will be a change of behaviour in the complexity when k ∼ √ N . Further and detailed investigations are needed to clarify this point.
After this digression on the SK model, let us come back to the computation of the LGS complexity that we have done in this paper. As we have already stressed, we do not use directly any of the previous definitions based on Σ k (N ). We have just argued that a LGS must correspond to a local solution of the cavity equations and we have counted the number of these solutions in the infinite N limit. This procedure does not give a clue on the finite N effects, and therefore on the definition of LGS on one given (finite N ) sample.
One possibility could be to write directly these equations on a finite lattice (these are nothing but TAP equations [37]) and to associate to each solution of these equation a LGS. However it is well known that in general global solutions of TAP equations for one given sample are rather rare [42]: the local analysis of the solution of the equations near to a given site of the lattice that we perform does not guarantee to provide a global solution to TAP equations. Probably the best one can hope is the following. We could define a quasi-solution of order α of TAP equations as a configuration of spins such that a fraction 1 − α of all TAP equations are satisfied and a fraction α of TAP equations are not satisfied. Our N → ∞ analysis considers on the same footing the true solutions of all TAP equations or the quasi-solutions where a vanishing fraction (for N → ∞) of the equations is not satisfied. So one could conjecture that the LGS correspond to quasi-solutions of order α, in the limit where α → 0. Notice that the α → 0 and N → ∞ limits may not commute and therefore one should be rather careful before making any precise statement on the number of TAP solutions.
In order to illustrate the previous point let us consider a simple example, a one dimensional chain, which is a particular example of a Bethe lattice with k = 1. Let us assume that the system is antiferromagnetic (i.e J(i, i + 1) = −1) and that periodic boundary conditions are used.
The TAP equations are given in terms of the local cavity fields: where we have used the obvious notation h ← (i) = h(i|i + 1) , h → (i) = h(i|i − 1). If N is even we have two degenerate GGS, i.e. σ(i) = (−1) i and σ(i) = (−1) i+1 On the contrary if N is odd we have N degenerate GGS with one defect in k, i.e. σ(i) = (−1) i for i < k ; σ(i) = (−1) i+1 for i ≥ k If we look at the TAP equations we see that for N even there are two solutions (corresponding to each of the two GGS) On the contrary there are no solutions of the TAP equations for odd N as an effect of the periodic boundary conditions. Obviously there are N quasi solutions of order 1, that luckily correspond to the N GGS.