Analytic solution of the two-star model with correlated degrees
Maíra Bolfe, Fernando L. Metz, Edgar Guzmán-González, Isaac Pérez Castillo
AAnalytic solution of the two-star model with correlated degrees
Maira Bolfe
Physics Department, Federal University of Santa Maria, 97105-900 Santa Maria, Brazil.
Fernando L. Metz ∗ Physics Institute, Federal University of Rio Grande do Sul, 91501-970 Porto Alegre, Brazil andLondon Mathematical Laboratory, 18 Margravine Gardens, London W6 8RH, United Kingdom
Edgar Guzm´an-Gonz´alez and Isaac P´erez Castillo † Departamento de F´ısica, Universidad Aut´onoma Metropolitana-Iztapalapa,San Rafael Atlixco 186, Ciudad de M´exico 09340, Mexico
Exponential random graphs are important to model the structure of real-world complex networks. Here wesolve the two-star model with degree-degree correlations in the sparse regime. The model constraints the averagecorrelation between the degrees of adjacent nodes (nearest neighbors) and between the degrees at the end-pointsof two-stars (next nearest neighbors). We compute exactly the network free energy and show that this modelundergoes a first-order transition to a condensed phase. For non-negative degree correlations between nextnearest neighbors, the degree distribution inside the condensed phase has a single peak at the largest degree,while for negative degree correlations between next nearest neighbors the condensed phase is characterized by abimodal degree distribution. We calculate the degree assortativities and show they are non-monotonic functionsof the model parameters, with a discontinuous behavior at the first-order transition. The first-order critical lineterminates at a second-order critical point, whose location in the phase diagram can be accurately determined.Our results can help to develop more detailed models of complex networks with correlated degrees.
I. INTRODUCTION
Random graphs constitute the main tool to model the com-plex behavior of large empirical networks observed in social,technological, and biological systems [1, 2]. In random graphmodels a network is typically represented by nodes that inter-act through edges. Random graph theory leads to importantinsights into the structure of networks as well as on the dy-namical processes occurring on them, such as the spreadingof diseases [3, 4], the stability of ecosystems to perturbations[5, 6], and the dynamics of sparsely connected neurons [7, 8].A fruitful approach to network modeling consists in mea-suring a set of observables in an empirical network and thenbuilding an ensemble of random graphs that matches thesefeatures in an average sense [9, 10]. The probability distribu-tion of graph configurations is derived by maximizing the net-work entropy subject to the constraints dictated by the empir-ical observations [11]. The resulting family of models, knownas exponential random graph (ERG) models, aims to repro-duce a set of empirical patterns while keeping other networkproperties entirely random. ERGs were introduced in the pi-oneering work of Holland and Leinhardt [12], and they soonbecame popular models in social network analysis [13–17].There are at least two main motivations to study ERGs.First, they serve as benchmark models to distinguish betweenrandom and non-random traits in the structure of real-worldnetworks [18–20]. Analytic solutions of ERGs give detailedinformation on the expected values of structural observablesand their fluctuations, which can be directly compared with ∗ [email protected] † [email protected] data from real-world networks. Second, ERG models mayexhibit degenerate configurations and phase transitions, i.e,abrupt changes in the macroscopic properties of the graphensemble. Discontinuous phase transitions can be a seriouslimitation in the generation of ERGs, since they prevent thatcertain configurations, with the desired structural features, aresampled. Analytic solutions of ERGs predict the existenceand location of phase transitions in the parameter space.ERG models with specific constraints have been widelystudied through numeric [21–23] and analytic [10, 11, 24–32]techniques from statistical mechanics. The two-star model isprobably the simplest of ERGs that undergoes a phase tran-sition [11, 25]. In this model, the graph ensemble is con-strained by the average number of edges and the average num-ber of two-stars (a two-star is a pair of edges that share a com-mon node). The two-star model has been originally solved inthe high-connectivity regime [11, 25], using mean-field tech-niques, and more recently in the sparse regime [29], whereeach node is adjacent to a finite number of others.Here we take the theory of ERGs one step further by solvingthe two-star model with degree-degree correlations [33, 34],which is an important structural feature of networks. In gen-eral, the degrees in real-world networks are not independent,but they can be positively or negatively correlated with eachother, as quantified by the Pearson correlation coefficient [33].Nodes with similar degrees have positive degree correlations,whereas nodes with highly distinct degrees have negative de-gree correlations. For instance, the degrees of adjacent ornearest neighbor nodes in social networks are positively cor-related [33, 34], but these correlations become negative forpairs of nodes connected through paths with more than oneedge [35, 36].Earlier works have focused on nearest neighbor degree cor-relations [33, 34] and their impact on dynamical processes a r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b on networks, such as the spreading of diseases [37, 38] andthe synchronization of coupled oscillators [39–41]. However,models that only take into account local properties [42, 43]do not reproduce certain global features of networks, such astheir community structure or the distribution of the shortestpath length. In fact, recent works [35, 36, 44–47] have shownthat long-range degree correlations between nodes separatedby more than one edge are important for the organization ofnetworks at a global level. For instance, results suggest thatthe fractal structure of scale-free networks requires long-rangedegree correlations between hubs [43]. Long-range degreecorrelations have been also observed in the airport transporta-tion network of the United States [35], transcriptional regula-tory networks [35], coauthorship networks [36], and the Twit-ter network [47]. Therefore, the solution of ERG models thatincorporate degree correlations in a systematic way representsa significant progress in network modeling.In this work we solve the two-star model with degree-degree correlations between nearest neighbors and betweennext nearest neighbors in the sparse regime. The free energy isexactly calculated thanks to the introduction of an upper cut-off in the degree sequence. We show that the phase diagram ofthe model exhibits a first-order critical line, surrounded by ametastable region, in which the graph sampling process mayget stuck in a local minimum of the free energy. The first-order transition separates a phase characterized by an approx-imate Poisson degree distribution from a condensed phase,where the degree distribution strongly depends on the degreecorrelations. We quantify the degree correlations through thedegree assortativity corresponding to nearest neighbor nodes[33, 34] and to next nearest neighbor nodes, located at theend-points of two-stars. When the degree assortativity of nextnearest neighbors is non-negative, the degree distribution inthe condensed phase is peaked at the maximum degree; whenthe assortativity of next nearest neighbors is negative, the con-densed phase is characterized by a bimodal degree distribu-tion. Both assortativities are non-monotonic functions of themodel parameters and exhibit a non-analytic behavior at thefirst-order transition. The main theoretical findings are wellcorroborated by Monte Carlo simulations.In the next section we introduce the generic framework ofERG models. In section III we define the main structural prop-erties of interest, including the two assortativity coefficients.Section IV explains how the model is analytically solved usingconventional techniques of statistical mechanics, and how thestructural properties follow from the free energy. The numer-ical results, obtained from the solutions of the saddle-pointequations, are discussed in section V, while in the last sectionwe present some final remarks. The appendix discusses thesymmetry properties of the order-parameter functions. II. EXPONENTIAL RANDOM GRAPH MODELS
A graph configuration of a network with N nodes can berepresented by a realization of the N × N adjacency matrix C . The entries of C fully encode the network topology, i.e.,the matrix element C ij tells whether there is an edge joining nodes i and j . We consider undirected and simple randomgraphs [1], which means that C is a symmetric matrix withall diagonal elements equal to zero. If there is an edge con-necting nodes i and j , then we set C ij = 1 , whereas C ij = 0 otherwise. The degree K i of a node iK i = N (cid:88) j =1 C ij (1)counts the number of edges attached to i , and the sequence K , . . . , K N contains important information about the net-work structure. In this work we consider random graph mod-els in which the maximum degree that may appear in a graphconfiguration is k max . The cutoff k max is a model parameter,independent of N , which can be freely adjusted. As we willsee below, the introduction of k max allows to compute exactlythe network properties in the limit N → ∞ .The probability P N ( C ) to observe a certain graph configu-ration C follows the Boltzmann-like form [11] P N ( C ) = e −H N ( C ) Z N N (cid:89) i =1 Θ( k max − K i ) , (2)where Θ( x ) = 1 if x ≥ , and Θ( x ) = 0 otherwise.The graph Hamiltonian H N ( C ) depends on the network con-straints, and Z N is the graph partition function Z N = (cid:88) C e −H N ( C ) N (cid:89) i =1 Θ( k max − K i ) . (3)The sum (cid:80) C runs over all possible realizations of the ad-jacency matrix. To study the stability of the graph configura-tions for N → ∞ , we need to compute the free energy density f = − lim N →∞ N ln Z N , (4)which plays the role of a generating function for the graphstructural properties. III. NETWORK OBSERVABLES
The topology of graphs sampled from P N ( C ) can be char-acterized by a set of structural observables. In this work weonly consider global observables, which are obtained by aver-aging a local quantity over the entire network.An important quantity to probe the network structure is theempirical degree distribution p k ( C ) = 1 N N (cid:88) i =1 δ k,K i , (5)which gives the probability that a randomly chosen node hasdegree k . The density of edges (cid:96) ( C ) and the density of two-stars s ( C ) are given by [11, 29] (cid:96) ( C ) = 12 N N (cid:88) ij =1 C ij = 12 N N (cid:88) i =1 K i , (6) s ( C ) = 12 N N (cid:88) ijn =1 (1 − δ in ) C ij C jn = 12 N N (cid:88) i =1 (cid:0) K i − K i (cid:1) . (7)A two-star (or path of length two) is a set with three differentnodes { i, j, k } such that C ij C jk = 1 .Degree correlations are commonly quantified by the degreeassortativity coefficient [33–36]. This is a global observabledefined as the Pearson correlation coefficient between the de-grees of two nodes. Here we characterize the degree-degreecorrelations by means of two assortativity parameters: thestandard assortativity A (1) ( C ) measures the degree correla-tions between adjacent nodes, while the assortativity A (2) ( C ) measures the degree correlations between nodes at the end-points of two-stars. In other words, A (1) ( A (2) ) quantifiesdegree correlations between nearest neighbors (next nearestneighbors).In what follows, the indexes k and l refer to degrees. For asingle graph instance, the assortativities are defined as A ( r ) ( C ) = (cid:80) ∞ kl =0 kl W ( r ) k,l ( C ) − (cid:104)(cid:80) ∞ k =0 kW ( r ) k ( C ) (cid:105) (cid:80) ∞ k =0 k W ( r ) k ( C ) − (cid:104)(cid:80) ∞ k =0 kW ( r ) k ( C ) (cid:105) , (8)with r = 1 , . The quantity W (1) k,l ( C ) = (cid:80) Nij =1 C ij δ k,K i δ l,K j (cid:80) Nij =1 C ij (9)is the probability that a randomly chosen edge joins two nodeswith degrees k and l , while W (2) k,l ( C ) = (cid:80) Nijn =1 (1 − δ in ) C ij C jn δ k,K i δ l,K n (cid:80) Nijn =1 (1 − δ in ) C ij C jn (10)is the probability that a randomly chosen two-star has degrees k and l at its end-points. The marginal distributions W ( r ) k ( C ) = ∞ (cid:88) l =0 W ( r ) k,l ( C ) ( r = 1 , (11)have the explicit forms W (1) k ( C ) = 12 N (cid:96) ( C ) N (cid:88) i =1 K i δ k,K i , (12) W (2) k ( C ) = 12 N s ( C ) N (cid:88) i =1 N (cid:88) j =1 C ij K j − K i δ k,K i . (13) Equation (12) shows that the contribution of a node to W (1) k isweighted according to its degree, while the weight of node i to the distribution W (2) k is determined by the number of edgesattached to the neighbors of i , except from the links comingfrom i itself. This is intuitive, as a node i with a large second-order degree (cid:80) Nj =1 C ij K j − K i is the end-point of a largenumber of two-stars.The assortativities A (1) and A (2) give the same type of sta-tistical information. Networks with statistically independentdegrees satisfy W ( r ) k,l ( C ) = W ( r ) k ( C ) W ( r ) l ( C ) , (14)and, consequently, A ( r ) ( C ) = 0 . Networks with A ( r ) ( C ) > are positively correlated or assortative, which means thatnodes connected through edges or two-stars are likely to havesimilar degrees. Finally, networks with A ( r ) ( C ) < are neg-atively correlated or disassortative, meaning that nodes withlarge degrees preferentially connect through edges or two-stars to nodes with small degrees.Equation (8) is not practical to calculate A ( r ) . To preparethe ground for the analytic computation of the assortativi-ties, let us derive more convenient expressions for A (1) and A (2) . Substituting Eqs. (9), (10), (12), and (13) in Eq. (8), werewrite A (1) and A (2) as follows A (1) ( C ) = Λ ( C ) − (cid:16) N(cid:96) ( C ) (cid:80) Ni =1 K i (cid:17) N(cid:96) ( C ) (cid:80) Ni =1 K i − (cid:16) N(cid:96) ( C ) (cid:80) Ni =1 K i (cid:17) , (15) A (2) ( C ) = χ ( C ) − [Σ ( C )] (cid:96) ( C ) s ( C ) Λ ( C ) − Ns ( C ) (cid:80) Ni =1 K i − [Σ ( C )] , (16)with Σ ( C ) = (cid:96) ( C ) s ( C ) Λ ( C ) − N s ( C ) N (cid:88) i =1 K i . (17)The object Λ qr ( C ) defines higher-order moments of nearestneighbor degrees Λ qr ( C ) = 12 N (cid:96) ( C ) N (cid:88) ij =1 C ij K qi K rj ( q, r ≥ , (18)and χ ( C ) is the correlation between next nearest neighbor de-grees χ ( C ) = 12 N s ( C ) N (cid:88) ijn =1 (1 − δ in ) C ij C jn K i K n . (19)Equations (15) and (16) hold for a single realization of C andthey show that A (1) and A (2) are ultimately given in terms ofmoments of the joint distribution of degrees at different pairsof nodes.In the limit N → ∞ , the fluctuations of intensive variablesvanish and a single realization of an intensive quantity coin-cides with its ensemble averaged value. Thus, we naturally as-sume that the assortativities and all other global observables ofinterest display such self-averaging behavior when N → ∞ ,and the assortativities become (cid:104) A (1) (cid:105) = (cid:104) Λ (cid:105) − (cid:16) (cid:104) K (cid:105) (cid:104) (cid:96) (cid:105) (cid:17) (cid:104) K (cid:105) (cid:104) (cid:96) (cid:105) − (cid:16) (cid:104) K (cid:105) (cid:104) (cid:96) (cid:105) (cid:17) , (20) (cid:104) A (2) (cid:105) = (cid:104) χ (cid:105) − (cid:104) (cid:104) (cid:96) (cid:105)(cid:104) s (cid:105) (cid:104) Λ (cid:105) − (cid:104) K (cid:105) (cid:104) s (cid:105) (cid:105) (cid:104) (cid:96) (cid:105)(cid:104) s (cid:105) (cid:104) Λ (cid:105) − (cid:104) K (cid:105) (cid:104) s (cid:105) − (cid:104) (cid:104) (cid:96) (cid:105)(cid:104) s (cid:105) (cid:104) Λ (cid:105) − (cid:104) K (cid:105) (cid:104) s (cid:105) (cid:105) , (21)where (cid:104)G(cid:105) denotes the ensemble average of an arbitrary ran-dom function G ( C ) for N → ∞(cid:104)G(cid:105) = lim N →∞ (cid:88) C G ( C ) P N ( C ) . (22)All ensemble averages in Eqs. (20) and (21) can be calculatedfrom the free energy f , which works as a generating functionfor the moments of degrees. IV. ANALYTIC SOLUTION OF THE TWO-STAR MODELWITH CORRELATED DEGREES
In this section we present the Hamiltonian of the two-starmodel with correlated degrees and we explain how to calcu-late, in the limit N → ∞ , the free energy f and the structuralobservables using standard tools from statistical mechanics. A. The Hamiltonian of the model
The ERG model is defined by the Hamiltonian H ( C ) = − Q (cid:88) r =1 α r N (cid:88) i =1 F r ( K i ) − γ N (cid:88) ij =1 C ij D ( K i , K j ) − β N (cid:88) ijk =1 (1 − δ ik ) C ij C jk K i K k + ln N (cid:88) i In this subsection we solve the model defined by Eq. (23).The aim is to calculate the free energy f in the limit N →∞ , from which ensemble averages of the network observablesreadily follow.The graph partition function reads Z N = (cid:89) i In the limit N → ∞ , the ensemble averages of the networkobservables, defined in section III, follow from the deriva-tives of f with respect to the model parameters. Let us il-lustrate this fact by deriving the analytic expression for thedegree distribution. If we set F ( l ) = δ k,l for arbitrary in-tegers ≤ k, l ≤ k max , then the ensemble average degreedistribution (cid:104) p k (cid:105) is determined from (cid:104) p k (cid:105) = − ∂f∂α , (40)which follows from Eqs. (2) and (4). From the explicit formof f , Eqs. (36) and (37), we get (cid:104) p k (cid:105) = k ! (cid:82) ∞−∞ Dx [ˆ ρ ∗ k ( x )] k e (cid:80) Qr =1 α r F r ( k ) (cid:80) k max l =0 1 l ! (cid:82) ∞−∞ Dx [ˆ ρ ∗ l ( x )] l e (cid:80) Qr =1 α r F r ( l ) . (41)This is a common strategy to calculate ensemble averages instatistical mechanics, namely, one performs the derivative ofthe free energy with respect to a parameter that is coupled toa certain observable in the Hamiltonian. Note that F ( k ) inEq. (41) is not necessarily given by F ( k ) = δ k,l . In otherwords, after the choice F ( k ) = δ k,l has served the purposeto obtain an expression for (cid:104) p k (cid:105) , we are free to choose F ( k ) as we please.Following an analogous procedure, the equations for the en-semble averages of all other observables introduced in sectionIII are obtained in the limit N → ∞ . The density of links andthe density of two-stars read (cid:104) (cid:96) (cid:105) = (cid:104) K (cid:105) , (42) (cid:104) s (cid:105) = 12 (cid:0) (cid:104) K (cid:105) − (cid:104) K (cid:105) (cid:1) , (43) where the moments (cid:104) K n (cid:105) ( n = 1 , , . . . ) of (cid:104) p k (cid:105) are deter-mined from (cid:104) K n (cid:105) = k max (cid:88) k =0 k n (cid:104) p k (cid:105) . (44)The moments (cid:104) Λ qr (cid:105) of the joint degree distribution at adjacentnodes read (cid:104) Λ qr (cid:105) = 12 (cid:104) (cid:96) (cid:105) k max (cid:88) kl =0 k q l r (cid:90) dxdx (cid:48) ρ ∗ k ( x ) ρ ∗ l ( x (cid:48) ) e W γ,β ( k,x ; l,x (cid:48) ) , (45)and the average degree correlation (cid:104) χ (cid:105) at the end-points oftwo-stars is given by (cid:104) χ (cid:105) = 12 (cid:104) s (cid:105) k max (cid:88) kl =0 (cid:90) dxdx (cid:48) ρ ∗ k ( x ) ρ ∗ l ( x (cid:48) ) × (cid:18) √ β xl − k (cid:19) e W γ,β ( k,x ; l,x (cid:48) ) . (46)Once we determine { ρ ∗ k ( x ) , ˆ ρ ∗ k ( x ) } from the solutions of thesaddle-point Eqs. (38) and (39), Eqs. (41-46) allow to com-pute the assortativities and characterize the network structurein the limit N → ∞ . In case Eqs. (38) and (39) have morethan a single solution, the structural observables are calculatedfrom the solution that corresponds to the global minimum of F [ ρ k , ˆ ρ k ] . V. RESULTS The exact equations derived in the previous section describeERGs with the generic Hamiltonian of Eq. (23) in the limit N → ∞ . In this section we solve these equations and studythe effect of degree correlations in the phase diagram of thetwo-star model. A. The two-star model without degree correlations In this section we present results for the two-star model inthe absence of degree correlations [9–11, 25, 29], where theaverage density of edges and the average density of two-starsare the only constraints. The model undergoes a discontinuoustransition as a function of the control parameters both in thedense regime [11, 25] and in the more realistic case of sparsenetworks [29]. In the latter case, reference [29] reports a dis-continuous behavior in the structural parameters, but the sta-bility of the macroscopic states and the corresponding phasediagram remain elusive. We complement the work of [29] byconstructing the full phase diagram of the two-star model inthe sparse regime.The Hamiltonian of the two-star model H ( C ) = − α N (cid:88) i =1 K i − α N (cid:88) i =1 K i + ln N (cid:88) i 76 78 8000 . α = 0 . k h p k i . . α µ FIG. 2. (a) Phase diagram of the two-star model with maximum de-gree k max = 80 (see Eq. (47)). The dashed black curve is the first-order critical line and the solid red curves delimit the metastable re-gion, within which the free energy has two minima. The inset showsthe two stable solutions of Eq. (48) as a function of α along thedashed curve. The two solutions for µ merge continuously at a crit-ical point, identified by the black dot in figure (a). Figures (b) and(c) show the degree distribution (green squares) for fixed α = 0 . , k max = 80 , and a value of α inside each phase. The blue crosses infigure (b) denote a Poisson distribution with the same mean degree. The behavior of the structural properties across the phasetransition is a subject of practical interest. Figure 3 showsthe density of links and the density of two-stars as a functionof α for different α . Both quantities are discontinuous atthe first-order transition. The discontinuity becomes gradu-ally smaller as we increase α , until it finally disappears atthe critical point, i.e., (cid:104) (cid:96) (cid:105) and (cid:104) s (cid:105) are continuous and mono-tonic functions of α provided α (cid:38) . . Since β = γ = 0 ,both assortativities are zero in this model. The theoretical re-sults of figure 3 are well corroborated by data obtained fromMonte Carlo simulations that sample graphs from the distri-bution of Eq. (2). Monte Carlo methods to generate ERGs arethoroughly discussed in [9]. . 004 0 . 008 0 . α h l i α = 1 . α = 1 . α = 1 . α = 1 . . 004 0 . 008 0 . , , , α h s i FIG. 3. Theoretical results (different line styles) for the average den-sity of edges (cid:104) (cid:96) (cid:105) and the average density of two-stars (cid:104) s (cid:105) as a functionof α for the two-star model (see Eq. (47)) with maximum degree k max = 80 and different α . The symbols are obtained from graphsamples generated through Monte Carlo simulations with α = 1 . and two system sizes: N = 200 (squares) and N = 800 (circles). B. Degree correlations between nearest neighbors In this subsection we analyze the role of nearest neigh-bor degree correlations on the phase diagram of the two-starmodel. We consider an ERG model that allows to tune thedensity of links and the degree correlations between adjacentnodes. The model is defined by the Hamiltonian H ( C ) = − α N (cid:88) i =1 K i − γ N (cid:88) i 001 0 . − − γf . 002 0 . − − γ FIG. 4. Free energy f as a function of γ for the two-star model withnearest neighbor degree correlations (see Eq. (51)), maximum degree k max = 30 , and α = 0 . . The non-analytic point of f marks thefirst-order critical point. The inset displays the two branches of thefree energy, each one obtained from a stable solution of the fixed-point Eq. (54). The metastable region is delimited by the dotted ver-tical lines. − . . . . . (a) α γ MetastableFirst-order FIG. 5. Phase diagram of the two-star model with nearest neighborcorrelated degrees (see Eq. (51)) and maximum degree k max = 30 .The dashed black curve marks the first-order critical line and thesolid red curves delimit the metastable region, within which the freeenergy has two minima. The phase diagram of the ERG model defined by Eq. (51),and obtained from the analysis of f , is shown in figure 5.The phase diagram has a first-order critical line, surroundedby a metastable region, which ends at a critical point. Thefirst-order critical line marks anew the condensation transi-tion, above which the degree distribution (cid:104) p k (cid:105) has a peak at k = k max . The profile of (cid:104) p k (cid:105) below and above the dashedline in figure 5 is qualitatively similar to, respectively, figures2-(b) and 2-(c). For γ > , adjacent nodes tend to have similardegrees, which strongly favors the formation of a regular ran-dom graph, driving the ERG model to the condensed phase.We do not observe a phase transition for γ < .The nearest neighbor assortativity (cid:104) A (1) (cid:105) is discontinuousat the first-order transition. Figure 6 shows (cid:104) A (1) (cid:105) as a func-tion of γ for different α . The discontinuity of (cid:104) A (1) (cid:105) becomessmaller as α increases, until it vanishes continuously at thecritical point located at α (cid:39) . . In fact, the two solutions − − − . − . − . (a) γ h A (1) i α = 0 . α = 0 . α = 1 . α = 1 . 20 0 . 001 0 . 002 0 . . . . . (b) γ h A (1) i FIG. 6. Theoretical results (different line styles) for the degree as-sortativity (cid:104) A (1) (cid:105) of adjacent nodes as a function of γ for the two-star model with nearest neighbor correlated degrees (see Eq. (51))and maximum degree k max = 30 . The symbols are results obtainedfrom sampling graphs using Monte Carlo simulations with α = 1 . , N = 1500 (circles), and N = 5000 (squares). for (cid:104) A (1) (cid:105) corresponding to each phase merge into a singlevalue as we approach the critical point along the dashed curvein figure 5. The assortativity (cid:104) A (1) (cid:105) is a non-monotonic func-tion of γ that vanishes when | γ | → ∞ . The latter propertycan be understood from the ground state configurations of theHamiltonian. For γ → −∞ , the minimum of Eq. (51) is at-tained when all degrees are zero; for γ → ∞ , Eq. (51) is min-imized for all degrees equal to k max . Since β = 0 , nodes atthe end-points of two-stars are uncorrelated and (cid:104) A (2) (cid:105) = 0 .Figure 6 also shows results generated through Monte Carlosimulations of finite random graphs, which confirms our the-oretical findings for N → ∞ . C. Degree correlations between next nearest neighbors In this subsection we focus on the competition betweennearest neighbor degree correlations and next nearest neigh- bor degree correlations. We consider the Hamiltonian H ( C ) = ln N (cid:88) i 006 0 0 . 006 0 . − . . . (a) γ h A (1) i β = 0 . β = 0 . β = 0 . − . − . 006 0 0 . 006 0 . . . . (b) γ h A (2) i FIG. 7. Degree assortativities (cid:104) A (1) (cid:105) and (cid:104) A (2) (cid:105) corresponding,respectively, to nearest neighbor nodes and next nearest neighborsnodes of the two-star model with degree-degree correlations (see Eq.(56)), β > , and maximum degree k max = 10 . First, we consider the ERG model of Eq. (56) for β > .In this regime the order-parameters ρ ∗ ( x ) , . . . , ρ ∗ k max ( x ) arereal-valued functions. Figure 7 shows that the degree assorta-tivities (cid:104) A (1) (cid:105) and (cid:104) A (2) (cid:105) have a discontinuous behavior as afunction of γ , and the model undergoes once more a first-ordertransition at γ = γ c ( β ) . For γ < γ c ( β ) , (cid:104) p k (cid:105) is approximatelygiven by a Poisson distribution, while for γ > γ c ( β ) the de-gree distribution has a single peak at k = k max . Figure 8 illus-trates the behavior of (cid:104) p k (cid:105) and ρ ∗ ( x ) , . . . , ρ ∗ k max ( x ) inside the0condensed phase. For γ → ∞ or β → ∞ , the degree distribu-tion converges to (cid:104) p k (cid:105) = δ k,k max and both assortativities arezero. Overall, figures 7 and 8 show that positive next nearestneighbor degree correlations do not change qualitatively thephase diagram.Let us now present results for β < , where ρ ∗ ( x ) , . . . , ρ ∗ k max ( x ) are complex-valued functions of x ∈ R .In the appendix, we demonstrate that Re ρ ∗ k ( x ) is an even func-tion and Im ρ ∗ k ( x ) is an odd function in the regime β < . Itis interesting to note that, for β < , we impose conflictingconstraints in the generation of graph samples, since negativevalues of β favor dissimilar degrees at the end-points of two-stars, leaving the degrees at the central nodes of two-stars ina frustrated situation. The appearance of these frustrated con-figurations should influence the graph structure. . (a) k h p k i γ = − . γ = − . 004 2 4 6 800 . . (b) xρ ∗ k ( x ) FIG. 8. Degree distribution (cid:104) p k (cid:105) and the order-parameter functions { ρ ∗ k ( x ) } of the two-star model with degree-degree correlations (seeEq. (56)), β = 0 . , maximum degree k max = 10 , and differentvalues of γ in the condensed phase. Figure (b) shows the functionalbehavior of ρ ∗ ( x ) (black solid line), ρ ∗ ( x ) (red dashed line), and ρ ∗ ( x ) (blue dot-dashed line) for γ = − . . The other components ρ ∗ ( x ) , . . . , ρ ∗ ( x ) are approximately zero. In figure 9 we present (cid:104) A (1) (cid:105) and (cid:104) A (2) (cid:105) as a function of γ for β < . The ERG model undergoes a first-order tran-sition at the critical point γ = γ c ( β ) . For γ < γ c ( β ) , thedegree correlations do not considerably affect (cid:104) p k (cid:105) , which iscloser to a Poisson distribution, similar to figure 2-(b). For γ > γ c ( β ) , negative degree correlations between next near-est neighbors have an important effect in the graph structureand (cid:104) p k (cid:105) can exhibit a bimodal shape, as illustrated in figure10. For increasing γ > γ c ( β ) , the weight (cid:104) p k max (cid:105) graduallyincreases until we attain (cid:104) p k (cid:105) = δ k,k max for γ → ∞ . Never-theless, the graph structure in the condensed phase for β < is qualitatively distinct from the regime β > , which is alsoattested by the functional behavior of the order-parameters,presented in figures 10-(b) and 10-(c) (as a comparison, seefigure 8-(b)). As shown by figure 9, the first-order transitiondisappears for β smaller than a certain threshold, which marksthe terminating point of the first-order critical line in the plane ( β, γ ) .The competition between nearest neighbor and next near-est neighbor degree correlations is summarized in figure 11,which depicts the phase diagram of the model defined by Eq.(56). The phase diagram exhibits a metastable region arounda first-order critical line γ c ( β ) that terminates at a negativevalue of β . For γ < γ c ( β ) , (cid:104) p k (cid:105) is closer to a Poisson distri-bution. For γ > γ c ( β ) and β > , we have (cid:104) A (2) (cid:105) > and . 02 0 . . . . (a) γ h A (1) i β = − . β = − . β = − . β = − . . 02 0 . − . − . − . (b) γ h A (2) i FIG. 9. Degree assortativities (cid:104) A (1) (cid:105) and (cid:104) A (2) (cid:105) corresponding,respectively, to nearest neighbor nodes and next nearest neighborsnodes of the two-star model with degree-degree correlations (see Eq.(56)), β < , and maximum degree k max = 10 . . . (a) k h p k i γ = 0 . γ = 0 . γ = 0 . − (b) − (c) x FIG. 10. Degree distribution (cid:104) p k (cid:105) and the order-parameter functions { ρ ∗ k ( x ) } of the two-star model with degree-degree correlations (seeEq. (56)), β = − . , maximum degree k max = 10 , and differentvalues of γ above the first-order transition. Figures (b) and (c) show,respectively, the functional behavior of Im ρ ∗ k ( x ) and Re ρ ∗ k ( x ) for γ = 0 . and three values of k : k = 10 (black solid lines), k = 9 (red dashed lines), and k = 8 (blue dot-dashed lines). The othercomponents ρ ∗ ( x ) , . . . , ρ ∗ ( x ) are approximately zero. the degree distribution has a single peak at k = k max . For γ > γ c ( β ) and β < , we have (cid:104) A (2) (cid:105) < and (cid:104) p k (cid:105) canexhibit two peaks, one of them located at k = k max , and anadditional peak at a smaller degree. VI. FINAL REMARKS In this paper we have solved the two-star model withdegree-degree correlations in the sparse regime. The modelallows to generate random graphs with prescribed degree cor-relations between adjacent nodes and between nodes at theend-points of two-stars. By introducing an upper cutoff in the1 − . 002 0 0 . . . 04 Bimodal β < β > βγ MetastableFirst-order FIG. 11. Phase diagram of the two-star model with degree-degreecorrelations (see Eq. (56)) and maximum degree k max = 10 . Thesymbols are theoretical results obtained from the solutions of Eqs.(37-39), and the solid lines are just a guide to the eye. The phasediagram has a first-order critical line surrounded by a metastable re-gion, where the free energy has two minima. Above the first-ordercritical line, the degree distribution is bimodal for β < , and it hasa single peak at the largest degree for β > . The degree distributionapproaches (cid:104) p k (cid:105) = δ k,k max as γ → ∞ . degree sequence, we have exactly calculated the network freeenergy, from which we derived complete phase diagrams andcharacterized the graph structure in the different phases.In terms of the degree distribution (cid:104) p k (cid:105) , the phase diagramof the model is characterized by three distinct regions. Thereis a phase where (cid:104) p k (cid:105) is approximately given by a Poisson dis-tribution, reminiscent from the structure of Erd¨os-R´enyi ran-dom graphs [1]. The phase diagram also exhibits a condensedphase, where the shape of (cid:104) p k (cid:105) strongly depends on the degreecorrelations. If the degree assortativities are non-negative in-side the condensed phase, then (cid:104) p k (cid:105) has a single peak at themaximum degree and the graph is approximately regular. Ifthe degree assortativity of next nearest neighbors is negativeinside the condensed phase, then (cid:104) p k (cid:105) is given by a bimodaldistribution, with one maximum at the largest degree and anadditional maximum at a smaller degree. While the Poissonand standard condensed phases appear even in the absence ofdegree-degree correlations, the existence of a bimodal degreedistribution is a genuine effect of negative degree correlationsbetween next nearest neighbor nodes. This result reveals theimportance of long-range degree correlations, beyond nearestneighbor nodes, in shaping the network structure.We have shown that the model undergoes a first-order tran-sition between the Poisson phase and the condensed phase.The first-order critical line is surrounded by a metastable re-gion, where the free energy has two minima, each one cor-responding to a phase. For combinations of model parame-ters inside the metastable region, algorithms to sample finitegraphs from this ERG model may get stuck in a local mini-mum of the free energy [9]. In addition, the jump of the struc-tural observables across the first-order critical line prevents usfrom generating graph samples with structural parameters ina certain range. Taken together, these features represent seri-ous limitations of the present ERG model as an effective toolto model real-world networks. The analytic solution of themodel for N → ∞ and the construction of its phase diagrams have a practical relevance, as these results allow to estimatethe metastable regions in the parameter space of finite graphs.The present paper constitutes a first step towards control-ling the generation of ERGs with correlated degrees. Overall,the results for the assortativities, the degree distribution, andthe phase diagrams allow to identify the regime of parameterswhere the model can be useful to reproduce certain propertiesof empirical networks. A drawback of the present model isthat the degree distribution in each phase does not bear anyresemblance to the broad degree distributions found in real-world networks. With the purpose of improving the model, itwould be interesting to solve it with a hard constraint in thedegree sequence or with a prescribed degree distribution [10].This work also opens the perspective to explore systematicallythe role of short-range and long-range degree correlations indynamical processes occurring on tree-like networks, since inthis case the equations for the dynamics are typically deter-mined only by the degree distribution [49]. Finally, we pointout that the free energy of ERG models can be mapped inthe cumulant generating function of certain structural observ-ables of Erd¨os-R´enyi random graphs [50–52] . Therefore, theresults of the present paper can be readily applied to studyanalytically the large deviations of higher-order topologicalproperties of Erd¨os-R´enyi random graphs in the limit N → ∞ [53]. ACKNOWLEDGEMENTS M.B. acknowledges a fellowship from CAPES/Brazil.F.L.M. and I.P.C. gratefully acknowledge London Mathemat-ical Laboratory for financial support. F.L.M. also acknowl-edges a fellowship from CNPq/Brazil Appendix A: Symmetry properties of the order-parameterfunctions for β < In this appendix we obtain the symmetry properties of theorder-parameters ρ ∗ ( x ) , . . . , ρ ∗ k max ( x ) under the transforma-tion x → − x . These properties allow to simplify the saddle-point Eqs. (38) and (39) and the computation of the structuralobservables introduced in section III.For β < , { ρ ∗ k ( x ) , ˆ ρ ∗ k ( x ) } are complex-valued functionsof x ∈ R . This can be seen from Eq. (38), which can bewritten for β < as ˆ ρ ∗ k ( x ) = k max (cid:88) l =0 U kl e i √ | β | xl (cid:90) ∞−∞ dx (cid:48) ρ ∗ l ( x (cid:48) ) e i √ | β | x (cid:48) k , (A1)where U kl = e γD ( k,l )+ | β | ( k + l ) . (A2)If we take the complex-conjugate ( . . . ) of Eq. (39) and makethe transformation x → − x , then we get ρ ∗ k ( − x ) = 1 T kk ! (cid:2) ˆ ρ ∗ k ( − x ) (cid:3) k − e − x + (cid:80) Qr =1 α r F r ( k ) , (A3)2with TT = k max (cid:88) k =0 k ! (cid:90) ∞−∞ dx (cid:2) ˆ ρ ∗ k ( − x ) (cid:3) k e − x + (cid:80) Qr =1 α r F r ( k ) . (A4)By taking the complex-conjugate of Eq. (A1), the function ˆ ρ ∗ k ( − x ) fulfills ˆ ρ ∗ k ( − x ) = k max (cid:88) l =0 U kl e i √ | β | xl (cid:90) ∞−∞ dx (cid:48) ρ ∗ l ( − x (cid:48) ) e i √ | β | x (cid:48) k . (A5)Since Eqs. (38) and (39) for { ρ ∗ k ( x ) , ˆ ρ ∗ k ( x ) } are the same asEqs. (A3) and (A5) for { ρ ∗ k ( − x ) , ˆ ρ ∗ k ( − x ) } , we conclude that ρ ∗ k ( x ) = ρ ∗ k ( − x ) (A6)and ˆ ρ ∗ k ( x ) = ˆ ρ ∗ k ( − x ) (A7)for arbitrary x . This implies that { Re ρ ∗ k ( x ) , Reˆ ρ ∗ k ( x ) } and { Im ρ ∗ k ( x ) , Imˆ ρ ∗ k ( x ) } are, respectively, even and odd func-tions of x .Let us use the above symmetry properties to simplify theorder-parameter equations. By setting ˆ ρ ∗ k ( x ) in polar form ˆ ρ ∗ k ( x ) = r k ( x ) e iϕ k ( x ) , (A8)with ϕ k ( x ) ∈ ( − π, π ] , and noting that r k ( − x ) = r k ( x ) (A9) and ϕ k ( − x ) = − ϕ k ( x ) , (A10)we readily obtain Im T = 0 from Eq. (A4). Hence Eq. (39)can be written as Re ρ ∗ k ( x ) = 1Re T kk ! [ r k ( x )] k − cos [( k − ϕ k ( x )] × e − x + (cid:80) Qr =1 α r F r ( k ) , (A11) Im ρ ∗ k ( x ) = 1Re T kk ! [ r k ( x )] k − sin [( k − ϕ k ( x )] × e − x + (cid:80) Qr =1 α r F r ( k ) . (A12)The above equations are coupled to Eq. (38), which can besimplified for β < using Eqs. (A6) and (A7) Reˆ ρ ∗ k ( x ) = 2 k max (cid:88) l =0 U kl cos (cid:16)(cid:112) | β | xl (cid:17) (cid:90) ∞ dx (cid:48) Y kl ( x (cid:48) ) , (A13) Imˆ ρ ∗ k ( x ) = 2 k max (cid:88) l =0 U kl sin (cid:16)(cid:112) | β | xl (cid:17) (cid:90) ∞ dx (cid:48) Y kl ( x (cid:48) ) , (A14)where Y kl ( x ) = Re ρ ∗ l ( x (cid:48) ) cos (cid:16)(cid:112) | β | x (cid:48) k (cid:17) − Im ρ ∗ l ( x (cid:48) ) sin (cid:16)(cid:112) | β | x (cid:48) k (cid:17) . (A15)The fixed-point functions { ρ ∗ k ( x ) } that solve Eqs. (A11-A14)determine the free energy f and all the structural parametersfor β < . Using the symmetry properties of Eqs. (A6) and(A7), it is straightforward to verify from Eq. (36) that f ∈ R for β < . [1] M. Newman, Networks: An Introduction (OUP Oxford, 2010).[2] S. N. Dorogovtsev and J. F. Mendes, Evolution of networks:From biological nets to the Internet and WWW (OUP Oxford,2013).[3] P. V. Mieghem, EPL (Europhysics Letters) , 48004 (2012).[4] A. V. Goltsev, S. N. Dorogovtsev, J. G. Oliveira, and J. F. F.Mendes, Phys. Rev. Lett. , 128702 (2012).[5] J. Grilli, T. Rogers, and S. Allesina, Nat. Commun. , 12031(2016).[6] I. Neri and F. L. Metz, Phys. Rev. Research , 033313 (2020).[7] N. Brunel, Comput. Neurosci. , 183 (2000).[8] S. Ostojic, Nat. Neurosci. , 594 (2014).[9] A. Coolen, A. Annibale, and E. Roberts, Generating RandomNetworks and Graphs (Oxford University Press, 2017).[10] G. Cimini, T. Squartini, F. Saracco, D. Garlaschelli,A. Gabrielli, and G. Caldarelli, Nat. Rev. Phys. , 58 (2019).[11] J. Park and M. E. J. Newman, Phys. Rev. E , 066117 (2004).[12] P. W. Holland and S. Leinhardt, Journal of the American Statis-tical Association , 33 (1981).[13] O. Frank and D. Strauss, Journal of the American Statistical Association , 832 (1986).[14] D. Strauss, SIAM Review , 513 (1986).[15] S. Wasserman and P. Pattison, Psychometrika , 401 (1996).[16] C. J. Anderson, S. Wasserman, and B. Crouch, Social Networks , 37 (1999).[17] T. A. Snijders, Annual Review of Sociology , 131 (2011).[18] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii,and U. Alon, Science , 824 (2002).[19] S. Maslov, K. Sneppen, and A. Zaliznyak, Physica A: Statisti-cal Mechanics and its Applications , 529 (2004).[20] T. Squartini and D. Garlaschelli, New Journal of Physics ,083001 (2011).[21] G. Palla, I. Der´enyi, I. Farkas, and T. Vicsek, Phys. Rev. E ,046117 (2004).[22] D. Jeong, M. Y. Choi, and H. Park, J. Phys. A , 9723 (2007).[23] A. Gorsky and O. Valba, Journal of Complex Networks (2020).[24] J. Berg and M. L¨assig, Phys. Rev. Lett. , 228701 (2002).[25] J. Park and M. E. J. Newman, Phys. Rev. E , 066146 (2004).[26] Z. Burda, J. Jurkiewicz, and A. Krzywicki, Phys. Rev. E , ,026106 (2004).[28] A. Annibale, A. C. C. Coolen, L. P. Fernandes, F. Fraternali,and J. Kleinjung, J. Phys. A , 485001 (2009).[29] A. Annibale and O. T. Courtney, J. Phys. A , 365001 (2015).[30] A. Coolen, Journal of Physics: Conference Series , 012022(2016).[31] F. A. L´opez and A. C. C. Coolen, J. Phys. A , 065002 (2020).[32] F. A. Lopez and A. C. Coolen, “Transitions in loopy randomgraphs with fixed degrees and arbitrary degree distributions,”(2020), arXiv:2008.11002 [cond-mat.dis-nn].[33] M. E. J. Newman, Phys. Rev. Lett. , 208701 (2002).[34] M. E. J. Newman, Phys. Rev. E , 026126 (2003).[35] M. Mayo, A. Abdelzaher, and P. Ghosh, Computational SocialNetworks , 4 (2015).[36] Y. Fujiki, T. Takaguchi, and K. Yakubo, Phys. Rev. E ,062308 (2018).[37] V. M. Egu´ıluz and K. Klemm, Phys. Rev. Lett. , 108701(2002).[38] M. Bogu˜n´a, R. Pastor-Satorras, and A. Vespignani, Phys. Rev.Lett. , 028701 (2003).[39] F. Sorrentino, M. di Bernardo, G. H. Cu´ellar, and S. Boccaletti,Physica D: Nonlinear Phenomena , 123 (2006).[40] I. Sendi˜na Nadal, I. Leyva, A. Navas, J. A. Villacorta-Atienza,J. A. Almendral, Z. Wang, and S. Boccaletti, Phys. Rev. E ,032811 (2015).[41] S. Jalan, A. Kumar, A. Zaikin, and J. Kurths, Phys. Rev. E , 062202 (2016).[42] C. Orsini, M. M. Dankulov, P. Colomer-de Sim´on, A. Ja-makovic, P. Mahadevan, A. Vahdat, K. E. Bassler, Z. Toroczkai,M. Bogun´a, G. Caldarelli, S. Fortunato, and D. Krioukov, Nat.Commun. , 8627 (2015).[43] Y. Fujiki, S. Mizutaka, and K. Yakubo, Eur. Phys. J. B , 126(2017).[44] D. Rybski, H. D. Rozenfeld, and J. P. Kropp, EPL (EurophysicsLetters) , 28002 (2010).[45] A. Arcagni, R. Grassi, S. Stefani, and A. Torriero, EuropeanJournal of Operational Research , 708 (2017).[46] S. Mizutaka and T. Hasegawa, Journal of Physics: Complexity , 035007 (2020).[47] Y. Fujiki and K. Yakubo, Phys. Rev. E , 032308 (2020).[48] J. Negele and H. Orland, Quantum Many-particle Systems , Ad-vanced Book Classics (Basic Books, 1988).[49] E. Guzm´an-Gonz´alez, I. P´erez Castillo, and F. L. Metz, Phys.Rev. E , 012133 (2020).[50] F. L. Metz and I. P´erez Castillo, Phys. Rev. Lett. , 104101(2016).[51] F. L. Metz and I. P´erez Castillo, Phys. Rev. E100