{'fraction_non_alphanumeric': 0.029598349809294, 'fraction_numerical': 0.009476920681871254, 'mean_word_length': 4.384994236613224, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The statistical mechanics of networks

We study the family of network models derived by requiring the expected properties of a graph ensemble to match a given set of measurements of a real-world network, while maximizing the entropy of the ensemble. Models of this type play the same role in the study of networks as is played by the Boltzmann distribution in classical statistical mechanics; they offer the best prediction of network properties subject to the constraints imposed by a given set of observations. We give exact solutions of models within this class that incorporate arbitrary degree distributions and arbitrary but independent edge probabilities. We also discuss some more complex examples with correlated edges that can be solved approximately or exactly by adapting various familiar methods, including mean-field theory, perturbation theory, and saddle-point expansions.


I. INTRODUCTION
The last few years have seen the publication of a large volume of work in the physics literature on networks of various kinds, particularly computer and information networks like the Internet and world wide web, biological networks such as food webs and metabolic networks, and social networks [1,2,3,4]. This work has been divided between empirical studies of the structure of particular networks and theoretical studies focused largely on the creation of mathematical and computational models. The construction of network models is the topic of this paper.
Models of networks can help us to understand the important features of network structure and the interplay of structure with processes that take place on networks, such as the flow of traffic on the Internet or the spread of a disease over a social network. Most network models studied in the physics community are of a practical sort. Typically one wishes to create a network that displays some feature or features observed in empirical studies. The principal approach is to list possible mechanisms that might be responsible for creating those features and then make a model incorporating some or all of those mechanisms. One then either examines the networks produced by the model for rewarding similarity to the realworld systems they are supposed to mimic, or uses them as a substrate for further modeling, for example of traffic flow or disease spread. Classic examples of models of this kind are the small-world model [5] and the many different preferential attachment models [6,7,8], which model network transitivity and power-law degree distributions respectively. However, there is another possible approach to the modeling of networks, which has been pursued comparatively little so far. An instructive analogy can be made here with theories of gases. There are (at least) two different general theories of the properties of gases. Kinetic theory explicitly models collections of individual atoms, their motions and collisions, and attempts to calculate overall properties of the resulting system from basic mechanical principles. Pressure, for instance, is calculated from the mean momentum transfered to the walls of a container by bombarding atoms. Kinetic theory is well motivated, easy to understand, and makes good sense to physicists and laymen alike. However, kinetic theory rapidly becomes complex and difficult to use if we attempt to make it realistic by the inclusion of accurate intermolecular potentials and similar features. In practice, kinetic theory models either make only rather rough and uncontrolled predictions, or they rely on large-scale computer simulation to achieve accuracy.
If one wants a good calculational tool for studying the properties of gases, therefore, one does not use kinetic theory. Instead, one uses statistical mechanics. Although certainly less intuitive, statistical mechanics is based on rigorous probabilistic arguments and gives accurate and reliable answers for an enormous range of problems, including many, such as problems concerning solids, for which kinetic theory is inapplicable. Equilibrium statistical mechanics provides a general framework for reasoning and a powerful calculational tool for very many problems in statistical physics.
Here we argue that the current commonly used models of networks are akin to kinetic theory. They posit plausible mechanisms or dynamics, and produce results in qualitative agreement with reality, at least in some respects. They are easy to understand and give us good physical insight. However, like kinetic theory, they do not make quantitatively accurate predictions and provide no overall framework for modeling, each network model instead concentrating on explaining one or a few features of the system of interest.
In this paper we discuss exponential random graphs, which are to networks as statistical mechanics is to the study of gases-a well-founded general theory with true predictive power. These advantages come at a price: exponential random graphs are both mathematically and conceptually sophisticated, and their understanding demands some effort of the reader. We believe this effort to be more than worthwhile, however. Theoretical techniques based on solid statistical foundations and capable of quantitative predictions have been of extraordinary value in the study of fluid, solid state, and other physical systems, and there is no reason to think they will be any less valuable for networks.
We are by no means the first authors to study exponential random graphs, although our approach is different from that taken by others. Exponential random graphs were first proposed in the early 1980s by Holland and Leinhardt [9], building on statistical foundations laid by Besag [10]. Substantial further developments were made by Frank and Strauss [11,12,13], and continued to be made by others throughout the 1990s [14,15]. In recent years a number of physicists, including ourselves, have made theoretical studies of specific cases [16,17,18,19,20,21]. Today, exponential random graphs are in common use within the statistics and social network analysis communities as a practical tool for modeling networks and several standard computer tools are available for simulating and manipulating them, including Prepstar, ERGM, and Siena [22].
In this paper we aim to do a number of things. First, we place exponential random graph models on a firm physical foundation, showing that they can be derived from first principles using maximum entropy arguments. In doing so, we argue that these models are not merely an ad hoc formulation studied primarily for their mathematical convenience, but a true and correct extension of the statistical mechanics of Boltzmann and Gibbs to the network world.
Second, we take an almost entirely analytic approach in our work, by contrast with the numerical simulations that form the core of most previous studies. We show that the analytic techniques of equilibrium statistical mechanics are ideally suited to the study of these models and can shed much light on their structure and behavior. Throughout the paper we give numerous examples of specific models that are solvable either exactly or approximately, including several that have a long history in network analysis. Nonetheless, the particular examples studied in this paper form only a tiny fraction of the possibilities offered by this class of models. There are many intriguing avenues for future research on exponential random graphs that are open for exploration, and we highlight a number of these throughout the paper.

II. EXPONENTIAL RANDOM GRAPHS
The typical scenario addressed in the creation of a network model is this: one has measurements of a number of network properties for a real-world network or networks, such as number of vertices or edges, vertex degrees, clustering coefficients, correlation functions, and so forth, and one wishes to make a model network that has the same or similar values of these properties. For instance, one might find that a network has a degree sequence with a power-law distribution and wish to create a model network that shows the same power law. Or one might measure a high clustering coefficient in a network and wish to build a model network with similarly high clustering.
Essentially all models considered in modern work, and indeed as far back as the 1950s and 1960s, have been ensemble models, meaning that a model is defined to be not a single network, but a probability distribution over many possible networks. We adopt this approach here as well. Our goal will be to choose a probability distribution such that networks that are a better fit to observed characteristics are accorded higher probability in the model.
Consider a set G of graphs. One can use any set G , but in most of the work described in this paper G will be the set of all simple graphs without self-loops on n vertices. (A simple graph is a graph having at most a single edge between any pair of vertices. A self-loop is an edge that connects a vertex to itself.) Certainly there are many other possible choices and we consider some of the others briefly in Sections III D and III F. The graphs can also be either directed or undirected and we consider both in this paper, although most of our time will be spent on the undirected case.
Suppose we have a collection of graph observables {x i }, i = 1 . . . r, that we have measured in empirical observation of some real-world network or networks of interest to us. We will, for the sake of generality, assume that we have an estimate x i of the expectation value of each observable. In practice it is often the case that we have only one measurement of an observable. For instance, we have only one Internet, and hence only one measurement of the clustering coefficient of the Internet. In that case, however, our best estimate of the expectation value of the clustering coefficient is simply equal to the one measurement that we have.
Let G ∈ G be a graph in our set of graphs and let P (G) be the probability of that graph within our ensemble. We would like to choose P (G) so that the expectation value of each of our graph observables {x i } within that distribution is equal to its observed value, but this is a vastly underdetermined problem in most cases; the number of degrees of freedom in the definition of the probability distribution is huge compared to the number of constraints imposed by our observations. Problems of this type however are commonplace in statistical physics and we know well how to deal with them. The best choice of probability distribution, in a sense that we will make precise in a moment, is the one that maximizes the Gibbs entropy subject to the constraints plus the normalization condition G P (G) = 1.
Here x i (G) is the value of x i in graph G.
Introducing Lagrange multipliers α, {θ i }, we then find that the maximum entropy is achieved for the distribution satisfying for all graphs G. This gives or equivalently where H(G) is the graph Hamiltonian and Z is the partition function Equations (6) to (8) define the exponential random graph model. The exponential random graph is the distribution over a specified set of graphs that maximizes the entropy subject to the known constraints. It is also the exact analogue for graphs of the Boltzmann distribution of a physical system over its microstates at finite temperature.
Using the exponential random graph model involves performing averages over the probability distribution (6). The expected value of any graph property x within the model is simply The exponential random graph, like all such maximum entropy ensembles, gives the best prediction of an unknown quantity x, given a set of known quantities, Eq. (2). In this precise sense, the exponential random graph is the best ensemble model we can construct for a network given a particular set of observations. In many cases we may not need to perform the sum (9); often we need only perform the partition function sum, Eq. (8), and the values of other sums can then be deduced by taking appropriate derivatives. Just as in conventional equilibrium statistical mechanics, however, performing even the partition function sum analytically may not be easy. Indeed in some cases it may not be possible at all, in which case one may have to turn to Monte Carlo simulation, to which the model lends itself admirably. As we show in this paper however, there are a variety of tools one can employ to get exact or approximate analytic solutions in cases of interest, including mean-field theory, algebraic transformations, and diagrammatic perturbation theory.

III. SIMPLE EXAMPLES
Before delving into the more complicated calculations, let us illustrate the use of exponential random graphs with some simple examples.

A. Random graphs
Consider first what is perhaps the simplest of exponential random graphs, at least for the case of fixed number of vertices n considered here.
Suppose we know only the expect number of edges m that our network should have. In that case the Hamiltonian takes the simple form We can think of the parameter θ as either a field coupling to the number of edges, or alternatively as an inverse temperature. Let us evaluate the partition function for this Hamiltonian for the case of an ensemble of simple undirected graphs on n vertices without self-loops. We define the adjacency matrix σ to be the symmetric n × n matrix with elements Then the number of edges is m = i<j σ ij , and the partition function is It is convenient to define the free energy which in this case is (Note that the free energy is extensive not in the number of vertices n, but in the number n 2 of pairs of vertices, since this is the number of degrees of freedom in the model.) Then, for instance, the expected number of edges in the model is Conventionally we re-express the parameter θ in terms of so that m = n 2 p. The probability P (G) of a graph in this ensemble can be written In other words, P (G) is simply the probability for a graph in which each of the n 2 possible edges appears with independent probability p.
This model is known as the Bernoulli random graph, or often just the random graph, and was introduced, in a completely different fashion, by Solomonoff and Rapoport [23] in 1951 and later famously studied by Erdős and Rényi [24,25]. Today it is one of the best studied of graph models, although, as many authors have pointed out, it is not a good model of most real-world networks [1,3,5]. One way in which its inadequacy shows, and one that has been emphasized heavily in networks research in the last few years, is its degree distribution. Since each edge in the model appears with independent probability p, the degree of a vertex, i.e., the number of edges attached to that vertex, follows a binomial distribution, or a Poisson distribution in the limit of large n. Most real-world networks however have degree distributions that are far from Poissonian, typically being highly right-skewed, with a small proportion of vertices having very high degree. Some of the most interesting networks, including the Internet and the world wide web, appear to have degree distributions that follow a power law [6,26,27]. In the next section we discuss what happens when we incorporate observations like these into our models.

B. Generalized random graphs
Suppose then that rather than just measuring the total number of edges in a network, we measure the degrees of all the vertices. Let us denote by k i the degree of vertex i. The complete set {k i } is called the degree sequence of the network. Note that we do not need to specify independently the number of edges m in the network, since m = 1 2 i k i for an undirected graph. The exponential random graph model appropriate to this set of observations is the model having Hamiltonian where we now have one parameter θ i for each vertex i. Noting that k i = j σ ij , this can also be written Then the partition function is and the free energy is More generally we could specify a Hamiltonian with a separate parameter Θ ij coupling to each edge [11]. Then This allows us for example to calculate the probability of occurrence p ij of an edge between vertices i and j: The model of Eq. (18) is the special case in which Θ ij = θ i +θ j and the normal (Bernoulli) random graph of Eq. (12) corresponds to the case in which the parameters Θ ij are all equal. Sometimes it is convenient to specify not a degree sequence but a probability distribution over vertex degrees. This can be achieved by specifying an equivalent distribution over the parameters θ i in (18). Let us define ρ(θ) dθ to be the probability that the parameter θ for a vertex lies in the range θ to θ + dθ. Then, averaging over the disorder so introduced, the free energy, Eq. (21), becomes The part of this free energy due to a single vertex with field parameter θ is and the expected degree of vertex i with field θ i is the derivative of this with respect to θ, evaluated at θ i : By a judicious choice of ρ(θ) we can then produce the desired degree distribution. (See also Sec. III E.) We studied this model in a previous paper [17], as a model for degree correlations in the Internet and other networks. We could alternatively specify a probability distribution ρ(Θ) for the parameters Θ ij in (22) that couple to individual edges. Or, taking the developments a step further, one could define joint distributions for the Θ ij on different edges, thereby introducing correlations of quite general kinds between the edges in the model. There are enormous possibilities to be explored in this regard, but we pass over them for now, our interests in the present paper lying in other directions.
One can calculate many other properties of our models. For example, for the model of Eq. (18), one can calculate the expectation value of any product of vertex degrees from an appropriate derivative of the partition function: Such derivatives are correlation functions of degrees within the model. Similarly, derivatives of the free energy give the connected correlation functions: and so forth. For instance, the two-vertex connected correlation is For the case of the Bernoulli random graph, which has all θ i equal, this gives k i k j c = p(1 − p) for i = j, where we have made use of Eq. (16). Thus the degrees of vertices in the random graph are in general positively correlated. One can understand this as an effect of the one edge that potentially connects the two vertices i and j. The presence or absence of this edge introduces a correlation between the two degrees. (For a sparse graph, in which p = O(n −1 ), the correlation disappears in the limit of large graph size.) In order to measure some quantities within exponential random graph models, it may be necessary to introduce additional terms into the Hamiltonian. For instance, to find the expectation value of the clustering coefficient C [5], one would like to evaluate which we can do by introducing an extra term linear in the clustering coefficient in the Hamiltonian. To measure clustering in the network of Eq. (18), for example, we could define Then Thus it is important, even in the simplest of cases, to be able to solve more general models, and much of the rest of the paper is devoted to the development of techniques to do this.

C. Directed graphs
Before we look at more complicated Hamiltonians, let us look briefly at what happens if we change the graph set G over which our sums are performed. The first case we examine is that of directed graphs. We define G to be the set of all simple loopless directed graphs, which is parameterized by the asymmetric adjacency matrix Thus, for instance, the Hamiltonian H = θm gives rise to a partition function and a corresponding free energy. The directed equivalent of the more general model of Eq. (18) in which we can control the degree of each vertex is a model that now has two separate parameters for each vertex, θ in i and θ out i , that couple to the in-and out-degrees: Then the partition function and free energy are From these we can calculate the expected in-and outdegree of a vertex: We note that i k in , as must be the case for all directed graphs, since every edge on such a graph must both start and end at exactly one vertex.
We can also define a probability distribution ρ(θ in , θ out ) for the fields on the vertices, and the developments generalize Eqs. (25)(26)(27) in a natural fashion.
We give a more complex example of a directed graph model in Section IV C 1, where we derive a solution to the reciprocity model of Holland and Leinhardt [9] using perturbative methods.

D. Fermionic and bosonic graphs
It will by now have occurred to many readers that results like Eqs. (21) and (27) bear a similarity to corresponding results from traditional statistical mechanics for systems of non-interacting fermions. We can look upon the edges in our networks as being like particles in a quantum gas and pairs of vertices as being like singleparticle states. Simple graphs then correspond to the case in which each single-particle state can be occupied by at most one particle, so it should come as no surprise that the results look similar to a system obeying the Pauli exclusion principle.
Not all networks need have only a single edge between any pair of vertices. Some can have multiple edges or multiedges. The world wide web is an example-there can be and frequently is more than one link from one page to another. The Internet, airline networks, metabolic networks, neural networks, citation networks, and collaboration networks are other examples of networks that can exhibit multiedges. There is no problem generalizing our exponential random graphs to this case and, as we might expect, it gives rise to a formalism that resembles the theory of bosons.
Let us define our set of graphs G to be the set of all undirected graphs with any number of edges between any pair of vertices (but still no self-edges, although there is no reason in principle why these cannot be included as well). Taking for example the Hamiltonian, Eq. (22), and generalizing the adjacency matrix, Eq. (11), so that σ ij is now equal to the number of edges between i and j, we have and The equivalent of the probability p ij of an edge appearing in the fermionic case is now the expected number of edges n ij between vertices i and j, which is given by Note that this quantity diverges if we allow Θ ij → 0, a phenomenon related to Bose-Einstein condensation in ordinary Bose gases.
For the special cases of Eqs. (18) and (10), we have and respectively. The connected correlation between the degrees of any two vertices in the latter case is for i = j. Thus the degrees are again positively correlated and the correlation diverges as θ → 0.

E. The sparse or classical limit
In most real-world networks the number of edges m is quite small. Typically m is of the same order as n, rather than being of order n 2 . Such graphs are said to be sparse. (One possible exception is food webs, which appear to be dense, having m = O(n 2 ) [28].) The probability p ij of an edge appearing between any particular vertex pair (i, j) is of order 1/n in such networks. Thus, for example, in the fermionic case of the network described by the Hamiltonian (22), Eq. (24) tells us that e Θij must be of order n in a sparse graph. The same is also true for the bosonic network of the previous section. This allows us to approximate many of our expressions by ignoring terms of order 1 by comparison with terms of order e Θij . We refer to such approximations as the "sparse limit" or the "classical limit," the latter by analogy with the corresponding phenomenon in quantum gases at low density.
In particular, the equivalent of Eq. (24) for either fermionic or bosonic graphs in the classical limit is p ij = e −Θij . For the case of Eq. (18), it is so that each edge appears with a probability that is a simple product of "fugacities" e −θi defined on each vertex. The classical limit of this model has been studied previously by a number of other authors [17,29,30,31,32], although again developed and justified in a different way from our presentation here; generally the edge probability (47) has been taken as an assumption, rather than a derived result.
For a given distribution ρ(θ) of θ, the expected degree of a vertex, Eq. (27), is which is simply proportional to e −θi . So we can produce any desired degree distribution by choosing the corresponding distribution for θ.

F. Fixed edge counts
Another possible choice of graph set G is the set of graphs with both a fixed number of vertices n and a fixed number of edges m. Models of this kind have been examined occasionally in the literature [16] and, if we once more adopt the view of the edges in a graph as particles, they can be considered to be the canonical ensemble of network models, where the variable edge-count models of previous sections are the grand canonical ensemble. As in conventional statistical mechanics, the grand ensemble is often simpler to work with than the canonical one, but progress can be made sometimes be made in the canonical case by performing the sum over all graphs regardless of edge count and introducing a Kronecker δ-symbol into the partition function to impose the edge constraint: where m is the desired number of edges. For instance, the fixed edge-count version of the generalized random graph, Eq. (22), would be one in which where we have made use of the integral representation for the δ-function The sum over graphs is now in the form of the partition function for the grand canonical version of the model, but with Θ ij → Θ ij + 2πiη, giving the field parameters an imaginary part. Thus, from Eq.
In general the integral cannot be done in closed form, which is why fixed edge-count graphs-and canonical ensembles in general-are avoided. The integral can in principle be carried out term by term for any finite n, but doing so is tantamount to performing the sum over all graphs with m edges explicitly, so there is little to be gained by the exercise. It is also possible to have a bosonic graph with a fixed number of edges-one would simply sum over the set of graphs that have m edges with any number of them being permitted to fall between any given pair of vertices.
We will not discuss further either fixed edge-counts or bosonic networks in this paper, concentrating instead on the grand canonical fermionic ones, which are more useful overall. However, essentially all of the results reported in the remainder of the paper can be generalized, with a little work, to these other cases if necessary.

IV. MORE COMPLEX HAMILTONIANS
Outside of the models described in the previous sections, and some minor variations on them, we know of few other exponential random graph models that are exactly solvable. (One exception is the reciprocity model of Holland and Leinhardt [9], for which we derive an exact solution in Sec. IV C 1.) To make further progress one must turn to approximate methods. There are (at least) three types of techniques that can yield approximate analytic solutions for exponential random graph models. The first and simplest is mean-field theory, which works well in many cases because of the intrinsically high dimensionality of network models; usually these models have an effective dimensionality that increases with the number of vertices n, so that the thermodynamic limit of n → ∞ also corresponds to the high dimension limit in which mean-field theory becomes accurate. Nonetheless, there are many quantities, such as those depending on fluctuations, about which mean-field theory says nothing, and for these other methods are needed. In some cases one can use non-perturbative approaches based on the Hubbard-Stratonovich transform or similar integral transforms, which are very effective and accurate but suitable only for models with Hamiltonians of specific forms polynomial in the adjacency matrix. More generally, one can use perturbation theory, which may involve larger approximations (although they are usually well controlled), but is applicable to Hamiltonians of essentially any form.
We discuss all of these approaches here. As an example of their application, we use one of the oldest and beststudied of exponential random graphs, the 2-star model. The Hamiltonian for the 2-star model is where m is the number of edges in the network and s is the number of "2-stars." A 2-star is two edges connected to a common vertex. (The minus sign in front of the parameter α is introduced for later convenience.) The quantities m and s can be rewritten in terms of the degree sequence thus: Substituting these expressions into Eq. (53), we can rewrite the Hamiltonian as where J = 1 2 (n−1)α and B = − 1 2 (θ+α). (The factor (n− 1) in the definition of J is also introduced for convenience later on.) Noticing once again that k i = j σ ij , where the variables σ ij are the elements of the adjacency matrix, we can also write We study the 2-star model in the fermionic case in which each vertex pair can be connected by at most a single edge, and within the grand canonical ensemble where the total number of edges is not fixed. Generalization to the other cases described above is of course possible, if not always easy.

A. Mean-field theory
The variables σ ij can be thought of as Ising spins residing on the edges of a fully connected graph, and hence the 2-star model can be thought of as an Ising model on the edge-dual graph of the fully connected graph [21]. (The edge-dual G * of a graph G is the graph in which each edge in G is replaced by a vertex in G * and two vertices in G * are connected by an edge if the corresponding edges in G share a vertex.) Using this equivalence, the mean-field theory of the 2-star model can be developed in exactly the same way as for the more familiar lattice-based Ising model.
We begin by writing out all terms in Eq. (56) that involve a particular spin σ ij : where we have explicitly taken account of all the ways in which σ ij can enter the first term in the Hamiltonian. (We have also dropped the term 2Jσ ij /(n−1) required to correctly count the terms diagonal in σ ij , since it vanishes in the large n limit.) Then, in classic mean-field fashion, we approximate the local field by its average: where, as before, p = σ is the mean probability of an edge between any pair of vertices, which is also called the connectance of the graph. Then H(σ ij ) = −(4Jp + 2B)σ ij , and we can write a self-consistency condition for p of the form Rearranging, this then gives us p = e 4Jp+2B 1 + e 4Jp+2B = 1 2 tanh(2Jp + B) + 1 .
For J ≤ 1 this equation has only one solution, but for J > 1 there may either be one solution or, if B is sufficiently close to −J, there may be three, of which the outer two are stable. Thus when B is close to −J we have a bifurcation at J c = 1, a continuous phase transition to a symmetry broken state with two phases, one of high density and one of low. We show in Fig. 1 a plot of the solution of (60) which displays clearly the characteristic hysteresis loop of the symmetry broken state.
Along the "symmetric line" B = −J there is always a solution p = 1 2 to Eq. (60) (although it may be unstable), and along this line we can think of p − 1 2 as a standard order parameter for the model which is zero in the highsymmetry phase and non-zero in the symmetry-broken phase. We can define a critical exponent β in the usual fashion by as J → 0 + , giving β = 1 2 , which is the usual Ising meanfield value and should come as no surprise, given the equivalence mentioned above between the 2-star model and the Ising model. One can define other critical exponents as well, which are also found to take Ising meanfield values. For instance, as we showed in [18], the variance χ of the connectance, which plays the role of a susceptibility, goes as χ ∼ |J − 1| −γ in the vicinity of the phase transition with γ = 1.

B. Non-perturbative approaches
We can go beyond the mean-field approximation of the previous section by making use of techniques borrowed from many-body theory. The developments of this section follow closely the lines of our previous paper on this topic [18], and rather than duplicate material needlessly, the reader is referred to that paper for details of the calculation. Here we merely summarize the important results.
The evaluation of the partition function for the 2-star model involves a sum of terms of the form e k 2 . The study of interacting quantum systems has taught us that such sums can be performed using the Hubbard-Stratonovich transform. We start by noting the well-known result for the Gaussian integral: Making the substitutions a → (n − 1)J and φ → φ i − k i /(n − 1), and rearranging, this becomes Then the partition function is where we have interchanged the order of sum and integral. The sum over graphs now has precisely the form of the partition function sum for the model of Eq. (18), and from Eq. (20) we can thus immediately write down the partition function where the quantity is called the effective Hamiltonian. Thus we have completed the partition function sum for the 2-star model, but at the expense of introducing the auxiliary fields {φ i } which must be integrated out to complete the calculation. The integral cannot, as far as we are aware, be evaluated exactly in closed form but, as we showed in [18], it can be evaluated approximately using a saddle-point expansion, with the result that the free energy of the 2-star model is given to leading order in the expansion by where is the position of the saddle point, i.e., the maximum of the Hamiltonian on the real-φ line. Note that Eq. (68) is identical to the mean-field equation, Eq. (60), for the connectance p of the 2-star model. Thus, φ 0 is the connectance of the model within the mean-field approximation and the saddle-point expansion, as is typically the case in such calculations, is an expansion about the mean-field solution.
From the free energy we can derive a number of quantities of interest. We showed in [18], for instance, that the variance of vertex degree in the model is given by which has a gradient discontinuity but no divergence at the phase transition. (This quantity is, by contrast, zero everywhere within mean-field theory.)

C. Perturbation theory
Exponential random graphs also lend themselves naturally to treatment using perturbation theory. Here we describe the simplest such theory, which is roughly equivalent to the high-temperature expansions of conventional thermal statistical mechanics. Expansions of this type have been examined previously by Burda et al. [19,20] for Strauss's model of a transitive network [12,13]. Here we develop the theory further for general exponential random graphs.
The fundamental idea of perturbation theory for random graphs is the same as for other perturbative methods: we expand about a solvable model in powers of the coupling parameters θ i in the Hamiltonian. We write the Hamiltonian for the full model in the form H = H 0 + H 1 , where H 0 is the Hamiltonian for the solvable model and H 1 takes whatever form is necessary to give the correct expression for H. Then the partition function is [19,20] where Z 0 = G e −H0 is the partition function for the unperturbed Hamiltonian, and . . . 0 indicates an ensemble average in the unperturbed model. The only case that has been investigated in any detail is the one where we expand around a random graph, H 0 = θm, so that the averages in Eq. (70) are averages in the ensemble of the random graph. (It is possible for θ to be zero, so this choice for H 0 does not place any restriction on the form of the overall Hamiltonian. If θ = 0 then the expansion is precisely equivalent to an ordinary high-temperature series.) However, for Hamiltonians H that give significant probability to networks substantially different from random graphs, the perturbation theory cannot be expected to give accurate answers at low order. In theory there is no reason why one could not expand about some other solvable case, although no such calculations have been done as far as we are aware. One obvious possibility, which we do not pursue here, is to expand around one of the generalized random graph forms, Eqs. (18) and (22).
Typically, to make progress with Eq. (70), we will expand the exponential in a power series of the form In practice, H 1 normally contains a coupling constant, such as the constant J in the 2-star model of Eq. (55), and thus our expression for the perturbed partition function is an expansion in powers of the coupling.
In this section, we apply the perturbation method to two example models. First, we study a simple model proposed about a quarter of a century ago by Holland and Leinhardt [9], which is exactly solvable by this method. Then we illustrate the application of the method to the 2-star model and compare its performance against the approximate saddle-point expansion results of the previous section.

The reciprocity model
Our first example of perturbation theory is a directed graph model. In the real world, many directed graphs display the phenomenon of reciprocity: a directed edge running from vertex A to vertex B predisposes the network to have an edge running from B to A as well. Put another way, the network has a higher fraction of vertex pairs that are joined in both directions than one would expect on the basis of chance ("mutual dyads" in the parlance of social network analysis). Behavior of this kind is seen, for example, in the world wide web, email networks, and neural and metabolic networks [33,34,35].
Holland and Leinhardt [9] have proposed a exponential random graph model of reciprocity, which we study here in a simplified version. As we now show, the perturbation expansion for this model can be written down to all orders and resummed to give an exact expression for the partition function.
The Hamiltonian for the model is where m is the total number of (directed) edges in the graph, and r is the number of vertex pairs with edges running between them in both directions. The unperturbed Hamiltonian H 0 is that of an undirected random graph (Sec. III C) with partition function given by Eq. (35) and each directed edge present with independent probability p = (e θ + 1) −1 . The perturbation H 1 can be written in terms of the adjacency matrix (34) as Then the perturbation series, Eq. (71), for the full Hamiltonian is with Thus the partition function is written as an expansion in powers of α whose coefficients are correlation functions of elements of the adjacency matrix, calculated within the ordinary random graph. If we can evaluate these correlation functions, at least up to some finite order, we can also evaluate the perturbed partition function. Since all edges are present or absent independently of one another in the random graph, the correlation functions factor: σ 12 σ 21 σ 34 σ 43 0 = σ 12 0 σ 21 0 σ 34 0 σ 43 0 = p 4 , (76) and so forth. The only exception is in cases where two or more of the elements σ ij being averaged are the same. In that case, noting that σ n ij = σ ij for any n, we have results like To evaluate expressions such as (75), therefore, we need to count the number of independent elements σ ij that appear in each term. This can be difficult for some models, but for the reciprocity model it is quite straightforward. The question we need to answer is this: if we choose k pairs of vertices (i, j) from the n Eqs. (74) and (75) where the function is the exponential generating function for a k,q . Now the number of ways of choosing k pairs such that all choices are made from a particular set of size q, but without the constraint that each pair in the set appear at least once, is just q k . Thus q m=1 q m a k,m = q k , or equivalently Multiplying by z k /k! and summing over k = 1 . . . ∞, this gives which immediately implies that g q (z) = (e z − 1) q , by induction on (81) with the initial condition g 1 (z) = ∞ k=1 z k a k,1 /k! = e z − 1. Substituting this result into Eq. (78) then gives us our solution: Or, making use of Eqs. (12) and (16) and F = − ln Z in the normal fashion. From these expressions we can, for instance, obtain the mean number of edges m and the mean number r of pairs of vertices connected by edges running both ways from A quantity of interest in directed networks is the reciprocity [34], which is the fraction of edges that are reciprocated. This quantity is found to be on the order of tens In Fig. 2, we show the reciprocity, along with the connectance of the network, as a function of α for the case p = 0.01. There is no phase transition or other unexpected behavior in this model: the measured properties are smooth functions of the independent parameters. Notice that there is a substantial range of values of α over which the connectance is low and the graph realistically sparse, but the reciprocity is still high, with values similar to those seen in real networks.

Example 2: The 2-star model
As our second example of the application of perturbation theory, we return to the 2-star model introduced at the beginning of Sec. IV. Unlike the case of the reciprocity model in the preceding section, perturbation theory does not lead to an exact solution of the 2-star model but, as we now show, we can get an approximate solution by studying the perturbation expansion to finite order-a different approximation from the saddle-point expansion of Sec. IV B.
We divide the Hamiltonian H = θm−αs into an unperturbed part H 0 = θm, which is the normal Bernoulli random graph, and a perturbation Hamiltonian H 1 = −αs. Then, following Eq. (71), the partition function for the full model is given by The number of 2-stars is and therefore Our strategy is to evaluate the series (86) up to some finite order in α to get an approximate solution for Z, but there is a problem. Each term in the series corresponds to states of the graph that have the corresponding number of 2-stars: the term in s 0 , for instance, counts the number of graphs that have a 2-star in any position in the graph. This is not enough for our purposes however. Realistic graphs will have not a finite number but a finite density of 2-stars in them, and the number of such graphs is counted by terms that appear at infinite order in the perturbation expansion in the limit n → ∞. So, without going to infinite order as we did in the reciprocity model, we are never going to get meaningful results from our expansion.
Similar problems appear in ordinary statistical mechanics and the solution is well known. Instead of expanding the partition function, we form an expansion for the free energy. We can write the free energy as where F 0 is the free energy of the unperturbed network and F 1 = − ln(Z/Z 0 ). Now we expand F 1 as a power series in α of the form where we have made use of the fact that F 1 = 0 when α = 0. Substituting into Z/Z 0 = e −F1 , we get and comparing terms with Eq. (86), we find and so forth. These are the cumulants of s within the ensemble defined by the unperturbed network. If we expand s in the form of Eq. (87) then they are connected correlations of elements of the adjacency matrix-"connected" because individual elements of the adjacency matrix are uncorrelated, so that all terms in the cumulants vanish unless they involve sets of 2-stars that share one or more edges. (Note that sharing a vertex, as in the more familiar spin models of traditional statistical mechanics, is not a sufficient condition for being connected. The fundamental degrees of freedom in a network are the edges.) We will proceed then as follows. We calculate the free energy F 1 in terms of connected correlations up to some finite order in α and from this we calculate the partition function Z = Z 0 e −F1 . Even though F 1 is known only to finite order, our expression for Z will include terms with all powers of the connected correlations in it, via the expansion of the exponential, and hence will include graphs with not only a finite number but a finite density of 2-stars. This idea, which will be routine for those familiar with conventional diagrammatic many-body theory, is entirely general and can be applied to any model, not just the 2-star model. In essence, the series given by e −F1 is a partial resummation to all orders of the partition function, including some but not all of the contributions to Z from disconnected correlations of arbitrarily high order.
Let us see how the calculation proceeds for the case of the 2-star model, to order α 3 , as above. The leading O(α) term in F 1 is simple: Since we are primarily interested in large networks, we can approximate this expression by its value to leading order in n, which is 1 2 n 3 p 2 . The second term, at order α 2 , is more complicated because there are several different ways in which two 2stars may combine to share one or more edges. In order to keep track of these different contributions, we make use of a diagrammatic representation similar to that employed by Burda et al. for Strauss's transitivity model [19]. Figure 3a shows the single diagram contributing to f 1 , which gives the result in Eq. (93). Figure 3b shows the three diagrams that contribute to f 2 . It is an assumption of our notation that each edge that appears in a diagram is distinct. Thus the third diagram in Fig. 3b, which represents the case in which the two 2-stars fall on top of one another, must be depicted separately, rather than being considered a special case of the first diagram. This turns out to be a good idea, since this term has a different functional form from the first diagram, and neither diagram is necessarily negligible by comparison with the other.
In general the basic "Feynman rules" for interpreting the diagrams are: 1. each edge contributes a factor of p; 2. each vertex contributes a factor of n; 3. the numerical multiplier is the number of distinct ways in which the diagram can be decomposed into overlapping 2-stars such that each edge occurs at least once, divided by the symmetry factor for the diagram. (The symmetry factor is the number of distinct permutations of the vertices that leave the diagram unchanged.) Then for the connected correlation functions one must subtract all other ways of composing lower order diagrams to make the given diagram, as in Eq. (92). To see how these rules work in practice, let us apply them to the first diagram in Fig. 3b. This diagram has four vertices and three edges, which gives a factor of n 4 p 3 , by the first two rules. The diagram can be decomposed into two 2-stars in 6 different ways, but the symmetry factor is also 6, so we end up with n 4 p 3 × 6/6 = n 4 p 3 . The contribution to the diagram from the term −f 2 1 in Eq. (92b) is −n 4 p 4 , so the final value of the diagram is n 4 (p 3 − p 4 ) to leading order in n. Proceeding in a similar fashion, the other diagrams of Fig. 3b contribute n 4 (p 3 − p 4 ) and 1 2 n 3 (p 2 − p 4 ), respectively. The diagrams for the O(α 3 ) term are shown in Fig. 3c, and are more complicated, but routine to evaluate using the rules above. The final expressions for the f s are: Note that we have retained the leading order terms in n separately at each order in p, since we have no knowledge a priori about the relative magnitude of n and p. In a sparse graph, we expect that p will be of order 1/n, in which case it may be possible to neglect some terms. Once we have the expansion of F 1 , it is straightforward to calculate statistical averages from derivatives of the free energy in the normal fashion. For example, the expected number of 2-stars in the network is given by And the expected number of edges is m = ∂F ∂θ = p(p − 1) ∂F 0 ∂p + ∂F 1 ∂p = 1 2 n 2 p + n 3 (1 − p)p 2 α 1 + 1 2 (1 + 6np − 8np 2 )α + 1 6 (1 + 21np + 58n 2 p 2 − 180n 2 p 3 + 129n 2 p 4 )α 2 . In Fig. 4, we show the connectance 2 m /n 2 and the density of 2-stars 2 s /n 3 calculated from the saddlepoint method of Sec. IV B and from the expressions above, at first, second, and third order. As the figure shows, the perturbation expansion agrees with the non-perturbative method at high and low values of J = 1 2 (n − 1)α, and markedly better for the third-order approximation than for the first-and second-order ones. However, in the region of the phase transition at J c = 1 the agreement is poor, as we would expect. In this region there will be large critical fluctuations and hence contributions to the free energy from large connected diagrams that are entirely missing from our series expansion. Presumably by extending the perturbation series we can derive successively more accurate answers in the critical region. We also note that the perturbation expansion gives results only for the sparse phase in the symmetry-broken region.
We have here studied in detail two examples of the treatment of exponential random graphs by perturbation theory (and another can be found in Ref. [19]). The techniques we have used, however, are entirely general and diagrammatic theories similar to these, with similarly simple "Feynman rules," can be derived for other examples as well.

V. CONCLUSIONS
In this paper we have discussed exponential random graphs, which in both a figurative and a quantitative sense play the role of a Boltzmann ensemble for the study of networks. Exponential random graphs are a formally well-founded framework for making predictions about the expected properties of networks given specific measurements of properties of those networks. We have shown in this paper how they can be derived in moderately rigorous fashion from maximum entropy assumptions about probability distributions over graph ensembles.
We have given many examples of particular calculations using exponential random graphs, starting with simple random graph models that have linear Hamiltonians, many of which have been presented previously by other authors, albeit it with rather different motivation. In most cases these linear models can be solved exactly, meaning that we can derive the partition function or equivalently the free energy of the graph ensemble exactly in the limit of large system size.
For nonlinear Hamiltonians it appears possible to find exact solutions only rarely, but we have been able to find approximation solutions in several cases using a number of different methods. Taking the particular example of the 2-star model, we have shown how its behavior can be understood using mean-field theory, perturbation theory, and non-perturbative methods based on the Hubbard-Stratonovich transform. We have also given one example, the reciprocity model of Holland and Leinhardt, that is exactly solvable by evaluating its perturbation expansion to all orders.
The results presented in this paper are only a tiny fraction of what can be done with exponential random graphs. There are many interesting challenges, both practical and mathematical, posed by this class of models. Exploration of the behavior and predictions of specific models as functions of their free parameters, development of other approximate solution methods, or expansion of those presented here, and the development of models to study network phenomena of particular interest, such as vertex-vertex correlations, effects of hidden variables, effects of degree distributions, and transitivity, are all excellent directions for further research. We hope to see some of these topics pursued in the near future.