The birth of geometry in exponential random graphs
TThe birth of geometry inexponential random graphs
Pawat Akara-pipattana a , Thiparat Chotibut a,b and Oleg Evnin a,c a Department of Physics, Faculty of Science, Chulalongkorn University, Bangkok, Thailand b Chula Intelligent and Complex Systems, Faculty of Science,Chulalongkorn University, Bangkok, Thailand c Theoretische Natuurkunde, Vrije Universiteit Brussel andThe International Solvay Institutes, Brussels, Belgium [email protected], [email protected], [email protected]
ABSTRACT
Inspired by the prospect of having discretized spaces emerge from random graphs, we con-struct a collection of simple and explicit exponential random graph models that enjoy, inan appropriate parameter regime, a roughly constant vertex degree and form very largenumbers of simple polygons (triangles or squares). The models avoid the collapse phenom-ena that plague naive graph Hamiltonians based on triangle or square counts. More thanthat, statistically significant numbers of other geometric primitives (small pieces of regu-lar lattices, cubes) emerge in our ensemble, even though they are not in any way explicitlypre-programmed into the formulation of the graph Hamiltonian, which only depends on prop-erties of paths of length 2. While much of our motivation comes from hopes to construct agraph-based theory of random geometry (Euclidean quantum gravity), our presentation iscompletely self-contained within the context of exponential random graph theory, and therange of potential applications is considerably more broad. a r X i v : . [ c ond - m a t . d i s - nn ] F e b Introduction
The idea that all of physics should be a manifestation of geometry has been at the heart ofthe theory of fundamental interaction over the last hundred years. But can geometry itselfarise from something that is more fundamental? In particular, starting with nothing but aset, can relations between its elements serve as precursors of geometric structures? One concrete realization of such abstract ideas is an attempt to derive discretized spacesfrom ensembles of random graphs. Extensive studies have been undertaken of random dis-retized spaces (simplicial complexes), with the aim of obtaining a Euclidean theory of quan-tum gravity [4, 5]. In recent works [6–8], an approach was formulated to deriving similarstructures from an ensemble of random graphs (the edges of simplicial complexes, after all,can be viewed as extremely constrained graphs). The random graphs considered in [6–8]are themselves rather constrained. For example, in [7], one deals with bipartite, 2 d -regulargraphs satisfying the hard-core constraint (no two squares may share more than one edge),while in [8] only the bipartiteness condition is dropped. Such constraints have to be rigidlyimposed at each vertex and each edge.We would like to develop less constrained models of appeal for studying the emergence ofgeometry. A context that we find particularly attractive is exponential random graph models(for textbook treatments, see [9, 10]). These models, which are natural from both physicaland information-theoretic viewpoints, assign Boltzmann-Gibbs exponential weights to allpossible ways of connecing a set of N points to each other, without a priori excluding anyspecific graphs. In accordance with the principle of maximum entropy [11], this probabilitymeasure maximizes the Shannon entropy (the lack of knowledge) subject to the constraintson the expectation values of selected observables, say, the total number of edges, triangles,etc.In the context of emergent geometry, one immediate question is then as follows: can oneformulate an exponential random graph model that gives rise to sparse graphs with smalldegree variation, and a large number of geometric primitives (say, triangles or squares, butwe shall also look at more extended geometric structures)? The last point is known to benontrivial: already the classic work [12] notes that attempting to specify the expectationvalue of the number of triangles in simple exponential graph models leads to the emergenceof very sparse (nearly empty) and very dense (nearly complete) graph configurations, ratherthan to ‘reasonable’ graphs with finite degrees and large numbers of triangles. A similarundesirable graph collapse phenomenon occurs in other exponential graph models whoseHamiltonians favor small clusters of particular shapes; see, for example, a description ofsuch phase transitions between extremely dense and extremely sparse phases in a simpletwo-star exponential random graph model in [13]. Some solutions to this graph collapseproblem have been proposed [7, 8, 14], but they involve controlling the connectivities of eachvertex individually, and such rigid constraints are what we would like to avoid. Anotherway to escape the collapse has been proposed in [15], which is in the form of a conventionalexponential graph model, but with a complicated graph Hamiltonian, and not adapted togeneration of graphs with approximately constant degrees. Some further related literature A notable research program that shares much of the same spirit, not directly related to our currentpursuits, is the theory of causal sets [1]. More broadly, considerable advocacy for the emergence of all lawsof physics from operations on discrete sets has been put forward in [2, 3].
1n random graphs viewed from a statistical physics perspective can be found in [16–20]; theproblem has also attracted attention of mathematicians working in probability theory [21,22].We shall then formulate and study two models that meet the above guidelines. In orderto attain an approximately constant vertex degree in the context of an exponential randomgraph model, we shall introduce a large negative two-star term in the graph Hamiltonian(this term can also be seen as the total number of paths of length 2). As explained in [19],such terms serve to reduce degree variance. We then have to design an additional termthat will boost the number of triangles or squares. For reasons explained above, naivetriangle or square counts do not work in this capacity, as they lead to graph collapse. Wetake inspiration from the hard-core constraint used in [7, 8] and replace this rigid constraintby a term in the graph Hamiltonian that favors configurations that meet the hard-corecondition. This can be done in a way adapted to the production of triangles or squares,which gives rise to two alternative, closely related models. These models, which we introducein section 2, justify our expectations and produce a large number of geometric primitives.Numerical results from Monte Carlo sampling of our graph ensembles are reported in section3. Besides the triangle and squares, we also observe small pieces of triangular lattices andcubes. These latter structures are emergent, as they do not explicitly appear in any way inthe graph Hamiltonian, which is entirely formulated in terms of properties of paths of length2. We conclude with brief comments on the analytic properties of our models in section 4,possible emergence of large-scale geometric features (as opposed to the small-scale geometricprimitives whose formation we see explicitly) in section 5, and a discussion of the prospectsin section 6.
It is customary in exponential random graph models [9, 10] to represent graphs as a set orrandom connections between N points (vertices) labelled by i = 1 , . . . , N . The connectionsare encoded in the adjacency matrix whose element c ij equal 1 if there is an edge betweenvertices i and j and zero otherwise. For undirected graphs without self-looping, c is asymmetric matrix with zero diagonal elements. One then defines a graph Hamiltonian H as a function of c and assigns Boltzmann weights proportional to e − H ( c ) to each graphconfiguration. The expectation value of an observable A ( c ) in this ensemble is then given by (cid:104) A (cid:105) = 1 Z (cid:88) { c } A ( c ) e − H ( c ) , (1)where the sum is over all possible adjacency matrices and Z ≡ (cid:80) { c } e − H ( c ) is the partitionfunction. In practice, as will be specified in the next section, one would use Monte Carlosimulations to sample this ensemble.The question is then how to choose the graph Hamiltonian. A good starting point is thetwo-star model that has been extensively studied in the literature [9, 10, 16, 17, 19]: H ∗ = − α (cid:88) ij c ij − β (cid:88) ijk c ij c jk . (2)2n the two-star model, the Lagrange multiplier α controls the number of edges in the graph,while the Lagrange multiplier β controls the number of ‘two-stars’, which is related to thetotal number of paths of length 2. Our first aim is to obtain graphs with a small ( N -independent) roughly constant vertex degree, in other words, graphs that are sparse and(approximately) regular. We can use α to control the vertex degree, while, as explainedin [19], a large negative β can suppress degree variance. This effect of β can be understoodintuitively as follows. Since the β -term in the Hamiltonian can be rewritten as (cid:80) i k i where k i is the degree of vertex i , a large negative β strongly suppresses high vertex degrees. Then,since the mean degree is kept at the desired value by adjusting α , it implies that degrees belowthe mean value are also strongly suppressed, leading to a small degree variance. Detailedanalysis can be found in [19].The two-star model by itself does not possess any significant geometric properties (forinstance, its triangle count is very low). We shall then use it as a point of departure forfurther modifications and look for Hamiltonians of the form H ( c ) = H ∗ ( c ) − σH mod ( c ) , (3)trying to find H mod that generates an ensemble with the desired properties. One couldnaively choose H mod as the total number of triangles in the graph, given by Tr( c ) /
6, butthis is well-known to fail since the early days of random graph simulations [12,18]. If one triesthis guess with β = 0, the triangle term, whose maximal value goes like O ( N ) is capableof overpowering in the thermodynamic limit the α -term which goes at most as O ( N ). As aresult, the model forms two phases, one very dense and one very sparse, and it is impossibleto tune the parameters to attain reasonable vertex degrees and high numbers of trianglesat the same time. If one tries to use a large negative β , as described above, to make thegraph approximately regular, it is possible to attain large numbers of triangles in such anapproximately regular graph, but at the cost of another terminal pathology: the graph burstsinto a collection of disconnected nearly complete graphs with just a few vertices each, whichhave been referred as ‘baby universes’ in [7, 8] following the spirit of the simplicial quantumgravity literature [4, 5]. Such dust-like configurations are of little interest insofar as one isconcerned with approaching large-volume discretized spaces. One can gain some intuitiveunderstanding of the above fracturing phenomenon by counting triangles in d -regular graphs.The maximal number of triangles formed around a given vertex in such a graph is d ( d − / N d ( d − / N/ ( d + 1)disconnected complete graphs with ( d + 1) vertices, and hence d ( d + 1)( d − / H mod replaced by the counts of other similar clusters(for example, squares), meet a similar fate. They tend to suffer condensation transitions ifthe vertex degree is unconstrained, and fracturing if the degree is constrained (either rigidly,or by a large negative β -term in the graph Hamiltonian). As an aside, counting formulas Our conventions differ from some of the literature, but are equivalent to it, and convenient for ourpurposes. The actual number of edges is (cid:80) i>j c ij , and the number of two-stars is (cid:80) i>k (cid:80) j c ij c jk , but theLagrangian multipliers for these quantities are related to α and β in (2) by a simple linear redefinition. hard-core condition . This hard-core conditionis crucial in the construction and says that no two squares (cycles of length 4) may possessmore than one edge in common. In particular, the hard-core condition allows subgraphs ofthe form , (4)but forbid subgraphs of the form . (5)We would like to avoid introducing such hard constraints (and keep nonzero, possibly small,probabilities for all conceivable graphs), but we can mimic this constraint in the context ofthe exponential random graph model (3) by choosing H mod that favors (4) but not (5). Thus,for every pair of vertices, we shall give a reward if there are exactly two paths of length 2between them. Formally, this is most conveniently expressed by first introducing the squareof the adjacency matrix q ij = (cid:88) k c ik c kj , (6)which simply counts the number of paths of length 2 between vertices i and j . We thenwrite H mod = H sq ≡ (cid:88) i>j δ ( q ij , , (7)where δ is the ordinary Kronecker symbol: δ ( n, m ) = 1 if n = m, and 0 otherwise . (8)Mnemonically, (7) is simply the number of entries 2 in the (off-diagonal) upper triangle ofthe matrix square of the adjacency matrix c . Configurations of the form (4-5) belong to theclass called k -independent two-paths in [15], though the way they are used in our action isdifferent from [15]. 4he model defined by (7) is adapted to producing squares while avoiding clumping/frac-turing, and we shall see in the next section that it is precisely what it does. Note, inparticular, that the small complete d -regular graphs that underlie fracturing as explained inthe passage under (3) are not favored by (7), while they would have been favored by a naive H mod counting the number of squares.Furthermore, one can straightforwardly modify (7) to make it adapted to production oftriangles rather than squares. Namely, we want a Hamiltonian that favors subgraphs of theform , (9)which contain two triangles, but does not favor subgraphs of the form , (10)which would have encouraged fracturing/clumping (more than two triangles share an edge).This is accomplished by choosing H mod = H tr ≡ (cid:88) i>j δ ( c ij q ij , , (11)This is just the number of entries 2 in the upper triangle of the Hadamard (entrywise)product of the adjacency matrix c and its ordinary matrix square q . Configurations of theform (9-10) belong to the class called k -triangles in [15], though the way they are used inour action is different from [15]. We now set out to explore the models defined by (1-3) and (7) or (11). While the graphHamiltonians are rather simple, and at some level are only slightly more complex than themuch-studied two-star model, the most viable approach is to analyze the correspondinggraph ensembles using computer simulations. We shall give some preliminary comments onthe analytic structure of the models in the next section.
General discussion of random graph simulations using Monte Carlo techniques can be foundin [10,25]. We use a simple sequence of stochastic flips that randomly select a pair of distinctvertices i and j , and then either erase the edge between them, if it exist, or add an edge, if5here is none to start with. The change of the adjacency matrix under such an ( ij )-flip canbe written as ∆ ij c mn = ( δ im δ jn + δ in δ jm )(1 − c ij ) . (12)As usual in Metropolis-type algorithms, this proposed ( ij )-flip is accepted with probabilitymin(1 , e H ( c ) − H ( c +∆ ij c ) ).For an efficient implementation of this algorithm, it is important to optimize the compu-tation of the variation of the graph Hamiltonian (3) under the flip. The two-star part (2) isvery simple, as it can be completely expressed through the vertex degree sequence k n = (cid:88) m c mn . (13)We can now rewrite (2) as H ∗ = − α (cid:88) i k i − β (cid:88) i k i . (14)The variation of the degree sequence is simply∆ ij k n = ( δ in + δ jn )(1 − c ij ) . (15)It then follows [10] from (15) that the variation of the two-star Hamiltonian is∆ ij H ∗ = − − c ij )( α + β ( k i + k j )) + 2 β. (16)What remains is to explain how to compute the variation of H mod . We start with thecase of the ‘square’ model (7), which is slightly simpler. Since H sq is expressed through q ,the square of the adjacency matrix, it suffices to write down the variation of q under the flip:∆ ij q mn = ( δ im c jn + δ jm c in + δ jn c im + δ in c jm )(1 − c ij ) + ( δ im δ in + δ jm δ jn ) . (17)All the changes in q thus occur in rows and columns number i and j . Without recomputingthe whole matrix, which would have been costly for large N , one can simply recomputethese two rows. Comparing the number of entries 2 in the old rows and the new ones givesthe variation H sq , and, together with the variation of H ∗ , can be converted into the flipacceptance probability. If the flip is accepted, the corresponding rows in q , which is storedas an explicit matrix, are updated. In this way, the number of operations required per flipgrows very moderately as O ( N ). One could in fact try to use the sparsity of the adjacencymatrix and its square in the low density regime that we focus on to further reduce the numberof operations per step to O (1), but this would require dealing with sparse matrix formatsthat have their own disadvantages, and we shall avoid it.A similar prescription can be given for the modified triangle model (11). The Hamiltonianis expressed through the products c mn q mn whose variation is given by∆ ij ( c mn q mn ) = q mn (∆ ij c mn ) + c mn (∆ ij q mn ) + (∆ ij c mn )(∆ ij q mn ) . (18)Note that ∆ ij c mn is only non-zero if m = i and n = j or m = j and n = i , in whichcase ∆ ij q nm is zero, hence the last term drops out. The term c mn ( δ im δ in + δ jm δ jn ) from the6xpansion of the second term in (18), upon the substitution of (17), also vanishes becausethe diagonal elements of c equal zero. Thus, (18) can be simplified to∆ ij ( c mn q mn ) = q mn ( δ im δ jn + δ in δ jm )(1 − c ij ) + c mn ( δ im c jn + δ jm c in + δ jn c im + δ in c jm )(1 − c ij ) . (19)As before, all the changes of the matrix ˜ q defined by ˜ q mn = c mn q mn are in rows and columnsnumber i and j . The way to proceed then is to store in memory the matrices c , q and ˜ q ,updating after each flip the relevant rows and columns, which requires O ( N ) operations perstep.As laid out among our motivations, we intend to obtain graphs of roughly constant smalldegree. We choose this reference degree to be 6, which for the triangle model is expectedto encourage formation of pieces of two-dimensional triangular lattice, and for the squaremodel, pieces of three-dimensional cubic lattice. The degree variance is kept low by takinga large negative value of the two-star parameter β . One then chooses σ as one wishes, but α , the Lagrange multiplier for the number of edges, has to be adjusted to make the meandegree close to 6. This is in practice implemented by running a sort of ‘adaptive’ MonteCarlo simulation, in which, every now and then, the current mean degree is measured and α is adjusted in proportion to how much it deviates from 6. Once this ‘adaptive’ simulationhas converged, one can take the resulting value of α and rerun a completely conventionalMonte Carlo simulation with α fixed at that value.It is a separate question how the parameters α , β , σ should depend on N for comparingthe simulations at different N . Our numerical experiments indicate that keeping β and σ fixed and choosing the value of α that brings the mean degree close to 6 results in approx-imately the same number of geometric features, to be defined below, in proportion to N .One expects that the value of α necessary for maintaining a small fixed degree varies loga-rithmically with N , as in the two-star model [19], but this variation is rather small in anycase for the range of N we consider (200 to 5000), and we do not keep track of it explicitly.(The value of α is thus simply chosen at each N for the given N -independent β and σ toascertain that the mean degree is close to 6.) We have performed simulations of the modified triangle model defined by (2), (3) and (11),and observed that the model stands up to our expectations. In an appropriate parameterregime that we shall specify shortly, the model forms a very large number of triangles andbigger geometric structures, without experiencing any clumping or fracturing that plaguenaive Strauss-type models.To quantify the geometric features of our random graphs, we introduce the trianglefraction η ∆ , given by the number of triangles [23, 24] n ∆ = Tr( c ) /
6, divided by the maximalnumber of triangles in a d -regular graph, given by N d ( d − /
6. Since we tune α and β tomake our graphs approximately 6-regular, N d ( d − / N , and hence η ∆ = n ∆ N . (20)7 .00% 0.00% -25 0.0 -50 -75 σβ Figure 1: Random graphs of the modified triangle model for a range of values of β and σ ,with α adjusted so that the mean degree is close to 6 (the actual values are within the range6 ± . , with the extra requirement that the degree of the central vertex is precisely six (which meansit has no connections to the vertices outside this subgraph), and the peripheral vertices do nothave any extra connections to each other besides the one displayed (while their connectionsto the rest of the graph are unrestricted). Since each graph vertex may lie at the center ofat most one hexagon as defined, it is natural to introduce the hexagon fraction relative tothe number of vertices in terms of the hexagon number n hex as η hex = n hex N . (21)With these preliminaries, we now set out to explore the parameter space of our model.A quick impression can be obtained from Fig. 1, where we have chosen a relatively smallnumber of vertices ( N = 200) to make the visualization transparent in a regular printout,8 Figure 2: Convergence of the mean degree, degree variance, and the fractions of trianglesand hexagons at N = 2000 , α = 177 , β = − , σ = 100, starting from an empty graph initialstate. Each curve is rescaled so that the final steady-state value is one. The horizontal axiscounts the elementary Monte Carlo steps in units of 10 million.though the picture is qualitatively similar for the higher values of N that we have explored(up to 5000). What one notices immediately in visual terms is the presence of two distinctregions. If σ is insufficiently large (below the diagonal), the graph is visually threadball-like,with very small triangle counts and no hexagons. Above the diagonal, on the other hand,the graphs have a crystalline appearance with a large number of triangles and appreciablenumber of hexagons. We emphasize that the observed number of triangles (a third, or more,of the maximal number possible at this vertex degree) is extremely large, expressing a high‘transitivity’ in the language of Strauss [12]. By contrast, in the simplest Erd˝os-R´enyi graphwith a fixed small ( N -independent) mean vertex degree, the fraction of triangles tends to 0for large graphs, while the probability of finding even a single hexagon is 0. For the rest ofour treatment, we shall focus on graphs with strong geometric features, above the diagonalin Fig. 1. Note that it would be interesting to explore the phase transition that presumablyseparates the two regions in Fig. 1, but it is a topic apart from our main objectives here.We comment on the convergence of our Monte Carlo simulations, exemplified by Fig. 2,now for a larger graph with N = 2000. One observes that the mean degree converges first,followed by the degree variance, while the geometric structures form an stabilize later. Themean degree and degree variance are defined as (cid:80) i k i /N and ( (cid:80) i k i /N − ( (cid:80) k i ) /N ) / ,while the triangle and hexagon fractions are defined as in (20) and (21). We have performedsimulations at various values of N and our rough estimate is that the number of Monte Carlosteps necessary for the number of hexagons to stabilize grows like O ( N . ). Graphs in this region visually appear rather similar to some of the graphs in [3], which are determin-istically generated. This is perhaps simply a reflection of the fact that triangle relations play a role in theformulations both here and in [3]. .0 0.2 0.4 0.6 0.8 1.05.765.845.926.006.08 (a) mean degree (b) degree variance (c) percentage of triangles (d) percentage of hexagons Figure 3: Outcomes of multiple simulations of the modified triangle model at N = 2000 , α =177 , β = − , σ = 100, with initial states in the form of Erd˝os-R´enyi graphs whose edgeoccupation probability is given on the horizontal axes. The green crosses indicate the meanof each vertical column.In the regime that interests us, namely, a large value of σ and a developed geometricstructure in the graphs, the model displays jamming behaviors, and the values of observablesat which the Monte Carlo evolution stabilizes vary somewhat between different runs with thesame values of thermodynamic parameters. This is perhaps not very surprising since, if thegraphs have a ‘crystalline’ appearance, they form in a process of spontaneous crystallization,where jammed behaviors often occur. Indeed, once a graph has formed with a large numberof configurations of the form (9), which our Hamiltonian favors, it becomes difficult to moveedges around without upsetting this order, so as to achieve even more optimal configurations.In practice, as seen from Fig. 2, the number of hexagons, for example, simply freezes aftera sufficiently long run, and stops fluctuating altogether on any time scales our numericalsimulations can reasonably access. In order to quantify the variations in these final jammedconfigurations and the failure to reach true equilibrium, we have prepared Fig. 3 whichsummarizes the outcomes of multiple runs with the same thermodynamic parameters, butstarting with initial configurations given by Erd˝os-R´enyi graphs of different densities. Onecan see that the mean degree converges to a predictable value within a narrow range, while thedegree variance, and the triangle (20) and hexagon (21) fractions display greater variation,10igure 4: A large graph obtained from simulations at N = 5000 , α = 177 , β = − , σ = 100.with the hexagon fraction being particularly variable. Irrespectively of that, all of the graphsthat come out of our simulations contain geometric primitives (triangles and hexagons inour case) at rates orders of magnitude higher than what one could hope to obtain fromexponential random graph models without the extra terms we have designed, which is whatwe see as a quantification of the success of our model.To conclude, we present in Fig. 4 a particular very large graph with N = 5000 , α =177 , β = − , σ = 100 that has come out of our simulations, in order to give the readers animpression of its visual qualities. 11 .3 The modified square model The modified square model defined by (2), (3) and (7) displays features very similar to whatwe have described above for the triangle model, with the role of triangles taken up by squares,and the role of hexagons, by cubes.The number of squares can be computed as [23, 24] n (cid:50) = (Tr( c ) + Tr( c ) − (cid:80) i k i ) / N/ ( d + 1) complete d -regular graphs with ( d + 1) points each, given by N d ( d − d − / N for d = 6, yielding our definition η (cid:50) = n (cid:50) N . (22)We furthermore consider the number of cubes defined as subgraphs of the form , in which the eight depicted vertices are connected to each other exactly as given, withoutany extra connections, but they may be connected to the rest of the graph in an arbitraryway. We then compute the number of cubes n cube and define the cube fraction in proportionto the total number of vertices η cube = n cube N . (23)The phenomenology of the modified square model is similar to the previous section. If σ isinsufficiently large, one obtains graphs reminiscent of pure two-star model simulations, withvery small numbers of geometric primitives. If σ is large, one obtains graphs with extremelylarge numbers of squares and cubes, constituting a significant fraction of the maximal numberpossible at this degree. In this regime, there are jamming phenomena, similar to what hasbeen described for the triangle model, in the sense that different Monte Carlo runs freezeat different numbers of geometric primitives, though for all runs the number of squares andcubes is very high, which represents a success of our model. We summarize the observationaldata of a particular set of runs in Fig. 5, which is analogous to Fig. 3 for the triangle model.Some differences in fine details are that the degree variance is higher (given the mean degreeand percentage of squares comparable to the mean degree and percentage of triangles in oursimulations of the triangle model), while the cube production is extremely high. While we have observed complex behaviors in our modified triangle and square models, suchas emergence of large counts of geometric primitives not explicitly coded into the graphHamiltonian, as well as moderately jammed dynamics in the high σ regime, we would liketo emphasize that the graph Hamiltonians given by (2-3) and either (7) or (11) are rathersimple in terms of their algebraic structure. In particular, if the modified square Hamiltonian(7), or the β -term of the two-star Hamiltonian (2), is expressed through the square of the12 .0 0.2 0.4 0.6 0.8 1.06.26.46.66.87.0 (a) mean degree (b) degree variance (c) percentage of squares (d) percentage of cubes Figure 5: Outcomes of multiple simulations of the modified square model at N = 2000 , α =80 , β = − , σ = 100, with initial states in the form of Erd˝os-R´enyi graphs whose edgeoccupation probability is given on the horizontal axes. The green crosses indicate the meanof each vertical column.adjacency matrix q ij , they become a sum of terms each of which only depends on one singleentry of q . The same holds for the α -term of the two-start Hamiltonian as a function of c ij .For the modified triangle Hamiltonian (11), the same holds with respect to the matrix whoseentries are c ij q ij .This relative simplicity of the dependence on the adjacency matrix allows for a convenientrewriting of the partition functions of our models in terms of integrals, rather than discretesums. A simpler version of the same procedure was applied in [17] to the two-star modelwithout any extra terms, resulting in a representation that elucidates the large N behavior ofthe model. We shall give a demonstration below how to implement this sort of transformationfor the modified square model. A similar (but more complicated) procedure must work for themodified triangle model. For the square model, in order to perform the discrete summationover the components of c , one must first introduce the square of c as an explicit summationvariable. This is done as follows.The partition function of the modified square model is Z ( α, β, σ ) = (cid:88) { c } exp (cid:104) − H ∗ ( c ) + σ (cid:88) i 00 800 16008910111213 Figure 6: Diameter values of the modified triangle model plotted as a function of the numberof vertices in a double-logarithmic presentation. At each graph size, 50 simulations have beenperformed. The area of the orange dots represents the number of hits this particular valueof diameter has received in the course of our simulations. The blue line is the best fit thatcorresponds to the diameter varying as N . .This quantity is expected to display various power laws depending on the scale (in relationto N ), and is connected to the Hausdorff dimension [28]. One may also consider the spectraldimension, which arises from analysis of diffusion processes on graphs [28]. To efficientlyemploy these notions, however, one must have at one’s disposal graphs much bigger thanwhat straightforward Monte Carlo simulations allow (we have considered graphs with afew thousand vertices). Another important quantity that connects graphs and manifolds,frequently seen in the recent literature, is the Ollivier curvature [6–8, 29–33].Perhaps the most basic quantifier of the geometric properties of graphs is the dependenceof the graph diameter (the length of the shortest path between a given pair of vertices,maximized over the choice of this pair) on N , the number of vertices in the graph. Thisquantity would grow like a power of N for regular lattices or regular simplicial complexes, butit is logarithmic in N for sparse non-geometric random graphs, as the simplest Erd˝os-R´enyimodel. Of course, distinguishing power laws from a logarithm is subtle for the limited rangeof N that we access, and furthermore diameter evaluation becomes computationally costlyfor large graphs. (It can be expressed through matrix multiplication and is therefore relatedto the numerical cost of multiplication of large matrices [39], a particular approach is theSeidel algorithm [40].) With all of these constraints, we nonetheless present a crude estimateof the diameter dependence on N , visible in Fig. 6. This fit suggests that the diameter varieslike N . ± . and thus displays a power-law dependence on N . Some related question have been considered in the random graph literature, such as the statistics ofdiameters of connected components [34,35] in the Erd˝os-R´enyi model, and the shortest path statistics [36–38]. Conclusions We have designed two exponential random graph models obtained by modifying the conven-tional two-star random graph Hamiltonian [17, 19] with either (7) or (11). In an appropriateparameter regime where the new terms are given substantial weight, the models displaypronounced geometric features, while avoiding the collapse or fracturing that plagues naivemodels with Hamiltonians based on subgraph counts, such as the Strauss model [12, 15, 17].We quantify the geometric features of our graphs by counting triangles, hexagons, squaresand cubes. These motifs appear in very significant numbers, a sizable fraction of the maxi-mal number possible in the relevant degree range. Such a behavior should be considered anemergent phenomenon, since the graph Hamiltonian explicitly depends only on properties ofpaths of length 2 and knows nothing about the bigger structures. By contrast, in the morefamiliar simpler models, as the Erd˝os-R´enyi or the two-star model, the probability of findingeven one hexagon or cube in a large graph of a finite ( N -independent) mean degree is simplyzero. In the range of graph sizes we can access, the diameter of our graphs appears to scaleas a power law of the number of vertices, similarly to lattices and discretized spaces, and incontrast to the logarithmic scaling of non-geometric exponential random graphs.There is a number of interesting properties of our models that would deserve furtherexploration. Fig. 1 suggest a phase transition separating geometric graphs that we havefocused upon here, and non-geometric graphs reminiscent of the simple two-star model. Itwould be very interesting to elucidate the nature of this transition, and, more ambitiously,use it to construct a continuum limit of our models and connect them to the topics ofrandom continuous geometry. All of these goals would require working with much biggergraphs, which is also the case for the topics of large-scale geometry outlined in the lastsection, and this remains a challenge. Similarly, it would be interesting to investigate furtherthe analytic representation of our models in terms of real fields on a complete graph, andconnect them to topics of 1 /N expansions. This would perhaps cast some further light onthe nature of jamming visible on Figs. 3 and 5, which suggest a rich landscape of metastableconfigurations in our models. Acknowledgments During the long course of formulation and development of this project, we have benefittedfrom conversations on related topics with Peter Grassberger, Christy Kelly, Dmitri Kri-oukov, Surachate Limkumnerd, Yong Lin and Bob Ziff. We especially thank Eytan Katzavand Luca Lionni for discussions about geometry of random graphs and comments on themanuscript. Research of O.E. has been supported by the CUniverse research promotioninitiative (CUAASC) at Chulalongkorn University. Research of T.C. has been supportedby the Program Management Unit for Human Resources and Institutional Development,Research and Innovation (grant number B05F630108). Numerical simulations have beenperformed using the computational platforms of the Extreme Condition Physics ResearchLaboratory and Chula Intelligent and Complex Systems, Chulalongkorn University. Graphvisualizations are prepared using Gephi software [41].17 eferences [1] S. Surya, The causal set approach to quantum gravity, Living Rev. Rel. (2019) 5arXiv:1903.11544 [gr-qc].[2] S. Wolfram, A new kind of science (Wolfram Media, 2002).[3] S. Wolfram, A class of models with the potential to represent fundamental physics, Compl. Syst. (2020) 107 arXiv:2004.08210 [cs.DM].[4] J. Ambjørn, B. Durhuus and Th. Jonsson, Quantum geometry: a statistical field theoryapproach (CUP, 2008).[5] J. Ambjørn, A. G¨orlich, J. Jurkiewicz and R. Loll, Nonperturbative quantum gravity .Phys. Rep., 519 (2012) 127 arXiv:1203.3591 [hep-th].[6] C. A. Trugenberger, Combinatorial quantum gravity: geometry from random bits, JHEP (2017) 045 arXiv:1610.05934 [hep-th].[7] C. Kelly and C. A. Trugenberger, Combinatorial quantum gravity: emergence of geomet-ric space from random graphs, J. Phys. Conf. Ser. (2019) 012016 arXiv:1811.12905[hep-th].[8] C. Kelly, C. A. Trugenberger and F. Biancalana, Self-assembly of geometric space fromrandom graphs, Class. Quant. Grav. (2019) 125012 arXiv:1901.09870 [gr-qc]; Con-vergence of combinatorial gravity, arXiv:2102.02356 [gr-qc].[9] M. Newman, Networks: an introduction (Oxford, 2010).[10] T. Coolen, A. Annibale and E. Roberts, Generating random networks and graphs (Oxford, 2017).[11] E. T. Jaynes, On the rationale of maximum-entropy methods, Proc. IEEE (1982) 939https://bayes.wustl.edu/etj/articles/rational.pdf[12] D. Strauss, On a general class of models for interaction , SIAM Rev. (1986) 513.[13] G. Bizhani, M. Paczuski and P. Grassberger, Discontinuous percolation transitionsin epidemic processes, surface depinning in random media, and Hamiltonian randomgraphs , Phys. Rev. E (2012) 011128 arXiv:1202.3136 [cond-mat.dis-nn].[14] D. Krioukov, Clustering implies geometry in networks , Phys. Rev. Lett. 116 (2016)208302 arXiv:1604.01575 [cond-mat.stat-mech].[15] T. A. B. Snijders, P. E. Pattison, G. L. Robins and M. S. Handcock, New specificationsfor exponential random graph models , Socio. Meth. (2006) 99.[16] Z. Burda, J. Jurkiewicz and A. Krzywicki, Network transitivity and matrix models , Phys.Rev. E (2004) 026106 arXiv:cond-mat/0310234.1817] J. Park and M. E. J. Newman, Solution of the two-star model of a network , Phys. Rev.E (2004) 066146 arXiv:cond-mat/0405457.[18] J. Park and M. E. J. Newman, Solution for the properties of a clustered network , Phys.Rev. E (2005) 026136 arXiv:cond-mat/0412579.[19] A. Annibale and O. T. Courtney, The two-star model: exact solution in the sparseregime and condensation transition , J. Phys. A (2015) 365001 arXiv:1504.06458[cond-mat.dis-nn].[20] A. Gorsky and O. Valba, Finite-size effects in exponential random graphs , J. Compl. Net. (2020) cnaa008 arXiv:1905.03336 [cond-mat.dis-nn]; Interacting thermofield doublesand critical behavior in random regular graphs, arXiv:2101.04072 [hep-th].[21] S. Chatterjee and P. Diaconis, Estimating and understanding exponential random graphmodels , Ann. Stat. (2013) 2428 arXiv:1102.2650 [math.PR].[22] C. Radin and M. Yin, Phase transitions in exponential random graphs , Ann. App. Prob. (2013) 2458 arXiv:1108.0649 [math.PR].[23] F. Harary and B. Manvel, On the number of cycles in a graph, Mat. ˇCas. Sloven. Akad.Vied (1971) 55.[24] E. W. Weisstein, Graph cycle at MathWorldhttp://mathworld.wolfram.com/GraphCycle.html.[25] T. A. B. Snijders, Markov chain Monte Carlo estimation of exponential random graphmodels, J. Soc. Struc. (2002).[26] J. Ambjørn, J. Jurkiewicz and R. Loll, Emergence of a 4-D world from causal quantumgravity, Phys. Rev. Lett. (2004) 131301 arXiv:hep-th/0404156.[27] L. Glaser and S. Steinhaus, Quantum gravity on the computer: impressions of a work-shop, Universe (2019) 35 arXiv:1811.12264 [gr-qc].[28] B. Durhuus, Hausdorff and spectral dimension of infinite random graphs, Acta Phys.Polon. B (2009) 3509.[29] Y. Ollivier, Ricci curvature of Markov chains on metric spaces, J. Func. Anal. (2009) 810 arXiv:math/0701886 [math.PR].[30] Y. Lin, L. Lu, and S.-T. Yau, Ricci curvature of graphs, Tohoku Math. J. (2011)605.[31] D. Cushing and S. Kamtue, Long-scale Ollivier-Ricci curvature of graphs, Anal. Geom.Metr. Sp. (2019) 22 arXiv:1801.10131 [math.CO].[32] P. van der Hoorn, W. J. Cunningham, G. Lippner, C. Trugenberger and D. Krioukov, Ollivier-Ricci curvature convergence in random geometric graphs, arXiv:2008.01209[math.PR]. 1933] P. van der Hoorn, G. Lippner, C. Trugenberger and D. Krioukov, Ollivier curvature ofrandom geometric graphs converges to Ricci curvature of their Riemannian manifolds, arXiv:2009.04306 [math.PR].[34] T. (cid:32)Luczak, Random trees and random graphs, Rand. Struct. Alg. (1998) 485.[35] A. K. Hartmann and M. M´ezard, Distribution of diameters for Erd˝os-R´enyi randomgraphs, Phys. Rev. E (2018) 032128 arXiv:1710.05680 [cond-mat.dis-nn].[36] E. Katzav, M. Nitzan, D. ben-Avraham, P. L. Krapivsky, R. K¨uhn, N. Ross and O. Bi-ham, Analytical results for the distribution of shortest path lengths in random networks, EPL (2015) 26006 arXiv:1504.00754 [cond-mat.dis-nn].[37] M. Nitzan, E. Katzav, R. K¨uhn and O. Biham, Distance distribution in configurationmodel networks, Phys. Rev. E (2016) 062309 arXiv:1603.04473 [cond-mat.dis-nn].[38] E. Katzav, O. Biham and A. K. Hartmann, The distribution of shortest path lengthsin subcritical Erd˝os-R´enyi networks, Phys. Rev. E (2018) 012301 arXiv:1806.05743[cond-mat.dis-nn].[39] F. R. K. Chung, Spectral graph theory ∼ fan/research/cb/ch3.pdf.[40] R. Seidel, On the all-pairs-shortest-path problem in unweighted undirected graphs,