-3.5ex 2.3ex Introduction

Many quantum and statistical systems simplify when the number of degrees of freedom becomes large, for example the number of internal degrees of freedom of fields at each spacetime point, or the size of matrices, etc. One is often talking about ‘large \(N\)’ limits [1], where \(N\) parametrizes the number of degrees of freedom. A ‘holy grail’ problem for such large \(N\) limits is the search for ‘master variables’ or ‘master fields’, such that, when the original problem is reformulated through these variables, \(N\) becomes a numerical parameter in the action (or in the statistical probability distribution), and changing \(N\) no longer affects the set of the degrees of freedom involved. As this parameter tends to infinity, the large \(N\) limit is often controlled by an explicit saddle point in terms of the master variables, which condense, as a result, to definite values.

The master variables can be of a different nature than the original variables. In the simplest case of the \(O(N)\) vector model whose fundamental variables \(\phi_i\) are vectors under the \(O(N)\) rotation, scalar fields constructed from these vectors condense to definite values at large \(N\) and play the role of master variables [1]. Similar ‘scalar condensation’ via different applications of the Hubbard-Stratonovich transformation occurs in multi-matrix models with a large number of matrices transformed into one another by rotations [2]–[4], or in computations of random tensor eigenvalue distributions in a large number of dimensions [5]–[7]. The situation can be much more exotic, however. For example, for large random matrices of size \(N\times N\), the resolvents, which are functions of an external spectral parameter, condense to definite values [8], [9]. Thus, while the original variables are numbers, the master variables are functions. This is especially visible in the corresponding treatments of sparse random matrix ensembles [10]–[15]: starting with random matrices, one arrives at explicit functional integrals in terms of the functional master variables that acquire a saddle point structure at large \(N\). Finally, in quantum gauge theories, loop variables make a prominent appearance at large \(N\) and are believed to condense to definite values [16]–[19]. Those are functionals on the loop space, while the original variable are fields in the physical spacetime.

Motivated by this diversity of large \(N\) limits and the master variables that control them, we would like to revisit here the two-star random graph model. This is probably the simplest graph model with nontrivial edge interactions and nontrivial thermodynamics, formulated as a Gibbs-like maximal entropy ensemble with controlled average number of edges and ‘2-stars’ (which one can think of as the number of paths of length 2). By contrast, if one only controls the number of edges, one ends up with the edges being filled or not filled randomly and independently, which is known as the Erdős-Rényi random graph. This model is of fundamental importance in the physics of networks [20], providing a point of departure for more sophisticated constructions, but its thermodynamics is trivial (this is for the case of distinguishable vertices; thermodynamics of unlabeled graphs with indistinguishable vertices is nontrivial even within the Gibbs ensemble that only controls the number of edges [21], [22]).

We thus turn to the two-star random graph on \(N\) vertices and look for a set of master variables that control the large \(N\) limit. One generally expects that the master variables are invariants of the symmetries of the theory, and here we must keep in mind that, unlike \(O(N)\) models, the matrix models defining random graph ensembles are typically invariant under vertex permutations, but not under any form of continuous rotations, leaving a much bigger set of invariants.

In the dense regime (a finite fraction of all possible edges occupied), a mean field solution of the two-star model was developed by Park and Newman [23]. This solution is simple, elegant and manifestly correct given the comparisons with numerics, but it does not bring under explicit control the \(1/N\) corrections, which is where the master variables proposed here lead to improvement. In the sparse regime (a finite average degree), the model was solved using auxiliary fields by Annibale and Courtney [24]. This solution, however, goes through functional integrals at the intermediate stages. The master variables proposed here will recover it in a more compact and elementary fashion. (We mention for additional perspective the mathematical works [25], [26] dealing with related topics, as well as some further relevant physics literature on the two-star model and its generalizations [27], [28].)

Our treatment will proceed with defining the relevant master variables in section [sec:secmaster], showing how they can be used to reproduce the Park-Newman mean field solution in the dense regime in section [sec:secPN], spelling out the computation of \(1/N\) corrections to this solution in section [sec:sec1N], followed by a very compact derivation of the large \(N\) solution in the sparse case in section [sec:secsparse]. We will conclude with a brief summary and discussion.

-3.5ex 2.3ex The two-star random graph and its auxiliary field representation

For analytic considerations, one starts by representing graphs as adjacency matrices. For each graph on \(N\) vertices, define a real symmetric zero-diagonal \(N\times N\) adjacency matrix \(A_{ij}\) whose entries equal 1 if there and edge between vertices \(i\) and \(j\), and 0 otherwise. The probability of each graph in the two-star model is then defined by \[\label{2stardef} P[A]=e^{-H[A]}/Z,\qquad Z\equiv\sum_A e^{-H[A]},\qquad H[A]\equiv \alpha\sum_{kl}{A_{kl}}+\beta\sum_{klm} A_{kl}A_{lm},\tag{1}\] where \(\sum_A\) implies summing over all possible adjacency matrices, that is, over all graph configurations. This is a typical Gibbs-like maximal entropy ensemble with two ‘fugacities’ \(\alpha\) and \(\beta\). Introducing the vertex degrees \(d_k\equiv \sum_l A_{kl}\), one can understand \(\sum_{kl}{A_{kl}}=\sum_k d_k\) as twice the number of edges, while \(\sum_{klm} A_{kl}A_{lm}=\sum_k d_k^2\). These are the two quantities whose expectation values are controlled by their thermodynamic conjugates \(\alpha\) and \(\beta\). By linear redefinitions of \(\alpha\) and \(\beta\), one can equivalently control instead the number of edges and the number of ‘two-star’ motifs (ordered triplets of distinct nodes with at least two edges connecting them), given by \(\sum_k d_k(d_k-1)\). While all of these definitions are physically equivalent, the conventions explicitly stated in (1 ) prove convenient for our subsequent analytic work.

To solve the ensemble (1 ) in the thermodynamic limit \(N\to\infty\), one needs to evaluate the partition function \(Z\). Direct summation over \(A\) is impossible due to the quadratic nonlinearity in \(H\). We start with a simple Hubbard-Stratonovich transform, as in [23]: \[e^{-\beta\sum_{km}A_{kl}A_{lm}}=\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty d\phi \,e^{-\phi^2 +2i\phi \sqrt{\beta}\sum_k A_{kl}}.\] We need to introduce one such variable \(\phi\) for each value of \(l=1..N\), obtaining \[Z= \frac{1}{\pi^{N/2}}\sum_{\mathbf{A}}\int_{-\infty}^\infty d\phi_1\ldots d\phi_N\,e^{-\sum_l\phi_l^2}e^{\sum_{k<l} A_{kl}[-2\alpha+2i(\phi_k+\phi_l). \sqrt{\beta}]},\] where we expressed the exponent through the independent entries \(A_{kl}\) with \(k<l\). The summation over \(A_{kl}\) with \(k<l\) then trivializes since the summand is factorized over the entries of the adjacency matrix. One thereby arrives at the following vector model in terms of \(\phi_l\): \[Z=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty d\phi_1\ldots d\phi_N\,e^{-\sum_k\phi_k^2}\prod_{k<l}\left(1+e^{-2\alpha+2i(\phi_k+\phi_l) \sqrt{\beta}}\right),\] which is often more convenient to represent as \[\label{vecmoddef} Z=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty d\phi_1\ldots d\phi_N\,e^{-\sum_k\phi_k^2\,+\,S[\phi]},\qquad S[\phi]\equiv\sum_{k<l}\log\left(1+e^{-2\alpha+2i(\phi_k+\phi_l) \sqrt{\beta}}\right).\tag{2}\]

One expects that such vector models simplify when \(N\) is large. This is what in fact happens, though the details depend crucially on the \(N\)-scalings of the thermodynamic parameters \(\alpha\) and \(\beta\) that determine whether the graph is in the sparse or dense regime (both of which are an option at large \(N\), defining thus different large \(N\) limits of the model).

A mean field solution of (2 ) in the dense regime was originally presented in [23]. The dense regime corresponds to \(\alpha=O(1)\), \(\beta=O(1/N)\) at \(N\to\infty\). In this regime, \(\phi_k\) condense to definite values of order \(\sqrt{N}\) (in the current conventions) plus fluctuations of order 1. The idea of [23] is then to expand around this condensation point as \(\phi_k=\sqrt{N}\phi_0+\varphi_k\) so that, with a Taylor expansion for \(S\), \[\label{phiTay} -\sum_k\phi_k^2+S[\phi]=S_0+\sum_k \phi_k^2 +\frac{1}{2}\sum_{kl}\frac{\partial^2 S}{\partial\phi_k\partial\phi_l}\varphi_k\varphi_l+\frac{1}{6}\sum_{klm}\frac{\partial^3 S}{\partial\phi_k\partial\phi_l\partial\phi_m}\varphi_l\varphi_l\varphi_m+\cdots.\tag{3}\] The mean field value \(\phi_0\) is chosen to ensure that the term linear in \(\varphi\) is absent. As an estimate for \(Z\) of (2 ), one then takes \(e^{S_0}\) times the Gaussian determinant from integrating the exponential the quadratic terms in (3 ), with all the higher order terms neglected.

All the indications are that the results of the above mean field analysis are exact at \(N\to\infty\), and they compare very well with the corresponding numerics [23]. And yet, why could the higher-order terms in (3 ) be neglected? Most naively, differentiating the \(\log\) in \(S\) produces negative powers of \(\phi\), and since \(\phi\) is of order \(\sqrt{N}\) at the point of expansion in (3 ), the higher-order terms come with manifest negative powers of \(N\), which is reassuring. This is, however, not yet the whole story. The number of variables grows with \(N\), and in the standard Feynman-like expansion while evaluating (2 ) with (3 ) substituted, one will have multiple sums over indices ranging from 1 to \(N\). This will introduce positive powers of \(N\) that will compete with the negative powers. (The ability of fluctuations of a large number of variables to compromise the \(1/N\) expansions is highlighted in a similar context in section 8.2 of [16].) More than that, if we repeatedly differentiate \(S\) with respect to the same component of \(\phi\), the resulting expression has \(N\) terms, contributing further positive factors of \(N\).

Due to all the above ingredients, the story with higher-order corrections to [23], which we briefly sketch in Appendix [sec:appPNexp], appears somewhat bewildering. While the naive evaluation within the Gaussian approximation produces the correct result at large \(N\) that can be verified numerically, the higher-order corrections come with a swarm of positive and negative powers of \(N\), and it is difficult to give a concise argument as to how exactly they are suppressed, and why they do not upset the Gaussian estimate (though it is known empirically that they do not). While one could attempt diagrammatic accounting of the powers of \(N\) in these corrections, somewhat in the spirit of random tensor considerations [29], [30], it is anything but easy.

Instead of attempting diagrammatic analysis of (2 ) with the expansion (3 ), we shall follow here a different strategy, introducing a further set of auxiliary variables, somewhat akin to those recently used for analyzing graphs with prescribed degree sequences [15], [31], [32]. In this way, conventional saddle points will emerge that control the large \(N\) limit of (2 ), where \(N\) becomes a numerical parameter in the ‘action,’ and the set of integration variables is \(N\)-independent. This will have other useful applications, besides elucidating the analytics behind the mean field solution of [23].

We first write (2 ) as \[Z=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty \left(\prod_k \frac{d\phi_k\,e^{-\phi_k^2}}{\sqrt{1+e^{-2\alpha+4i\phi_k \sqrt{\beta}}}}\right) \prod_{k,l=1}^N\exp\left[\frac{1}{2}\log\left(1+e^{-2\alpha+2i(\phi_k+\phi_l) \sqrt{\beta}}\right)\right],\] and then expand the logarithm as \(\log(1+x)=x-x^2/2+x^3/3+\cdots\) to obtain \[\begin{align} Z&=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty \left(\prod_k \frac{d\phi_k\,e^{-\phi_k^2}}{\sqrt{1+e^{-2\alpha+4i\phi_k \sqrt{\beta}}}}\right) \prod_{k,l=1}^N\exp\left[-\sum_{J=1}^\infty\frac{(-1)^J}{2J} e^{J(-2\alpha+2i(\phi_k+\phi_l) \sqrt{\beta})}\right]\nonumber\\ &=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty \left(\prod_k \frac{d\phi_k\,e^{-\phi_k^2}}{\sqrt{1+e^{-2\alpha+4i\phi_k \sqrt{\beta}}}}\right) \exp\left[-\sum_{J=1}^\infty\frac{(-1)^Je^{-2\alpha J}}{2J} \Big(\sum_k e^{2iJ\phi_k \sqrt{\beta}}\Big)^{\!2\,}\right]. \end{align}\] An attractive feature of the last expression is that it can be completely factorized over \(\phi_k\) at the cost of introducing one more Hubbard-Stratonovich transformation with respect to the variables \(x_J\) that couple to \(\sum_k e^{2iJ\phi_k \sqrt{\beta}}\). In this way, one gets \[\label{xJHS} Z=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty \left(\prod_{J=1}^\infty \frac{dx_J}{\sqrt{2\pi J}} \,e^{- x_J^2/2J}\right)\prod_k\int_{-\infty}^\infty d\phi_k \frac{e^{-\phi_k^2}\exp\left[\sum_{J=1}^\infty i^{J+1}e^{-J\alpha}e^{2iJ \phi_k \sqrt{\beta}}x_J/J\right]}{\sqrt{1+e^{-2\alpha+4i\phi_k \sqrt{\beta}}}},\tag{4}\] or more suggestively, since all the \(\phi_k\) integrals are identical to each other: \[\begin{align} \tag{5} &Z=\int_{-\infty}^\infty \left( \prod_{J=1}^\infty\frac{dx_J}{\sqrt{2\pi J}} \,e^{- x_J^2/2J}\right)\, e^{NS[x]},\\ &e^{S[x]}=\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty d\phi \frac{e^{-\phi^2}}{\sqrt{1+e^{-2\alpha+4 i\phi\sqrt{\beta}}}}\exp\left[\sum_{J=1}^\infty\frac{i^{J+1}}{J} e^{-J\alpha}x_Je^{2iJ\phi\sqrt{\beta}}\right].\tag{6} \end{align}\] The integral in the last line gives an ‘effective action’ \(S[x]\), in terms of which one has to deal with a saddle point problem at large \(N\), with \(N\) playing no more of a role than providing a large saddle point parameter.

-3.5ex 2.3ex The dense regime at leading order

We now turn to the dense regime of [23], where \(\alpha\sim 1\) and \(\beta\sim 1/N\), which we write as \[\beta=\frac{B}{N}.\] These scalings are necessary to ensure a nontrivial infinite \(N\) limit with a finite mean connectivity (the fraction of all possible edges that are filled).

For the first-pass treatment in this section, as in [23], we only keep the extensive part of the free energy, that is, contributions to \((\log Z)/N\) that survive at \(N\to\infty\). For that, we have compute \(e^S\) up to the order \(O(1)\) at \(N\to\infty\), so that \(e^{NS}\) is computed correctly up to factors that stay finite at \(N\to\infty\). In this way, we reproduce the solution of the model at the precision level of [23], while the subleading corrections will be discussed in the next section.

For constructing a saddle point estimate of (5 6 ), it is important to identify the \(N\)-scaling of the variables responsible for the dominant contribution. In view of the mean field picture of [23] that has been verified by comparisons with numerics, the field \(\phi\) condenses to values of order \(\sqrt{N}\). Since \(x_J\) are Hubbard-Stratonovich conjugates of \(\sum_k e^{iJ\phi_k \sqrt{\beta}}\), they are expected to condense at values of order \(N\), since the sum over \(k\) consists of \(N\) identical terms of order 1. We shall see that this regime \(\phi\sim\sqrt{N}\), \(x_J\sim N\) is indeed consistent and produces a saddle point estimate of \(Z\) with controlled corrections.

Motivated by the above picture, we write \(\phi=i\sqrt{BN}\phi_0+\varphi\). We will choose \(\phi_0\) to ensure that the integral over \(\varphi\) is dominated by \(\varphi\) of order 1. Then, neglecting all terms suppressed by \(1/N\), \[\begin{align} e^{S[x]}=&\frac{e^{NB\phi_0^2}}{\sqrt{\pi\left(1+e^{-2\alpha- 4B\phi_0}\right)}}\int_{-\infty}^\infty d\varphi\,e^{-\varphi^2-2i\sqrt{BN}\phi_0\varphi}\nonumber\\ &\times\exp\left[\sum_{J=1}^\infty\frac{i^{J+1}}{J} e^{-J(\alpha+2B\phi_0)}x_J\left(1+2iJ\varphi\sqrt{B/N}-2J^2\varphi^2B/N\right)\right]. \label{action} \end{align}\tag{7}\] We then choose \(\phi_0\) to force the terms of order \(O(\sqrt{N})\) that are linear in \(\varphi\) cancel out, which ensures that only values of \(\varphi\) of order 1 contribute to the integral: \[\label{phi0eq} \phi_0=\sum_{J=1}^\infty {i^{J+1}} e^{-J(\alpha+2B\phi_0)}\frac{x_J}{N}.\tag{8}\] One is then left with an elementary Gaussian integral over \(\varphi\) that yields \[e^{S[x]}=\frac{e^{NB\phi_0^2}}{\sqrt{1+e^{-2\alpha-4 B\phi_0}}}\exp\left[\sum_{J=1}^\infty\frac{i^{J+1}}{J} e^{-J(\alpha+2B\phi_0)}x_J\right]\left(1+{2B}\sum_{J=1}^\infty {i^{J+1}} Je^{-J(\alpha+2B\phi_0)}\frac{x_J}{N}\right)^{-1/2}.\] The corrections to this formula are of order \(1/N\), which will give at most multiplicative contributions of order \(O(1)\) in \(e^{NS}\), and those cannot contribute to the extensive part of the free energy.

From the above, we write \(S\) as \[\begin{align} \tag{9} &S=NS_{\mathrm{main}}+S_{\mathrm{corr}},\qquad S_{\mathrm{main}}\equiv B\phi_0^2+\sum_{J=1}^\infty\frac{i^{J+1}}{J} e^{-J(\alpha+2B\phi_0)}\frac{x_J}{N}, \\ &S_{\mathrm{corr}}\equiv -\frac{1}{2}\log\left(1+{2B}\sum_{J=1}^\infty {i^{J+1}} Je^{-J(\alpha+2B\phi_0)}\frac{x_J}{N}\right)-\frac{1}{2}\log\left(1+e^{-2\alpha-4 B\phi_0}\right).\tag{10} \end{align}\] where the leading terms of order \(N\) have been separated out explicitly, and \(\phi_0\) is a function of \(x_J\) implicitly determined by (8 ). To construct a saddle point estimate of (5 ), we look for stationary points of \(NS[x]-\sum_Jx_J^2/2J\), determined at leading order in \(N\) by \[\label{sadd95est} \frac{\partial}{\partial x_K}\left(N^2S_{\mathrm{main}}[x]-\sum_J \frac{x_J^2}{2J}\right)=0.\tag{11}\] With the dominant terms in \(S\) collected in (9 ), this yields \[x_K/K=2{N^2B}\phi_0\frac{\partial\phi_0}{\partial x_K}+N\frac{i^{K+1}}{K}e^{-K(\alpha+2B\phi_0)}-2BN\frac{\partial\phi_0}{\partial x_K}\sum_{J=1}^\infty i^{J+1} e^{-J(\alpha+2B\phi_0)}x_J.\] The terms with \(\partial\phi_0/\partial x_K\) cancel out in view of (8 ), leaving the following saddle point configuration: \[\label{Xsddl} x_K=X_K\equiv N i^{K+1} e^{-K(\alpha+2B\phi_0)}.\tag{12}\] Note that \(X_K\) are consistently of order \(N\). At this saddle point, (8 ) becomes an explicit equation for \(\phi_0\) \[\label{phi0atsadd} \phi_0=\sum_{J=1}^\infty {(-1)^{J+1}} e^{-J(2\alpha+4B\phi_0)}=\frac{1}{e^{2\alpha+4B\phi_0}+1},\tag{13}\] as in the mean field solution of [23]. If we expand (5 ) around the saddle point configuration \(X_K\), the Gaussian integral over fluctuations is manifestly of order \(O(1)\) at large \(N\), as we shall see more explicitly in the next section. Then, nonvanishing contributions to the free energy per vertex may only come from the saddle point exponential evaluated at \(X_K\), that is \[\frac{\log Z}{N}=-\sum_J \frac{X_J^2}{2JN} + NS_{\mathrm{main}}[X]+S_{\mathrm{corr}}[X].\] From (9 ), (10 ) and (12 ), this is evaluated as \[\label{PNres} \frac{\log Z}{N}={NB}\phi_0^2+\frac{N-1}{2}\log\left(1+e^{-2\alpha-4B\phi_0}\right)-\frac{1}{2}\log\left[1+{2B}\phi_0(1-\phi_0)\right].\tag{14}\] where we have used the evident formulas \[\sum_{J=1}^\infty \frac{(-1)^{J+1}}{J}z^J=\log(1+z),\qquad \sum_{J=1}^\infty {(-1)^{J+1}}z^J=\frac{z}{1+z},\qquad \sum_{J=1}^\infty J{(-1)^{J+1}}z^J=\frac{z}{(1+z)^2},\] and also equation (8 ) to simplify the terms of order 1. The expression is identical to the result of [23] up to a change of notation: what we call \(\alpha\) here is \(-B\) in [23], and what we call \(B\) is \(-JN/(N-1)\) in [23].

-3.5ex 2.3ex The dense regime: subleading corrections

To keep the story clean, we shall take the saddle point values deduced in the previous section as an input for identifying the relevant point of expansion, and restart the derivation independently departing, once again, from the partition function (5 6 ). We will see that this results in an explicit controllable \(1/N\) expansion that we will process in a matter-of-fact manner.

In (5 6 ) with \(\beta=B/N\), we introduce the following variable redefinitions motivated by the previous section: \[\begin{align} \tag{15} &\phi \equiv i\sqrt{BN}\phi_0 + \varphi,\qquad \phi_0\left(e^{2\alpha+4B\phi_0}+1\right)\equiv 1,\\ &x_K \equiv X_K + \chi_K, \qquad X_K\equiv N i^{K+1} e^{-K(\alpha+2B\phi_0)}. \tag{16} \end{align}\] To avoid clutter in the sums, we define the shorthand \[c_J\equiv\frac{X_J}{JN}=\frac{i^{J+1}}{J} e^{-J(\alpha+2B\phi_0)}.\]

In the previous section, we computed \(S\) up to the order \(O(1)\), and then the free energy is computed up to the order \(O(N)\), since \(S\) enters the expression for \(Z\) as \(e^{NS}\) and one has to take a logarithm of \(Z\). In order to compute the next correction, which is the nonextensive contribution to the free energy of order \(O(1)\), one correspondingly has to compute \(S\) up to the order \(O(1/N)\). To do so, with the above notation, we rewrite (6 ) as \[e^{S[\chi]}=e^{NB\phi_0^2}\int_{-\infty}^\infty d\varphi\,\frac{e^{-\varphi^2-2i\sqrt{BN}\phi_0\varphi}}{\sqrt{\pi\left(1+e^{-2\alpha- 4B\phi_0 }e^{4i\varphi\sqrt{B/N}}\right)}}\exp\left[\sum_{J=1}^\infty c_J(X_J + \chi_J)e^{2iJ\varphi\sqrt{B/N}}\right].\] We expand the denominator in the integrand as \[\left(1+ae^{b/\sqrt{N}}\right)^{-1/2} = (1+a)^{-1/2}\left(1-\frac{a}{2(1+a)}\frac{b}{\sqrt{N}}+\frac{a^2-2a}{8(1+a)^2}\frac{b^2}{N}+\cdots\right),\] and, keeping in mind that \(X_J\sim N\), expand the integrand up to order \(1/N\). Importantly, terms involving \(\varphi \sqrt{N}\) in the exponent cancel out because of our choice of \(\phi_0\) and \(X_J\), and one is left with the following Gaussian integral (where we omit odd powers of \(\varphi\) that would, in any case, integrate to 0): \[\begin{align} e^{S[\chi]}=&\frac{e^{NB\phi_0^2+N\sum_J c_J^2J+\sum_J c_J\chi_J}}{\sqrt{\pi\left(1+e^{-2\alpha- 4B\phi_0}\right)}} \int_{-\infty}^\infty d\varphi \;e^{-\varphi^2\left(1 + {2B}\sum_J c_J^2J^3\right)}\Bigg\{ 1 + O(N^{-2})\\ &+\frac{1}{N}\Bigg[ \frac{2B\varphi^2(2e^{2\alpha+ 4B\phi_0}-1)}{(1+e^{2\alpha+ 4B\phi_0})^2}- 2B\varphi^2\left({\textstyle\sum_J}c_JJ\chi_J\right)^2 - 2B\varphi^2{\textstyle\sum_J}c_JJ^2 \chi_J \\ & + \frac{4B\varphi^2}{1+e^{2\alpha+ 4B\phi_0}}{\textstyle\sum_J}c_JJ\chi_J+ \frac{2B^2\varphi^4}{3}{\textstyle\sum_J}c_J^2 J^5+\frac{8B^2\varphi^4}{3}\textstyle\sum_J c_J J \chi_J \textstyle\sum_Kc_K^2 K^4\\ & -\frac{8B^2\varphi^4}{3(1+e^{2\alpha+ 4B\phi_0})}\,{\textstyle\sum_J}c_J^2 J^4 - \frac{8B^3\varphi^6}{9}\left({\textstyle\sum_J}c_J^2 J^4\right)^2\Bigg]\Bigg\}. \end{align} \label{phintexpand}\tag{17}\] We will often encounter the following sums, for which we introduce explicit notation: \[\label{sigdef} \sigma_n \equiv\sum_{J=1}^{\infty} c_J^2 J^{n+1}= \sum \frac{(-1)^{J+1}}{J}\left(\frac{\phi_0}{1-\phi_0}\right)^{\!\! J}J^n = \left(x\frac{d}{dx}\right)^{\!\! n}\log(1+x)\bigg|_{x=\frac{\phi_0}{1-\phi_0}}.\tag{18}\] Evaluating the \(\varphi\)-integral in (17 ) then yields \[\label{eS} \begin{align} e^{S[\chi]}=&\frac{e^{NB\phi_0^2+N\sigma_0+\sum_J c_J\chi_J}}{\sqrt{(1+e^{-2\alpha- 4B\phi_0})(1+2B\sigma_2)}} \Bigg\{ 1+ O(N^{-2})\\ &+\frac{1}{N}\Bigg[ \frac{B(2e^{2\alpha+ 4B\phi_0}-1)}{(1+e^{2\alpha+ 4B\phi_0})^2(1+2B\sigma_2)}- \frac{B\left({\textstyle\sum_J}c_JJ\chi_J\right)^2}{1+2B\sigma_2} - \frac{B\sum_J c_J J^2 \chi_J}{1+2B\sigma_2} \\ &+ \frac{2B\sum_J c_J J\chi_J}{(1+e^{2\alpha+ 4B\phi_0})(1+2B\sigma_2)}+ \frac{B^2\sigma_4}{2(1+2B\sigma_2)^2}+\frac{2B^2\sigma_3{\textstyle\sum_J}c_JJ\chi_J}{(1+2B\sigma_2)^2} \\ & -\frac{2B^2\sigma_3}{(1+e^{2\alpha+ 4B\phi_0})(1+2B\sigma_2)^2} - \frac{5B^3\sigma_3^2}{3(1+2B\sigma_2)^3}\Bigg]\Bigg\}. \end{align}\tag{19}\] Raising this expression to the power of \(N\), so as to substitute it in (5 ), and making use of \((1+x/N)^N \sim e^x\) at large \(N\), we arrive at \[\begin{align} e^{NS[\chi]}=&\frac{e^{N^2B\phi_0^2+N^2\sigma_0+N\sum_J c_J\chi_J}}{\left[(1+e^{-2\alpha- 4B\phi_0})(1+2B\sigma_2)\right]^{N/2}} \,\,\exp\Bigg[\frac{B(2e^{2\alpha+ 4B\phi_0}-1)}{(1+e^{2\alpha+ 4B\phi_0})^2(1+2B\sigma_2)}\\ &+ \frac{B^2\sigma_4}{2(1+2B\sigma_2)^2} -\frac{2B^2\sigma_3}{(1+e^{2\alpha+ 4B\phi_0})(1+2B\sigma_2)^2} - \frac{5B^3\sigma_3^2}{3(1+2B\sigma_2)^3}\nonumber\\ &+\frac{2B^2\sigma_3{\textstyle\sum_J}c_JJ\chi_J}{(1+2B\sigma_2)^2}+ \frac{2B\sum_J c_J J\chi_J}{(1+e^{2\alpha+ 4B\phi_0})(1+2B\sigma_2)}- \frac{B\sum_J c_J J^2 \chi_J}{1+2B\sigma_2}+ - \frac{B\left({\textstyle\sum_J}c_JJ\chi_J\right)^2}{1+2B\sigma_2}\Bigg],\nonumber \end{align}\] where all possible corrections to this expression are suppressed by powers of \(1/N\). Evaluating \(\sigma_n\) in terms of \(\phi_0\) from (18 ) – the explicit formulas are tabulated in Appendix [sec:appsigma] – and also keeping in mind that \(e^{2\alpha+4B\phi_0} = (1-\phi_0)/{\phi_0}\) and \(1+e^{-2\alpha-4B\phi_0} = {1}/(1-\phi_0)\), we obtain \[\begin{align} e^{NS[\chi]}=&\frac{e^{N^2B\phi_0^2}(1-\phi_0)^{-N^2+N/2}e^{N\sum_J c_J\chi_J}}{[1+2B\phi_0(1-\phi_0)]^{N/2}} \,\exp\Bigg[\frac{B\phi_0(2-3\phi_0)}{1 + {2B}\phi_0(1-\phi_0)} - \frac{2B^2\phi_0^2(1-\phi_0)(1-2\phi_0)}{[1 + {2B}\phi_0(1-\phi_0)]^2}\nonumber\\ & + \frac{B^2 \phi_0(1-\phi_0)(1-6\phi_0+6\phi_0^2)}{2[1 + {2B}\phi_0(1-\phi_0)]^2} - \frac{5B^3\phi_0^2(1-\phi_0)^2(1-2\phi_0)^2}{3[1 + {2B}\phi_0(1-\phi_0)]^3}\nonumber\\ &+\frac{2B\phi_0[1+B(1-\phi_0)]{\textstyle\sum_J}c_JJ \chi_J}{[1 + {2B}\phi_0(1-\phi_0)]^2} - \frac{B\left(\sum_J c_J J\chi_J\right)^2}{1 + {2B}\phi_0(1-\phi_0)} - \frac{B\,\sum_J c_J J^2 \chi_J}{1 + {2B}\phi_0(1-\phi_0)}\Bigg]. \label{eNS} \end{align}\tag{20}\] From this, \[\begin{align} \label{Zintexponent} &-\sum_{J}\frac{x_J^2}{2J} + NS[x] = N^2B\phi_0^2 - \frac{N(N-1)}{2}\log(1-\phi_0) -\frac{N}{2}\log[1+2B \phi_0(1-\phi_0)]\\ &+\frac{B\phi_0(2-3\phi_0)}{1 + {2B}\phi_0(1-\phi_0)} +\frac{B^2\phi_0(1-\phi_0)(1-10\phi_0+14\phi_0^2) }{2[1+2B \phi_0(1-\phi_0)]^{2}}- \frac{5B^3\phi_0^2(1-\phi_0)^2(1-2\phi_0)^2}{3[1 + {2B}\phi_0(1-\phi_0)]^3}\nonumber\\ & +\frac{2B\phi_0[1+B(1-\phi_0)]{\textstyle\sum_J}c_JJ \chi_J}{[1 + {2B}\phi_0(1-\phi_0)]^2} - \frac{B\,\sum_J c_J J^2 \chi_J}{1 + {2B}\phi_0(1-\phi_0)}- \frac{B\left({\textstyle\sum_J}c_JJ\chi_J\right)^2}{1 + {2B}\phi_0(1-\phi_0)}.\nonumber \end{align}\tag{21}\] Importantly, the term involving \(N\sum_J c_J\chi_J\) in the first line of (20 ) cancels out in (21 ). This leaves in (5 ) a Gaussian integral over \(\chi_J\) where all integration variables are only allowed to take values of order 1, leading to an explicit finite result. This reflects the usefulness of our parametrization in (15 16 ). Had we chosen a different expansion point, large terms of the form \(N\chi_J\) would have compromised the usefulness of evaluating the Gaussian integrals at this stage, as they would induce rearrangements in the \(1/N\) corrections.

To evaluate (5 ), it remains to integrate \(e^{-\sum_{J}{x_J^2}/{2J} + NS[x]}\) over \(\chi_J\). To do so, we employ the following Hubbard-Stratonovich transformation for the exponential of the last term in (21 ): \[\label{HSy} \exp\left[-\frac{B\left({\textstyle\sum_J}c_JJ\chi_J\right)^2}{1+2B \phi_0(1-\phi_0)}\right] = \int_{-\infty}^\infty dy \frac{e^{-y^2}}{\sqrt\pi}\exp\left[2yi\sqrt{\frac{B}{1+2B \phi_0(1-\phi_0)}}\,\,{\textstyle\sum_J}c_JJ\chi_J\right].\tag{22}\] It is now straightforward, even if somewhat laborious, to evaluate the remaining Gaussian integral over \(\chi_J\) where the quantities \(\sigma_n\) defined in (18 ) appear once again, and can thereafter be explicitly expressed through \(\phi_0\). The result is \[\begin{align} &\prod_{J=1}^\infty\int_{-\infty}^\infty \frac{d\chi_J\,e^{-{\chi_J^2}/{2J}}}{\sqrt{2\pi J}}\exp{\textstyle\left[c_J J\left( \frac{2B\phi_0[1+B(1-\phi_0)]}{[1 + {2B}\phi_0(1-\phi_0)]^2} - \frac{JB}{1 + {2B}\phi_0(1-\phi_0)}+ 2yi \sqrt{\frac{B}{1+2B \phi_0(1-\phi_0)}}\right)\chi_J \right]}\nonumber\\ &=\exp\Big[{\textstyle\frac{B^2\phi_0(1-\phi_0)(1-10\phi_0+18\phi_0^2)+8B^3\phi_0^3(1-\phi_0)^2(5\phi_0-2)+8B^4\phi_0^4(1-\phi_0)^3(3\phi_0-1)}{2[1 + {2B}\phi_0(1-\phi_0)]^4}}\label{chiint}\\ &{\textstyle +2yi \left(\frac{B}{1+2B \phi_0(1-\phi_0)}\right)^{3/2}\frac{\phi_0(1-\phi_0)(4\phi_0-1)+4B\phi_0^3(1-\phi_0)^2}{1 + {2B}\phi_0(1-\phi_0)} -\frac{2y^2B\phi_0(1-\phi_0)}{1+2B \phi_0(1-\phi_0)}}\Big].\nonumber \end{align}\tag{23}\] We have once again used the expressions for sums of the form \(\sum_J c_J^2 J^{n+1}\) given in Appendix [sec:appsigma]. To undo the Hubbard-Stratonovich transformation in (22 ), we need to multiply the last line with \(e^{-y^2}/\sqrt{\pi}\) and integrate over \(y\). This yields: \[\begin{align} &\int_{-\infty}^\infty \frac{dy}{\sqrt\pi}\exp\Bigg[{\textstyle -y^2\Bigg(1+\frac{2B\sigma_2}{1+2B \phi_0(1-\phi_0)}\Bigg)+2yi \left(\frac{B}{1+2B \phi_0(1-\phi_0)}\right)^{3/2}\frac{\phi_0(1-\phi_0)(4\phi_0-1)+4B^2\phi_0^3(1-\phi_0)^2}{1 + {2B}\phi_0(1-\phi_0)}}\Big]\nonumber\\ &=\exp\left\{-\frac{B^3\phi_0^2(1-\phi_0)^2[4\phi_0-1+4B\phi_0^2(1-\phi_0)]^2}{[1+2B \phi_0(1-\phi_0)]^4[1+4B \phi_0(1-\phi_0)]}\right\}\left(\frac{1+4B\phi_0(1-\phi_0)}{1+2B \phi_0(1-\phi_0)}\right)^{-1/2}. \label{yintres} \end{align}\tag{24}\] Finally, we need to gather all the contributions to \(Z\) of (5 ), consisting of the \(\chi\)-independent part of (21 ), the second \(y\)-independent line of (23 ), and (24 ). This leaves for the free energy \[\begin{align} \label{Forder1} &-F \equiv \log Z = N^2B\phi_0^2 - \frac{N(N-1)}{2} \log(1-\phi_0)-\frac{N}{2}\log\left[1+2\Delta\right]\\\ &+\frac{2\Delta-B\phi_0^2}{1 + 2\Delta} -\frac{5B\Delta^2(1-2\phi_0)^2}{3(1 + 2\Delta)^3}-\frac{B\Delta^2[4\phi_0(1+\Delta)-1]^2}{(1+2\Delta)^4(1+4\Delta)} - \frac{1}{2}\log\frac{1+4\Delta}{1+2\Delta}\nonumber\\ &+\frac{B\Delta\big[1-10\phi_0+16\phi_0^2+2\Delta(1-14\phi_0+24\phi_0^2)+2\Delta^2(1-12\phi_0+20\phi_0^2)\big]}{(1 + 2\Delta)^4}.\nonumber \end{align}\tag{25}\] with \(\Delta\equiv B \phi_0(1-\phi_0)\) introduced for compactness. The first line represents the mean-field result of [23], while the subsequent two lines give the first nontrivial correction that the novel representations developed here allowed us to derive.

As the computations leading to (25 ) are rather convoluted, it is important to implement some independent checks. One such check is provided by the particle-hole duality respected by the two-star model: in the original partition function (1 ), one can change the summation variable from the adjacency matrix \(A\) to the ‘inverted’ adjacency matrix, where all filled edges are replaced by empty ones and vice versa. This manifestly relates partition functions, and hence free energies, at two different values of \(\alpha\) at any given \(N\). As (25 ) is expressed through \(\phi_0(\alpha,B)\) rather than \(\alpha\), one has to deal with the corresponding transformation for \(\phi_0\) that maps it to \(1-\phi_0+O(1/N)\). The subleading \(1/N\) corrections in this transformation mix the different orders of \(N\) in (25 ) and since the different contributions must cancel out in the end, this provides a nontrivial check of the subleading corrections in the last two lines of (25 ). Our formula passes this test. We provide an explicit implementation in Appendix [sec:appdual].

Further validation of (25 ) is provided by comparisons with numerical Monte Carlo sampling of the two-star model. We follow a strategy similar to what we have previously employed in [33] for more complicated related models, reviewed in Appendix [sec:appMC]. In our comparisons between analytics and numerics, we focus on the degree variance. In terms of the free energy \(F\equiv -\log Z\), the averages of the mean degree \(\langle k\rangle \equiv \sum_j \langle d_j\rangle/N=\langle d_1\rangle\) and mean square degree \(\langle k^2\rangle \equiv \sum_j \langle d_j^2\rangle/N=\langle d_1^2\rangle\) are expressed as \[\left<k\right> = \frac{1}{N}\frac{\partial F}{\partial\alpha},\qquad \left<k^2\right> = \frac{\partial F}{\partial B}.\] To compute these derivatives, one also needs the derivatives of \(\phi_0\) with respect to \(\alpha\) and \(B\): \[\frac{\partial\phi_0}{\partial\alpha} = -\frac{2\phi_0(1-\phi_0)}{1 + 4B\phi_0(1-\phi_0)},\qquad \frac{\partial\phi_0}{\partial B} = -\frac{4\phi_0^2(1-\phi_0)}{1 + 4B\phi_0(1-\phi_0)}.\] Then, the variance is \[V\equiv\left<k^2\right> -\left<k\right>=\frac{\partial F}{\partial B} - \left(\frac{1}{N}\frac{\partial F}{\partial\alpha}\right)^2.\] We separate out the contribution due to the extensive part of free energy, given by the first line of (25 ): \[-F_0 \equiv N^2B\phi_0^2 - \frac{N(N-1)}{2} \log(1-\phi_0)-\frac{N}{2}\log\left[1+2B \phi_0(1-\phi_0)\right].\] The contribution to the degree variance \(V\) coming from \(F_0\) is \[\label{F0var} \frac{\partial F_0}{\partial B} - \left(\frac{1}{N}\frac{\partial F_0}{\partial\alpha}\right)^2 = N\frac{\phi_0(1-\phi_0)}{1+2\Delta} - \frac{\phi_0^2[2B\phi_0(\phi_0-2)+2B+1]^2}{(1+2\Delta)^2(1+4\Delta)^2},\tag{26}\] with \(\Delta=B\phi_0(1-\phi_0)\). The leading order variance \[\label{defV0} V_0\equiv N\frac{\phi_0(1-\phi_0)}{1+2\Delta}\tag{27}\] matches the one derived in [23]. We are keeping the subleading order in (26 ) as it will combine with the higher-order corrections. The remaining nonextensive piece of the free energy (25 ) is \[\begin{align} -F_{corr} \equiv&\frac{2\Delta-B\phi_0^2}{1 + 2\Delta} -\frac{5B\Delta^2(1-2\phi_0)^2}{3(1 + 2\Delta)^3}-\frac{B\Delta^2[4\phi_0(1+\Delta)-1]^2}{(1+2\Delta)^4(1+4\Delta)} - \frac{1}{2}\log\frac{1+4\Delta}{1+2\Delta}\nonumber\\ &+\frac{B\Delta\big[1-10\phi_0+16\phi_0^2+2\Delta(1-14\phi_0+24\phi_0^2)+2\Delta^2(1-12\phi_0+20\phi_0^2)\big]}{(1 + 2\Delta)^4}. \label{eq:F95corr} \end{align}\tag{28}\] Since this expression is bulky, it is more practical to use symbolic computation software for handling its derivatives, and we rely on SymPy [34] for this purpose, following by evaluating the resulting expressions for each set of parameters. To verify (25 ) numerically, we then extract point-by-point the degree variance from Monte Carlo simulations, with \(V_0\) subtracted, and compare it with \[\label{defVcorr} V_{corr}\equiv \frac{\partial F_{corr}}{\partial B} - \left(\frac{1}{N}\frac{\partial F_{corr}}{\partial\alpha}\right)^2- \frac{\phi_0^2[2B\phi_0(\phi_0-2)+2B+1]^2}{(1+2\Delta)^2(1+4\Delta)^2}.\tag{29}\]

In Fig. 1, we present the results of this comparison, showing excellent agreement. We opt for the moderate number of vertices \(N=200\) as it allows us to keep the higher order \(1/N\) corrections to the analytics sufficiently small, but at the same time the Monte Carlo equilibration for such moderately sized graphs happens sufficiently fast to attain a very high precision. This high precision is needed since we are measuring a contribution that is small in comparison with the leading order variance \(V_0\).

We have been dealing above with the dominant \(1/N\) correction to the mean field result of [23]. If one is to compute further \(1/N\) corrections to the free energy, the process is completely algorithmic, even if it will become more and more laborious with each next order. Indeed, higher order correction will manifest themselves as contributions of order \(1/N^2\) and smaller in the formula for \(e^S\) on top of what is written explicitly in (19 ). All of these contributions are furthermore polynomial in \(\chi_J\). Then, when writing \(e^{-\sum_{J}{x_J^2}/{2J} + NS[x]}\) in (5 ), these corrections will give contributions suppressed by \(1/N\) or more in the exponent, all of which can be re-expanded as additive contributions polynomial in \(\chi_J\). Thus, one will end up with the same Gaussian integrals as in our treatment above, except that they will have extra polynomial insertions in terms of \(\chi_J\) suppressed by powers of \(N\). All such Gaussian integrals with polynomial insertions can be evaluated using the standard formulas, leaving behind an explicit series in powers of \(1/N\) that will correct our result in (25 ).

-3.5ex 2.3ex The sparse regime

We now turn to the sparse regime, \(2\alpha=\log( N/c)\), \(\beta=O(1)\). In this case, the \(J\)-sum in (6 ) works very differently from the dense case considered above, since \(e^{-J\alpha}\) turns into powers of \(c/N\), introducing stronger and stronger \(1/N\) suppression in contributions from higher \(J\). This leads to considerable simplifications.

One can selfconsistently assume \(x_J\sim \sqrt{N}\), and then, because of the suppression of the contributions from higher \(x_J\) mentioned above, \(S\) will only depend on \(x_1\) at large \(N\). Writing \(x_1=-x\,\sqrt{N/c\,}\), \[\label{sparseZS} Z= \sqrt{ \frac{N}{2\pi}}\int_{-\infty}^\infty dx \,e^{N[S(x,\beta)-x^2/2c]},\qquad S(x,\beta)\equiv\log \frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty d\phi\, \exp\left[-\phi^2+x\, e^{2i\phi\sqrt{\beta}}\right].\tag{30}\] The saddle point is at \(x=X\) satisfying \[\label{sddlsparse} X=c\,\frac{\partial S}{\partial x}\Bigg|_{x=X}.\tag{31}\] The free energy per vertex \(f\) satisfies \[-f\equiv \frac{\log Z}{N}=S(X,\beta)-X^2/2c,\] where \(X\) is the solution of (31 ). It is convenient to represent \(e^{S(x)}\) as a series \[e^{S(x,\beta)}=\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty d\phi\, e^{-\phi^2}\sum_{d=0}^\infty\frac{x^d}{d!} e^{2id\phi\sqrt{\beta}}=\sum_{d=0}^\infty\frac{x^d}{d!} e^{-\beta d^2}.\] Thus, \[\label{Xser} X=c\,\frac{\sum_{d=1}^\infty X^{d-1}e^{-\beta d^2}/(d-1)!}{\sum_{d=0}^\infty X^d e^{-\beta d^2}/d!}.\tag{32}\] The mean degree is \[\langle k \rangle\equiv \frac{\partial f}{\partial\alpha}=\frac{dc}{d\alpha}\frac{\partial f}{\partial c}=-2c\frac{\partial f}{\partial c}=\frac{X^2}{c}+2c\frac{\partial X}{\partial c}\frac{\partial}{\partial X}\left[S(X,\beta)-X^2/2c\right].\] The last term vanishes by the saddle point equation, giving \[\label{kXc} \langle k \rangle=\frac{X^2}{c}.\tag{33}\] At the same time, \[\langle k^2 \rangle\equiv \frac{\partial f}{\partial\beta}=-\frac{\partial S(X,\beta)}{\partial\beta}-\frac{\partial X}{\partial\beta}\frac{\partial}{\partial X}\left[S(X,\beta)-X^2/2c\right].\] Again, due to the saddle point equation, \[\label{k2eq} \langle k^2 \rangle=-\frac{\partial S(X,\beta)}{\partial\beta}=\frac{\sum_{d=1}^\infty d^2 X^d e^{-\beta d^2}/d!}{\sum_{d=0}^\infty X^d e^{-\beta d^2}/d!}.\tag{34}\] Expressing \(X\) from (33 ) and substituting it into (32 ), multiplied by \(X\), and into (34 ), we obtain: \[\begin{align} \tag{35} \langle k \rangle&=\frac{\sum_{d=1}^\infty d\,c^{d/2}\,\langle k \rangle^{d/2}e^{-\beta d^2}/d!}{\sum_{d=0}^\infty c^{d/2}\,\langle k \rangle^{d/2} e^{-\beta d^2}/d!},\\ \langle k^2 \rangle&=\frac{\sum_{d=1}^\infty d^2 c^{d/2}\,\langle k \rangle^{d/2}e^{-\beta d^2}/d!}{\sum_{d=0}^\infty c^{d/2}\,\langle k \rangle^{d/2} e^{-\beta d^2}/d!}.\tag{36} \end{align}\] This is identical to the Annibale-Courtney solution of [24], but now derived using only elementary integrals, and in a few lines rather than a few pages. At any given \(c\) and \(\beta\), (35 ) has to be solved as a nonlinear equation for \(\langle k\rangle\), whereupon \(\langle k^2\rangle\) is computed directly from (36 ).

-3.5ex 2.3ex Conclusions

We have revisited solutions of the two-star random graph model in the thermodynamic limit, in both dense and sparse regime. Previous approaches relied on the variables \(\phi_k\) dating back to [23], which are Hubbard-Stratonovich conjugates of the vertex degrees. We pointed out that considerable empowerment of the formalism results from introducing further variables \(x_J\) as in (4 ). These variables can be thought of as Hubbard-Stratonovich conjugates of \(\sum_k e^{2iJ\phi_k\sqrt{\beta}}/N\). The usage of such exponentials of the fields is ubiquitous in conformal field theories [35] – see [36] for a recent discussion – but to the best of our knowledge they have not appeared previously in studies of random graphs.

In the dense regime, it is known that the variables \(\phi_k\) condense themselves to definite values at large \(N\). It is predictable that all variables \(x_J\) condense to definite values as well, with small fluctuations. Our representation provides for effective control over these fluctuations, with a straightforward bookkeeping of \(1/N\) factors, which is challenging in the language of \(\phi_k\). We thus manage to compute the nonextensive part of the free energy in the dense regime. In the sparse regime, the situation is more peculiar: the variables \(\phi_k\) do not condense – the representation for \(S\) in (30 ) does not have a large \(N\) saddle point structure – but the variables \(x_J\) do, leading to a simple saddle point calculation of the free energy. This calculation is phrased entirely in terms of elementary one-dimensional integrals.

Besides elucidating the analytic structure of the two-star model, we hope that the techniques presented here will usefully transfer to further related considerations. A number of interesting extensions of the two-star model can be seen in the literature [37]–[41], and they provide an attractive avenue for exploration in this regard.

-3.5ex 2.3ex *Acknowledgments

We thank Sergei Nechaev for valuable discussions on closely related topics. OE is supported by the C2F program at Chulalongkorn University and by NSRF via grant number B41G680029.

-3.5ex 2.3ex Naive expansion around the Park-Newman solution

To illustrate the power of our formalism, we would like to briefly contrast it with naive expansion around the Park-Newman solution [23]. In our notation, the representation for the partition function used for obtaining that solution is given by (2 ), which we reproduce here for convenience \[\label{vecmoddefapp} Z=\frac{1}{\pi^{N/2}}\int_{-\infty}^\infty d\phi_1\ldots d\phi_N\,e^{-\sum_k\phi_k^2\,+\,S[\phi]},\qquad S[\phi]\equiv\sum_{k<l}\log\left(1+e^{-2\alpha+2i(\phi_k+\phi_l) \sqrt{B/N}}\right).\tag{37}\]

To derive the mean field solution, one assumes that the dominant configuration is site-independent, \(\phi_k=i\sqrt{BN}\phi_0\), where \(\phi_0\) satisfies the same equation we derived in section [sec:secPN]: \[\label{phi0app} \phi_0=\frac{1}{e^{2\alpha+4B\phi_0}+1}.\tag{38}\] One then expands as \(\phi_k=i\sqrt{BN}\phi_0+\varphi_k\) and writes \[-\sum_k\phi_k^2+S[\phi]=S_0+\sum_k \varphi_k^2 +\frac{1}{2}\sum_{kl}\frac{\partial^2 S}{\partial\phi_k\partial\phi_l}\varphi_k\varphi_l+\frac{1}{6}\sum_{klm}\frac{\partial^3 S}{\partial\phi_k\partial\phi_l\partial\phi_m}\varphi_l\varphi_l\varphi_m+\cdots.\] Equation (38 ) ensures that there are no terms linear in \(\varphi_k\) in this expansion. Then, the exponential of all of these terms starting with the cubic one should be expanded in a Taylor series, so that a Gaussian integral over \(\varphi_k\) emerges from (37 ) with an infinite series of polynomial insertions.

To get a sense of how this expansion works, we write out the general formulas for the derivative tensors \[\frac{\partial^nS}{\partial\phi_{j_1}\dots\partial\phi_{j_n}} = \sum_{k<l} \frac{d^n}{dx^n}\log\left(1+e^x\right)\bigg|_{x={-2\alpha+2i(\phi_k+\phi_l) \sqrt{B/N}} }\prod_{m=1}^n\left(2i\sqrt{\frac{B}{N}}(\delta_{kj_m} + \delta_{lj_m})\right).\] Notice the Kronecker symbols involving only the two indices \(k\), \(l\) and one external index \(j_m\), so that the derivative tensors vanish unless there are at most 2 distinct values present among the indices \(j_m\). The first few orders evaluated at the saddle point are \[\begin{align} \frac{\partial^2 S}{\partial\phi_k\partial\phi_l}&= \begin{cases} -\frac{4B}{N} \phi_0 (1-\phi_0) & \text{for } i \neq j,\\ -4B\frac{(N-1)}{N} \phi_0 (1-\phi_0) & \text{for } i = j; \end{cases}\\ \frac{\partial^3 S}{\partial\phi_k\partial\phi_l\partial\phi_m}&= \begin{cases} 0\text{with more than 2 distinct values in \{k,l,m\}} ,\\ 8i\;\left(\frac{B}{N}\right)^{3/2} \phi_0 (1-\phi_0)(1-2\phi_0)\text{with 2 distinct values in \{k,l,m\}},\\ 8i\;\left(\frac{B}{N}\right)^{3/2}(N-1)\phi_0 (1-\phi_0)(1-2\phi_0)\text{all indices are the same;} \end{cases}\\ \frac{\partial^4 S}{\partial\phi_k\partial\phi_l\partial\phi_m\partial\phi_n}&= \begin{cases} 0\text{with more than 2 distinct values in \{k,l,m,n\},} \\ 16\; \frac{B^2}{N^2} \phi_0(1-\phi_0)(1-6\phi_0+6\phi_0^2)\text{with 2 distinct values in \{k,l,m,n\},} \\ 16\; \frac{B^2}{N^2} (N-1) \phi_0(1-\phi_0)(1-6\phi_0+6\phi_0^2)\text{all indices are the same.} \end{cases} \end{align}\] One can in principle build the standard Feynman perturbation theory from this expansion, with the inverse of the matrix \({\partial^2 S}/{\partial\phi_k\partial\phi_l}\) used as the propagator, and the higher-order terms in the \(\varphi\)-expansion used as vertices. However, the expressions for the diagrams will involve summations over indices taking values from \(1\) to \(N\), contributing positive powers of \(N\) that will compete with the negative powers seen in the denominators of the derivative formulas above. Furthermore, the powers of \(N\), both in the propagator and the vertices, depend on the number of coincident indices. This leads to complicated accounting of the powers of \(N\). By contrast, with the variables \(x_J\) and the representation (5 6 ) used in the main text, \(N\) is a purely numerical parameter in the action, and the powers of \(N\) are recovered straightforwardly, essentially, by means of Taylor expansions.

-3.5ex 2.3ex Some useful formulas

While the expressions below can be straightforwardly derived from (18 ), they are used extensively in our manipulations and we summarize them for the convenience of the reader: \[\begin{align} &\sigma_0 = \sum_{J=1}^{\infty} c_J^2 J = -\log(1-\phi_0),\\ &\sigma_1 = \sum_{J=1}^{\infty} c_J^2 J^2= \frac{1}{1+e^{2\alpha+4B\phi_0}} = \phi_0 ,\\ &\sigma_2 = \sum_{J=1}^{\infty} c_J^2 J^3 = \phi_0(1-\phi_0),\\ &\sigma_3 = \sum_{J=1}^{\infty} c_J^2 J^4 = \phi_0(1-\phi_0)(1-2\phi_0), \\ &\sigma_4 = \sum_{J=1}^{\infty} c_J^2 J^5 = \phi_0(1-\phi_0)(1-6\phi_0+6\phi_0^2). \end{align}\]

-3.5ex 2.3ex Particle-hole duality

The partition function of the two-star model is defined as \[Z\equiv\sum_A e^{-H[A]},\qquad H[A]\equiv \alpha\sum_{kl}{A_{kl}}+\beta\sum_{klm} A_{kl}A_{lm},\] We want to change the summation variable as \(A_{jk} = 1- \tilde{A}_{jk} -\delta_{jk}\), which precisely corresponds to flipping the state of all the edges (filled to empty, empty to filled). Then, \[H[A(\tilde{A})] = \left[-\alpha- 2\beta(N-1)\right] \sum_{kl}{\tilde{A}_{kl}}+\beta\sum_{klm} \tilde{A}_{kl}\tilde{A}_{lm} + \alpha N(N-1) + \beta(N^3 - 2N^2 + N).\] This means that under this transformation, the free energy \(F\equiv -\log Z\) changes as (substitute \(\beta = B/N\)) \[F(\alpha,B) = \alpha N(N-1) + B(N-1)^2 +F\left(-\alpha-2B +\frac{2B}{N},B\right).\] Since we express our results as functions of \(\phi_0(\alpha,B)\) defined by (15 ), rather than as functions of \(\alpha\) and \(B\), it is handy to correspondingly recast the transformation \(\alpha\rightarrow \tilde{\alpha}\equiv -\alpha-2B +\frac{2B}{N}\) as \(\phi_0 \rightarrow \tilde{\phi_0}\). To do so, we first express \(\alpha\) through \(\phi_0\) as \[\alpha= \frac{1}{2}\log\frac{1-\phi_0}{\phi_0}-2B\phi_0,\] enact the transformation \(\alpha\rightarrow \tilde{\alpha}\), and then write the equation defining \(\tilde{\phi}_0\) through \(\tilde{\alpha}\): \[\tilde{\phi}_0 = \frac{1}{1+e^{2\tilde{\alpha}+ 4B\tilde{\phi}_0}} = \frac{1}{1+ \frac{\phi_0}{1-\phi_0}e^{-4B(1-\phi_0)+4B/N}e^{4B\tilde{\phi_0}}}.\] This equation can be solved in terms of \(1/N\) expansion as \[\tilde{\phi}_0 \equiv 1-\phi_0 - \frac{1}{N} \frac{4B\phi_0(1-\phi_0)}{1+4B\phi_0(1-\phi_0)} - \frac{1}{N^2}\frac{8\phi_0B^2(1-\phi_0)(1-2\phi_0)}{[1+4B\phi_0(1-\phi_0)]^3}+\cdots. \label{phi095trans}\tag{39}\] With this transformation, we write \[F(\phi_0,B) = N(N-1)\left[\frac{1}{2}\log\frac{1-\phi_0}{\phi_0}-2B\phi_0\right] + B(N-1)^2 + F\left(\tilde{\phi}_0(\phi_0,B), B\right).\] This relation is respected by our result (25 ) up to order \(O(1)\), providing a validation of our derivations.

-3.5ex 2.3ex Monte Carlo sampling

We sample graphs from the two-star ensemble using a particular version of Metropolis Monte Carlo algorithm that we have previously employed in [33] for related models . We summarize it below. General discussion can be found in [42].

We start with an adjacency matrix of an Erdős-Rényi graph generated by randomly connecting any two vertices with probability \(1/2\). Then, for each Monte Carlo step, we propose an update move which flips the current edge (filled to empty, empty to filled) for a randomly picked pair of vertices \((i,j)\). The change in the adjacency matrix can be formally written as \[\label{eq:varA} \Delta_{ij}A_{mn} = (\delta_{im}\delta_{jn}+\delta_{in}\delta_{jm})(1-2A_{ij}).\tag{40}\] This move is then accepted with probability \(\text{min}\left(1,e^{H(A) - H(A+\Delta_{ij}A)}\right)\) as per the usual Metropolis algorithm.

Computing the full Hamiltonian takes a number of operations of order \(\mathcal{O}(N^2)\), which can be time-consuming when one has to take a large number of samples from a big system. One can optimize the numerical process by computing only the change of the Hamiltonian from each proposed update move. For this, one needs to compute the change in the total number of degrees \(\sum_{k}d_k = \sum_{kl}A_{kl}\) and two-stars \(\sum_{k}d_k^2 = \sum_{klm}A_{kl}A_{lm}\) of the graph. From 40 , the change in degree of node \(m\) is \[\Delta_{ij} d_m = (\delta_{im}+\delta_{jm})(1-2A_{ij}).\] It follows that \[\Delta_{ij}\sum_m d_m = 2(1-2A_{ij}),\qquad \Delta_{ij}\sum_{m}d^2_{m} = 2\left[(1-2A_{ij})(d_i + d_j) +1\right].\] The change in the Hamiltonian is thus \[\Delta_{ij}H(A)\equiv = H(A+\Delta_{ij}A)-H(A) = -2(1-2A_{ij})(\alpha+ \beta(d_i+d_j)+\beta).\] This allows one to compute the acceptance probability replying on a number of operations of order 1.

We compute the mean degree \(\left<k\right>\) and the mean degree squared \(\left<k^2\right>\) by averaging over \(5\times10^5\) Monte Carlo sample points, which are taken at intervals of \(N^2/10\) elementary Monte Carlo steps with the update rule described above. The long interval between the samples ensures that a significant part of the graph gets updated by the stochastic evolution before each next sample is taken.

References↩︎