Grand Canonical Approach to an Interacting Network

We consider a network composed of an arbitrary number of directed links. We employ a grand canonical partition function to study the statistical averages of the network in equilibrium. The Hamiltonian is composed of two parts: a “free” Hamiltonian H0 attributing a constant energy E to each link, and an interacting Hamiltonian Hint involving terms quadratic in the number of links. A Gaussian integration leads to a reformulated Hamiltonian, where now the number of links appears linearly. The reformulated Hamiltonian allows obtaining the exact behavior in limiting cases. At high temperatures the system reproduces the behavior of the free model, while at low temperatures the thermodynamic behavior is obtained by using a renormalized chemical potential, μeff = μ + λ, where λ is the strength of the interaction. We also resort to a mean field approximation, describing accurately the system over the entire range of all dynamical parameters. A detailed Monte-Carlo simulation verifies our theoretical expectations. We indicate that our model may serve as a prototype model to address a number of different systems.


Introduction
We are used to wondering first about particles or states and then about their interactions.We are first to ask about what it is and afterwards how it is.In most of our models in physics, or other branches of science, a system is considered as a collection of objects, with interactions among the objects introduced at a later stage.The examples are abundant (Ising model, lattice QCD, etc).It is not unjustified to consider that what lies beneath this general approach is a form of logic, Aristotelian logic.Aristotle's logic primes the issue of existence (or nonexistence) of a property attributed to an object.Boole gave an algebraic form to this type of logic, used as an instrument set theory (membership or non-membership of an object to a set).Most of our constructions in mathematics and physics rely precisely on set theory and Boolean logic.
Yet, there are signs that other forms of logic may be present.Quantum mechanics especially defy the classical logic.Neumann and Birkhoff noticed, as early as 1936, that the quantum really needed a non-Boolean logical syntax.We are led to reorient our thinking and consider that things have no meaning in themselves, and that only the correlations among them are real.Recently the relational logic of Peirce [1] [2] has been adopted as a consistent framework to prime correlations and it has been shown that this form of logic may serve as the conceptual foundation of quantum mechanics and string theory [3].Within the Peircean logic, relation is the fundamental, primary constituent, and everything else is expressed in terms of relations.Relation R lm is a generic term indicating a binary relation between two individual terms.Examples for possible R lm : 1) it may stand for the motion from initial point m to the final point l; 2) it may represent the transition from the m state to the l state; 3) it constitutes a proof, where starting from the initial theorem m we prove theorem l; 4) it stands for the internet connection from computer m to computer l.
Multiple relations among a number of sites lead to a network, characterized by its own topology or geometry.There are highly original propositions where space-time itself is considered as a network.For Wheeler, pregeometry, the stage preceding geometry, is the result of a calculus of propositions-relations [4] [5].Penrose suggested to view space-time as a network of relations-links carrying angular momenta [6] [7].A graph model, endowed with a Hamiltonian, may also lead to geometrogenesis [8] [9].Extensive studies indicate the usefulness of causal dynamical triangulation for the study of quantum space-time [10], while the connection of complex networks to de Sitter universe has been outlined [11].On the other hand, starting from different queries, networks have been employed for the study and analysis of quantum walks, internet, neural functioning of the brain, cognitive science, language acquisition, and biological systems [12]- [20].
A relation R lm may hold valid, or not valid, thus receiving two values, "yes" or "no".It resembles a fermion, which, upon measurement of its spin, is revealed as "spin up" or "spin down".We expect then that a network of relations will bear similarities to a network of fermions.In a network composed of relations-links, each individual link costs an energy E. The entire system has an average energy determined by the temperature T (inverse β).The number of valid occupied links is another parameter and is influenced by the chemical potential μ.Finally an interaction among relations, with strength determined by a parameter λ, gives rise to complexity and unexpected patterns for our relational network.We assume that the entire system, after multiple interactions, reaches a state of equilibrium and thus we are entitled to use methods of statistical mechanics for its study.In the next section we present our model and the available exact results.In the third section we present approximate solutions, involving Gaussian integration or the mean field approximation.In the fourth section a detailed Monte Carlo simulation is carried out.We may notice the excellent agreement between the Monte Carlo simulation and the mean field approximation.In the last section we present our conclusions and we mention directions for future work.

The Model
Essentially all network models considered in modern work have been ensemble models, meaning that a model is defined to be not a single network, but a probability distribution over many possible networks.The best choice of probability distribution is that from statistical physics that maximizes the Gibbs entropy.We adopt this approach here as well, where we discuss the equilibrium statistics of networks.The partition function Z can be defined as a sum over all graphs with a fixed number of sites and a fixed number of relations-links.Alternatively one can define the grand partition function Z G with an arbitrary number of links adjusted by a suitable chemical potential.

Description
After the introduction of the frame we are working in, let us describe our interacting network.We consider a loopless network (there is no link connecting a site to itself), with a fixed number of sites N, parameterized by a M × N matrix 1 if there is a link from to 0 otherwise The number of links is , , with the maximum number of links being ( ) . Suppose that rather than measuring the total number of links in a network, we measure the degrees of all the sites.For each site i we define the outgoing links out i k and the entering links Since each link appears twice, once as an "out" link and the other as an "in" link, we have In that way, the actual number of links is given by ( ) The Hamiltonian for our model is with H 0 representing the "free" Hamiltonian and H int providing an interaction among the different links.Within our model 0 , where E is the energy cost for an individual link.
For the interaction term we opted . 1 The interaction term resembles the two-star model [21] [22].The two-star model is an undirected network where the interaction terms couples all links attached to the same site.In our case we are studying a directed network, where we distinguish the incoming links from the outgoing links and the interaction couples the incoming flux to the outgoing flux for every single site.This type of interaction favors (for positive λ) the flow of a signal, the flow of information, the flow of a current, the spread of a virus, the propagation of a fire etc.

λ =
Consider first the "free" Hamiltonian (we set λ = 0).The grand partition function is We can evaluate then the average number of links where P 0 is the probability of occupation of a single link ( ) Notice the similarity of our "free" model to the Fermi-Dirac gas model.As T → ∞ ( ) β → , all links are equally occupied with probability 1 2 .As 0 all links are unoccupied, while for E µ < all links are occupied with unit probability.The average number of occupied links as a function of the temperature is shown in Figure 1, to be compared later with the corresponding quality of the complete model.

λ ≠
The full Hamiltonian is not amenable to exact analytic results, the main problem being the appearance of quadratic terms.To circumvent this problem we resort to a gaussian integration.Let us concentrate on the individual term exp 1 . By defining the column vector we obtain with A the 2 × 2 matrix Using well known Gaussian matrix integrals we may rewrite ( ) where now , in out k k appear as linear terms.Including all sites we obtain the full grand partition function where the integration measure is Summing over all possible configurations amounts to include the possible values of R ij (1 and 0).We end up with

Approximate Solutions
The expression above, Equation (18), allows to extract the behavior of Z G in limiting cases 1) When T → ∞ , or equivalently 0 β → , the integration measure, Equation ( 17), forces all , ja jb x x to re- main close to zero.We obtain then where ( ) We conclude that at high temperature, since , the interaction is irrelevant and our interacting model is reproduced by the results offered by the free model (with λ = 0).
2) In the case 0  dominates over unity and the integration for a single site is simplified to Including all sites we arrive at Thus at low temperatures the presence of the interaction amounts to a renormalization of the chemical potential The net result is that the interaction allows more links to be located under the Fermi energy with a high probability being occupied.

Mean-Field Approximation
For a network with a big number of links, we may resort to a mean field approximation, applicable to temperature.The grand canonical partition function may be rewritten as For a large number of links we may approximate in k and out k by their mean values ( ) ( ) where P is the occupation probability of a single link.
The above expression is similar to the grand canonical partition function of the free model, with μ replaced by 2 P µ λ

+
. We conclude that the occupation probability of a single link satisfies the equation ( ) We may attempt a numerical solution of the above equation.Let us define The solution(s) of Equation ( 28), we call it P * , is found as the intersection of the two curves f 1 and f 2 ( ) ( )

Monte Carlo Algorithm
We performed Monte Carlo simulations of the above described system and compared them with our exact results.Monte Carlo simulations are usually performed for thermal systems using the canonical ensemble and the well known "Metropolis" algorithm.The grand canonical ensemble is much more difficult to simulate and not routinely used in the literature [23].In the classic "Metropolis" method system configurations are generated from a previous state using a transition probability which depends on the energy difference between the initial and final states.The sequence of states produced follows a time ordered path, but the time in this case is "Monte Carlo time" and should not be confused with actual time.It is simply an effort to mimic a dynamic process that will lead the system to equilibrium so that statistical averages can be computed.Consider a discrete state-continuous time system.If ( ) n P t denotes the probability that a system is found at state n at time t and n m W → denotes the transition rate from state n to state m then the dynamics of the system are described by the following master equations (one for each of the discrete system states) In equilibrium ( ) Thus, one obtains the "detailed balance" condition: For the canonical ensemble the probability to find the system at the n-th state is: From Equation (37) and Equation (38) one deduces that the ratio of the transition probabilities depends only on the energy difference between the states and not on the unknown Z function.
( ) ( ) The calculation of thermodynamical statistical averages would be trivial if we knew the partition function Z in Equation (38).Since this is, however, generally unknown in Monte Carlo simulations we start from a random system configuration and allow it to evolve with suitable chosen transition probabilities.Any choice of transition probabilities that satisfies Equation (38) is acceptable.The standard "Metropolis" algorithm consists in choosing ( ) The "Metropolis" algorithm briefly procceeds as follows: 1) Choose a random initial state m.
2) Choose a new state n from the "neighboring" states to state m.
3) Calculate the energy difference ΔE between the two states.4) Draw a random number r uniformly distributed between 0 and 1. 5) Use Equation (39) to decide whether to accept the new state, i.e. accept if The above equations are valid for systems with a fixed number of "particles" N p at constant temperature T.An alternative way of studying physical systems is by using the grand canonical ensemble.In this case the number of "particles" of the system may vary i.e. the system can exchange particles with its environment.Thus, the independent variables are T and the chemical potential μ.For the grand-canonical ensemble the probability to find the system at the n-th state is: ( ) where Z G is the grand partition function, μ is the chemical potential and N p is the number of "particles" of the system.One may think that simulating the grand canonical ensemble is a trivial extension of the above algorithm and that it suffices to use Equation (41) in the place of Equation (38), while calculating the acceptance probabilities.This is true in principle; however there are practical complications that should be considered.Consider for simplicity a spin system, i.e. a system where each particle can be found in one of two states (spin up or down) and that can exchange particles with its environment.
A rough sketch of the proper algorithmic process that should be used for a grand canonical simulation is as follows: 1) Choose at random a state i.This state has energy E i and n particles.
2) Find all the "neighboring" states of i.This includes all the states accessible from i if we insert one particle, all states accessible from i if we remove one particle and all states resulting from a random spin flip (a rearrangement of the internal degrees of freedom).Select randomly one of these states.
3) Calculate the energy new E and the number of particles new n of this trail state.4) Calculate the acceptance ratio starting from Equation (41) instead of Equation (38).Accept or reject the trial state.
5) Continue until the system reaches a steady state.In practice, Step 2 of finding the "neighboring" states can become rather complicated and "heuristic" tricks have to be implemented to simplify the process.
Our network system is a system with variable number of links L. We performed grand canonical Monte Carlo simulations of this system using the temperature T and chemical potential μ as control parameters using the following algorithm.
1) We start from a network of N nodes and we randomly insert L 0 links out of the possible ( ) 2) We calculate the energy of the system H and the number of links L.
3) We draw a random number x between 0 and 1.If 1 3 x < we try to transpose a link.If 1 3 2 3 x < < we try to remove a link.And if 2 3 x > we try to insert a link in the system.4) In case we decide to transpose a link then we generate a new trial configuration with just one link in a different "place".We accept the new configuration with probability a correction term that accounts for the difference in the number of available states when the number of links decreases by one.7) We continue the simulation for a number of Monte Carlo steps which ensures that the system has reached equilibrium.
8) We repeat the calculations for several initial system realizations.

Results
In Figure 1 we plot the mean number of links L as a function of β for the case 0 λ = , after the system has reached the equilibrium state.The exact analytic results (solid lines) and the Monte Carlo simulations (points) are in excellent agreement.We observe that for small β ( ) T → ∞ , the mean number of links tends to max 1 2 L , while for large β ( ) 0 T → there are two limits: for E µ < , the mean number of links tends to max L and for E > the mean number of links tends to zero.We notice here that the results are identical to the corresponding results from traditional statistical mechanics for systems of non-interacting fermions (Fermi-Dirac distribution).
In and the bifurcation takes place at the predicted value 2 βλ = .In Figure 3 we display characteristic pictures of the network for both the free and the interacting case.We notice that in the free model there are very few connecting links.In the interaction model on the other hand, we obtain a very dense network.
What is really striking in our findings is the impressive agreement between the mean field approximation, Equation (28), and the Monte Carlo simulation we carried out (see all the figures).We attribute this agreement to the large number N of the sites involved.Our findings imply a sort of "equivalence" among all sites and therefore we do not expect the creation of hubs in the network.

Conclusions
We have studied a dynamical directional network model with a variable number of links.In most of similar studies  the links are added or removed following a given prescription.In our model the number of links is affected by the presence of an interaction term in the Hamiltonian.The whole network pattern is sensitive upon the values of the strength of interaction λ, the inverse temperature β and the chemical potential μ.Though the model is not simple, we manage to obtain a semi-analytic description.At high temperatures the dynamics of the network is described by the free model, while at low temperatures the determining parameter is the renormalized chemical potential eff µ µ λ = + .It is rather gratifying the nice agreement between the theoretical expectations and the Monte-Carlo simulation.
We consider our model as a prototype model, able to address different situations.Variation of the assumptions and the parameters involved may lead to distinct patterns.For example: 1) We assume the same energy for all links.We may examine the possibility that links connected to specific sites have a reduced energy.Then we anticipate the emergence of hubs.Such a structure may resemble the internet structure.
2) There is no geometry in the present model.We may envisage that the sites belong to a fixed geometry, for example a lattice.The interaction term in the Hamiltonian may stand for the propagation of a germ from a human organ to another, or the propagation of a fire in a forest.
3) We consider λ as a constant.A space-dependent λ (including negative values) will give rise to a cluster of multiple networks, or a "multiverse".
Work along these lines is in progress.

Figure 1 .
Figure 1.Mean number of links L vs. β for λ = 0 and for values of the energy E below and above the chemical potential μ.The points are Monte Carlo simulation results for systems of N = 50 nodes, the solid line represents the analytic results and the error bars represent the standard deviation of L.

the energy of the new configuration. 5 )
In the case of insertion we generate a new configuration with one additionally randomly placed link and accept the new configuration with probability that accounts for the difference in the number of available states when the number of links increases by one.6)Links are removed with probability ( )

Figure 2
we present the results for the interacting case, with 2 λ = , 2 µ = , for different energies, as a function of β.Distinctive features emerge.First at high temperatures ( ) in agreement with our findings in section 2.3, that at high temperatures the interaction term is not important and the results are reproduced by the grand partition function of the free model.As β increases (temperature decreases) the results are sensitive upon the precise value of the energy.In Figure 2(a) the energy is smaller than the effective chemical potential.The probability starts to rise with increasing β (as expected) and at a critical β value the bifurcation takes place.The original solution becomes unstable and two new branch-solutions emerge.In Figure 2(b) on the other hand the energy is bigger than the effective chemical potential.As β increases, the probabiity decreases and again at a critical β value the bifurcation occurs.In Figure 2(c) the energy is equal to the effective chemical potential.The probability remains constant at the value 1 2

Figure 2 .
Figure 2.For λ = 2 and μ = 2.The solid line represents mean-field approximation results and the points are Monte Carlo simulation results for systems of N = 50 nodes.
When the slope of f 2 is bigger than the slope of f 1 at P * , then f 2 rises above f 1 for