Small worlds and clustering in spatial networks

PHYSICAL REVIEW RESEARCH 2, 023040 (2020) Small worlds and clustering in spatial networks Marián Boguñá ,1,2,*Dmitri Krioukov,3,4,†Pedro Almagro ,1,2and M. Ángeles Serrano1,2,5,‡ 1Departament de Física de la Matèria Condensada, Universitat de Barcelona, Martí i Franquès 1, E-08028 Barcelona, Spain 2Universitat de Barcelona Institute of Complex Systems (UBICS), Universitat de Barcelona, Barcelona, Spain 3Network Science Institute, Northeastern University, 177 Huntington avenue, Boston, Massachusetts 022115, USA 4Department of Physics, Department of Mathematics, Department of Electrical & Computer Engineering, Northeastern University, 110 Forsyth Street, 111 Dana Research Center, Boston, Massachusetts 02115, USA 5Institució Catalana de Recerca i Estudis Avançats (ICREA), Passeig Lluís Companys 23, E-08010 Barcelona, Spain (Received 20 August 2019; accepted 12 March 2020; published 14 April 2020) Networks with underlying metric spaces attract increasing research attention in network science, statistical physics, applied mathematics, computer science, sociology, and other fields. This attention is further amplified by the current surge of activity in graph embedding. In the vast realm of spatial network models, only a few reproduce even the most basic properties of real-world networks. Here, we focus on three such properties—sparsity, small worldness, and clustering—and identify the general subclass of spatial homogeneous and heterogeneous network models that are sparse small worlds and that have nonzero clustering in the thermodynamic limit. We rely on the maximum entropy approach in which network links correspond to noninteracting fermions whose energy depends on spatial distances between nodes. DOI: 10.1103/PhysRevResearch.2.023040 I. INTRODUCTION In spatial networks, nodes are positioned in a geometric space, and the distances between them in the space affect their linking probability in the network [1]. In real-world systems, such spaces can be explicit/physical, as in geographically embedded networks [2,3] or in the Ising model [4–6]or epidemic models with long-range interactions [7,8]. Yet these spaces can be also hidden/latent. Latent similarity spaces have been employed for nearly a century to model homophily in social networks, for instance [9,10]: the closer the two people are in a virtual similarity space, the more similar they are, the more likely they know each other [11]. Another field where the space can be virtual are graph embeddings in computer science and machine learning, with applications including network compression, visualization, and node labeling [12,13]. In models of spatial networks, the space is usually explicit. Perhaps the simplest spatial network model is that of random geometric graphs that have been extensively studied in mathematics and physics since the early 1960s [14–18]. In random geometric graphs, nodes are positioned in a space randomly using a point process, usually a Poisson point process, and two nodes are linked in the graph if the distance between them in the space is less than a fixed threshold. If the intensity of *[email protected] †[email protected] ‡[email protected] Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license. Further distribution of this work must maintain attribution to the author(s) and the published article’s title, journal citation, and DOI. the point process does not depend on the graph size n, then the resulting graphs are sparse and have nonzero clustering in the thermodynamic n→∞limit, thus sharing these two properties with many real-world complex networks [19,20]. Yet many real-world networks are also heterogeneous small worlds, while random geometric graphs are homogeneous large worlds. This mismatch was resolved in Refs. [21,22] where a class of models of spatial networks that are sparse heterogeneous small worlds with nonzero clustering was introduced. Networks in these models have some additional properties commonly observed in real-world networks, such as selfsimilarity [21,23] and community structure [24–26]. Yet the following question remains: under what conditions do spatial network models have the mentioned properties of real-world networks? Here, we first focus on just three properties: (i) sparsity, (ii) small worldness, and (iii) nonzero clustering. The main contribution of our work is to show that spatial networks in Rdhave all these three properties at once if and only if the probability pij of connection between nodes iand jscales with the distance xij between them in Rdas pij ∼x−β ij with β∈(d,2d). In other words, the condition pij ∼x−β ij with β∈(d,2d) is not only a sufficient but also the necessary condition to have maximally random ensembles with realistic topological properties. In particular, we show that, among all possible choices of link energy Eas a function of distance x in a general ensemble of random geometric graphs, the only one that satisfies the properties stated above is the one where E∝ln xand so p∼x−βwith d<β<2d, leading to sparse, small-world, and clustered networks. We then add (iv) heterogeneity to the list of the requirements, and show that βmust be within the same range (d,2d) if the variance of the degree distribution is finite. 2643-1564/2020/2(2)/023040(10) 023040-1 Published by the American Physical Society BOGUÑÁ, KRIOUKOV, ALMAGRO, AND SERRANO PHYSICAL REVIEW RESEARCH 2, 023040 (2020) If it is infinite, however, e.g., if the degree distribution is a power law with exponent γ∈(2,3), then the networks are always ultrasmall worlds, and any β>dsatisfies all the four requirements. Finally, we show that if we also want to suppress nonstructural degree correlations, then the unique shape of the connection probability in the heterogeneous case is as in Refs. [21,22]: pij ∼(κiκj)β/dx−β ij , where κi,κj are the expected degrees of nodes i,j. This unique solution corresponds to the S1/H2model defined in Refs. [21,22]. We thus conclude that the S1/H2model plays the same role for random geometric graphs as the configuration model does for random graphs. Since random geometric graphs provide the simplest explanation for the emergence of clustering in complex networks, our models should be considered as the appropriate null models of clustering and degree distributions observed in many real-world networks. To obtain these results, we take a statistical physics stance in which we interpret spatial network models as probabilistic mixtures of grand canonical ensembles that maximize ensemble entropy under certain constraints, and are thus statistically unbiased. We call these mixtures hyper-grand-canonical ensembles, as some of their parameters are random. II. SETTINGS AND NOTATIONS We consider a very general class of spatial network models. The space is any compact homogeneous and isotropic Riemannian manifold of dimension dand volume n, and with no boundaries. We require the curvature of the manifold to go to zero at n→∞. That is, the space is locally the Euclidean space, and it is exactly the Euclidean space Rd in the thermodynamic limit. Examples are the d-sphere or d-torus of size growing with nsuch that its volume is n.Any growing compact d-dimensional hyperbolic manifold with no boundaries is also fine. On such a manifold we sprinkle n points uniformly at random according to the manifold metric. These points are thus the binomial point process of rate 1 on the manifold, and they form the node set of a random graph. Conditioned on node coordinates on the manifold, nodes iand jare connected independently with probabilities pij =p(xij), where xij is the distance between iand jon the manifold. By aij we denote the adjacency matrix of these random graphs: conditioned on node coordinates, aijsare independent Bernoulli random variables with success rates pij. These graphs are known as soft random geometric graphs [27,28]. We interpret these random graph ensembles as probabilistic mixtures of grand canonical ensembles that maximize ensemble entropy under the constraints that the average number of particles and average energy are fixed to given values. Particles are edges aij here, and their energies εij depend on distances xij:εij =f(xij). The connection probability function p(x) then takes the familiar Fermi-Dirac form, see Appendix A, p(xij)=1 eβ(f(xij)−μ)+1.(1) The Lagrange multipliers corresponding to the number-of- particles and energy constraints are the chemical potential μand inverse temperature β⩾0, as usual. We assume that neither f(x) nor βdepend on n,butμ, and consequently the absolute activity λ=eβμ, can depend on nas they usually do in statistical physics. Since energies εij are not fixed as in grand-canonical ensembles but are random instead, we call this ensemble a hyper-grand-canonical ensemble. We require our networks to be always sparse, meaning that the expected average degree in them is fixed to a finite positive constant kfor any network size n. We call a network model a small world if the average hop distance of shortest paths in the model networks grows slower than any polynomial of n. In particular, average distances growing as any polynomial of ln nin a model would render the model a small world. The model is also an ultrasmall world if the average distance grows slower than any polynomial of ln n. If a model is not a small world, we call it a large world. By clustering we mean the average local clustering coefficient. Symbol “∼”inan∼bn or a(x)∼b(x) means that an/bnor a(x)/b(x) converge to a finite positive constant at n→∞or x→∞, respectively. III. HOMOGENEOUS SPATIAL NETWORKS In any network model satisfying the settings above, the degree distribution is homogeneous: in the thermodynamic limit n→∞ it converges to the Poisson distribution with the mean equal to the average degree in the network [29,30]. By network homogeneity we mean here not only degree homogeneity, but also all the consequences of the manifest invariance of these ensembles with respect to the group of isometries of the manifold. In particular, the expected values of any graph property of any two nodes in these random graphs are the same. For instance, not only the expected degree, but also the expected clustering of any two nodes is the same and equal to the average clustering in the network. The main question is under what conditions these networks are small worlds and have nonzero clustering in the thermodynamic limit. The results are summarized in Table I. Intuitively, they are easy to comprehend. If f(x) grows too fast with x,so that p(x) decays too fast, then the network does not have sufficiently long links that are needed for small worldness. The network is thus necessarily a large world. On the other hand, if f(x) grows too slow with x, so slow that with an n-independent p(x) the average degree diverges, then to have an n-independent average degree the absolute activity λmust depend on nand go to zero at n→∞, meaning that pijsgoto zero as well. But since clustering scales with nthesameway as pijs do (recall that average clustering is the probability that two random neighbors of a random node are connected), it is zero at n→∞. Luckily, there exists a sweet spot at which the rate of growth of f(x) is not too fast and not too slow, so that the networks are small worlds and have nonzero clustering at the same time, the second regime in Table I. To show that this sweet spot (or range indeed) is as shown in Table I, we first observe that the average degree in our graphs is k=ki= j aij∼n1/d xd−1p(x)dx ∼λn1/d xd−1e−βf(x)dx =λIn.(2) 023040-2 SMALL WORLDS AND CLUSTERING IN SPATIAL NETWORKS PHYSICAL REVIEW RESEARCH 2, 023040 (2020) TABLE I. The result summary for homogeneous networks. Small world: yes/no: the networks are small/large worlds. Clustering: yes/no: the networks have nonzero/zero clustering in the thermodynamic (n→∞) limit. ER: Erd˝ os-Rényi random graphs [31]. RGG: sharp random geometric graphs in Rd[14]. Parameters: βthe inverse temperature, dthe space dimension, linf =lim infx→∞ f(x)/ln x, lsup =lim supx→∞ f(x)/ln x,wheref(x) is the energy function: εij =f(xij), where xij is the distance between nodes iand jin the space. Note that lsup =0 corresponds to f(x) growing slower than logarithmically, in which case the networks are in the first regime for any value of β<∞. Note that linf =∞corresponds to f(x)growing faster than logarithmically, in which case the networks are in the last regime for any value of β>0. Note that f(x)=cln x+o(ln x) corresponds to linf =lsup =c. The cases with β∈[d/lsup,d/linf], β∈[2d/lsup,2d/linf], and lsup ⩾2linf require further details about the specific shape of f(x) to classify the network into one of the three shown classes. Parameter regime Small world Clustering β→0 (ER) Yes No β<d/lsup d/linf <β<2d/lsup Yes Yes β>2d/linf No Yes β→∞(RGG) This is because the number of nodes at distances [x,x+dx] from a given node iin Rdis proportional to xd−1dx, node i is connected to each of those nodes with probability p(x), and we integrate up to the space diameter, which is ∼n1/d.The lower integration limit is any positive constant. If the integral Indiverges with n, then λmust depend on nand go to zero at n→∞ to yield an n-independent k above. But if λtends to zero, then pij ∼λe−βεij tends to zero as well, and so does clustering. The integral Indiverges if the monotonic function f(x) does not grow sufficiently fast. In particular, Indiverges if lsup =lim supx→∞ f(x)/ln x<d/β, the first regime in Table I. On the other hand, if Inconverges— in particular, it does so if linf =lim infx→∞ f(x)/ln x>d/β, the second row in Table I—then λis strictly positive, and so are pijs and clustering. Turning to small worldness now, the network is a small world only if it contains links connecting nodes located at distances xij of the order of the space diameter ∼n1/d,aswell as at all other smaller distances. Let l(x) be the distribution of link lengths x, defined as distances between linked nodes in the space. Observe that l(x)∼xd−1p(x)∼xd−1e−βf(x).If f(x)=cln x, then l(x)∼x−δwith exponent δ=cβ−d+ 1. Since networks are sparse, there are ∼nlinks. The expected maximum value among nsamples from a power law with exponent δis ∼nξwith ξ=1/(δ−1) [32]. The network is a small world only if this expected maximum link length is larger than the space diameter ∼n1/d, which implies ξ>1/d or β<2d/c, cf. the second row in Table Iwith lsup =c.If f(x) grows faster than logarithmically, linf =∞, then l(x) decays faster than a power law, ξ=0, and there are no long links at all, so that our networks are necessarily large worlds, the last regime in Table I. This logic is about the necessary conditions for small worldness, but they have been proven to be also sufficient 104105 size n 5 10 15 20 Average shortest path lenght ls β=0.5d β=1.1d β=1.3d β=1.5d β=1.7d β=1.9d β=1.98d 104105 size n 101 102 Average shortest path lenght ls β=2.0d β=2.2d β=2.4d β=2.6d β=2.8d β=3.0d 0.6 0.8 1 1.2 1.4 1.6 1.8 2 β/d 0.0 0.5 1.0 1.5 2.0 2.5 3.0 2 2.2 2.4 2.6 2.8 3 β/d 0.20 0.25 0.30 0.35 0.40 0.45 0.50 Exponent of the fit ls=a(ln n)bExponent of the fit ls=anb (a) (b) (d) (c) FIG. 1. The average shortest path length lsin the homogeneous spatial networks as a function of the network size n. This function is measured for different values of inverse temperature βin the connection probability (1) with f(x)=ln xused to generate random networks on the d=2-dimensional sphere of area n. The average degree in all these networks is fixed to k=10 by the appropriate choice of the chemical potential μ, and the results are averaged over 10 random network realizations for each data point. The functions ls(n) are then fit with alnbnfor β<2din (a), and with anbfor β⩾2din (b). The exponents bof these fits as functions of β/d are shown in the bottom panels (c) and (d). The dashed red line is ∼n1/2, the distance scaling in the two-dimensional sharp RGGs corresponding to β→∞. [33–35], and we confirm all the results above in simulations in Fig. 1(small worldness) and Fig. 2(clustering). Figure 1 shows that the average shortest path length lsscales with the network size nas ls∼lnbnif β<2d, and as ls∼nbif β⩾2d. In the small world regime β<2d, the exponent bin ls∼lnbnis close to 1 for any β<d, while for β∈(d,2d) it is a growing function of βthat appears not to diverge but to approach some finite maximum value as βapproaches 2d. In the large-world regime β⩾2d, exponent bin ls∼nb is also growing function of βranging in values from some minimum value at β=2dthat does not appear to be zero, to its theoretical maximum b=1/dat zero temperature β→∞ corresponding to sharp RGGs. The nature of the small-to- large world phase transition at β=2dappears to be an interesting open question [34]. Network sizes that are sufficiently large to provide any hints regarding whether this transition is continuous or discontinuous, can certainly not be reached in simulations. However, the results in Fig. 1suggest that the transition is discontinuous since the continuous transition would yield small-world b→∞at β→2d−and large-world b→0atβ→2d+. IV. HETEROGENEOUS SPATIAL NETWORKS Instead of the chemical potential μ, the Lagrange multiplier that fixes the expected average degree in the homogeneous ensemble, in the heterogeneous ensemble we have n Lagrange multipliers αithat fix the expected degree ki= jaijof each individual node to a desired value κi.The 023040-3 BOGUÑÁ, KRIOUKOV, ALMAGRO, AND SERRANO PHYSICAL REVIEW RESEARCH 2, 023040 (2020) 103104105106 size n 10-5 10-4 10-3 10-2 10-1 100 Clustering C(n) 103104105106 size n 10-5 10-4 10-3 10-2 10-1 100 Clustering C(n) β=0.5d β=0.75d β=1.1d β=1.5d β=2.0d β=3.0d 103104105106 size n 10-5 10-4 10-3 10-2 10-1 100 Clustering C(n) 0123 β/d 0.0 0.2 0.4 0.6 Clustering, C Homogeneous γ=3.5 γ=2.5 Homogeneous Heterogeneous, γ=3.5 Heterogeneous, γ=2.5 (a) (b) (c) (d) FIG. 2. The average local clustering coefficient in homogeneous and heterogeneous spatial networks as a function of the network size nfor different values of inverse temperature βand power-law exponent γ, (a)–(c). The space is the d=2-sphere of area n.The average degree is fixed to k=10. The results are averaged over 10 random network realizations for each data point. The bottom right panel (d) shows the average clustering as a function of β/dfor the largest network size n=3×105. The simulations confirm the continuous transition from the limiting zero to nonzero clustering at β=d. relations between κiand αiare documented in Appendix C. Here, we assume that the parameters κiare hyperparameters, meaning they are random and sampled from a fixed distribution ρ(κ), in which case we have the same hyper-grand- canonical ensemble as in the homogeneous case above, except that the connection probability changes from (1)to p(xij,α i,αj)=1 eβf(xij)+αi+αj+1.(3) The degree distribution in this ensemble converges to the mixed Poisson distribution P(k)=(1/k!)κκke−κρ(κ)dκ whose shape follows the shape of ρ(κ)[29,30]. This type of heterogeneous spatial network models were first introduced in Ref. [21], and many other similar classes of models have been defined and studied since then [36–38]. The qualitative behavior of clustering—zero versus nonzero in the thermodynamic limit—is exactly the same in these heterogeneous models as in the homogeneous one. Indeed, the expression (2) for the average degree kchanges to k∼n1/d xd−1dxα,α p(x,α,α)ρ(α)dαρ(α)dα ∼α e−αρ(α)dα2n1/d xd−1e−βf(x)dx =ˆ λIn,(4) where ˆ λ=e−α2, and ρ(α) is the distribution of Lagrange multipliers determined by the distribution of expected degrees ρ(κ). Following exactly the same reasoning as in the homogeneous case, albeit applied to ˆ λIninstead of λIn,we thus conclude that clustering is zero or nonzero at n→∞ TABLE II. The result summary for the heterogeneous networks with f(x)=ln xand Pareto ρ(κ)asinRef.[21]. The abbreviations are: HSCM: the hypersoft configuration model [39]; RHG: sharp random hyperbolic graphs in Hd+1[22]; USW: ultrasmall worlds; SW: small worlds; LW: large worlds; ZC: zero clustering at n→∞; PC: positive clustering at n→∞.Iff(x) grows slower or faster than logarithmically, then the networks are in the first and last rows, respectively. The γ=∞case is the homogeneous ensemble in Table I. Parameter regime 2 <γ <3γ>3, γ=∞ β→0(HSCM) USW,ZC SW,ZC β<d d<β<2dUSW, PC SW, PC β>2dUSW, PC LW, PC β→∞(RHG) depending on whether Indiverges or converges. For f(x)= ln xfor example, this means that the situation is exactly the same as in the homogeneous case: the clustering is zero if β<dand nonzero if β>d. Turning to small worldness, we assume henceforth that f(x)=ln x. We do so not only to simplify the discussion, but also because we prove in Appendix Bthat f(x)=ln x is unique in the sense that this is the only possible form of f(x) that does not induce any degree correlations other than the structural ones [32]. We also assume that the distribution ρ(κ) of expected degrees κis the Pareto distribution ρ(κ)=(γ−1)κγ−1 0κ−γ,where κ⩾κ0>0 and γ>2. (5) We note that the networks defined by (3,5) with f(x)=ln x were introduced in Ref. [21] and are equivalent to random hyperbolic graphs [22]. The calculation of the link length distribution l(x)inthis case yields l(x)∼x−δwith δ=β−d+1ifβ<d(γ−1), and δ=d(γ−2) +1 otherwise (details are given in Appendix D). Following the same logic behind the necessary conditions for small worldness as in the homogeneous case, which says that the networks can be small worlds only if ξ=1/(δ−1) >1/d, we conclude that small worlds are possible if β<2dor γ<3, or both. The networks are necessarily large worlds if β>2dand γ>3. A more detailed analysis proves that these necessary conditions for small worldness are also sufficient [40–42]. In fact, the qualitative clustering/small-worldness yes/no diagram for any γ>3is exactly the same as in Table Ifor the homogeneous ensemble with f(x)=ln xand linf =lsup =1. If γ<3, then our networks are worlds that are not only small but also ultrasmall, regardless of the value of β[43,44]. Table II collects all the results, and Figs. 2,3, and 4 confirm them in simulations. One sees in Fig. 3that if γ>3, then the simulation results are qualitatively similar to the homogeneous case, except that they are noisier, and the values of exponent bare significantly smaller. If the network size is small, such small values of bcan be deceiving, making these large worlds appear as small worlds. Figure 4shows the 023040-4 SMALL WORLDS AND CLUSTERING IN SPATIAL NETWORKS PHYSICAL REVIEW RESEARCH 2, 023040 (2020) 104105 size n 3 4 5 6 7 8 9 10 Average shortest path lenght ls β=0.5d β=1.1d β=1.3d β=1.5d β=1.7d β=1.9d 104105 size n 101 Average shortest path lenght ls β=2.0d β=2.4d β=3.0d β=4.0d β=6.0d β=10d 0.60.811.21.41.61.82 β/d 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 246810 β/d 0.10 0.15 0.20 0.25 0.30 0.35 Exponent of the fit ls=a(ln n)bExponent of the fit ls=anb (a) (b) (c) (d) FIG. 3. The average shortest path length lsin the heterogeneous spatial networks as a function of the network size n, (a)–(b). The settings are the same as in Fig. 1, except that the networks are heterogeneous (3) with Pareto ρ(κ)(5)andγ=3.5. (c) and (d) show the results of the same fittings as in Fig. 1. average shortest path length in the heterogeneous case with γ=2.5 as a function of n. However, due to the relatively small network sizes in the simulations, it is not possible to detect the ln ln nultrasmall world distance scaling that holds for any γ<3 and any β. However, one can see that the average shortest path length grows slower than ln nfor any β. 104105 size n 0.30 0.32 0.34 0.36 0.38 0.40 ls/ln(n) β=1.1d β=1.5d β=2.0d β=2.5d β=3.0d β=10d β=100d Heterogeneous γ=2.5 FIG. 4. The average shortest path length lsin the heterogeneous spatial networks with γ=2.5 as a function of the network size n. This function is measured for different values of the inverse temperature βin the connection probability Eq. (C9) with Pareto ρ(κ) with γ=2.5 used to generate random networks on the d= 2-dimensional sphere of area n. The average degree in all these networks is fixed to k=10 by the appropriate choice of the chemical potential ˆμ, and the results are averaged over 10 random network realizations for each data point. The functions ls(n)are divided by ln nto highlight the ultra-small-world distance scaling. We finally remark that the homogeneous ensemble is the γ→∞limit of the heterogeneous one, because at γ→∞ the Pareto distribution ρ(κ) becomes the degenerate distribution δ(κ−κ0), so that ρ(α)→δ(α+βμ/2) recovering (1) from (3). In the infinite temperature β→0 limit, the connection probability (3) is equal to 1/(eαi+αj+1), which is the connection probability in the hyper-grand-canonical or hypersoft configuration model that defines the unique ensemble of unbiased random graphs whose entropy is maximized across all graphs with a given degree distribution [39]. In the opposite zero temperature β→∞limit, the ensemble is equivalent to random hyperbolic graphs with a sharp connectivity threshold [22]. Finally, the γ→∞,β →0 limit is ER. V. OUTLOOK In summary, in spatial networks that are either homogeneous or have a finite degree distribution variance, the decay of the connection probability function with distance xin a space of dimension dmust be between ∼x−dand ∼x−2d to yield sparse small worlds with nonzero clustering. If the degree distribution variance is infinite though, then the spatial networks are ultrasmall worlds with any connection probability, and they have nonzero clustering if this probability decays with xfaster than x−d. Small worldness is linked to link energy and the distribution of link lengths. Networks are small worlds if they contain links of all lengths up to the space diameter. Clustering is dictated by the integrability of the connection probability function. If it is not integrable, then it must decay with the network size nto let the network be sparse, but then clustering is zero. This is directly related to the important notion of projectivity [45,46]: if the connection probability depends on n, then the network model is not projective, leading to nonlocal effects that cannot be present in any real-world network [39,47,48]. We thus see that any realistic model of sparse spatial networks must necessarily have nonzero clustering. As a final comment, we have presented spatial network models as hyper-grand-canonical ensembles, probabilistic mixtures of grand-canonical ones. In the latter ensembles, the constraints under which the ensemble entropy is maximized are clear: the average energy and the average number of particles in the ensemble, which fix the average link length and average degree, or a sequence of expected degrees, respectively. What remains unclear is under what constraints the considered hyper-grand-canonical ensembles are entropy maximizers. Are these constraints similar to the grand-canonical ones, or are they completely different, perhaps related to the expected number of triangles in the network [49]? In other words, what are the unbiased maximum entropy spatial network models for sparse heterogeneous small worlds with nonzero clustering? ACKNOWLEDGMENTS We acknowledge support from a James S. McDonnell Foundation Scholar Award in Complex Systems; the ICREA Academia prize, funded by the Generalitat de Catalunya; Ministerio de Economía y Competitividad of Spain Project No. FIS2016-76830-C2-2-P (AEI/FEDER, UE); the project Mapping Big Data Systems: embedding large complex networks 023040-5 BOGUÑÁ, KRIOUKOV, ALMAGRO, AND SERRANO PHYSICAL REVIEW RESEARCH 2, 023040 (2020) in low-dimensional hidden metric spaces - Ayudas Fundación BBVA a Equipos de Investigación Científica 2017; Generalitat de Catalunya Grant No. 2017SGR1064; the NSF Grant No. IIS-1741355; and ARO Grants No. W911NF-16-1-0391 and No. W911NF-17-1-0491. APPENDIX A: SPATIAL NETWORKS AS HYPER-GRAND-CANONICAL ENSEMBLES Let G{A;P(A)}be an ensemble of networks with adjacency matrices A={aij}and probability measure P(A). Let also Fl(A) be an arbitrary set of network functions. The canonical ensemble of random graphs that maximize the Gibbs entropy S=− A P(A)lnP(A)(A1) under the constraints that the ensemble averages of Fl, Fl= A Fl(A)P(A),(A2) are fixed to some values ¯ Fl, is given by the Boltzmann/Gibbs distribution P(A)=e−lαlFl(A) Z,(A3) where Z= A e−lαlFl(A)(A4) is the partition function, and αlthe Lagrange multipliers coupled to the constraints FL= ¯ Fl. The values of ¯ Fldetermine the values of αl[50]. The Gibbs distribution is known as an exponential family distribution in statistics, so that such canonical ensembles are called exponential random graphs there [51]. This distribution is known to be the unique unbiased distribution: it is proven that given the constraints, this is the unique distribution that encodes all the information contained in the constraints, and more importantly, it does not encode any other information [52–54]. Note that node pairs enumerate the n 2particle states {i,j}, i<j, that particles—that is, links—can occupy. If the graphs are simple and unweighted, then particles are fermions: there can be either zero or one particle at any particle state. If state {i,j}is occupied, then aij =1, and aij =0 otherwise. Different system states then corresponds to different networks A, and the number of particles in a system state Ais M(A)= i<j aij.(A5) Suppose now that nodes in these networks are nfixed points in any Riemannian manifold. The coordinates of these points define the distance matrix X={xij}between them on the manifold. Given any function f(x), we call εij =f(xij) the energy of the particle state {i,j}. The energy of the system state Ais then E(A)= i<j εijaij = i<j f(xij)aij.(A6) Consider now the canonical ensemble defined by just two functions F0(A)=E(A),(A7) F1(A)=M(A),(A8) and two constraints E= ¯ E,(A9) M= ¯ M,(A10) where ¯ E,¯ Mare given real numbers. Note that this ensemble is a vanilla grand-canonical ensemble in statistical physics that maximizes ensemble entropy under the average energy and number of particles constraints. The latter constraint fixes the average number of links and consequently the average degree, while the former constraint fixes the average link length. Denoting the Lagrange multipliers by α0=β, (A11) α1=−βμ, (A12) one can check [50] that the ensemble distribution (A3) can be written as P(A)=i<jeβ(μ−εij)aij Z = i<j paij ij (1 −pij)1−aij,(A13) where pij =1 eβ(εij−μ)+1(A14) is the connection probability that takes the standard Fermi- Dirac form, and where the values of the chemical potential μ and the inverse temperature βdetermine the average degree and link length, respectively. In heterogeneous networks, instead of one average degree constraint (A10), we have nper-node constraints ki= j aij=κi,(A15) where κiis any given sequence of expected degrees. Consequently, instead of one Lagrange multiplier (A12), we have nLagrange multipliers αi. The values of κidetermine the values of αivia the set of equations derived below for a particular case of interest. One can check that the homogeneous Fermi-Dirac connection probability (A14) changes in the heterogeneous case to pij =1 eβεij+αi+αj+1.(A16) Suppose now that the coordinates of nodes are no longer fixed on the manifold, but that they are random, e.g., that they are a binomial or Poisson point process. Then distances xij are random as well, and so are energies εij. Similarly, in the heterogeneous case, suppose that κiare no longer fixed either, but also random, e.g., sampled from a fixed distribution ρ(κ). The Lagrange multipliers αiare then also random. In both the homogeneous and heterogeneous cases 023040-6 SMALL WORLDS AND CLUSTERING IN SPATIAL NETWORKS PHYSICAL REVIEW RESEARCH 2, 023040 (2020) the connection probabilities pij are now random. Yet since links are still established with the same, albeit random, Fermi- Dirac connection probabilities (A14), the resulting ensembles with random coordinates (and random expected degrees) are probabilistic mixtures of the grand canonical ensembles defined above. We call these mixtures hyper-grand-canonical ensembles since energies εij (and degrees κi) are no longer parameters but hyperparameters instead. These hyper-grand- canonical mixtures are conceptually no different from how a (grand) canonical ensemble is itself a probabilistic mixture of microcanonical ensembles [55]. APPENDIX B: NO DEGREE CORRELATIONS ⇔f(x)=cln x Here we show that the energy function f(x)=cln xis the unique one leading to the absence of correlations of expected degrees in the considered ensembles in the thermodynamic limit, so that all degree correlations (if any) are only structural [32]. Formally, the absence of expected degree correlations means that the probability distribution P(κ|κ) of expected degrees κof nodes to which a random node of expected degree κis connected, does not depend on κ. Instead of expected degrees κ,κand distribution P(κ|κ), it is more convenient to work with the corresponding Lagrange multipliers α,α and distribution P(α|α), the latter also independent of αif there are no expected degree correlations. Using results from Ref. [29], P(α|α)=ρ(α)F(α+α) dαρ(α)F(α+α ),(B1) where ρ(α) is the distribution of αdefined by ρ(κ)giventhe relations between κs and αs as documented in the subsequent section, and function Fis defined to be F(α+α)=∞ 0 xd−1dx 1+eβf(x)+α+α.(B2) To find under which conditions P(α|α) is independent of α, we differentiate (B1) with respect to αand equate the result to zero to obtain F(α+α) F(α+α)=ρ(α )F(α+α)dα ρ(α)F(α+α )dα .(B3) Since the right-hand side of this equation does not depend on α, function Fis of the form F(x)=aebx, with aand bsome constants. Define q(x)≡ef(x)and z≡e−(α+α)/β to rewrite the uncorrelatedness condition as ∞ 0 xd−1dx 1+(q(x)/z)β=azbβ.(B4) That is, the network is uncorrelated at the level of hidden variables α,αwhenever Eq. (B4) holds for any value of z∈R+, with aand bsome constants. 1. If f(x)=cln x, then Eq. (B4)holds We first notice that the energy function f(x)=cln xis a sufficient condition for uncorrelatedness, since then Eq. (B4) trivially holds with bβ=d/cand a=∞ 0 td−1dt 1+tcβ.(B5) 2. If Eq. (B4) holds, then f(x)=cln x We are next to prove that in the small-world regime where f(x)/ln x∈(linf,lsup)forx>X,linf = lim infx→∞ f(x)/ln x,lsup =lim supx→∞ f(x)/ln x, and some constant X>0, the assumption that Eq. (B4) holds implies that f(x)=cln x∀x∈R+. We consider two cases. a. Case linf =lsup =c In this case, function f(x)/ln xhas a limit, lim x→∞ f(x) ln x=c,(B6) and since q(x)=ef(x),wehave lim x→∞ q(x) xc=1.(B7) If we define ˆq(x)≡q(x)/xc,Eq.(B4) can be written as zd c∞ 0 td−1dt 1+tcβˆqz1 ctβ=azbβ.(B8) If this equation holds for all values of z∈R+, then the integral in it must be a power of zfor any zincluding z1. Let us split this integral in two: xc(ε) z1/c 0 td−1dt 1+tcβˆq(z1 ct)β+∞ xc(ε) z1/c td−1dt 1+tcβˆq(z1 ct)β,(B9) where xc(ε) is such that for any x>xc(ε) we have that |ˆq(x)−1|<ε. We thus see that function ˆqis bounded in the integration domain of the second integral in Eq. (B9). This implies that lim z→∞ ∞ xc(ε) z1/c td−1dt 1+tcβˆq(z1 ct)β=∞ 0 td−1dt 1+tcβ,(B10) which is a constant independent of z. At the same time, the limit z→∞of the first integral in Eq. (B9) is zero because the domain of integration goes to zero and the integrand does not diverge at zero. Combining all these observations with Eq. (B8) we conclude that bβ=d/cand that a=∞ 0 td−1dt 1+tcβ[ˆq(z1 ct)]β=∞ 0 td−1dt 1+tcβ∀z∈R+.(B11) This is possible only if ˆq(x)=1, and hence f(x)=cln x. b. Case linf = lsup Let us assume now that linf and lsup are both positive and finite but not necessarily equal. The condition for uncorrelatedness in Eq. (B4) implies that there must exist a value of b 023040-7 BOGUÑÁ, KRIOUKOV, ALMAGRO, AND SERRANO PHYSICAL REVIEW RESEARCH 2, 023040 (2020) such that the limit lim z→∞ z−bβ∞ 0 xd−1dx 1+(q(x)/z)β(B12) exists. However, if f(x)/ln xis squeezed between linf and lsup at x1, then this integral is squeezed between z d lsup and z d linf at z1, and the limit does not exist, so that we arrive at a contradiction. We thus conclude that the only possibility is that linf =lsup =c, so that f(x)=cln x. Finally we remark that ccan be always set to 1 by a proper choice of energy units. APPENDIX C: RELATIONS BETWEEN κAND α Here, we derive these relations for the heterogeneous hyper-grand-canonical ensemble with the energy function f(x)=ln x, the Poisson point process of intensity 1 in Rd, and any distribution of expected degrees ρ(κ). The cases with β>dand β<dmust be considered separately. 1. Case β>d In this case, thanks to the integrability of the connection probability with respect to the spatial distance, we can work directly in the thermodynamic limit in Rd. Since the space is homogeneous we assume without loss of generality that a node with variable αis at the origin, and we want to calculate its expected degree κ. It is convenient to work in spherical coordinates in Rd, in which the volume element is dVRd=rd−1dr dVSd−1,(C1) where dVSd−1is the volume element on the unit (d−1) sphere whose volume is Sd−1=dVSd−1=2πd 2 d 2.(C2) The expected degree κof our node is then κ=Sd−1ρ(α)dα∞ 0 rd−1dr 1+reα+α ββ.(C3) Changing variables we simplify this to κ=Sd−1e−αd βe−αd β∞ 0 td−1dt 1+tβ =Sd−1e−αd βe−αd βπ βsin dπ β .(C4) By taking the average of Eq. (C4), we find the relation between the term e−αd βand the average degree k, which plugged again in Eq. (C4) leads to the final relation between αand κ: α=− β dln κ+1 2ln ˆμ,(C5) ˆμ= βsin dπ β πSd−1k= βd 2sin dπ β 2π1+d 2k .(C6) This result implies that the edge-state energy εij and chemical potential μin the ensemble are given by εij =ln xij (κiκj)1 d,(C7) μ=1 dln ˆμ. (C8) The connection probability can then be written as p(xij,κ i,κj)=1 1+xij (ˆμκiκj)1 dβ∼(ˆμκiκj)β/dx−β ij .(C9) 2. Case β<d In this case, the connection probability is not integrable with respect to distance, so that we have to take the finite size effects into account. This implies that the answer depends on a particular choice of the manifold family. Yet we remind that our general settings are such that for any nthe manifold volume is n, so that n=VdRd,(C10) where Ris the linear size of the manifold, while Vdis its volume at R=1. For example, if the manifold is a dtorus, then Ris its side length and Vd=1. If it is a dsphere, then R is its radius, and Vdis the volume Sdof the unit dsphere: Vd=Sd=2πd+1 2 d+1 2.(C11) We consider the case with the dsphere for concreteness. Since the dsphere is homogeneous, we assume without loss of generality that a node with variable αis at its north pole. The volume element on the dsphere with θthe polar angle is dVSd=sind−1dθdVSd−1,(C12) so that the expected degree of our node is κ=Sd−1Rdρ(α)dαπ 0 sind−1θdθ 1+Rθeα+α ββ.(C13) For R1 and β<dthe integral in the last expression can be approximated as π 0 sind−1θdθ 1+Rθeα+α ββ∼πd−β d−βR−βe−(α+α).(C14) Using this expression in Eq. (C13), we conclude that for n1 the relation between αand the expected degree κis given by α=− ln κ+1 2ln ˆμ,(C15) where ˆμ=d−β πd−β 1 Sd−1kSd n1−β d =(d−β)d 2 2π3d 2−βkn1−β d2πd+1 2 d+1 21−β d .(C16) 023040-8 SMALL WORLDS AND CLUSTERING IN SPATIAL NETWORKS PHYSICAL REVIEW RESEARCH 2, 023040 (2020) The connection probability is then p(xij,κ i,κj)=1 1+xβ ij ˆμκiκj ∼ˆμκiκjx−β ij .(C17) This connection probability depends on nand tends to zero as ∼1/n1−β/dsince so does ˆμ. We also note that it cannot be written as a Fermi-Dirac distribution function, meaning that the energy of edges cannot be defined in the case β<d. APPENDIX D: LINK LENGTH DISTRIBUTION l(x) Here, we calculate the tail of this distribution in the heterogeneous hyper-grand-canonical ensemble with the energy function f(x)=ln x, the Poisson point process of intensity 1 in Rd, and Pareto ρ(κ)∼κ−γwith γ>2. If β<d, then according to (C17), the distribution of link lengths is given by l(x)∼ p(x,κ,κ)xd−1ρ(κ)dκρ(κ)dκ ∼x−β+d−1κρ(κ)dκ2 ∼x−(β−d+1),(D1) since the Pareto distribution has a finite mean if γ>2. If β>d, then according to (C9), the expression for l(x) becomes l(x)∼x−β+d−1κβ/dρ(κ)dκ2 ∼x−β+d−1κ−γ+β/ddκ2 ,(D2) where the κintegral is finite if β<d(γ−1), so that l(x)is still ∼x−(β−d+1). If β>d(γ−1), then the κintegral in the last equation is infinite, so that slightly more care is needed to derive the scaling of l(x) with x. Specifically, for large βthe exact expression for the connection probability p(x,κ,κ)in(C9) can be approximated by 1 for x/(ˆμκκ)1/d<1, and by 0 otherwise, in which case we get l(x)∼xd−1xd ρ(κ)dκρ(κ)dκ∼x−[d(γ−2)+1].(D3) Collecting all the cases, l(x)∼x−δ,where δ=β−d+1,if β<d(γ−1) d(γ−2) +1,if β>d(γ−1).(D4) Note that the link length distribution l(x) is well-defined in the n→∞limit only if δ>1, i.e., if β>dcorresponding to the nonzero clustering regime. [1] M. Barthélemy, Phys. Rep. 499,1(2011). [2] V. Latora and M. Marchiori, Phys. A: Stat. Mech. Appl. 314, 109 (2002). [3]R.Guimera,S.Mossa,A.Turtschi,andL.a.N.Amaral,Proc. Natl. Acad. Sci. USA 102,7794 (2005). [4] M. Gitterman, J. Phys. A: Math. Gen. 33,8373 (2000). [5] J. Viana Lopes, Y. G. Pogorelov, J. M. B. Lopes dos Santos, and R. Toral, Phys. Rev. E 70,026112 (2004). [6] T. N. Teles, Fernanda P. da C. Benetti, R. Pakter, and Y. Levin, Phys. Rev. Lett. 109,230601 (2012). [7] D. Mollison, in Proceedings of the Sixth Berkeley Symposium on Mathematical Statistics and Probability, Vol. 3: Probability Theory,editedbyL.M.LeCam,J.Neyman,andE.L.Scott (University of California Press, 1972), p. 579. [8] D. Mollison, J. Roy. Stat. Soc. B 39,283 (1977). [9] A. P. Sorokin, Social Mobility (Harper, New York, 1927). [10] G. Majone, Qual. Quant. 6,239 (1972). [11] M. McPherson, L. Smith-Lovin, and J. M. Cook, Ann. Rev. Soc. 27,415 (2001). [12] A. Grover and J. Leskovec, in Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, San Francisco 2016 (Association for Computing Machinery, New York, 2016), p. 855. [13] P. Goyal and E. Ferrara, Knowledge-Based Syst. 151,78 (2018). [14] E. N. Gilbert, J. Soc. Ind. Appl. Math. 9,533 (1961). [15] M. Penrose, Random Geometric Graphs (Oxford University Press, Oxford, 2003). [16] J. Dall and M. Christensen, Phys. Rev. E 66,016121 (2002). [17] J. P. Coon, C. P. Dettmann, and O. Georgiou, Phys. Rev. E 97, 042319 (2018). [18] K. Kosmidis, S. Havlin, and A. Bunde, EPL 82,48005 (2008). [19] A.-L. Barabási, Network Science (Cambridge University Press, Cambridge, 2016). [20] M. E. J. Newman, Networks (Oxford University Press, Oxford, 2018). [21] M. Á. Serrano, D. Krioukov, and M. Boguñá, Phys.Rev.Lett. 100,078701 (2008). [22] D. Krioukov, F. Papadopoulos, M. Kitsak, A. Vahdat, and M. Boguñá, Phys.Rev.E82,036106 (2010). [23] G. García-Pérez, M. Boguñá, and M. Á. Serrano, Nature Phys. 14,583 (2018). [24] K. Zuev, M. Boguñá, G. Bianconi, and D. Krioukov, Sci. Rep. 5,9421 (2015). [25] A. Muscoloni and C. V. Cannistraci, New J. Phys. 20,052002 (2018). [26] G. García-Pérez, M. Á. Serrano, and M. Boguñá, J. Stat. Phys. 173,775 (2018). [27] M. Penrose, Ann. Appl. Prob. 26,986 (2016). [28] C. P. Dettmann and O. Georgiou, Phys.Rev.E93,032313 (2016). [29] M. Boguñá and R. Pastor-Satorras, Phys. Rev. E 68,36112 (2003). [30] B. Bollobás, S. Janson, and O. Riordan, Random Struct. Algor. 31,3(2007). [31] P. Erdös and A. Rényi, Publ. Math. Debrecen 6, 290 (1959). 023040-9