Mean-field theory of vector spin models on networks with arbitrary degree distributions

Understanding the relationship between the heterogeneous structure of complex networks and cooperative phenomena occurring on them remains a key problem in network science. Mean-field theories of spin models on networks constitute a fundamental tool to tackle this problem and a cornerstone of statistical physics, with an impressive number of applications in condensed matter, biology, and computer science. In this work we derive the mean-field equations for the equilibrium behavior of vector spin models on high-connectivity random networks with an arbitrary degree distribution and with randomly weighted links. We demonstrate that the high-connectivity limit of spin models on networks is not universal in that it depends on the full degree distribution. Such nonuniversal behavior is akin to a remarkable mechanism that leads to the breakdown of the central limit theorem when applied to the distribution of effective local fields. Traditional mean-field theories on fully-connected models, such as the Curie-Weiss, the Kuramoto, and the Sherrington-Kirkpatrick model, are only valid if the network degree distribution is highly concentrated around its mean degree. We obtain a series of results that highlight the importance of degree fluctuations to the phase diagram of mean-field spin models by focusing on the Kuramoto model of synchronization and on the Sherrington-Kirkpatrick model of spin-glasses. Numerical simulations corroborate our theoretical findings and provide compelling evidence that the present mean-field theory describes an intermediate regime of connectivity, in which the average degree $c$ scales as a power $c \propto N^{b}$ ($b<1$) of the total number $N \gg 1$ of spins. Our findings put forward a novel class of spin models that incorporate the effects of degree fluctuations and, at the same time, are amenable to exact analytic solutions.


I. INTRODUCTION
Several man-made and natural complex systems are represented by networks of nodes joined by links [1]. The study of the interplay between the structure of complex networks and dynamical processes on top of them has grown into a major research field [2,3], with many applications in physics, biology, information theory, and technology. As opposing to the homogeneous structures typically studied in solid state physics, such as Bravais and Bethe lattices [4], the striking feature of complex networks is the existence of strong local fluctuations in their structure. This heterogeneous character is responsible for most of the nontrivial dynamical properties of networked systems [2].
Spin models on networks describe systems formed by a large number of state variables, represented by scalars or continuous vectors, which are coupled through the links of networks [3,5,6]. The study of such models is of utmost importance for at least two main reasons. First, they are minimal models to address the impact of heterogeneous structures on the cooperative behavior of a large number of interacting degrees of freedom. Second, seemingly unrelated problems across disciplines can be cast in terms of the unifying framework of spin models on random networks [7,8]. Models of scalar spins on networks have a vast number of applications in a variety of research fields, such as opinion dynamics [9,10], models of socio-economic phenomena [11], artificial neural networks [12][13][14], agent-based models of the market behavior [15][16][17], dynamics of biological neural networks [18][19][20], * fmetzfmetz@gmail.com information theory and computer science [8], sparse randommatrix theory [21,22], and the stability of large dynamical systems [23,24]. Models of vector spins on networks are relevant for the study of synchronization phenomena [25][26][27][28], random lasers [29][30][31], vector spin-glasses [32][33][34], and the collective dynamics of swarms [35][36][37][38].
Mean-field theories stand out as one of the most celebrated tools in physics and the natural starting point to address the collective behavior of many interacting spins. The heart of the mean-field approach is the assumption that all spins are statistically equivalent, in the sense that each spin experiences an effective random field drawn from the same distribution. By virtue of that, the original problem of many interacting elements is replaced by a problem of a single spin coupled to an effective field. Paradigmatic examples of meanfield theories are derived from spin models on fully-connected networks, such as the Curie-Weiss model of ferromagnetism [39], the Kuramoto model of synchronization [25,27], and the Sherrington-Kirkpatrick model of spin-glasses [40,41].
Fully-connected mean-field theories are expected to provide an universal description of spin models on highconnectivity networks, for which the mean degree is infinitely large. In fact, it seems sensible to argue that spin models on networks gradually flow to their fully-connected behavior as the mean degree increases, since the detailed structure of a network should become irrelevant in the high-connectivity limit. Although this intuitive argument has been confirmed for a few network models [42][43][44], understanding the impact of degree fluctuations on the high-connectivity limit of spin models remains a key open question. In this work we provide a comprehensive solution to this problem. Surprisingly, we find that the high-connectivity limit of spin models on networks is not universal, since it depends on the full degree dis-arXiv:2110.11153v2 [cond-mat.dis-nn] 9 Feb 2022 tribution. It follows that traditional, fully-connected meanfield theories do not generally predict the macroscopic behavior of spin models on high-connectivity networks.
Such nonuniversal character is intimately related to the breakdown of the central limit theorem as applied to the distribution of effective local fields. To illustrate this point, let us consider the equilibrium behavior of Ising spins on random networks at inverse temperature β [7,44,45]. The local magnetization m i = tanh(βh i,eff ) at node i is determined by the effective field where ∂ i is the set of neighbors of i, J ij is the random coupling strength between spins at nodes i and j, and m (i) j is the local magnetization at node j in the absence of i. The network degrees are random variables with average c. Figure 1 depicts simulation results for two different degree distributions. For random graphs with a Poisson degree distribution, the central limit theorem holds and the distribution of effective fields converges to a Gaussian distribution when c → ∞, which corresponds to the fully-connected mean-field behavior [41]. For networks with an exponential degree distribution, the central limit theorem fails and the effective fields are no longer Gaussian in the high-connectivity limit. The breakdown of the central limit theorem is caused by the strong fluctuations of the random number of summands in Eq. (1), which is nothing more than the degree of node i. This compelling mechanism for the failure of the central limit theorem has been studied for more than seventy years in probability theory [46][47][48], but its evident importance for spin models on networks has so far eluded a careful analysis. In this paper we fill this gap and derive a novel family of mean-field theories that emerge from such breakdown of the central limit theorem.

A. Main results
The central result of this work is a set of equations for the equilibrium behavior of spin models on high-connectivity networks with an arbitrary degree distribution. The spins are Ising variables or continuous vectors with finite dimension, while the random pairwise interactions between spins follow an arbitrary distribution. The high-connectivity limit of the model is cast in terms of an effective problem of a single spin, whose configurations follow the Boltzmann distribution with an effective energy given by Eq. (43). The analytic expression for the distribution of the effective energy, Eq. (44), is one of the main outcomes of this work.
The remarkable consequence of Eq. (44) is that, even in the high-connectivity limit, the behavior of spin models on graphs is not universal, but strongly dependent on the degree distribution. Networks for which the high-connectivity limit depends on the degree distribution are called heterogeneous networks. In contrast, the behavior of spin models on the socalled homogeneous networks is universal, i.e., independent of the degree distribution and consistent with the behavior of fully-connected models. The analytic results of section IV are very general and they can be applied to a variety of specific models. We illustrate the effects of degree fluctuations on the mean-field behavior of spin models by focusing on two examples: the Kuramoto model of synchronization and the Ising model of spin-glasses. We obtain the complete phase diagrams of both models [see figures 4 and 6].
Degree fluctuations dramatically affect the distribution of effective fields and the magnetization inside the ferromagnetic phase. We show that the distribution of effective local fields of Ising spin models exhibits a long tail for large fields and a divergence at zero field [see Eq. (85)], in contrast to the Gaussian [41] or the Dirac-δ [39] distributions of effective fields in fully-connected models. In the case of the Kuramoto model, the magnetization or phase coherence displays a singular point that separates two distinct regions in the ferromagnetic phase [see figure 3], each one identified by a different behavior of the magnetization for strong couplings. In addition, we show that degree fluctuations have an ambiguous effect in the ordered phase. Although the ferromagnetic phase of the Kuramoto model expands over the entire phase diagram when the variance of the degree distribution diverges, the magnetization becomes arbitrarily small. This result sheds light on a key paradigm of network science, namely the idea that collective behavior is improved as the critical point vanishes with the increase of degree fluctuations [2,49]. Our results provide compelling evidence that in the thermodynamic limit of networks with diverging degree variance the only possible emerging phase is ultimately the paramagnetic (asynchronous) one. Last but not least, we find that the shape of the phase distribution of the Kuramoto model on heterogeneous networks fluctuates from site to site [see figure 5], which exposes a lack of correspondence between local and global ensemble averages, in sharp contrast to the behavior on fullyconnected networks.
The results of section IV B generalize the replicasymmetric mean-field theory and the phase diagram of Ising spin-glasses [40,41] to the case of heterogeneous networks. We explain how to derive the Almeida-Thouless line [50], which bounds the region of the phase diagram where the replica-symmetric theory is unstable and our results are no longer exact. When the variance of the degree distribution diverges, the spin-glass phase as well as the replica-symmetry breaking region are confined to an arbitrarily small sector of the phase diagram. Interestingly, the low-temperature portion of the Almeida-Thouless line exhibits a non-monotonic behavior as a function of the degree fluctuations.
Throughout the paper, we present simulation results for spin models on random networks with a finite number N of nodes. Besides confirming our theoretical findings, the simulations reveal that the heterogeneous mean-field theory derived in this work is valid in the limit N → ∞ when the mean degree scales as c ∝ N b (0 < b < 1). This intermediate regime of connectivity lies between the sparse (b = 0) and the fully-connected (b = 1) regimes. Our work uncovers the nonuniversal behavior of spin models on heterogeneous networks in such intermediate scaling of c.
The paper is organized as follows. In the next section we define the Hamiltonian of the vector spin model on an ensemble of random networks with an arbitrary degree distribution. Section III presents the cavity or message-passing equations for the equilibrium behavior of spin models on networks with finite mean degree, as derived in previous works [51,52]. This section also features the distributional version of the cavity equations. Section IV is the core of the paper. Initially, we thoroughly explain how to calculate the high-connectivity limit of the distributional cavity equations, from which the heterogeneous mean-field theory emerges. The explicit results for the effect of degree fluctuations are presented in two subsections. Subsection IV A is focused on the Kuramoto model with ferromagnetic couplings, while subsection IV B presents results for Ising spin models with random pairwise interactions. We further discuss our findings and main conclusions in section V. The paper contains two appendices. Appendix A explains in detail how to derive the Almeida-Thouless line for Ising spin models on heterogeneous networks, while appendix B presents some details of the numerical simulations.

II. THE MODEL SET-UP
We consider D-dimensional spins σ 1 , . . . , σ N , with σ i = (σ i,1 σ i,2 . . . σ i,D ) T , placed on the nodes of a simple random graph [53]. Each spin σ i is a continuous vector that identifies a point on the surface of the D-dimensional hypersphere R D with unit radius. The probability to observe a global spin configuration {σ} = (σ 1 , . . . , σ N ) in thermal equilibrium follows the Boltzmann distribution where β = 1/T is the inverse temperature, H(σ) is the Hamiltonian, and Z is the partition function The shorthand notation R D dσ i denotes an integral over all possible configurations of σ i such that σ 2 i = 1. We study a generic family of spin models invariant under orthogonal transformations and defined by the Hamiltonian [54] where {c ij } N i,j=1 are the elements of the graph adjacency matrix C, and {J ij } N i,j=1 are the coupling strengths between the spins. The set-up specified by Eq. (4) comprises a broad class of traditional models of spins randomly coupled through the links of a network. The Ising model, the XY or Kuramoto model with identical oscillators, and the classical Heisenberg model on a random network are obtained by setting, respectively, D = 1, D = 2, and D = 3.
The binary matrix elements c ij ∈ {0, 1}, with c ii = 0 ∀ i, specify the network topology: if c ij = c ji = 1, there is an undirected edge between nodes i and j, whereas c ij = 0 otherwise. The coordination number or degree k i = N j=1 c ij gives the number of nodes connected to i. The adjacency random matrix C is generated according to the configuration model of networks [55][56][57], in which a graph instance is uniformly drawn from the set of all possible simple graphs with a prescribed degree distribution [55][56][57] The average degree c, defined as is independent of the total number of nodes N in the graph. The configuration model of networks, in which p k is specified at the outset, is the ideal setting to explore how p k impacts the macroscopic behavior of the spin model of Eq. (4) in the highconnectivity limit c → ∞. We assume that p k is an arbitrary degree distribution with finite moments as long as c < ∞.
Although the main results presented below are valid for any p k , we will be particularly interested on random graphs with a negative binomial degree distribution [58] parametrized by c and α > 0. The parameter α is related to the variance σ 2 b of p (b) k as follows The distribution p k becomes identical to the exponential distribution for α = 1, and it converges to the Poisson distribution when α → ∞. In the limit c → ∞, the relative variance σ 2 b /c 2 is solely governed by α and, as we will show below, the mean-field theory depends on the degree distribution only through α, which enables to study the role of the degree fluctuations in a clear-cut way. In the context of complex networks, the negative binomial distribution finds applications in models of weighted random graphs [59], in studies of the spread of infectious diseases on networks [60,61], and as the empirical degree distribution of real-world contact networks [62].
The variable J ij ∈ R in Eq. (4) denotes the strength of the mutual interaction between σ i and σ j . The coupling strengths {J ij } N i,j=1 are independently and identically distributed random variables drawn from an arbitrary distribution p J with mean K 0 /c and standard deviation K 1 / √ c. Note that the ensemble of random graphs is fully specified by p k and p J .

III. THE DISTRIBUTIONAL VERSION OF THE MESSAGE-PASSING EQUATIONS
In this section we obtain a set of distributional equations for the single-site marginals of the spins when c is finite. These distributional equations, originally derived in references [51,52] and further studied in [42,43], build on the message-passing or cavity equations for network models with a local tree-like structure, such as the configuration model [63]. The mean-field equations for heterogeneous networks follow from the limit c → ∞ of the distributional equations. The probability density to observe a configuration σ i at an arbitrary node i follows from the Boltzmann distribution Thanks to the approximate local tree-like structure of the graph when N 1 [56,63], the local marginals {p i (σ i )} i=1,...,N on a single graph instance fulfill [51] where ∂ i is the set of nodes adjacent to node i, Z i is the normalization factor of p i (σ i ), and p (i) l (σ l ) is the distribution of σ l on the so-called cavity graph G (i) , obtained from the original graph G by deleting node i and all its edges. The distributions p (i) j (σ j ) solve the following self-consistency equations l (σ l ), (11) in which ∂ j \ i is the set of nodes adjacent to j except for i ∈ ∂ j , and Z (i) j is the normalization factor of p (i) j (σ j ). After solving the message-passing Eqs. (11), the marginals on the original graph are reconstructed from Eqs. (10).

B. The distributional equations for an ensemble of networks
Equations (10) and (11) yield an approximation for the local marginals of an ensemble of spin models defined on a single graph instance with N 1 nodes. These equations become asymptotically exact for K 1 = 0 when N → ∞ [44,45], due to the existence of a vanishing fraction of short cycles in the graph. Nevertheless, if the coupling strengths are random (K 1 > 0), Eqs. (10) and (11) exhibit a large number of solutions at low temperatures [44,45], which reflects the existence of many local minima in the free-energy [7,64]. In this regime of parameters, Eqs. (10) and (11) only provide an approximate description of the system. Inasmuch as p i (σ i ) and p (i) j (σ j ) are (positive) random functions whose shape fluctuates along the nodes of a network, it is instrumental to adopt an ensemble viewpoint and work with the functional probability densities of p i (σ i ) and p (i) j (σ j ), defined respectively as and with δ F [f (σ)] representing the functional Dirac-δ defined over the space of all possible functions f (σ).
Let us obtain the distributional equations for the functionals W [p] and R[p]. Equations (10) and (11) can be intuitively rewritten as where and Since the left and right hand sides of Eqs. (14) and (15) where Dh is a functional integration measure, while yields the functional distribution of i (σ i ) by setting n = 0 or n = 1, respectively. Equations (18)  The central quantity F n [h] is determined, in the limit N → ∞, by the degree distribution p k , the distribution of coupling strengths p J , and the functional distribution R[p]. As we will show below, F n [h] simplifies in the high-connectivity limit c → ∞.

C. Local and global ensemble averages
The distributions W [p] and R [p] are key quantities in the study of the macroscopic behavior of spin models on networks as they allow to compute ensemble averages of the spins and derive equations for the order-parameters [51,52]. Due to the fluctuations in the local structure of a random network, it is important to clearly distinguish between local and global ensemble averages. Let O(σ i ) be an observable defined in terms of σ i . The local ensemble average of O(σ i ) on the original graph is defined as while the local average on the cavity graph reads Local averages only depend on the random functions p i (σ) and p (j) i (σ) defined at node i. By taking the averages of Eqs. (21) and (22) with respect to the functional distributions of p i (σ) and p For instance, the choice O(σ) = σ yields the global magnetization m on the original graph, and the global magnetization M on the cavity graph, Notably, the present formalism also gives access to the fluctuations of the random function p i (σ) along the different nodes. Instead of working directly with the functional distribution W [p], whose domain is infinite dimensional, it is sensible to study the moments of p i (σ) for fixed σ In particular, the variance quantifies the spread of the functional shape of the local marginals p i (σ) around their average p(σ) W .

IV. THE HETEROGENEOUS MEAN-FIELD THEORY
The heterogeneous mean-field theory describes spin models on high-connectivity random graphs characterized by strong degree fluctuations, in contrast to the standard, homogeneous mean-field theory, whose predictions are limited to random graphs with vanishing degree fluctuations.
The heterogeneous mean-field equations are derived by taking the c → ∞ limit of Eqs. (18) and (19) for arbitrary p k . The core of the calculation lies in the high-connectivity limit of the distribution F n [h]. Before performing this calculation in full generality, let us consider the behavior of the random variable H i (σ i ) for large c in the particular case of Ising spins, in order to gain insight into the physical meaning of H i (σ i ) and its distribution F 0 [h]. Thanks to the scaling of the moments of J il with c, the coupling strengths become very small for c 1 and it suffices to expand Eq. (16) in powers of J il up to O(J 2 il ), leading to the following expression with σ p (i) l the local magnetization at node l in the absence of i ∈ ∂ l . Equation (29) reveals that H i (σ i ) can be regarded as the interaction energy of the spin σ i with an effective random field, specified by the local magnetizations of the neighboring spins in the absence of σ i . The physical meaning of the random variable H (j) i (σ i ) is completely analogous. Equation (29) involves two sums over a large number of independent random variables. Hence it is tempting to invoke the central limit theorem and argue that H i (σ i ) follows a Gaussian distribution, which is in fact correct for Erdös-Rényi and regular random graphs [42,43]. However, the total number of terms in the sums of Eq. (29) is itself a random variable, whose fluctuations are controlled by the degree distribution p k . As we will show below, the central limit theorem breaks down when p k is not sufficiently concentrated around its mean value c [48], which leads to a family of heterogeneous meanfield equations that explicitly depend on p k .
With the purpose of calculating the c → ∞ limit of F n [h] for arbitrary p k , it is convenient to introduce the characteristic functional from which F n [h] is determined via the inverse Fourier transform where S D is the total number of single-spin states and noting that the above expression factorizes in terms of a product over the neighborhood of a single node, it is straightforward to show that G n [t] assumes the compact form in the limit c → ∞. The functional η[t] is given by and ν(g) is the high-connectivity limit of the distribution of rescaled degrees k 1 /c, . . . , k N /c The assumption that the distribution of rescaled degrees attains a well-defined limit ν(g) for c → ∞ underlies the validity of Eq. (33). The next step is to calculate η[t] from Eq. (34). Since the coupling strengths become very weak for c 1, we expand Eq. (34) in powers of J r and in turn compute the limit c → ∞, finding the D × D connected correlation matrix χ between the components of a spin and the D × D correlation matrix C between the components of the local magnetization The order-parameters of Eqs. (37)(38)(39) are defined in terms of the distribution R[p] of the single-site marginals on the cavity graph. The final step is to perform the inverse Fourier transform in Eq. (31). Nonetheless, the exponential in G n [t] has a quadratic term in t(σ), as can be noticed by inserting Eq. (36) in Eq. (33). This quadratic term, given by is linearized by a multivariate Gaussian integral over z = (z 1 z 2 . . . z D ) T with the Gaussian measure and dz = D α=1 dz α . By carrying out the above linearization, we arrive at the appealing expression for where H eff (σ|g, z) is the effective Hamiltonian of a single spin Inserting Eq. (42) in Eq. (31), the functional integral over t(σ) is immediately performed, leading to Equations (43) and (44) constitute the central result of this paper, as they provide the analytic expression for the highconnectivity limit of the distribution of the coupling energy between a single spin and the effective field coming from its neighborhood. These equations are valid for the generic spin model defined by Eq. (4). Remarkably, the distribution of H eff is generally not Gaussian, but it explicitly depends on the distribution ν(g) of rescaled degrees, which means the model retains information about the graph structure even when c → ∞. The non-Gaussian nature of H eff is a direct consequence of the breakdown of the central limit theorem, as discussed in the beginning of this section. For random graph models in which ν(g) = δ(g − 1), such as regular and Erdős-Rényi graphs [1], the central limit theorem holds and H eff follows a Gaussian distribution. The strength of the degree fluctuations for c → ∞ is quantified by the variance of ν(g) In the present context, we say that a network model is homogeneous if ∆ 2 ν = 0, whereas an heterogeneous network is characterized by ∆ 2 ν > 0. Plugging Eq (44) back into Eqs. (18) and (19), we find from which a closed set of equations for the order-parameters is derived from Eqs. (37)(38)(39). Equation (46) summarizes the gist of the mean-field approach: the macroscopic behavior of an infinitely large system is reduced to an effective problem of a single spin, in which the single-site states are sampled from a Boltzmann distribution with effective energy H eff (σ|g, z). Equations (46) and (47), together with the order-parameter Eqs. (37)(38)(39), define what we call the heterogeneous mean-field theory. It is unfeasible to solve the generic order-parameter equations that follow from Eqs. (46) and (47). Thus, in the sequel we explore the role of degree fluctuations in specific models.

A. Ferromagnetic couplings
In this section we address the impact of heterogeneous degrees on the behavior of ferromagnetic spin models. For K 1 = 0, the coupling energy is given by and the distribution of the effective field gK 0 M on a spin σ follows from ν(g). Due to the orthogonal invariance of the Hamiltonian, Eq. (4), the magnetization m is determined by its absolute value |m|. For heterogeneous networks we get where I a (x) is the modified Bessel function of the first kind. The order-parameter M is the magnetization on the ensemble of cavity graphs, and its absolute value solves the selfconsistency equation Equations (49) and (50) generalize the fully-connected meanfield equations for ferromagnetic models with D-dimensional spins [32] to the case of heterogeneous networks. In particular, for Ising spins (D = 1) we have I 1 2 (x)/I − 1 2 (x) = tanh(x), and the scalar magnetization m is determined from which generalizes the Curie-Weiss mean-field equations [39] to heterogeneous random graphs. Both parameters, |M | and |m|, quantify the coherence of vector spins, but in slightly different ways: |m| is the centroid of the spin configuration {σ} on the surface of the Ddimensional hypersphere, while |M | is the average of the spins weighted by the rescaled degree g. For homogeneous random graphs, identified by ν(g) = δ(g − 1), Eqs. (49)  from fully-connected graphs. For heterogeneous networks, however, these parameters do not coincide, since highly connected nodes (hubs) contribute more to the integral in Eq. (50). Theoretical approaches for the Kuramoto model (D = 2) on heterogeneous networks have often characterized the synchronization phase transition solely in terms of |M | [27]. While this is appropriate for homogeneous networks, it leads to discrepancies between theory and simulations for the behavior of |m| when the network has strongly fluctuating degrees, such as in scale-free networks [65]. Here we do not face this problem of choosing suitable observables beforehand, since both quantities, |m| and |M |, and the relationship between them emerge naturally from the high-connectivity limit of local tree-like networks, without any additional assumption on the network topology.
The solutions of Eqs. (49) and (50) allow to study the effect of degree fluctuations on the phase diagram of the ferromag-netic spin model. In the limit T → ∞, we find the solution M = m = 0 and the system lies in the paramagnetic phase. By expanding Eq. (50) up to O(|M |), we conclude that |M | has a nontrivial solution provided where ∆ 2 ν is the heterogeneity parameter, given by Eq. (45). The above expression defines a line in the parameter space at which the model undergoes a continuous phase transition between the paramagnetic and the ferromagnetic phase. Equation (53) shows that the size of the ferromagnetic phase increases as a function of the heterogeneity parameter. Indeed, the critical value of βK 0 vanishes as ∆ 2 ν → ∞ for any finite D, which is in line with previous results for the critical temperature of the Ising model on scale-free networks [49,66].
In the rest of this section, we present explicit results for random networks with the negative binomial degree distribution. Substituting Eq. (7) in Eq. (35), we calculate the corresponding distribution of rescaled degrees where α > 0 is related to ∆ 2 ν as follows Equation (54) allows us to interpolate continuously between homogeneous networks (α → ∞) and strongly heterogeneous networks (α → 0) by varying a single parameter α that controls the degree fluctuations in the high-connectivity limit. Figure 2 illustrates the effect of α on the behavior of |m| as a function of the coupling strength K 0 for the Ising (D = 1) and the Kuramoto model (D = 2). The ferromagnetic or synchronized phase, identified by |m| > 0, appears through a continuous phase transition at the critical coupling strength predicted by Eq. (53). Moreover, figure 2 compares the theoretical results, derived from the numerical solutions of Eqs. (49) and (50), with data from numerical simulations of the dynamics of each model. The agreement between theory and simulations is excellent in both cases, confirming the exactness of Eqs. (49) and (50).
Even though the critical value of βK 0 vanishes for α → 0, which could suggest the absence of a paramagnetic phase in the strongly heterogeneous regime, we note that |m| becomes gradually smaller inside the ordered phase as α is reduced. To better understand how the heterogeneity parameter α impacts the ordered phase of spin models, we present in figure 3 the behavior of |m| as a function of α for the Kuramoto model. The main outcome is that, for βK 0 1, deeply in the synchronized phase, |m| has the functional form In the sector α > 1, we are able to show analytically that δ(α) = 1 and f (α) = (50) for βK 0 1. In the range 0 < α ≤ 1, there are strong degree fluctuations and the perturbative expansion fails, leading to a divergence in |m|. In spite of that, figure 3 shows that the results for different βK 0 collapse onto a single curve by choosing δ(α) = α in the sector 0 < α ≤ 1, confirming the functional form put forward in Eq. (56). This numerical procedure does not give the explicit form of f (α), but it clearly shows that lim α→0 f (α) = 1.
On the whole, Eq. (56) leads to two interesting conclusions. First, the magnetization |m| vanishes inside the ordered phase as α → 0 (∆ ν → ∞). This is a surprising finding given that highly heterogeneous networks have been regarded as optimal structures for the emergence of synchronization. Indeed, the vanishing of the critical coupling for the Kuramoto model on heterogeneous networks has been reported in previous studies [27,[67][68][69][70][71], but it has remained unclear how |m| behaves when the variance of the degree distribution diverges. Our results confirm that, in the limit c → ∞, the critical coupling does vanish as ∆ ν → ∞, but they also reveal that |m| goes to zero, suggesting that highly heterogeneous networks do not sustain any level of synchronous oscillation in such limit. The second interesting conclusion is that the point α = 1 separates two different regions inside the synchronous phase, each one marked by a particular behavior of the magnetization |m| for βK 0 1. For 0 < α ≤ 1, degree fluctuations have a significant impact in the ordered phase and the magnetization is given by 1 − |m| ∼ (βK 0 ) −α . For α > 1, degree fluctuations are immaterial and the magnetization behaves as 1 − |m| ∼ (βK 0 ) −1 , regardless of the heterogeneity parameter α. Interestingly, the function f (α) displays a cusp at α = 1, reflecting the non-analytic behavior of |m| as a function of α. Although Eq. (56) and figure 3 are results for the Kuramoto model (D = 2), we expect that the above two conclusions hold for arbitrary D < ∞.
The main results discussed up to now are summarized in figure 4, in which we present the phase diagram of the Kuramoto model on heterogeneous networks. Notice that for α → ∞ the transition line approaches βK 0 = 2, which is the critical point of fully connected networks with identical stochastic oscillators [72,73].
We end this section by discussing another consequence of heterogeneous degrees, namely the lack of correspondence between local and global ensemble averages. We present results for the Kuramoto model (D = 2), but the main conclusions should be once more valid for any D < ∞. The state of a spin σ i in D = 2 is fully specified by the phase φ i ∈ (−π, π], distributed according to the local marginal P i (φ i ). Combining Eqs. (27), (46), and (48), it is straightforward to write the analytic expression for the moments of the random function P i (φ i ) where we have set M = (|M | 0) T without loosing any generality. The above equation yields all moments of the functional distribution W [P ] of the single-site marginals P i (φ i ) over the entire network. In particular, we are interested in the mean µ P (φ) and the variance σ 2 Figure 5 illustrates the effect of the heterogeneity parameter α on the mean µ P (φ) inside the synchronized phase of the Kuramoto model. The different panels are for the same value of the magnetization |m|. More importantly, the shaded area in figure 5 represents the dispersion [µ P (φ)−σ P (φ), µ P (φ)+ σ P (φ)] around the mean µ P (φ), quantifying the fluctuations of the random function P i (φ i ) from site to site. Each set of gray points in figure 5 is the histogram of a single-site phase, which weights the amount of time that a randomly chosen spin σ i has a certain orientation φ i . Each one of these single-site histograms is constructed from a single run of the numerical simulation by storing the values of a given φ i over a long time interval after the system has reached equilibrium.
For homogeneous networks (α → ∞), the variance σ 2 P (φ) is zero and the distribution W [P ] is peaked at the mean value µ P (φ). Accordingly, each single-site histogram coincides with µ P (φ), confirming that all spins are identical, in the sense that the fluctuations of the individual phases are described by the same distribution P (φ) W . Thus, local and global ensemble averages become equivalent in the highconnectivity limit of ferromagnetic spin models on homogeneous networks, such as regular and Erdős-Rényi random graphs. In contrast, we observe that in the case of heterogeneous networks (α < ∞), the size of the shaded area in figure  5 becomes gradually larger as α decreases, confirming that the functional form of P i (φ i ) fluctuates from site to site, and the statistical properties of the spins remain different from each other, even in the high-connectivity limit.
Although the stochastic Kuramoto model on heterogeneous networks has been extensively studied in previous works (see, e.g. [27,73,74]), little attention has been paid to the characterization of site to site fluctuations. It is thus noteworthy that the theory presented in Sec. IV allows us to study not only global observables of the system, but also fluctuations at the level of individual nodes.

B. Random couplings
In this subsection we discuss the effect of heterogeneous degrees in mean-field spin models with random coupling strengths. We restrict ourselves to the simplest case of Ising spins (D = 1), with σ i ∈ {−1, 1}, in which the weights p i (σ) of a single-site marginal are encoded in the two-dimensional vector p i = (p i (1) p i (−1)) T .
The mean-field equations for the high-connectivity limit of Ising spins coupled through the random links of heterogeneous networks follow from Eqs. (43), (46), and (47). The distributions of single-site marginals read where the distributions ω 0 (h) and ω 1 (h) follow from by setting n = 0, 1. The function P g (z) is the normal probability density The order-parameters M and Q, appearing in Eq. (61), are defined as with d p = dp(1)dp(−1). The local magnetization σ p is computed by replacing the integral in Eq. (22) by a discrete sum over σ ∈ {1, −1}. The quantities M and Q are, respectively, the global magnetization and the Edwards-Anderson order-parameter [40] on the ensemble of cavity graphs. Inserting Eq. (60) in the above definitions and integrating over p, we arrive at the self-consistency equations The global magnetization and the Edwards-Anderson order-parameter follow from the distribution W ( p) of marginals on the original graph. Substituting Eq. (59) in the above definitions and integrating over p, we get Equations (65), (66), (69) and (70) generalize the replicasymmetric mean-field equations of Ising spin-glass models with random interactions to the case of heterogeneous networks. For homogeneous networks, characterized by ν(g) = δ(g − 1), we obtain ω 0 (h) = ω 1 (h) and the above equations reduce to the replica-symmetric equations of the Sherrington-Kirkpatrick model [40,41]. The solutions of the selfconsistency Eqs. (65) and (66) allow us to address the impact of heterogeneous degrees on the phase diagram of the system. According to Eq. (59), the local marginal p i (σ) at node i is a Boltzmann factor parametrized by the scalar effective field h i,eff [64], given in terms of the local magnetization σ pi as follows where p i is drawn from W ( p). Combining Eqs. (29) and (14), one can show that h i,eff is also given by Eq. (1). The object ω 0 (h) is nothing more than the empirical distribution of effective fields which also gives information about the fluctuations of the weights p i throughout the network. The random variable h i,eff is defined for a single realization of the random graph and its distribution P eff (h) is determined by the fluctuations in the graph structure. According to Eq. (61), ω n (h) is determined only by the probability densities P g (z) and ν(g), which enables to derive the general formula for arbitrary n. Note that Eq. (73) is valid for an arbitrary distribution ν(g) of rescaled degrees. If the network is homogeneous (ν(g) = δ(g −1)), then ω 0 (h) is a Gaussian distribution with mean K 0 m and standard deviation K 1 √ q EA [40,41].
If the network is heterogeneous, then ω 0 (h) is not Gaussian, which stems from the breakdown of the central limit theorem due to the large variance of the random number of summands in Eq. (1). The average p(σ) W in the present formalism corresponds with the replica-symmetric ansatz employed in the study of spin models on random graphs [64]. Before discussing the phase diagram for K 1 > 0, we point out that the heterogeneous mean-field theory is not exact in the whole phase diagram. In the context of the replica method [7], we say that equations (65), (66), (69) and (70) are the replica symmetric (RS) solution of the model. In the case of the Sherrington-Kirkpatrick model [40,41], obtained for the choice ν(g) = δ(g − 1), the RS solution becomes unstable at low temperatures owing to the existence of an exponentially large number of metastable states [7]. The sector of the phase diagram where the replica symmetric theory fails is bounded by the so-called Almeida-Thouless (AT) line [50], which is determined by the eigenvalues of the Hessian of the free energy.
It is thus important to establish the limit of validity of the heterogeneous mean-field theory and ask how degree fluctuations impact the location of the AT line. Instead of tackling this problem by analyzing the Hessian of the free-energy, here we follow an alternative strategy [75][76][77], based on the number of solutions of the cavity equations (11). By defining {p where the joint distribution R 12 ( p 1 , p 2 ) of weights is defined as follows (75) When the RS theory is stable, Eqs. (11) exhibit a single solution, the distribution R 12 ( p 1 , p 2 ) becomes diagonal and we simply have ρ = Q, with Q given by Eq. (66). When the RS theory is unstable, the cavity Eqs. (11) have a large number of fixed-point solutions resulting from the existence of a large number of extrema in the free-energy. To calculate the AT line, it suffices to consider the stability of the RS theory under a perturbation of ρ. By plugging ρ = Q+δρ in Eq. (74) and expanding its right hand side up to O(δρ), we conclude that the RS theory is unstable provided where ω 2 (h) is defined through Eq. (73). Equation (77) is the condition for the breakdown of the RS theory, generalizing the AT line of Ising spin models with random couplings to heterogeneous networks. For homogeneous networks with ν(g) = δ(g − 1), we recover the AT line for the Sherrington-Kirkpatrick model [50]. The details involved in the derivation of Eq. (77) are explained in appendix A. We remark that, although we have limited ourselves to Ising spins, the current approach to obtain the AT line can be also applied to Ddimensional spins. Let us discuss our results for the phase diagram of the Ising model on heterogeneous networks with K 1 > 0. Equations (65) and (66) exhibit the paramagnetic solution M = Q = 0 at sufficiently high temperatures. The model undergoes a second-order phase transition from the paramagnetic to the ferromagnetic phase, characterized by M = 0 and Q = 0, at the critical temperature derived from an expansion of Eqs. (65) and (66) for |M | 1 and Q 1. As long as K 1 > 0, the model has a spinglass phase, in which M = 0 and Q = 0. By setting M = 0 and expanding Eq. (66) for Q 1, we find the critical temperature that marks the second-order transition between the paramagnetic and the spin-glass phase.
Equations (78) and (79) yield the continuous phase transitions for an arbitrary degree distribution ν(g). The AT line as well as the continuous transition between the ferromagnetic and the spin-glass phase are calculated by numerically solving the order-parameter Eqs. (65) and (66). As before, we present results for networks with a negative binomial degree distribution, in which ν(g) is given by Eq. (54). The function ω n (h) (n = 0, 1, 2) is the key quantity that determines the orderparameters and the AT line. Inserting Eq. (54) in Eq. (73) and integrating over g [78], we get the analytic expression where and K a (x) is the modified Bessel function of the second kind with order a [78]. In the limit α → ∞, ω n (h) converges to a Gaussian distribution with mean K 0 M and variance K 2 1 Q, independently of the index n = 0, 1, 2. Figure 6 depicts the phase diagram of the Ising spin-glass model on networks with a negative binomial degree distribution for different values of the heterogeneity parameter α. The phase diagram of the Sherrington-Kirkpatrick model is recovered for α → ∞ [40,41]. The different critical lines in figure  6 meet at the point which serves as a useful guide to understand the effect of α on the phase diagram. From Eq. (82) we note that, in the limit α → 0, the ferromagnetic phase essentially expands over the entire phase diagram, while the spin-glass phase is confined to an arbitrary small region around K 0 /K 1 0. The critical temperature between the paramagnetic and the spin-glass phase diverges as 1/ √ α for α → 0. According to figure 6, the location of the AT line shows that the RS theory is unstable throughout the spin-glass phase and in the low-temperature sector of the ferromagnetic phase, for any value of α. Although decreasing values of α seem to gradually stabilize the replica-symmetric ferromagnetic phase for lower temperatures, the impact of the network heterogeneity on the AT line is fundamentally different in the regime T → 0. In fact, the inset in figure 6 demonstrates that the critical value of K 0 marking the AT transition for T → 0 is a nonmonotonic function of α, with a maximum around α 0.5. Moreover, the inset suggests that, for strong heterogeneous networks with α → 0, the RS theory may become stable in the low-temperature sector of the ferromagnetic phase.
We end this section by studying how the network heterogeneity impacts the distribution P eff (h) of effective fields. The full analytic expression for P eff (h) is obtained directly from the distribution ω 0 (h) (see Eq. (72)), defined by means of Eq. (61). We point out that P eff (h) is always given by P eff (h) = δ(h) in the paramagnetic phase.
Let us first consider P eff (h) for networks with ferromagnetic couplings (K 0 > 0 and K 1 = 0). The distribution of effective fields for homogeneous networks has the typical form P eff (h) = δ(h − K 0 m) of the Curie-Weiss model, while we obtain for heterogeneous networks with arbitrary degree distribution. The above equation holds for the choice M > 0, and the symbol Θ(h) represents the Heaviside step function. For a negative binomial degree distribution, the above expression becomes We see that strong degree fluctuations lead to striking modifications in P eff (h) when compared to homogeneous networks. The first interesting aspect concerns the behavior of P eff (h) for large fields. For α < 1, P eff (h) exhibits a power-law regime P eff (h) ∼ h α−1 (h 1) up to the point h = O(K 0 M/α), above which it decays exponentially fast. Therefore, the power-law decay of P eff (h) persists over an arbitrarily large range of h as α → 0. The second interesting effect of the network heterogeneity occurs in the behavior of P eff (h) around h = 0. In the limit h → 0, P eff (h) converges to a constant if α ≥ 1, whereas it diverges as P eff (h) ∼ h α−1 if α < 1. Thus, strong degree fluctuations induce the appearance of a substantial fraction of vanishing effective fields, which has a detrimental effect on the ordered phase.
In the sequel, we discuss how the randomness in the coupling strengths impacts P eff (h). For homogeneous networks with K 1 > 0, the central limit theorem can be applied to Eq. (1) and P eff (h) is a Gaussian distribution with mean K 0 m and variance K 2 1 q EA [40]. The analytic expression of P eff (h) for networks with negative binomial degrees and K 1 > 0 is derived by setting n = 0 in Eq. (80). Similarly to the above results for K 1 = 0, P eff (h) features a crossover between a power-law and an exponential decay for sufficiently large |h|. The main difference with respect to K 1 = 0 appears around h = 0. In the regime |h| → 0, P eff (h) converges to a constant if α > 1/2, whereas it displays the following asymptotic forms if α ≤ 1/2 Thus, the interplay between random coupling strengths and heterogeneous degrees leads to a logarithmic divergence in P eff (h), a feature that appears as well in the spectral density of heterogeneous networks [79]. We further note that fluctuations in the coupling strengths mitigate the divergence around h = 0, in the sense that the power-law exponent in Eq. (85) is smaller when compared to the case K 1 = 0. The comparison between numerical simulations for the effective fields and the analytic expressions for P eff (h) in the cases α → ∞ and α = 1 are shown in figure 1. The agreement between the theory and simulations is excellent, confirming the exactness of our theoretical findings for sufficiently high temperatures. The heterogeneous mean-field theory is derived by taking the high-connectivity limit c → ∞ after the limit N → ∞. As such, one expects that our theory is valid for high-connectivity networks in which the density of edges c/N goes to zero for N → ∞. In fact, figure 7 presents numerical simulation results for the global magnetization and the Edwards-Anderson order-parameter of spin models on heterogeneous networks with average degree c = √ N and with the negative binomial degree distribution of Eq. (7). For homogeneous networks (α → ∞), such as regular and Erdős-Rényi random graphs, the order-parameters flow as N → ∞ to their fully-connected values, obtained from the solution of the Sherrington-Kirkpatrick model. For heterogeneous networks, in which α is finite, the order-parameters converge, in the limit N → ∞, to the results of the heterogeneous mean-field theory derived in the present paper. Figure 7 confirms that the mean-field theory presented in this work describes spin models on heterogeneous networks where c scales as c ∝ N b , with 0 < b < 1. This regime of connectivity lies between sparse networks (b = 0) and diluted networks (b = 1) [80][81][82][83].

V. DISCUSSION AND CONCLUSIONS
Mean-field theories are formidable tools to study the macroscopic properties of spin models on networks. The most well-studied family of mean-field theories are realized on fully-connected architectures, in which a given spin interacts with all others. As a network becomes more densely connected, it is natural to expect that local fluctuations in the network structure are gradually washed out, and the macroscopic properties of the underlying spin model converge to those of a fully-connected system. In this work we have shown that this is generally not the case. We have derived a novel class of exact mean-field equations that explicitly depend on the degree distribution and that apply to the high-connectivity limit of heterogeneous networks. Fully-connected mean-field theories, in contrast, are limited to homogeneous networks, for which the degree distribution is peaked at its mean value.
Put differently, we have proven that the high-connectivity limit of spin models on networks is nonuniversal, as it depends on the full degree distribution. On the other hand, the universality with respect to the randomness in the coupling strengths is robust to degree fluctuations. In fact, the heterogeneous mean-field equations only depend on the first two moments of the distribution of couplings. The nonuniversal behavior of the heterogeneous mean-field theory is accompanied by the failure of the central limit theorem for the effective fields, caused by large fluctuations in the number of summands in Eq. (1). Apart from a few exceptions [79,84,85], the consequences of this interesting mechanism for the breakdown of the central limit theorem to statistical physics remain largely unexplored. We have illustrated its crucial role for the equilibrium of spin models, but one can envisage the far-reaching importance of this mechanism for a variety of problems on networks, such as the nonequilibrium dynamics of spin models [86,87], the storage capacity of neural networks [14], and the stability of large dynamical systems [23].
We have presented several results that highlight the importance of degree fluctuations to the high-connectivity limit of spin models, with particular emphasis on the mean-field theory of synchronization and of Ising spin-glasses. Although the heterogeneous mean-field equations are valid for any degree distribution, most of the explicit results have been derived for a negative binomial degree distribution, which allows to smoothly interpolate between homogeneous and heterogeneous networks by changing the variance of the degree distribution.
Degree fluctuations have a conspicuous effect on the ferromagnetic phase. As the variance of the degree distribution increases, the ferromagnetic phase gradually covers the whole phase diagram of the Kuramoto model, suggesting that the oscillators synchronize at any coupling strength, as argued in several previous works [27,[67][68][69][70][71]. However, we have shown that not only the critical coupling becomes arbitrarily small as the network heterogeneity diverges, but also the magnetization in the synchronous phase. Therefore, our results of Sec. IV A show that instead of facilitating synchronization, degree fluctuations actually hamper the formation of an ordered synchronous component, and ultimately inhibit the emergence of any coherent behavior when the variance of the degrees is infinitely large.
Moreover, the magnetization of the Kuramoto model displays a cusp that separates the ferromagnetic phase in two qualitatively distinct regions, each one characterized by a different impact of degree fluctuations on the magnetization. Such non-analytic behavior is not exclusive of the magnetization, but it seems to be a generic feature of the macroscopic behavior. Indeed, the effective field distribution of Ising spin-glass models on heterogeneous networks exhibits either a power-law or a logarithmic divergence at zero field, which contrasts with the Gaussian effective fields of fullyconnected models. This divergence reflects the substantial fraction of spins with a vanishing local magnetization inside the ferromagnetic phase.
As another genuine feature of strong degree fluctuations in the ferromagnetic phase, we have found that the shape of the single-site phase distribution of the Kuramoto model fluctuates from site to site, in sharp contrast to fully-connected models, for which the functional form of the local phase distribution is fixed. These results follow from the mean and the variance of the single-site phase distribution. In fact, the path-integral formalism developed in this paper gives access to all moments of the marginal distribution of single-site configurations. This is particularly relevant in the case of vector spins, for which the functional distribution of local marginals is not parametrized in terms of a finite number of fields. Thus, when compared to other mean-field techniques for coupled oscillators [27], the path-integral formalism has the technical advantage of allowing the calculation of dynamical properties of vector spin models at the level of individual nodes. It would be interesting to investigate the individual phase fluctuations for three-dimensional spins [38], a problem that can give insights into the synchronous dynamics of swarms and flocks in three dimensions.
We have studied the interplay between random coupling strengths and degree fluctuations in Ising spin models. Figure 6 generalizes the replica-symmetric phase diagram of the fully-connected Sherrington-Kirkpatrick model [41] to heterogeneous networks. In this case, the heterogeneous meanfield theory is not exact on the entire phase diagram, but it becomes unstable at temperatures below the Almeida-Thouless line [50]. In spite of that, the results of figure 6 allow to conclude that degree fluctuations promote the ferromagnetic phase in detriment of the spin glass phase. In addition, the low-temperature sector of the Almeida-Thouless line exhibits a non-monotonic behavior as the degree fluctuations increase. In view of this distinctive behavior and given that the replica symmetry breaking theory on sparse networks is a notorious difficult problem [44,45,88], it would be very interesting to employ the replica symmetry breaking machinery [89,90] and derive the exact version of the heterogeneous mean-field theory. This would give an alternative, simpler route to exploit how the complex picture describing the spin-glass phase [7] is modified by the presence of network heterogeneities.
The heterogeneous mean-field theory has been derived by taking the high-connectivity limit c → ∞ after the thermodynamic limit N → ∞, where c is the average degree and N is the total number of spins. An implicit assumption in this procedure is that the local tree-like structure of the network is preserved and hence the distributional cavity equations remain valid. This can be only achieved if the density of links c/N approaches zero for c → ∞, which suggests that the heterogeneous mean-field theory should apply when c ∝ N b (0 < b < 1). We have confirmed this conjecture by means of numerical simulations. The connectivity regime c ∝ N b (0 < b < 1) lies between sparse (c = O(1)) and diluted (c = O(N )) networks [80][81][82][83]. Even though this intermediate connectivity range, called the extreme diluted regime, has been known for a long time in the field of neural networks [91,92], it has been studied only in the case of homogeneous networks, for which degree fluctuations are unimportant. Here we have examined the extreme diluted regime of heterogeneous networks and unveiled its nonuniversal features. In the diluted regime, degree fluctuations are of O(N 0 ) and the macroscopic behavior of spin models is captured by the fullyconnected mean-field equations. The universality properties of spin models in the three different regimes of connectivity are summarized in  On the whole, our findings demonstrate that network heterogeneities, here expressed by degree fluctuations, play a surprisingly important role in the high-connectivity limit of spin models. Other types of topological features, such as modular structure and the presence of loops [93][94][95], should as well play a fundamental role in the high-connectivity behavior. While spin models on sparse networks pose many technical challenges [44,45,88], the mean-field theory of fullyconnected models has a simpler formal structure [7], at the cost of completely neglecting the network structure. Our work puts forward a novel family of mean-field theories that explicitly takes into consideration the network structure and, at the same time, are simple enough that they can be thoroughly analyzed. This framework opens the door to the development of similar mean-field theories for other complex systems.
completely analogous to those discussed in section IV, with the main difference that the two copies of the system become correlated after taking the ensemble average. The resulting expression for R 12 ( p 1 , p 2 ) reads where D u is the Gaussian bivariate measure of u = (u 1 u 2 ) T with A denoting the 2 × 2 matrix The order-parameters M 1 , M 2 , Q 1 , Q 2 , and ρ are determined from R 12 ( p 1 , p 2 ) as follows with the local average σ p defined by Eq. (22). The orderparameters M and Q are duplicated for the simple reason that we are dealing with two replicas of the original system. The quantity ρ measures the correlation between the local magnetizations of each replica. If the RS solution is stable, then R 12 ( p 1 , p 2 ) simplifies to the diagonal form of Eq. (76) and the order-parameters fulfill M 1 = M 2 = M , Q 1 = Q 2 = Q, and ρ = Q, where M and Q follow from the solutions of Eqs. (65) and (66). The correlation order-parameter ρ fulfills the self-consistent equation derived by substituting Eq. (A1) in Eq. (A7). In order to analyze the stability of Eq. (A8) around the RS solution, it is convenient to make an orthogonal change of integration variables. Let T be the 2 × 2 orthogonal matrix that diagonalizes The matrix T and the eigenvalues of A have, respectively, the explicit forms where we have defined The dependency of some quantities with respect to the orderparameters has been omitted in order to simplify the notation. By making the orthogonal change of integration variables u = T z, with z = (z 1 z 2 ) T , we can rewrite Eq. (A8) as ×tanh β gK 0 M 1 + K 1 √ g t 11 λ + z 1 + t 12 λ − z 2 ×tanh β gK 0 M 2 + K 1 √ g t 21 λ + z 1 + t 22 λ − z 2 , after rescaling the components of z as z 1 → λ + z 1 and z 2 → λ − z 2 . The normal probability density P g (z) is defined by Eq. (62). Equation (A11) is the proper starting point to analyze the stability of the RS theory. The fixed-point solutions describing the macroscopic states of the system composed of two replicas are given in terms of five order-parameters: M 1 , M 2 , Q 1 , Q 2 , and ρ. In order to probe the stability of the RS solution, it suffices to consider fluctuations solely in the direction of ρ. Thus, by setting M 1 = M 2 = M , Q 1 = Q 2 = Q, and ρ = Q + δρ in Eq. (A11), and then expanding its right-hand side up to O(δρ), we obtain the following equation for the perturbation δρ δρ = δρ β 2 K 2 1 ∞ 0 dg ν(g) g 2 ∞ −∞ dz P g (z) × sech 4 β gK 0 M + K 1 gQz .
If the coefficient of δρ in the above equation is larger than one, then the RS fixed-point solution is unstable, since the iteration of the above linear equation leads to the growth of the perturbation δρ. After a simple change of integration variables, this condition becomes identical to Eq. (77).
For D = 2, the spins are two-dimensional vectors with unit norm that rotate on a plane. As such, we can simulate the spin system as a set of stochastic Kuramoto oscillators [26,28], in which the state of a node i at time t is fully specified by the phase φ i (t) ∈ (−π, π]. The stochastic dynamics of the model is dictated by the following set of equations [27] where ξ i (t) is a Gaussian white noise that satisfies with T being the temperature. For all numerical results shown in this paper, we numerically integrate the stochastic equations (B3) with the Heun's method with time step dt = 0.01. The complex magnetization m c (t) is calculated as the average of the phasors that rotate in the complex unit circle, that is m c (t) = R(t)e iψ(t) = 1 N where 0 ≤ R(t) ≤ 1 measures the level of phase-coherence in the system, and ψ(t) is the average phase. In this context, the absolute value |m(t)| of the vector magnetization m(t) in section IV A is obtained as |m(t)| = R(t). In particular, the long-time behavior of the magnetization |m(t)| in Fig