{'fraction_non_alphanumeric': 0.0524144068825092, 'fraction_numerical': 0.017314149327795684, 'mean_word_length': 3.9618627057406663, 'pattern_counts': {'":': 0, '<': 46, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Effective field theory for models defined over small-world networks. First and second order phase transitions

We present an effective field theory method to analyze, in a very general way, models defined over small-world networks. Even if the exactness of the method is limited to the paramagnetic regions and to some special limits, it provides, yielding a clear and immediate (also in terms of calculation) physical insight, the exact critical behavior and the exact critical surfaces and percolation thresholds. The underlying structure of the non random part of the model, i.e., the set of spins filling up a given lattice L_0 of dimension d_0 and interacting through a fixed coupling J_0, is exactly taken into account. When J_0\geq 0, the small-world effect gives rise, as is known, to a second-order phase transition that takes place independently of the dimension d_0 and of the added random connectivity c. When J_0<0, a different and novel scenario emerges in which, besides a spin glass transition, multiple first- and second-order phase transitions may take place. As immediate analytical applications we analyze the Viana-Bray model (d_0=0), the one dimensional chain (d_0=1), and the spherical model for arbitrary d_0.


I. INTRODUCTION
Since the very beginning of the pioneer work by Watts and Strogatz [1], the interest toward small-world networks -an interplay between random and regular networks -has been growing "exponentially". Mainly, there are two reasons that have caused such a "diffusion".
The first reason is due to the topological properties of the small-world network. In synthesis, if N is the size of the system, for any finite probability p of rewiring, or for any finite added random connectivity c (the two cases correspond to slightly different procedures for building a small-world network) one has: a "short-distance behavior", implying that the average shortest distance between two arbitrarily chosen sites grows as l(N ) ∼ log(N ), as in random networks, and a large clustering coefficient, C(N ) ∼ O(1), as in regular lattices. The interplay between these two features makes small-world networks representative of many realistic situations including social networks, communications networks, chemical reactions networks, protein networks, neuronal networks, etc.
The second reason is due the fact that, in models defined over small-world networks, despite the presence of an underlying finite dimensional structure -a lattice L 0 of dimension d 0 < ∞ -the existence of short-cut bonds makes such models mean-field like and -hopefully -exactly solvable. However, even if such a claim sounds intuitively correct, the complexity of these models turns out to be in general quite high and, compared to numerical works, there are still few exact results for small-world networks [2]- [9] (for the percolation and synchronization problem see [10] and [11]).
In particular, for d 0 > 1 a mean-field critical behavior is expected and has been also supported by Monte Carlo (MC) simulations [12]. Some natural questions then arise. Are we able to prove analytically such a be-havior? If for example d 0 = 2, does the mean-field critical behavior hold for any situation? Yet, does the correlation length diverge at the critical temperature?
Furthermore, even if for d 0 = 1 an exact analytical treatment has been developed at the level of replica symmetry breaking (RSB) [5] and one step replica symmetry breaking (1RSB) [6], the calculations are quite involved and the solutions of the coupled equations to evaluate the order parameters require a certain numerical work, which becomes rapidly hard in the 1RSB case. In any case, even if these methods are able to give in principle exact results at any temperature, they are not in general suitable to provide a clear simple and immediate physical picture of the model, even within some approximations. The main problem in fact resides in the presence of short loops: as soon as d 0 > 1 these loops cannot be neglected and the "traditional" cavity and replica methods seem hardly applicable. In particular, we are not able to predict what happens, for example, if we set J 0 negative. Should we still expect a second-order phase transition? And what about the phase diagram?
In this paper, we present a general method to study random Ising models defined on small-world graphs built up by adding a random connectivity c over an underlying arbitrary lattice L 0 having dimension d 0 . We will show that this method -in a very simple and physically sound way -provides an answer to the above questions as well as to many others.
Roughly speaking, as an effective field theory the method generalizes the Curie-Weiss mean-field equation m = tanh(βJm) to take into account the presence of the short-range couplings J 0 besides the long-range ones J. As we will show, the magnetization m of the model defined over the small-world network, shortly the random model, behaves as the magnetization m 0 of the model defined over L 0 , shortly the unperturbed model, but im-mersed in an effective external field to be determined self-consistently. Even if the exactness of this method is limited to the paramagnetic regions (P), it provides the exact critical behavior and the exact critical surfaces, as well as simple qualitative good estimates of the correlation functions in the ferromagnetic (F) and spin glass regions (SG). Furthermore, in unfrustrated systems, the method becomes exact at any temperature in the two limits c → 0 + and c → ∞.
The consequences of such a general result are remarkable from both the theoretical and the practical point of view. Once the explicit form of the magnetization of the unperturbed model, m 0 = m 0 (βJ 0 , βh), is known, analytically or numerically, as a function of the couplings J 0 and of the external field h, we get an approximation to the full solution of the random model, which is analytical or numerical, respectively, and becomes exact in the P region. If we do not have m 0 = m 0 (βJ 0 , βh) but we know at least some of its properties, we can still use these properties to derive certain exact relations and the critical behavior of the random model.
In the first part of the paper, after presenting the selfconsistent equations, we focus on their application for a general study of the critical surfaces and of the critical behavior. In the second part, we apply the method to study models of interest which can be solved analytically (and very easily) as for them we know m 0 (βJ 0 , βh): the Viana-Bray model (which can be seen as a d 0 = 0 dimensional small-world model), the one-dimensional chain small-world model, and the spherical small-world model in arbitrary d 0 dimension.
We stress that the critical surfaces as well as the correlation functions in the P region provided by the present method are exact and not based on any special ansatz as the replica-symmetry and the tree-like ansatz. We prove in particular that: independently of the added random connectivity c, of the underlying dimension d 0 , of the structure of the underlying lattice L 0 , and of the nature of the phase transition present in the unperturbed model (if any), for J 0 ≥ 0 we always have a second-order phase transition with the classical mean-field critical indices but with a finite correlation length if calculated along the "Euclidean distance" defined in L 0 ; on the other hand, for J 0 < 0 we show that, as soon as c is sufficiently large, there exist at least two critical temperatures which, depending on the behavior of χ 0 (βJ 0 , βh) -the susceptibility of the unperturbed system -correspond to first-or second-order phase transitions. This phenomenon will be explicitly shown in the example of the one-dimensional small-world model. Note that, as it will result from the detailed analysis of the self-consistent equation (Sec. IIIB), in any case, the critical behavior of the unperturbed model -if any -can never influence the behavior of the random model.
The paper is organized as follows. In Sec. II we introduce the class of small-world networks over which we define the random Ising models, stressing some important differences concerning the definition of the correla-tion functions with respect to those usually considered in "ordinary" random models. In Sec. III we present our method: in Sec. IIIA we provide the self-consistent equations and their relations with physical correlation functions, in Sec. IIIB we analyze the stability of the solutions of the self consistent equations and the critical surface and behavior of the system. We separate the Sec. IIIB in the sub-cases J 0 ≥ 0 and J 0 < 0. In Sec. IIIC we discuss the limits of the method. In Sec. IIID we study the stability between the F and the SG phases and the phase diagram. Finally, in Sec. IIIE we mention how to generalize the method to cases with more different shortrange couplings J 0 , and to analyze possible disordered antiferromagnetic systems. In Secs. IV, V and VI the theory is applied to the three above mentioned example cases. The successive Secs. VII, VIII and IX are devoted to the derivation of the method. The starting point of the proof is given in Sec. VII and is based on a general mapping between a random model and a non random one [13]- [15] suitably adapted to the present case. The selfconsistent equations are then easily derived in Sec. VIII. Note that, apart from the equations concerning the stability between the P-F and the P-SG transitions, which are derived in Sec. IX, the derivations of the equations presented in Sec. IIIB are mostly left to the reader, since they can be easily obtained by standard arguments of statistical mechanics using the Landau free energy ψ(m) that we provide and that is derived in Sec. VIII too. Finally, in Sec. X we draw some conclusions. In Appendix A we generalize the method to inhomogeneous external fields to make clear the subtle behavior of the correlation functions in small-world models.

II. RANDOM ISING MODELS ON SMALL-WORLD NETWORKS
The family of models we shall consider are random Ising models constructed by super-imposing random graphs with finite average connectivity c onto some given lattice L 0 whose set of bonds (i, j) and dimension will be indicated by Γ 0 and d 0 , respectively. Given an Ising model (the unperturbed model ) of N spins coupled over L 0 through a coupling J 0 with Hamiltonian and given an ensemble C of unconstrained random graphs c, c ∈ C, whose bonds are determined by the adjacency matrix elements c i,j = 0, 1, we define the corresponding small-world model as described by the following Hamiltonian the free energy F and the averages O l , with l = 1, 2, being defined in the usual (quenched) way as (β = 1/T ) and where Z c,J is the partition function of the quenched system O c,J the Boltzmann-average of the quenched system ( O depends on the given realization of the J's and of c: O = O c;J ; for shortness we later will omit to write these dependencies) and dP (J) and P (c) are two product measures given in terms of two normalized measures dµ(J i,j ) ≥ 0 and p(c i,j ) ≥ 0, respectively: The variables c i,j ∈ {0, 1} specify whether a "longrange" bond between the sites i and j is present (c i,j = 1) or absent (c i,j = 0), whereas the J i,j 's are the random variables of the given bond (i, j). For the J i,j 's we will not assume any particular distribution, while, to be specific, for the c i,j 's we shall consider the following distribution This choice leads in the thermodynamic limit N → ∞ to a number of long range connections per site distributed according to a Poisson law with mean c > 0 (so that in average there are in total cN/2 bonds). Note however that the main results we report in the next section are easily generalizable to any case in which Eq. (8) holds, or holds only in the thermodynamic limit due a sufficiently small number of constrains among the matrix elements c i,j . When we will need to be specific, for the J i,j 's we will assume either the distribution or (11) to consider ferromagnetism or glassy phases, respectively. In Eq. (11) p ∈ [0, 1]. The quantities of major interest are the averages, and the quadratic averages, of the correlation functions which for shortness will be indicated by C (1) and C (2) . For example, the following are non connected correlation functions of order k: where k ≥ 1 and the indices i 1 , . . . , i k are supposed all different. For shortness we will keep on to use the symbols C (1) and C (2) also for the connected correlation function since, as we shall see in the next section, they obey to the same rules of transformations. We point out that the set of indices i 1 , . . . , i k is fixed along the process of the two averages. This implies in particular that, if we consider the spin with index i and the spin with index j, their distance remains undefined, or more precisely, the only meaningful distance between i and j, is the distance defined over L 0 , i.e., the Euclidean distance between i and j, which we will indicate as ||i − j|| 0 . Therefore, throughout this paper, it must be kept in mind that, for example, C (1) (||i − j|| 0 ) = σ i σ j is very different from the correlation function G (1) (l) of two points at a fixed distance l, l being here the distance defined over both L 0 and the random graph c, i.e., the minimum number of bonds to join two points among both the bonds of Γ 0 and the bonds of the random graph c. In fact, if, for J 0 = 0, one considers all the possible realizations of the Poisson graph, and then all the possible distances l between two given points i and j, one has where here P N (l) is the probability that, in the system with N spins, the shortest path between the vertices i and j has length l. If we now use G (1) (l) ∼ (tanh(βJ)) l [16] (in the P region holds the equality) and the fact that the average of l with respect to P N (l) is of the order log(N ), we see that the two point connected correlation function (14) goes to 0 in the thermodynamic limit. Similarly, all the connected correlation functions defined in this way are zero in this limit. Note however, that this independence of the variables holds only if J 0 = 0. This discussion will be more deeply analyzed along the proof by using another point of view, based on mapping the random model to a suitable fully connected model.

A. The self-consistent equations
Depending on the temperature T, and on the parameters of the probability distributions, dµ(·) and p(·), the random model may stably stay either in the paramagnetic (P), in the ferromagnetic (F), or in the spin-glass (SG) phase. In our approach for the F and SG phases there are two natural order parameters that will be indicated by m (F) and m (SG) . Similarly, for any correlation function, quadratic or not, there are two natural quantities indicated by C (F) and C (SG) , and that in turn will be calculated in terms of m (F) and m (SG) , respectively. To avoid confusion, it should be kept in mind that in our approach, for any observable O, there are -in principle -always two solutions that we label as F and SG, but, as we shall discuss in Sec. IIID, for any temperature, only one of the two solutions is stable and useful in the thermodynamic limit.
In the following, we will use the label 0 to specify that we are referring to the unperturbed model with Hamiltonian (1). Note that all the equations presented in this paper have meaning and usefulness also for sufficiently large but finite size N . For shortness we shall omit to write the dependence on N .
Let m 0 (βJ 0 , βh) be the stable magnetization of the unperturbed model with coupling J 0 and in the presence of a uniform external field h at inverse temperature β. Then, the order parameters m (Σ) , Σ=F,SG, satisfy the following self-consistent decoupled equations where the effective couplings J (F) , J (SG) , J and βJ (SG) 0 = tanh −1 (tanh 2 (βJ 0 )).
Note that |J . For the correlation functions C (Σ) , Σ=F,SG, for sufficiently large N we have where C 0 (βJ 0 , βh) is the correlation function of the unperturbed (non random) model.
Concerning the free energy density f we have where l = 1, 2 for Σ=F,SG, respectively, and the non trivial free energy term L (Σ) is given by f 0 (βJ 0 , βh) being the free energy density in the thermodynamic limit of the unperturbed model with coupling J 0 and in the presence of an external field h, at inverse temperature β. For given β, among all the possible solutions of Eqs. (15), in the thermodynamic limit, for both Σ=F and Σ=SG, the true solutionm (Σ) , or leading solution, is the one that minimizes L (Σ) : Finally, let k be the order of a given correlation function C (1) or C (2) . The averages and the quadratic averages over the disorder, C (1) and C (2) , are related to C (F) and C (SG) , as follows and From Eqs. (27) and (28) for k = 1, we note that the Edward-Anderson order parameter [17] C (2) = σ 2 = q EA is equal to (C (SG) ) 2 = (m (SG) ) 2 only in the SG phase, whereas in the F phase we have q EA = (m (F) ) 2 . Therefore, since m (SG) = m (F) , m (SG) is not equal to √ q EA ; in our approach m (SG) represents a sort of a spin glass order parameter [18].
The localization and the reciprocal stability between the F and SG phases will be discussed in Sec. IIID. Note however that, at least for lattices L 0 having only loops of even length, the stable P region is always that corresponding to a P-F phase diagram, so that in the P region the correlation functions must be calculated only through Eqs. (24) and (27).
As an immediate consequence of Eq. (15) we get the susceptibilityχ (Σ) of the random model: (29) whereχ 0 stands for the susceptibility χ 0 of the unperturbed model divided by β (we will adopt throughout this dimensionless definition of the susceptibility): and similarly for the random model. For the case Σ =F without disorder (dµ(J ′ ) = δ(J ′ − J)dJ ′ ), Eq. (29) was already derived in [7] by series expansion techniques at zero field (h = 0) in the P region (where m = 0).
Another remarkable consequence of our theory comes from Eq. (20). We see in fact that in the thermodynamic limit any correlation function of the random model fits with the correlation function of the unperturbed model but immersed in an effective field that is exactly zero in the P region and zero external field (h = 0). In other words, in terms of correlation functions, in the P region, the random model and the unperturbed model are indistinguishable (modulo the transformation J 0 → J (SG) 0 for Σ =SG). Note however that this assertion holds only for a given correlation function calculated in the thermodynamic limit. In fact, the corrective O(1/N ) term appearing in the rhs of Eq. (20) cannot be neglected when we sum the correlation functions over all the sites i ∈ L 0 , as to calculate the susceptibility; yet it is just this corrective O(1/N ) term that gives rise to the singularities in the random model.
More precisely, for the two point connected correlation functionχ where l = 1, 2 for Σ = F, SG, respectively, if we havẽ where the dependence on N inχ i,j andχ 0 are understood. Eq. (33) clarifies the structure of the correlation functions in small-world models. In the rhs we have two terms: the former is a distance-dependent short-range term whose finite correlation length, for T = T being the critical temperature of the unperturbed model with coupling J (Σ) 0 ), makes it normalizable, the latter is instead a distance-independent long-range term which turns out to be normalizable thanks to the 1/N factor. Once summed, both the terms give a finite contribution to the susceptibility. It is immediate to verify that by summingχ (Σ) i,j over all the indices i, j ∈ L 0 and dividing by N we get back -as it must be -Eq. (29). Eq. (33) will be derived in Appendix A where we generalize the theory to a non homogeneous external field.
For what follows, we need to rewrite the non trivial part of the free energy density L (Σ) (m) as where the introduced term ψ (Σ) plays the role of a Landau free energy density and is responsible for the critical behavior of the system. Around m = 0, up to terms O(h 2 ) and O(m 3 h), ψ (Σ) (m) can be expanded as follows where finally, the last term ∆ (βf 0 ) βJ 0 , βJ (Σ) m is defined implicitly to render Eqs. (35) and (36) exact, but terms O(h 2 ) and O(m 3 h), explicitly: We recall that the k−2-th derivative ofχ 0 βJ (Σ) 0 , βh with respect to the second argument, calculated at h = 0, gives the total sum of all the k-th cumulants normalized to N : stands for the cumulant, or connected correlation function, of order k of the unperturbed model, σ i1 σ i2 . Note that, apart from the sign, these terms are proportional to the Binder cumulants [19] (which are all zero above T c0 for k > 2) only for N finite. In the thermodynamic limit the terms b (Σ) , c (Σ) , . . ., in general are non zero and take into account the large deviations of the block-spin distribution functions from the Gaussian distribution.
Let T be the critical temperatures, if any, of the random model and let t (Σ) be the corresponding reduced temperatures: Here, the term "critical temperature", stands for any temperature where some singularity shows up. However, if we limit ourselves to consider only the critical temperatures crossing which the system passes from a P region to a non P region, from Eq. (34) it is easy to see that, independently on the sign of J 0 and on the nature of the phase transition, we have the important inequalities where we have introduced β (Σ) c0 , the inverse critical temperature of the unperturbed model with coupling J (Σ) 0 and zero external field. If more than one critical temperature is present in the unperturbed model, β c0 is the value corresponding to the smallest value of these critical temperatures (highest in terms of β). Formally we set β (Σ) c0 = ∞ if no phase transition is present in the unperturbed model. A consequence of Eq. (40) is that, in studying the critical behavior of the system for h = 0, we can puth (Σ) = 0. Throughout this paper, we shall reserve the name critical temperature of the unperturbed model as a P-F critical temperature through which the magnetization m 0 (βJ 0 , 0) passes from a zero to a non zero value, continuously or not. This implies, in particular, that for J 0 < 0 we have -formally -β c0 = ∞.
In this paper we shall study only the order parameters m (F) and m (SG) , whereas we will give only few remarks on how to generalize the method for possible antiferromagnetic order parameters. We point out however that the existence of possible antiferromagnetic transitions of the unperturbed model does not affect the results we present in this paper.
It is convenient to distinguish the cases J 0 ≥ 0 and J 0 < 0, since they give rise to two strictly different scenarios.
1. The case J0 ≥ 0 In this case βJ (Σ)χ 0 βJ 0 , βJ (Σ) m (Σ) + βh is an increasing function of β for β < β can be determined by developing -for h = 0 -Eqs. (15) for small m (F) and m (SG) , respectively, which, in terms ofχ 0 gives the following exact equatioñ where the constrain β c0 excludes other possible spurious solutions that may appear when d 0 ≥ 2 (since in this case β (Σ) c0 may be finite). The critical behavior of the system can be derived by developing Eqs. (15) for small fields. Alternatively, one can study the critical behavior by analyzing the Landau free energy density ψ (Σ) (m (Σ) ) given by Eq. (36).
In the following we will suppose that for J 0 > 0, b (Σ) be positive. We have checked this hypothesis in all the models we have until now considered and that will be analyzed in Secs. IV, V and VI. Furthermore, even if the sign of c (Σ) cannot be in general a priori established, for the convexity of the function f 0 with respect to βh, the sum of the six-th term with ∆ (βf 0 ) βJ in Eq. (36) must go necessarily to +∞ for m (Σ) → ∞. In conclusion, when J 0 ≥ 0, for the critical behavior of the system, the only relevant parameters of so that the critical behavior can be immediately derived as in the Landau theory for the so called m 4 model [20]. On noting that it is convenient to define so that we have Note that, due to the fact that J 0 ≥ 0, A (Σ) > 0, and, as already mentioned, b (Σ) ≥ 0 as well. By using Eq. (47) for β < β c0 and near β (Σ) c , we see that the minimumm (Σ) of ψ (Σ) , i.e., the solution of Eq. (15) near the critical temperature, is given bȳ Similarly, we can write general formulas for the susceptibility and the equation of state. We havẽ Finally, on using Eqs. (36) and (48) we get that the specific heat C (Σ) has the following finite jump discontinuity at β where C (Σ) c is the continuous part of the specific heat corresponding to the part of the free energy density without ψ (Σ) .
Hence, as a very general result, independently of the structure of the underlying graph L 0 and its dimension d 0 , independently of the nature of the phase transition present in unperturbed model (if any), and independently of the added random connectivity c, provided positive, we recover that the random model has always a meanfield critical behavior with a second-order phase transition with the classical exponents β = 1/2, γ = γ ′ = 1, δ = 3 and α = α ′ = 0, and certain constant coefficients depending on the susceptibilityχ 0 and its derivatives calculated at β = β (Σ) c and external field h = 0. Note however, that the correlation length of the system calculated along the distance of L 0 , || · || 0 , remains finite also at β (Σ) c . In fact, from Eq. (20), for the two point correlation function at distance r If we now assume for C 0 (βJ 0 , 0; r) the following general Ornstein-Zernike form f 0 (r) = f 0 (βJ 0 ; r) being a smooth function of r (which has not to be confused with the free energy density), and ξ 0 = ξ 0 (βJ 0 ) the correlation length, which is supposed to diverge only at β c0 (if any), on comparing Eqs. (52) and (53) we have (notice that, as explained in Sec. IIIA, at least for lattices L 0 having only loops of even length, the physical correlation function is only that corresponding to Σ =F, i.e., C (1) = C (F) ) where and Therefore, due to the inequalities (43), we see that The knowledge of C 0 (βJ 0 , βh; r) also for h = 0 would allows us to find the general expression for c . However, since C 0 (βJ 0 , βh; r) has no critical behavior for h = 0, it follows that C (Σ) (r) cannot have a critical behavior for β > β (Σ) c either (and then also for β → β (Σ) c from the right). This result is consistent with [2].

The case J0 < 0
In this case J is no longer a monotonic function of β. However, it is easy to see that that for β = 0 and β → ∞, this function goes to 0. Therefore, for a sufficiently large connectivity c, from Eq. (34) we see that there may appear at least two regions where the paramagnetic solution m (F) = 0 is stable, separated by a third region in which a non zero solution is instead stable. However the situation is even more complicated since, unlike the case J 0 ≥ 0, the non monotonicity of βJ (F)χ 0 βJ 0 , βJ (F) m (F) reflects also in the fact that the self-consistent Eq. (15) for Σ =F may have more solutions of the kind ±m (F) , ±m ′ (F ), . . . which are still stable with respect to the stability condition (34), for h = 0. We face in fact here the problem to compare more stable solutions. According to Eq. (23), in the thermodynamic limit, among all the possible stable solutions, onlym (F) , the solution that minimizes L (F) , survives, whereas the not leading ones play the role of metastable states. This kind of scenario, which includes also finite jump discontinuities, has been besides observed in the context of small-world neural networks in [21] where we even observe some analogy in the used formalism, at least for the simplest case of one binary pattern. From Eqs. (38) and (39) we see that the signs of the Landau coefficients a (Σ) , b (Σ) , c (Σ) , . . ., are functions of β and J 0 only. Given J 0 < 0, the most important quantity that features the non monotonicity of The equation for β (F) * , as a function of J 0 , defines a point where b (F) = 0. If J 0 < 0, the most general equation for a generic critical temperature is no longer given by Eq. (44). In fact, in general, a critical temperature now is any temperature where the stable and leading solutionm (F) may have a singular behavior, also with finite jumps between two non zero values.
There are some simplification when for the Landau coefficient c (F) , we have c (F) > 0, or at least c (F) > 0 out of the P region In this situation in fact, from Eq. (36) we see that a (F) , b (F) and c (F) are the only relevant terms for the critical behavior of the system and -for small values ofm (F) -we can again apply the Landau theory, this time for the so called m 6 model [20]. In such a case, for the solutionm (F) we havē whereas with a (F) ≥ 0 establishes a line of first-order transitions through whichm (F) changes discontinuously from zero to The point a (F) = b (F) = 0 is a tricritical point where the second and first-order transition lines meet. If we approach the tricritical point along the line b (F) = 0 we get the critical indices α = 1/2, α ′ = 0, β = 1/4, γ = γ ′ = 1 and δ = 5. However, this critical behavior along the line b (F) = 0 has not a great practical interest since from Eq. (38) we see that it is not possible to keep b (F) constant and zero as the temperature varies. Finally, we point out that, even if c (F) > 0, when the transition is of the first-order, Eqs. (59) and (61) hold only for b (F) , and then a (F) , sufficiently small, since only in such a case the finite discontinuity ofm (F) is small and then the truncation of the Landau free energy term ψ (F) to a finite order meaningful. Note that this question implies also that we cannot establish a simple and general rule to determine the critical temperature of a first-order phase transition (we will return soon on this point).
When c (F) < 0, the Landau theory of the m 6 model cannot be of course applied. However, as in the case J 0 > 0, even if the sign of c (F) cannot be a priori established, for the convexity of the function f 0 with respect to βh, the sum of the six-th term with ∆ (βf 0 ) βJ (36) must go necessarily to +∞ for m (F) → ∞ and a qualitative similar behavior of the m 6 model is expected. In general, when J 0 < 0, the exact results are limited to the following ones.
From now on, if not otherwise explicitly said, we shall reserve the name critical temperature, whose inverse value of β we still indicate with β (F) c , to any temperature on the boundary of a P region (through whichm (F) passes from 0 to a non zero value, continuously or not). For each critical temperature, depending on the value of β (F) * , we have three possible scenario of phase transitions: Note that, according to our definition of critical temperature, the critical behavior described by Eqs. (59) -(62) represents a particular case of the general scenario expressed by Eqs. (63). We see also that, in general, when b (F) ≤ 0, approaching the tricritical point, for the critical exponent β we have β ≤ 1/4.
In the case in which β (F) c corresponds to a second-order phase transition, or in the case in which a (F) < 0 out of the P region (at least immediately near the critical temperature), β (F) c can be exactly calculated by Eq. (44). When we are not in such cases, the only exact way to determine the critical temperature is to find the full solution form (F) which consists in looking numerically for all the possible solutions of Eq. (15) and -among those satisfying the stability condition (34) -selecting the one that gives the minimum value of L (F) .

C. Level of accuracy of the method
In the P region, Eqs. (15)(16)(17)(18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28)(29)(30)(31)(32)(33) are exact, whereas in the other regions provide an effective approximation whose level of accuracy depends on the details of the model. In particular, in the absence of frustration the method becomes exact at any temperature in two important limits: in the limit c → 0 + , in the case of second-order phase transitions, due to a simple continuity argument; and in the limit c → ∞, due to the fact that in this case the system becomes a suitable fully connected model exactly described by the self-consistent equations (15) (of course, when c → ∞, to have a finite critical temperature one has to renormalize the average of the coupling by c).
However, for any c > 0, off of the P region and infinitely near the critical temperature, Eqs. (15)(16)(17)(18)(19)(20) are able to give the exact critical behavior in the sense of the critical indices and, in the limit of low temperatures, Eqs. (15)(16)(17)(18)(19) provide the exact percolation threshold. In general, as for the SK model, which can be seen as a particular model with J 0 = 0, the level of accuracy is better for the F phase rather than for the SG one and this is particularly true for the free energy density f (Σ) , Eq. (21). In fact, though the derivatives of f (Σ) are expected to give a good qualitative and partly also a quantitative description of the system, f (SG) itself can give wrong results when the SG phase at low temperatures is considered. We warn the reader that in a model with J 0 = 0, and a symmetrical distribution dµ(J i,j ) with varianceJ, the method gives a ground state energy per site u (SG) , which grows with c as u (SG) ∼ −Jc, whereas the correct result is expected to be u (SG) ∼ −J √ c [22]. As a consequence, in the SK model, in the limit β → ∞, the method gives a completely wrong result with an infinite energy. We stress however that the order parameters m (F) and m (SG) , and then also the correlation functions, by construction, are exact in the zero temperature limit.

D. Phase diagram
The physical inverse critical temperature β c of the random model is in general a non single-value function of X: β c = β c (X), where X represents symbolically the parameters of the probability dµ for the couplings J i,j , and the parameter c, the average connectivity (which is also a parameter of the probability distribution of the short-cut bonds). The parameters of dµ can be expressed through the moments of dµ, and as they vary the probability dµ changes. For example, if dµ is a Gaussian distribution, as in the SK model, there are only two parameters given by the first and second moment. A concrete example for the one dimensional small-world model will be shown in Figs. (17) and (18).
In the thermodynamic limit, only one of the two solutions with label F or SG survives, and it is the solution having minimum free energy. In principle, were our method exact at all temperatures, we were able to derive exactly all the phase diagram. However, in our method, the solution with label F or SG are exact only in their own P region, i.e., the region where m (F) = 0 or m (SG) = 0, respectively. Unfortunately, according to what we have seen in Sec. IIIC, whereas the solution with label F is still a good approximation also out of the P region, in frustrated model (where the variance of dµ is large if compared to its first moment) the free energy of the solution with label SG becomes completely wrong at low temperatures. Therefore, we are not able to give in general the exact boundary between the solution with label F and the solution with label SG, and in particular we are not able to give the physical frontier F/SG. However, within some limitations which we now prescribe, we are able to give the exact critical surface, i.e., the boundary with the P phase, establishing which one -in the thermodynamic limit -of the two critical boundaries, P-F or P-SG, is stable (we will use here the more common expression "stable" instead of the expression "leading"), and to localize some regions of the phase diagram for which we can say exactly whether the stable solution is P, F, or SG. We will prove the stability of these solutions in Sec. IX. When for a region we are not able to discriminate between the solution with label F and the solution with label SG and they are both out of their own P region, we will indicate such a region with the symbol "SG and/or F" (stressing in this way that in this region there may be also mixed phases and re-entrance phenomena).
In Sec. IX we prove that there are four possible kind of phase diagrams that may occur according to the cases The four kind of possible phase diagrams are schematically depicted in Figs. 1-4 in the plane (T, X).
As we have seen in Sec. IIIB1, if J 0 ≥ 0, for both the solution with label F and SG, we have one -and only onecritical temperature. In the following, to avoid confusion, it should be kept in mind the distinction between the physical β c = β c (X) and β (Σ) c , with Σ =F or SG. β c satisfies the following rules.
Case (1): If d 0 < 2 and J 0 is a finite range coupling, or else d 0 = ∞ at least in a broad sense (see [15]), β c (X) is a single-value function of X, and we have A schematic representation of this case is given in Fig. 1.
Notice in particular that the second line of Eq. (65) does not exclude that β c (X) might be a non single-value function of X. A schematic representation of this case is given in Fig. 2.

J0 < 0
As we have seen in Sec. IIIB2, if J 0 < 0, for a sufficiently large connectivity c, the solution with label F has at least two separated P regions corresponding to two critical temperatures. Here we assume that the underlying lattice L 0 has only loops of even length so that, for example, triangular lattices are here excluded. Let us suppose to have for the solution with label F only two critical temperatures (the minimum number, if J 0 < 0), and let be In general we have the following scenario.
Case (3): If d 0 < 2 and J 0 is a finite range coupling, or d 0 = ∞ in a broad sense (see [15]), β c2 (X) is a singlevalue function of X and satisfies Eq. (64). The other critical inverse temperature β c1 (X) is instead: either a two-value function of X and we have or where ∄ in Eq. (68) means that if β there is no stable boundary with the P region. A schematic representation of this case is given in Fig. 3.
Case (4): If 2 ≤ d 0 < ∞, β c2 satisfies Eq. (65); whereas for β c1 we have either where in Eq. (69) we have introduced the symbol SG> to indicate that in general the stable P-SG surface is above (or below in terms of temperatures) the surface coming from the solution with label SG: β . Notice that, similarly to the case (3), we cannot exclude that β c1 in Eq. (70) be a non single-value function of X, as well as β (SG>) c in Eq. (69). A schematic representation of this case is given in Fig. 4.
If more than two critical temperatures are present, the above scheme generalizes straightforwardly.
Keeping our definition for the introduced symbol "SG and/or F", we stress that: in all the fours cases the phases F and "SG and/or F" are exactly localized; in the cases (1) and (3) the phases P and SG are exactly localized; in the cases (2) and (4) the SG phase is always limited below (in terms of temperatures) by the unstable P-SG surface coming from the solution with label SG (indicated as P-SG unst in Figs. 2 and 4). Finally, we stress that -under the hypothesis that L 0 has only loops of even length -the stable P regions correspond always to the solution with label F. For 2 ≤ d 0 < ∞, from the second line of Eqs. (65) and (69) and from Eq. (70), we see that the method is not able to give the complete information about the P-SG boundary since we have only inequalities, not equalities. Furthermore, in these regions of the phase diagram the physical critical temperature in general may be a non single-value function of X. On the other hand, we have the important information that in these equations the inequalities between the physical β c and β (SG) c (the solution with label SG) are always strict. As a consequence, we see that, when 2 ≤ d 0 < ∞, in these regions the SG "magnetization" m (SG) will always have a finite jump discontinuity in crossing the surface given by β c . In other words, along such a branch of the critical surface corresponding to the second line of Eqs. (65) and (69) and Eq. (70), we have a first-order phase transition, independently of the fact that the phase transition corresponding to the β (SG) c surface is second-order, and independently of the sign of J 0 .

E. Generalizations
The generalization to the cases in which the unperturbed model has an Hamiltonian H 0 involving couplings depending on the bond b ∈ Γ 0 is straightforward. In this case we have just to substitute everywhere in the formulae (15) 0b }. However, the critical behavior will be in general different and more complicated than that depicted in the Subsections IIIB1 and IIIB2. In particular, even in the case in which all the couplings J (Σ) 0b are positive, we cannot assume that the Landau coefficient b (Σ) be positive so that, even in such a case, first-order phase transitions are in principle possible, as has been seen via Monte Carlo simulations in directed small-world models [23].
As anticipated, our method can be generalized also to study possible antiferromagnetic phase transitions in the random model. There can be two kind of sources of antiferromagnetism: one due to a negative coupling J 0 in the unperturbed model, the other due to random shortcuts J i,j having a measure dµ with a negative average.
In the first case, if for example the sublattice L 0 is bipartite into two sublattices L   (4): J0 < 0 and 2 ≤ d0 < ∞. T and X as in Fig. 1 introduce the site-dependent solution m 0i to find correspondingly in a set of coupled equations (at most N ), the effective fields m i of the random model.
In the second case, following [24] we consider a lattice L 0 which is composed of, say, p sublattices L (ν) 0 , ν = 1, . . . , p. Then, we build up the random model with the rule that any shortcut may connect only sites belonging to two different sublattices. Hence, as already done in [13] for the generalized SK model, we introduce p effective fields m (ν) which satisfy a system of p self-consistent equations involving the p fields m (ν) 0 and calculated in the p external fields J (F) m (ν) (note that here the symbol F stresses only the fact that the effective coupling must be calculated through Eq. (16)). As an immediate example, let us consider the Viana-Bray model [25]. It can be seen as the simplest smallworld model in which N spins with no short-range couplings (here J 0 = 0) are randomly connected by longrange connections J (possibly also random). Note that formally here L 0 has dimension d 0 = 0. Since J 0 = 0, for the unperturbed model we have It is interesting to check that the first and second derivatives ofχ 0 in h = 0 are null and negative, respectively. In fact we have and Applying these results to Eqs. (15)(16)(17)(18)(19) we get immediately the self-consistent equations for the F and the SG magnetizations and the Viana-Bray critical surface On choosing for dµ a measure having average and variance scaling as O(1/c), for c ∝ N , we recover the equations for the Sherrington-Kirkpatrick model [26] already derived in this form in [13] and [14]. In these papers Eqs.
(74-77) were derived by mapping the Viana-Bray model and, similarly, the Sherrington-Kirkpatrick model to the non random fully connected Ising model. In this sense it should be also clear that, at least for β ≤ β c and zero external field, in the thermodynamic limit, the connected correlation functions (of order k greater than 1) in the Sherrington-Kirkpatrick and in the Viana-Bray model are exactly zero. In fact, in the thermodynamic limit, the non random fully connected model can be exactly reduced to a model of non interacting spins immersed in an effective medium so that among any two spins there is no correlation. Such a result is due to the fact that, in these models, all the N spins interact through the same coupling J/N , no matter how far apart they are, and the net effect of this is that in the thermodynamic limit the system becomes equivalent to a collection of N non interacting spins seeing only an effective external field (the medium) like in Eq. (71) with βh replaced by βJm.
For the measure (10) our approximated Eq. (74) can be compared with the exact known equation that can be derived by using the Bethe-Peierls or the replica approach and is given by (see for example [27] and references therein) where the effective field H is determined by the following integral equation In the limit β → ∞, Eqs. (74) and (75) give the following size (normalized to 1) of the giant connected component These equations are not exact, however they succeed in giving the exact percolation threshold c = 1. In fact, concerning the equation (78) for the F phase, the exact equation for m (F) is (see for example [27] and references therein) which, in terms of the function tanh, becomes so that Eqs. (78) and (80) are equivalent at the order O(m (F) ). We see also that, as stated in the Sec. IIIC, Eqs. (78) and (80) become equal in the limits c → 0 and c → ∞.

B. Gas of Dimers
Let us consider for L 0 a set of 2N spins coupled through a coupling J 0 two by two. The expression "gas of dimers" stresses the fact that the dimers, i.e. the couples of spins, do not interact each other. As a consequence, the free energy, the magnetization, and the susceptibility of the unperturbed model can be immediately calculated. We have − βf 0 (βJ 0 , βh) = 1 2 log 2e βJ0 cosh(2βh) + 2e −βJ0 , m 0 (βJ 0 , βh) = e βJ0 sinh(2βh) e βJ0 cosh(2βh) + e −βJ0 , Let us calculate also the second derivative ofχ 0 . From we get We note that, as expected, the second derivative ofχ 0 in h = 0, for J 0 ≥ 0 is always negative, whereas, for J 0 < 0 it becomes positive as soon as β|J 0 | > log( √ 2). By using the above equations, from Sec. III we get immediately the following self-consistent equation for the magnetizations and -at least for J 0 ≥ 0 -the equation for the critical temperature As it will be clear soon, this model lies between the Viana-Bray model and the more complex d 0 = 1 dimensional chain small-world model, which will be analyzed in detail in the next section. Our major interest in this simpler gas of dimers small-world model is related to the fact that, in spite of its simplicity and d 0 =0 dimensionality -since the second derivative ofχ 0 may be positive when J 0 is negative -according to the general result of Sec. IIIB, it is able to give rise to also multiple first-and second-order phase transitions.

V. SMALL WORLD IN d0=1 DIMENSION
In this section we will analyze the case in which L 0 is the d 0 =1-dimensional chain with periodic boundary conditions (p.b.c.). The corresponding small-world model with Hamiltonian (2) in zero field has already been analyzed in [5] by using the replica method. Here we will recover the results found in [5] for β c and will provide the self-consistent equations for the magnetizations m (F) and m (SG) whose solution, as expected, turns out to be in good agreement with the corresponding solutions found in [5] for c small and large (the latter when the frustration is relatively small). It will be however rather evident how much the two methods differ in terms of simplicity and intuitive meaning. Furthermore, we will derive also an explicit expression for the two-points connected correlation function which, to the best of our knowledge, had not been published yet. Finally, we will analyze in the detail the completely novel scenario for the case J 0 < 0 which, as mentioned, produces multiple first-and secondorder phase transitions.
In order to apply the method of Sec.III we have to solve the one dimensional Ising model with p.b.c. immersed in an external field. The solution of this non random model is easy and well known (see for example [28]). If we indicate with λ 1 and λ 2 the two eigenvalues coming from the transfer matrix method, one has λ 1,2 = e βJ0 cosh(βh) ± e 2βJ0 sinh 2 (βh) + e −2βJ0 1 2 , from which it follows that, for the free energy density, the magnetization and the two-points connected correlation function, we have − βf 0 (βJ 0 , βh) = log (λ 1 ) , m 0 (βJ 0 , βh) = e βJ0 sinh(βh) e 2βJ0 sinh 2 (βh) + e −2βJ0 1 2 , where the phase ϕ is defined by and ||i − j|| 0 is the (euclidean) distance between i and j. Let us calculateχ 0 and its first and second derivatives. From Eq. (83) we havẽ From Eq. (88) we see that for J 0 > 0 and any β we have, for sufficiently small h, ∂ 2 ∂(βh) 2χ 0 (βJ 0 , βh) < cosh(βh) whereas for J 0 < 0 we have (92) For J 0 ≥ 0 the above equation gives the exact P-F and P-SG critical temperatures in agreement with [5]. When J 0 < 0 -unless the transition be of second-order -Eq. (92) for Σ =F does not signal a phase transition. In general, as J 0 < 0 the P-F critical temperature must be determined by looking at all the stable solutions m (F) of the self-consistent equation (91) 23)). In all these examples we have chosen the measure (10). Figs. 5 and 6 concern two cases with J 0 > 0 so that one and only one second-order phase transition is present. The input data of these two cases are the same as those analyzed numerically in [5] (note that in the model considered in [5], the long range coupling J is divided by c). As already stated in Sec. III, the self-consistent equations become exact in the limit c → 0, for second-order phase transitions, and in the limit c → ∞. Therefore, for the magnetization, by comparison with [5], in Figs. 5 and 6, where c is relatively small and big, respectively, we see good agreement  also below the critical temperature. Figs. 7-14 concern eight cases with J 0 < 0. In these figures we plot also the line y = 1, to make evident when the stability conditions for the solutions m (F) = 0 and m (F) = 0, which are given byχ 0 β (F) J , are satisfied, respectively. As explained above, the critical behavior and the localization of the critical temperatures is more complicated when J 0 < 0. In particular, given |J 0 |, if c is not sufficiently high the solution m (F) = 0 remains stable at all temperatures and if it is also a leading solution, no phase transition occurs. Let us consider Eq. (92). For J 0 < 0 the lhs of this equation has some maximum at a finite valueβ given bȳ where r def = |J 0 |/J, and we have introduced Hence, we see that a sufficient condition for the solution m (F) = 0 to become unstable is that be Note that the above represents only a condition for the instability of the solution m (F) = 0, but the true solution is the one that is both stable and leading. In fact, when J 0 < 0, a phase transition in general may be present also when Eq. (93) is not satisfied and, correspondingly the possible critical temperatures will be not determined by Eq. (92). In Fig. 7 we report a case with J = 1, J 0 = −1 (r = 1) and a relatively low value of c, c = 5, so that no phase transition is present. Similarly, in Fig. 8 we report again a case in which no phase transition is present due to the  fact that here r is relatively big, r = 1.1. It is interesting to observe that for r = 1 Eq. (93) requires a value of c greater than the limit value c = 3 + 2 √ 2 = 5.8284.... In both Figs. 9 and 10 we report a case in which Eq. (93) is still not satisfied, but nevertheless two first-order phase transitions are present. In both Figs. 11 and 12 we have one first-and one second-order phase transition. In both Figs. 13 and 14 we have two second-order phase transitions. As anticipated in Sec. IIIB, we note that in Figs. 9-12, i.e., the cases in which there is at least one firstorder phase transition, there are always regions where both the solutions, m (F) = 0 and m (F) = 0, are simultaneously stable but only one solution is leading (the one drawn); whereas in Figs. 13 and 14, as it was in Figs. 5 and 6, since we have only second-order phase transitions, the stability condition turns out to be a necessary and sufficient condition for determining the leading solution and the critical temperature can be determined also by Eq. (92) with Σ =F. Σ =SG. Note that, unlike the P-F critical surface, the P-SG critical surface does not depend on the parameter p entering in Eq. (11). For the reciprocal stability between the P-F and the P-SG critical surface we remind the reader to the general rules of Sec. IIID (see cases (1) and (3)) which, for J 0 ≥ 0, reduce to the results reported in Sec. 6.1 of the Ref. [5]. Here we stress just that, if J 0 ≥ 0, for p ≤ 0.5, only the P-SG transition is possible. However, when J 0 < 0 and c is not sufficiently large, the SG phase may be the only stable phase even when p = 1. In fact, although when J 0 < 0 the solution m (F) may have two P-F critical temperatures, in general, if the P-SG temperature is between these, we cannot exclude that the solution m (SG) starts to be the leading solution at sufficiently low temperatures. In Figs. 17 and 18  and c = 0.5, respectively, the P-F transition is always second-order. We plan to investigate in more detail the phase diagram in future works.

VI. SMALL-WORLD SPHERICAL MODEL IN ARBITRARY DIMENSION d0
In this section we will analyze the case in which the unperturbed model is the spherical model built up over a d 0 -dimensional lattice L 0 (see [28] and references therein) [31]. In this case the σ's are continuous "spin" variables ranging in the interval (−∞, ∞) subjected to the sole constrain i∈L0 σ 2 i = N , however our theorems and formalism can be applied as well and give results that, within the same limitations prescribed in Sec. III, are exact.
Following [28], for the unperturbed model we have − βf 0 (βJ 0 , βh) = 1 2 log π βJ 0 + φ (βJ 0 , βh,z) , where andz =z(βJ 0 , βh) is the (unique) solution of the equation ∂ z φ (βJ 0 , βh, z) = 0: from which it follows the equation for m 0 The derivative g ′ can in turn be expressed as J 0 (it) being the usual Bessel function whose behavior for large t is given by The critical behavior of the unperturbed system depends on the values of g ′ (z) and g ′′ (z) near z = 0. It turns out that for d 0 ≤ 2 one has g ′ (0) = ∞ and there is no spontaneous magnetization, whereas for d 0 > 2 one has g ′ (0) < ∞ and at h = 0 the unperturbed system undergoes a second-order phase transition with magnetization given by Eq. (96) which, for β above β c0 , becomes where the inverse critical temperature β c0 is given by Furthermore, it turns out that for d 0 ≤ 4 one has g ′′ (0) = ∞, whereas for d 0 > 4 one has g ′′ (0) < ∞. This reflects on the critical exponents α, γ and δ, which take the classical mean-field values only for d 0 > 4. According to Sec. III, to solve the random modelfor simplicity -at zero external field, we have to perform the effective substitutions βJ 0 → βJ

A. Random Models defined on Quenched Graphs
Let us consider the following random model. Given a graph g, which can be determined through the adjacency matrix for shortness also indicate by g = {g b }, with g b = 0, 1, b being a bond, let us indicate with Γ g the set of the bonds b of g and let us define over Γ g the Hamiltonian where J b is the random coupling at the bond b, and σ i b , σ i b are the Ising variables at the end-points of b. The free energy F and the physics are defined as in Sec. II by Eqs. (3)- (6): where dP ({J b }) is a product measure over all the possible bonds b given in terms of normalized measures dµ b ≥ 0 (we are considering a general measure dµ b allowing also for a possible dependence on the bonds) where Γ full stands for the set of bonds of the fully connected graph. As in Sec. II, we will indicate a generic correlation function, connected or not, by C with understood indices i 1 , . . . , i k all different, see Eqs. (12) and (13). In the following, given an arbitrary vertex i of g, we will consider as first neighbors j of i only those vertices for which dµ i,j (J i,j )J i,j or dµ i,j (J i,j )J 2 i,j are at least O(1/N ). Note that we can always neglect couplings having lower averages. We will indicate with D(Γ g ) the average number of first neighbors of the graph g. For a d-dimensional lattice, D(Γ g ) = 2d − 1, for a Bethe lattice of coordination number q, D(Γ g ) = q − 1, and for long range models, D(Γ g ) ∝ N . We will exploit in particular the fact that D( Given a random model defined trough Eqs. (101-104), we define, on the same set of bonds Γ g , its related Ising model trough the following Ising Hamiltonian where the Ising couplings J have non random values In the following a suffix I over quantities such as H I , F I , f I , g I , etc. . . , or J c , etc. . . , will be referred to the related Ising system with Hamiltonian (105).
We can always split the free energy of the random system with N spins as follows φ being the high temperature part of the free energy. Let ϕ be the density of φ in the thermodynamic limit Let us indicate with ϕ I the high temperature part of the free energy density of the related Ising model defined through Eqs. (105-107). As is known, ϕ I can be expressed in terms of the quantities z b = tanh(βJ ) and z i = tanh(βh i ), i.e., the parameters of the high temperature expansion: The related Ising model is defined by a set of, typically few, independent couplings {J Asymptotically, at sufficiently high dimensions D(Γ g ), the critical inverse temperature of the spin glass model β c is given by and in the paramagnetic phase for D(Γ g ) > 2 the following mapping holds and and Σ=F/AF or SG, for l=1 or 2, respectively. In the limit D(Γ g ) → ∞ and h i = 0, i = 1, . . . , N , Eqs. (111-117), give the exact free energy and correlation functions in the paramagnetic phase (P); the exact critical paramagnetic-spin glass (P-SG), β (SG) c , and paramagnetic-F/AF (P-F/AF), β (F/AF) c , surfaces, whose reciprocal stability depends on which of the two ones has higher temperature. In the case of a measure dµ not depending on the bond b, the suffix F and AF stand for ferromagnetic and antiferromagnetic, respectively. In the general case, such a distinction is possible only in the positive and negative sectors in the space of the parameters of the probability distribution, { dµ b J b ≥ 0} and { dµ b J b < 0}, respectively, whereas, for the other sectors, we use the symbol F/AF only to stress that the transition is not P-SG.
It is not difficult to see that, when the measure dµ does not depend on the specific bond b, i.e., if dµ b ≡ dµ b ′ ∀b, b ′ ∈ Γ g , in the P region Eqs. (111-117) lead to the following exact limit for ϕ and C [15] lim therefore, the basic role of Eqs. (114-117), is to show how, in the limit D(Γ g ) → ∞, ϕ and C approach zero and which are their singularities. In particular this proves that for all the (random) infinite dimensional models and any disorder non bond-dependent, the critical exponent α ′ for the specific heat has the mean-field classical value, α ′ = 0, and that the correlation functions (with different indices) above the critical temperature are exactly zero. We point out however that, when the measure dµ b depends explicitly on the bond b, Eq. (118) in general does not hold [32]. In fact, when the measure dµ b is bond-dependent, the symmetry expressed by Eq. (118) is broken since the bonds are no longer equivalent. As we will see in the next section, in small-world models with an underlying lattice L 0 having d 0 < 2, even if Eq. (118) may still holds for ϕ, the symmetry is broken for C since the direction(s) of the axis(es) of L 0 is(are) now favored direction(s). Yet, if 2 ≤ d 0 < ∞, the symmetry (118) for ϕ is broken as well.
The analytic continuation of Eqs. (114-117) to β > β c and/or for h = 0 provide certain estimations which are expected to be qualitatively good. In general such estimations are not exact, and this is particularly evident for the free energy density of the SG phase. However, the analytic continuation for the other quantities gives a good qualitative result and provide the exact critical behavior (in the sense of the critical indices) and the exact percolation threshold.
For models defined over graphs satisfying a weak definition of infinite dimensionality, as happens on a Bethe lattice with coordination number q > 2, a more general mapping has been established [15]. In this case, all the above equations -along the critical surface (at least)still hold exactly in the thermodynamic limit, where we can set effectively D(Γ g ) = ∞. However, for the aims of this paper we do not need here to consider this generalization of the mapping.
We have yet to make an important comment about Eqs. (25), (26) and (28), concerning the evaluation of a correlation function in the SG phase here for a random system with J 0 = 0 (for the moment being). In fact Eq. (116), for both a normal and a quadratic correlation function C (1) or C (2) , has a factor 1/2 not entering in the physical Eqs. (25), (26) and (28). The difference is just due to an artefact of the mapping that separates the Gibbs state into two pure states [18] not only in the F case, but also in the SG case. In fact, let us consider the correlation functions of order k = 1, that is, C (1) = σ 1 and C (2) = q EA = σ 1 2 . We see that, for C (1) , Eq. (116) in the SG phase gives C (1) = m (SG) /2. On the other hand, for any non zero solution m (SG) of the self-consistent Eq. (15), there exists another solution −m (SG) , and both the solutions have 1/2 probability to be realized in the random model. Since the SG phase is expected to be the phase characterized by having q EA = 0 and σ 1 = 0, we see that if we introduce both the solutions m (SG) and −m (SG) , we get σ 1 = 0 in the SG phase. Similarly, for C (2) , Eq. (116) in the SG phase gives C (2) = (m (SG) ) 2 /2, which at zero temperature gives 1/2, whereas a completely frozen state with q EA = 1 is expected. Again, we recover the expected physical q EA by using both the solutions m (SG) and −m (SG) . Repeating a similar argument for any correlation function of order k, and recalling that for k even (odd) the correlation function is an even (odd) function of the external magnetic field h, we arrive at Eqs. (25), (26) and (28).

B. Random Models defined on Unconstrained Random Graphs
Let us consider now more general random models in which the source of the randomness comes from both the randomness of the couplings and the randomness of the graph. Given an ensemble of graphs g ∈ G distributed with some distribution P (g), let us define The free energy F and the physics are now given by and similarly for O l , l = 1, 2. Here Z g ({J b }) is the partition function of the quenched system onto the graph realization g with bonds in Γ g For unconstrained random graphs, or for random graphs having a number of constrains that grows sufficiently slowly with N , the probability P (g), for large N , factorizes as In such a case we can exploit the mapping we have previously seen for models over quenched graphs as follows. Let us define the effective couplingJ b : correspondingly, since the random variables J b and g b are independent we have with the sum rule As a consequence, if we define the following global measure we see that the mapping of the previous section can be applied as we had a single effective graph Γ p given by where Z p is the partition function of the model with Hamiltonian H p given by VIII. DERIVATION OF THE SELF-CONSISTENT EQUATIONS By using the above results, we are now able to derive easily Eqs. (15)(16)(17)(18)(19)(20)(21)(22)(23). Sometimes to indicate a bond b we will use the symbol (i, j), or more shortly ij.
It is convenient to look formally at the coupling J 0 also as a random coupling with distribution Let us rewrite explicitly the Hamiltonian (2) as follows and let us introduce the random variables Taking into account that the random variable J 0 + c ij J ij , up to terms O(1/N ), is distributed according to dµ 0 (J 0 ), the independent random variables J ′ b and g ′ b have distributions dµ ′ b and p ′ b respectively given by where the measures dµ and p are those of the model introduced in Sec. II. As a consequence, Eq. (122) can be cast in the form of Eq. (120) with the measure Finally, since p b (g b ) = 0 for any b ∈ Γ full , we have also and due to the fact that D(Γ full ) ∝ N , in the thermodynamic limit the mapping becomes exact.
According to Eqs. (105-107), the related Ising model of our small-world model has the following Hamiltonian with two free couplings: J (I) 0 , for Γ 0 , and J (I) , for Γ full After solving this Ising (I) model the mapping allows us to come back to the random model by performing simultaneously for any b ∈ Γ full the reverse substitutions tanh βJ where l = 1, 2 for Σ = F or SG solution, respectively. Since the couplings J (I) 0 and J (I) are arbitrary, we find it convenient to renormalize J (I) as J (I) /N and at the end of the calculation to put again J (I) instead of J (I) /N . Note that for the mapping nothing changes if we do not make this substitution; the choice to use J (I) /N instead of J (I) is merely due to a formal convenience, since in this way the calculations are presented in a more standard and physically understandable form. In fact, according to Eqs. (123) and (126) what matters after solving the related Ising model with J (I) /N instead of J (I) is that, once for Σ=F and once for Σ=SG, we perform -simultaneously in the two couplings -the following reverse mapping transformations (l = 1, 2 for Σ = F or SG, respectively): for (i, j) / ∈ Γ 0 , and tanh βJ for (i, j) ∈ Γ 0 . Explicitly, by applying Eqs. (123), (9) and (121), the transformations (127) and (128) become, respectively and where we have made use of the definitions (16)- (19) introduced in Sec. III. Let us now solve the related Ising model. We have to evaluate the following partition function In the following we will suppose that J (I) (and then J (Σ) ) is positive. The derivation for J (I) (and then J (F) ) negative differs from the other derivation just for a rotation of π/2 in the complex m-plane, and leads to the same result one can obtain by analytically continue the equations derived for J (I) > 0 to the region J (I) < 0. By using the Gaussian transformation we can rewrite Z I as where c N is a normalization constant and, in the exponent of Eq. (131), we have again neglected terms of order O(1). For finite N we can exchange the integral and the sum over the σ's. By using the definition of the unperturbed model with Hamiltonian H 0 , Eq. (1), whose free energy density, for given βJ 0 and βh, is indicated with f 0 (βJ 0 , βh), we arrive at so that, if the stability condition is satisfied, in the thermodynamic limit we arrive at the following expression for the free energy density f I of the related Ising model Similarly, any correlation function C I of the related Ising model is given in terms of the correlation function C 0 of the unperturbed model by the following relation Of course, the saddle point solution m sp represents the magnetization of the related Ising model, as can be checked directly by deriving Eq. (135) with respect to βh and by using Eq. (134).
If the saddle point equation (134) has more stable solutions, the "true" free energy and the "true" observable of the related Ising model will be given by Eqs. (135)  As stressed in Sec. IIIB, for the unperturbed model we use the expression "critical temperature" for any temperature where the magnetization m 0 at zero external field passes from 0 to a non zero value, continuously or not. Note that, as a consequence, if J (I) 0 < 0, we have formally β (I) c0 = ∞, independently from the fact that some antiferromagnetic order may be not zero.
Let us start to make the obvious observation that a necessary condition for the related Ising model to have a phase transition at h = 0 and for a finite temperature, is the existence of some paramagnetic region P I where m I = 0. We see from the saddle point equation (134) that, for h = 0, a necessary condition for m I = 0 to be a solution is that be β ≤ β (I) c0 for any β in P I , from which we get also β . In a few lines we will see however that the inequality must be strict if β (I) c0 is finite, which in particular excludes the case J 0 < 0 (for which the inequality to be proved is trivial).
Let us suppose for the moment that be β The critical behavior of the related Ising model can be studied by expanding Eq. (134) for small fields. However, we find it more convenient to expand L(m) in series around m = 0 since in this way everything can be cast in the standard formalism of the Landau theory of phase transitions. From Eq. (133), taking into account that the functionχ 0 (βJ 0 , βh) is an even function of βh, we have the following general expression valid for any m, β and small h L(m) = βf 0 βJ where we have introduced the Landau free energy density ψ(m) given by where a = 1 − βJ (I)χ 0 βJ where ϕ 0 (βJ 0 , βh) = −βf 0 βJ log cosh(βJ is the high temperature part of the free energy density of the unperturbed model with coupling J (Σ) 0 and external field h. There are some important properties for the function ϕ 0 (βJ 0 , 0): it is a monotonic increasing function of βJ 0 ; if the lattice L 0 has only loops of even length, ϕ 0 (βJ 0 , 0) is an even function of βJ 0 ; furthermore, if d 0 < 2 and the coupling-range is finite, or if d 0 = ∞ at least in a wide sense [15], in the thermodynamic limit we have ϕ 0 (βJ 0 , 0) = 0; if instead 2 ≤ d 0 < ∞, ϕ 0 (βJ 0 , 0) = 0. We see here therefore what anticipated in Sec. VIIA: when J 0 = 0, the symmetry among the random couplings is broken and for d 0 sufficiently high this reflects in a non zero ϕ (Σ) also in the P region.
Next we prove Eqs. (64)-(70). To this aim we have to calculate Eq. (150) at the leading solutionm (Σ) and to compare ϕ (F) and ϕ (SG) . Note that the term in the square parenthesis of Eq. (150) is non negative sincē m (Σ) is the absolute minimum of L (Σ) . We recall that for critical temperature we mean here any temperature lying on the boundary P-F or P-SG, so thatm (Σ) | β = 0 for any β in the P region.
A. J0 ≥ 0 If J 0 ≥ 0, for both the solution with label F and SG, we have only one second-order phase transition so that m (F) = 0 andm (SG) = 0, respectively, are the stable and leading solutions even on the boundary with the P region.
Let us suppose β turbed magnetization m 0 (βJ 0 , βh) is known analytically, i.e., the small-world models in dimension d 0 = 0, 1, ∞, corresponding to an ensemble of non interacting units (spins, dimers, etc...), the one dimensional chain, and the spherical model, respectively. In particular, we have studied in detail the small-world model defined over the one dimensional chain with positive and negative shortrange couplings showing explicitly how, in the second case, multicritical points with first-and second-order phase transition arise. Finally, the small-world spherical model -an exact solvable model (in our approach) with continuous spin variables -has provided us with an interesting case study to explore what happens as d 0 changes continuously from 0 to ∞. As expected from general grounds, unlike the non random version of the model, the small-world model presents always a finite temperature phase transition, even in the limit d 0 → 0 + . This latter result -besides to be consistent with what we have found in the d 0 = 0-dimensional discrete models -has a simple physical explanation in our approach.
In fact, it consists on mapping the small-world model (a random model) to a corresponding non random model (no long-range bonds) but immersed in an effective uniform external field which is active as soon as the added random connectivity c is not zero (see Eqs. (15)- (17)). Many interesting variants of the above models can be considered and are still analytically solvable by our approach (see also the generalizations considered in Sec. IIIE). However, our approach can be also applied numerically to study more complex small-world models for which the corresponding unperturbed model is not analytically available [29]. In fact, the numerical complexity in solving such small-world models is comparable to that required in solving a non random model immersed in a uniform external field.
Models defined on complex small-world networks [30] are an interesting subject of future work.  In this appendix we prove Eq. (33) calculating the O(1/N ) correction responsible for the divergence of the susceptibility of the random system at T c . To this aim we firstly need to generalize our method to an arbitrary external field. Let us consider again a fully connected model having -as done in Sec. VIII -long-range couplings J (for brevity we will here omit the label I) and short-range couplings J 0 but now immersed in an arbitrary (non homogeneous) external field {h n }, where n = 1, . . . , N . After using the Gaussian transformation we have the following partition function: hence, by using we see that if the stability condition 1 − βJχ 0 (βJ 0 , {βJm sp + βh n }) > 0, is satisfied, in the thermodynamic limit we arrive at the following expression for the free energy density f of the related Ising model immersed in an arbitrary external field βf = β 2 Jm 2 + βf 0 (βJ 0 , {βJm + βh n }) On the other hand, by derivation with respect to βh i and by using Eq. (A7), it is immediate to verify that m i def = σ i = m 0i (βJ 0 , {βJm sp + βh n }), and then also (from now on for brevity on we omit the symbol sp ) We want now to calculate the correlation functions. From Eq. (A11), by deriving with respect to βh j we haveχ where for brevity we have omitted the argument inχ 0;l,j andχ 0;i,k , which is the same ofχ 0 appearing in the denominator. If we now come back to choice a uniform external field h n = h, n = 1, . . . , N , we can use translational invariance and for the related Ising model (fully connected) we obtain the following correlation functioñ [χ 0 (βJ 0 , βJm + βh)] 2 1 − βJχ 0 (βJ 0 , βJm + βh) +χ 0;i,j (βJ 0 , βJm + βh). (A16) Finally, by performing the mapping substitutions (129) and (130) we arrive at Eq. (33). Similarly, any correlation function C of the related Ising model will be given by a similar formula with the leading term C 0 plus a correction O(1/N ) becoming important only near T c .