Statistical mechanics of multiedge networks

I. INTRODUCTION

The increasing and unprecedented quality and quantity of available data coming from very different areas is boosting the field of complex networks. Interdisciplinary science demands new efforts and new tools, addressed not only to develop more efficient computational strategies to analyze incoming data but also to span a theoretical framework where both more accurate and more tractable models can provide predictions closer to the reality one wants to face. In this context, a standard approach consists of representing in a graph the complex structure of interactions among the elements of a given system. Statistical mechanics is an extraordinary framework where such a complex structure can be appropriately modeled and with this aim a large amount of studies has appeared in recent years .

The simplest representation of a network assumes the existence of nodes and edges. The edges do not need to be necessarily symmetric and they can provide information about the relative influence (interaction) of a node onto another, i.e., they can be directed and have a certain intensity. However, the early studies in the field were essentially focused on a binary projection of the network on a graph where only the existence of an edge and its distribution were required to determine some of its properties, without considering the effect of weights. In this way one could compute the probability for two arbitrary sites (nodes) to be connected through an edge. It is well known that such probability keeps certain analogies with occupation numbers in quantum statistics, in particular, with fermionic systems . Further developments have extended these results to richer and more complex structures such as directed and weighted graphs finding analogy with bosonic systems .

A successful and complete description of any system amenable to be represented as a complex network through statistical mechanics requires an appropriate characterization regarding the detailed features of the microscopic components forming the network. Up to now, the actions represented by an edge have always been considered as indistinguishable events, but in some cases, the data represented by edges and their weights can be generated by independent, identifiable actions. This circumstance needs to be taken into account whenever using any statistical technique applied to complex networks. Its impact being important as can be seen in analogy to statistical mechanics where the distinguishability of particles leads to different descriptions in terms of Fermi-Dirac, Bose-Einstein or Maxwell-Boltzmann statistics.

Processes generated by single agents represent single events, for instance, a trip between two locations, a call between two cities or communications at a given time. A particularly clear example where this distinguishability of weight units is crucial are transportation networks , where human flows between different locations are used to build the so-called Origin Destination matrices which collect global information about their mobility or allocation. A naive approach based on a standard weighted description of a network is not satisfactory for this problem as it was already pointed out by Wilson who mapped transport systems to statistical physics using an entropy maximization approach. Hence, modern complex networks theory calls for the introduction of new general models taking into account the possible distinguishability of the edges forming the considered network.

The present work addresses this issue by presenting a general statistical mechanics approach to networks created from distinguishable single-unit events that can be grouped in edges, in analogy to the allocation of particles in discrete energy levels of classical statistical mechanics systems. In such networks, the edge weights encode information about independent processes performed by different agents at (possibly) different times, and hence those are necessarily represented by integer numbers and constitute quantized, distinguishable entities. The difference between this representation, which leads to what we call multiedge networks, and the already considered weighted ones is clear as schematically shown in Fig. . Multiedge networks assume the existence of a minimal weight of unit value (and hence a quantization) representing an indissociable event. In some cases, groups of these events connect the same pair of nodes, hence forming an edge with multiple connections, which has a different nature from a weighted one where no quantization is imposed. If one aims to apply the well-known ensemble theory used in statistical physics to networks , one needs to take into account this subtle yet important difference, which has profound implications on the associated statistics. The choice of one or the other representation will thus depend on the problem at hand and makes a big difference in terms of collective behavior of the whole network as will be shown along this paper. This study fills the existing gap in the field by introducing a complete and necessary discussion on distinguishability that is to be added to the already existing models and thus complements the complex networks body of theory.

FIG. 1

The distinguishability issue: Out of a set of three distinguishable nodes and three distinguishable (directed) events, one can see that the number of configurations (described in terms of edge names) compatible with a given occupation number sequence is degenerate, since all permutations of events generate *different* microstates.

Many examples of systems studied using network theory can be understood as multiedge networks. All those networks composed by the actions of independent agents (mobility networks , call networks ) and in general those whose edges are formed by the sum of independent actions developed during a certain time fall within this category . The present work deals with the challenge of introducing null models that provide exact analytical expectations for the main network observables under some given fixed constraints. This, in turn, allows one to assess whether any observed network characteristic from real data is due to the imposed constraints or to the contrary, to some unexpected phenomena worth exploring.

The paper is structured as follows: Sec. introduces the concept of multiedge networks together with some basic nomenclature and definitions used. Sections introduce the general derivation of the considered ensembles using different sets of constraints (mainly linear constraints on occupation numbers, binary-projection constraints, and both). In each section, some important examples are explicitly derived and formulas for the state statistics and the ensemble entropies are obtained. Finally some concluding remarks are given discussing the work done. The Appendices contain details, additional discussion, and some mathematical developments on mentioned examples in the text.

II. ENSEMBLE APPROACH TO MULTIEDGE NETWORKS: NOMENCLATURE AND DEFINITIONS

A multiedge network is a collection of nodes that may be connected by none, one or more than one element (event) taken from a set of independent, distinguishable entities. Such an object admits a coarse-grained representation in terms of a multiedge matrix T. The matrix entries $𝑡 𝑖 𝑗$ are bounded integer stochastic variables $𝑡 𝑖 𝑗 \leq 𝑡 (max) 𝑖 𝑗 = 𝑇$ denoting the number of events joining nodes $𝑖$ and $𝑗$ and $𝑇 = \sum 𝑖 𝑗 𝑡 𝑖 𝑗$ represents the total number of events. In the context of statistical mechanics they play an analog role to occupation numbers. The matrix $𝐓$ need not necessarily be symmetric so $𝑡 𝑖 𝑗 \neq 𝑡 𝑗 𝑖$ in general. Let $𝑁$ be the number of nodes of the network which generate $𝑁 (𝑁 - 1)$ possible states where individual events can be allocated. $𝐸 = \sum 𝑖 \neq 𝑗 𝛩 (𝑡 𝑖 𝑗)$ represents the number of occupied states (regardless of their occupation number as long as $𝑡 𝑖 𝑗 > 0$ ) in the network, representing $𝛩 (𝑥)$ the Heaviside step function,

𝛩 (𝑥) = {0 𝑥 = 0 1 𝑥 \geq 1 (𝑥 \in 𝒩) .

Such a function ensures that $𝐸$ accounts for the number of existing connections regardless of the number of events contained in each entry of the matrix, therefore creating a binary projection of the network.

Our goal is to construct an ensemble framework that allows one to treat multiedge networks with any given set of constraints $⃗ 𝐶 = ⃗ 𝐶 ({𝑡 𝑖 𝑗})$ defined in terms of the variables of the system. The constraints define a macrostate and restrict the available phase space to all the possible graphs compatible with such constraints. In this context, we assume two main starting hypotheses. First, all the configurations (microstates) compatible with the observed constraints have the same a priori probability of appearance. Secondly, we assume that $𝑇, 𝑁 ≫ 1$ being $𝑁$ given (this means that the topological structure of the network, the number of available states, does not change) which allows a statistical treatment of the problem. We further assume that the distribution of occupation numbers is stationary and defines in turn a probability that fixes our thermodynamic limit. This probability $𝑝 (𝑡 𝑖 𝑗)$ indicates the asymptotic (relative) distribution of occupation numbers:

𝑝 * 𝑖 𝑗 = ⟨ 𝑡 𝑖 𝑗 \sum 𝑡 𝑖 𝑗 ⟩ \to 𝑝 (𝑡 𝑖 𝑗) as 𝑇 \to \infty,

(1)

where $⟨ \dots ⟩$ denotes an average performed over the ensemble considered. Under these considerations it is possible to establish a complete mapping between the problem at hand and a classical statistical mechanics problem. The events correspond to distinguishable particles occupying any of the $𝑁 (𝑁 - 1)$ available energy states [the methodology is laid for directed graphs, but can be easily adapted to undirected ones, yielding $𝑁 (𝑁 - 1) / 2$ available states]. The main difference between the proposed system and the ones studied by classical statistical equilibrium mechanics regards the constraints used: They will not be in general of extensive, global nature but local at the level of nodes. The standard procedure, already used by different authors, consists of finding an expression for the probability of obtaining a macrostate defined by its multiedge adjacency matrix $𝐓$ ,

𝒫 (𝐓) = 𝒫 ({𝑡 𝑖 𝑗}),

and then maximizing its associated entropy,

𝑆 S H = - \sum 𝛤 𝒫 (𝐓) l n 𝒫 (𝐓),

(2)

where the sum runs over all possible configurations of events (microstates) on the accessible phase space $𝛤$ of the given ensemble. Since distinguishability plays an important role in our work, it is important to note that if one wishes to compute the sum over possible values of the occupation numbers ${𝑡 𝑖 𝑗}$ (our canonical variables) rather than over configurations, a degeneracy term needs to be added to expression (see Fig. ),

𝒟 ({𝑡 𝑖 𝑗}) = (\sum 𝑖 𝑗 𝑡 𝑖 𝑗)! \prod 𝑖 𝑗 𝑡 𝑖 𝑗! .

(3)

In the remainder of the paper, we shall work on the occupation-number space $𝛺$ (which is a coarsened representation of the $𝛤$ configurational space), including the degeneracy term and hence considering an expression for the entropy of the form,

𝑆 S H = - \sum 𝛤 𝒫 ({𝑡 𝑖 𝑗}) l n 𝒫 ({𝑡 𝑖 𝑗}) = - \sum 𝛺 ({𝑡 𝑖 𝑗} 𝒟 ({𝑡 𝑖 𝑗}) 𝒫 ({𝑡 𝑖 𝑗}) l n (𝒟 ({𝑡 𝑖 𝑗}) 𝒫 ({𝑡 𝑖 𝑗})) = - \sum 𝛺 ({𝑡 𝑖 𝑗}) 𝑃 ({𝑡 𝑖 𝑗}) l n 𝑃 ({𝑡 𝑖 𝑗}),

(4)

where we changed the notation $𝒫 \to 𝑃$ to denote the counting over occupation numbers in the (coarse-grained) $𝛺$ space rather than event configurations in $𝛤$ space.

A final comment deserves the seminal work by Wilson on transport theory . He followed a similar scheme and found expressions for the expected number of events connecting two arbitrary nodes according to some considered constraints. However, in the present paper we go further; we present a modern methodology able to handle more complicated constraints as compared to those analyzed by Wilson and also consider constraints which affect simultaneously the distribution of occupation numbers and its binary projection on the graph, for instance, those affecting the strength (weighted degree ) of a node and its degree.

III. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS ON THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$

In this work we proceed following a microcanonical scenario. In this ensemble, all the configurations are equally probable and the constraints are considered “hard,” i.e., the phase space accessible always fulfils strictly the constraints. Thus, we must only maximize the expression $l n 𝛺 ({𝑡 𝑖 𝑗})$ with respect to $𝑡 𝑖 𝑗$ which allows one to find the most probable value of some relevant observables, i.e., their statistically expected value in the ensemble of the maximally random graph which strictly fulfills the constraints:

max {𝛺' ({𝑡 𝑖 𝑗}) | ⃗ 𝐶 ({𝑡 𝑖 𝑗}) = ⃗ 𝐶} \equiv max {𝛺 ({𝑡 𝑖 𝑗})} .

(5)

We follow the same methodology developed in seminal works by Bianconi and write down the volume of the $𝛺$ space introducing auxiliary fields $ℎ 𝑖 𝑗$ ,

𝛺 = \sum {𝑐} 𝑒 \sum {𝑖, 𝑗} ℎ 𝑖 𝑗 𝑡 𝑖 𝑗 𝑄 \prod 𝑞 𝛿 (𝐶 𝑞 - 𝐶 𝑞 ({𝑡 𝑖 𝑗})) = \sum {𝑐} \int 𝑒 \sum {𝑖, 𝑗} ℎ 𝑖 𝑗 𝑡 𝑖 𝑗 𝑄 \prod 𝑞 𝐷 [𝜃 𝑞] 𝑒 𝜃 𝑞 (𝐶 𝑞 ({𝑡 𝑖 𝑗} - 𝐶 𝑞)),

(6)

where ${𝑐}$ denotes a sum over microscopic configurations ( $𝛤$ space), $𝑄$ is the total number of constraints, and $⃗ 𝜃 𝑞$ are the related Lagrange multipliers. Note also that an integral representation of the Kronecker delta has been used. The introduction of auxiliary fields $ℎ 𝑖 𝑗$ allows one to recover all the central moments of the distribution of $𝑡 𝑖 𝑗$ . In fact, one can see that is closely related to the cumulant generating function of $𝑡 𝑖 𝑗$ and hence all its cumulants can be recovered by differentiation. In particular,

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝜕 ℎ 𝑖 𝑗 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 0 \forall 𝑖, 𝑗, 𝜎 2 𝑡 𝑖 𝑗 = 𝜕 2 ℎ 𝑖 𝑗 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 0 \forall 𝑖, 𝑗, 𝜎 2 𝑡 𝑖 𝑗, 𝑡 𝑘 𝑙 = 𝜕 2 ℎ 𝑖 𝑗 ℎ 𝑘 𝑙 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 0 \forall 𝑖, 𝑗 .

(7)

If we wish to perform the sum over occupation numbers ( $𝛺$ space) rather than over configurations of the system, we need to take into account the degeneration given in ,

𝛺 = \sum {𝑡' 𝑖 𝑗} \int 𝑇! \prod 𝑖 𝑗 𝑡' 𝑖 𝑗! 𝑒 \sum 𝑖 𝑗 ℎ 𝑖 𝑗 𝑡' 𝑖 𝑗 \prod 𝑞 𝐷 [𝜃 𝑞] 𝑒 - 𝜃 𝑞 𝐶 𝑞 𝑒 𝜃 𝑞 𝐶 𝑞 ({𝑡' 𝑖 𝑗}) .

(8)

For this ensemble and for linear constraints on occupation numbers, one can write

𝐶 𝑞 ({𝑡 𝑖 𝑗}) = \sum 𝑖 𝑗 𝑐 (𝑖 𝑗) 𝑞 𝑡 𝑖 𝑗,

(9)

being $𝑐 (𝑖 𝑗) 𝑞$ a quantity that usually depends on a “property” of the edge between nodes $𝑖$ and $𝑗$ (a distance, for instance) or be a real number ( $𝑐 (𝑖 𝑗) 𝑞 = 𝛿 𝑖, 𝑞$ for outgoing strength sequence constraints, for example, as we shall see). In such a situation, the sum over ${𝑡 𝑖 𝑗}$ sequences such that $\sum 𝑡 𝑖 𝑗 = 𝑇$ can be exactly performed inside the integral in yielding

𝛺 = \int (\prod 𝑞 𝐷 [𝜃 𝑞]) 𝑒 - \sum 𝑞 𝜃 𝑞 𝐶 𝑞 𝑒 𝑇 l n \sum 𝑖 𝑗 e x p {(ℎ 𝑖 𝑗 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞)} = \int (\prod 𝑞 𝐷 [𝜃 𝑞]) 𝑒 𝑓 ({𝜃 𝑞, 𝐶 𝑞}, {ℎ 𝑖 𝑗}) .

(10)

The occupation number statistics can be shown to have multinomial nature (see Appendix ),

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑇 𝑝 𝑖 𝑗; 𝜎 2 𝑡 𝑖 𝑗 = 𝑇 𝑝 𝑖 𝑗 (1 - 𝑝 𝑖 𝑗), 𝜎 2 𝑡 𝑖 𝑗, 𝑡 𝑘 𝑙 = - 𝑇 𝑝 𝑖 𝑗 𝑝 𝑘 𝑙 .

(11)

This fact assures that the relative fluctuations of occupation numbers vanish in the thermodynamic limit ( $𝑇 \to \infty$ ) since

𝜎 𝑡 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ = \sqrt 1 - 𝑝 𝑖 𝑗 𝑝 𝑖 𝑗 𝑇 \to 0, as 𝑇 \to \infty .

We identified $𝑝 𝑖 𝑗$ as the probability for an individual event to be assigned to state $𝑖 𝑗$ . Explicitly,

𝑝 𝑖 𝑗 = 𝑒 \sum 𝑞 𝑐 𝑖 𝑗 𝑞 𝜃 𝑞 \sum 𝑖 𝑗 𝑝 𝑖 𝑗 .

(12)

Concerning the entropy of the graphs in this ensemble, the integral in can be approximated to first order by using the steepest descent methods,

𝑆 B G = l n 𝛺 | ℎ 𝑖 𝑗 \forall 𝑖, 𝑗 = 0 ≃ - 𝑓 * = - \sum 𝑞 𝐶 𝑞 𝜃 * 𝑞 + 𝑇 l n 𝑒 \sum 𝑞 𝜃 * 𝑞 \sum 𝑖 𝑗 𝑐 (𝑖 𝑗) 𝑞,

(13)

where $𝜃 * 𝑞$ are the solutions of the saddle point equations given by

𝜕 𝜃 𝑞 𝑓 | ℎ 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝐶 𝑞 = 𝐶 𝑞 (⟨ 𝑡 𝑖 𝑗 ⟩) = 𝑇 \sum 𝑖, 𝑗 𝑐 (𝑖 𝑗) 𝑞 𝑝 𝑖 𝑗 .

(14)

Merging one obtains the event-specific entropy of a given graph in this ensemble:

𝑆 B G 𝑇 = - (\sum 𝑞 𝜃 𝑞 \sum 𝑖, 𝑗 𝑐 (𝑖 𝑗) 𝑞 𝑝 𝑖 𝑗 - l n \sum 𝑖 𝑗 𝑒 \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞) = - (\sum 𝑖 𝑗 𝑝 𝑖 𝑗 \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞 - l n \sum 𝑖 𝑗 𝑒 \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞) = - \sum 𝑖 𝑗 𝑝 𝑖 𝑗 l n 𝑝 𝑖 𝑗 = 𝑆 S H 𝑇,

(15)

which has a final Shannon entropy form, equivalent to the Boltzmann-Gibbs entropy.

If one is able to exactly solve the saddle point equations obtained from the steepest descent approximation, the full distribution of occupation numbers is recovered. Despite being in a microcanonical framework, one could consider a canonical ensemble where the constraints are fulfilled on average, hence $𝐶 𝑞 (⟨ 𝑡 𝑖 𝑗 ⟩) = ⟨ 𝐶 𝑞 (𝑡 𝑖 𝑗) ⟩$ . Having proven that the partition function has a multinomial-cumulant form over the occupation numbers and considering only linear soft constraints, the requirement that those constraints need to be fulfilled only on average in the sampling over the phase space is automatically satisfied. Moreover, the constraints have vanishing relative fluctuations in the thermodynamic limit,

⟨ 𝐶 𝑞 ({𝑡 𝑖 𝑗}) ⟩ = \sum 𝑖, 𝑗 𝑐 (𝑖 𝑗) 𝑞 ⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑇 \sum 𝑖, 𝑗 𝑐 (𝑖 𝑗) 𝑞 𝑝 𝑖 𝑗, 𝜎 2 𝐶 𝑞 = \sum 𝑖, 𝑗, 𝑘, 𝑙 𝑐 (𝑖 𝑗) 𝑞 𝑐 (𝑘, 𝑙) 𝑞 𝜎 2 𝑖 𝑗, 𝑘 𝑙 \propto 𝑇, 𝜎 2 𝐶 𝑞 ⟨ 𝐶 𝑞 ⟩ 2 \to 1 𝑇 \to 0 as 𝑇 \to \infty,

where we have used the properties of the multinomial distribution presented in . Although the theoretical basis for the generation of graphs in different ensembles is introduced in this paper (see Appendix ), the challenges for the exact and efficient generation of such ensembles will be shortly tackled and presented in future work.

In the following, we shall consider some explicit cases of linear constraints on $𝑡 𝑖 𝑗$ .

A. No constraints

This is the simplest case where we have a single hard constraint (apart from the number of nodes $𝑁$ ) $𝑇 = \sum 𝑖 𝑗 𝑡 𝑖 𝑗$ . Therefore reads

𝛺 = \sum \sum 𝑡 𝑖 𝑗 = 𝑇 𝑇! \prod 𝑖 𝑗 𝑡 𝑖 𝑗! \prod 𝑖 𝑗 (𝑒 ℎ 𝑖 𝑗) 𝑡 𝑖 𝑗 = (\sum 𝑖 𝑗 𝑒 ℎ 𝑖 𝑗) 𝑇 .

In this case it is straightforward to determine the average value of the occupation numbers (in this case $‾ 𝑡$ denotes an average over a single graph realization):

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑇 (𝑁 (𝑁 - 1)) \equiv ‾ 𝑡 \equiv 𝑇 𝑝 \forall 𝑖, 𝑗 .

(16)

Therefore, all the occupation numbers have constant probability $𝑝 𝑖 𝑗 = 𝑝 \forall 𝑖, 𝑗$ of being chosen per event sorted. The result is according to what intuition would tell us: Under no constraints events are equally distributed among levels which reminds the high-temperature regime in classical systems where there is an arbitrary large (but finite) number of energy levels. It is also possible to compute the covariances on occupation numbers,

𝜎 2 𝑡 𝑖 𝑗, 𝑡 𝑘 𝑙 = {- 𝑇 𝑝 2 𝑖 𝑗 \neq 𝑘 𝑙 𝑇 𝑝 (1 - 𝑝) 𝑖 𝑗 = 𝑘 𝑙 .

Other typical network magnitudes of interest such as the strengths (both incoming $𝑠 (i n) 𝑗 = \sum 𝑖 𝑡 𝑖 𝑗$ and outgoing $𝑠 (o u t) 𝑖 = \sum 𝑗 𝑡 𝑖 𝑗$ ), which will be extensively used in this paper, can also be computed. They are random integer variables with fixed mean $⟨ 𝑠 𝑖 ⟩ \equiv ‾ 𝑠 = (𝑁 - 1) 𝑝 𝑇 \equiv 𝑝 𝑠 𝑇$ and variance $𝜎 2 𝑠 = 𝑇 𝑝 𝑠 (1 - 𝑝 𝑠)$ $\forall 𝑖 \in [1, 𝑁]$ which are on average equal for each node and have also a multinomial character.

Concerning the Boltzmann-Gibbs entropy of the ensemble,

𝑆 B G = 𝑇 l n (𝑁 (𝑁 - 1)) = - 𝑇 l n 𝑝 = 𝑆 S H .

(17)

Which recovers a Shannon form over events as expected, $𝑆 S H = - \sum 𝑗 𝑖 𝑝 l n 𝑝$ .

B. Fixed average event cost $‾ 𝑐 = 𝐶 / 𝑇$

If additionally to the number of events $𝑇$ , we consider a cost matrix (symmetric, dense, and positive definite) $𝐃 = {𝑑 𝑖 𝑗}$ and fix the total cost $𝐶 𝑇 = \sum 𝑑 𝑖 𝑗 𝑡 𝑖 𝑗$ , we trivially get

𝛺 = \int 𝐷 [𝜃] e x p {- 𝜃 𝐶 𝑇 + 𝑇 l n \sum 𝑒 ℎ 𝑖 𝑗 + 𝜃 𝑑 𝑖 𝑗} = \int 𝐷 [𝜃] e x p {𝑐 (𝜃, {ℎ 𝑖 𝑗}, {𝑡 𝑖 𝑗})},

which leads to the saddle point equation,

𝜕 𝜃 𝑐 | ℎ 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝐶 𝑇 𝑇 \equiv ‾ 𝑐 = \sum 𝑑 𝑖 𝑗 𝑒 𝜃 𝑑 𝑖 𝑗 \sum 𝑒 𝜃 𝑑 𝑖 𝑗,

and finally,

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑇 𝑒 𝜃 𝑑 𝑖 𝑗 \sum 𝑒 𝜃 𝑑 𝑖 𝑗 = 𝑇 𝑝 (𝜃, 𝑑 𝑖 𝑗), 𝜎 2 𝑡 𝑖 𝑗 = 𝑇 𝑝 (𝜃 𝑑 𝑖 𝑗) (1 - 𝑝 (𝜃, 𝑑 𝑖 𝑗)) .

And this leads to a weighted version of the Waxman graph . This reasoning can be extended to study the interesting case where the distribution of costs is also fixed, which is of particular interest in the field of O-D matrices used to analyze mobility, where the mobility of users using certain types of transports is assumed to follow particular statistical forms (). This case is analyzed in detail and solved in Appendix .

For any of the cases involving cost matrices, especially those related with distances, it is very important to remark that the allocation of occupation numbers is not independent in each state and hence it is not true that the probability $𝑊 (𝑡 𝑖 𝑗, 𝑑 𝑖 𝑗, 𝜃)$ of having a state occupied by $𝑡 𝑖 𝑗$ events at distance $𝑑 𝑖 𝑗$ is $𝑊 (𝑡 𝑖 𝑗, 𝑑 𝑖 𝑗) = 𝐾 𝑓 (𝜃, 𝑑 𝑖 𝑗)$ . Rather $𝑊 (𝑡 𝑖 𝑗)$ represents a conditional probability of observing $𝑡 𝑖 𝑗$ events (trips in this scenario) at distance $𝑑 𝑖 𝑗$ given the distance matrix $𝐃$ and the rest of the constraints (included in $𝐾$ ). This means in particular that if a deference function $𝑓 (𝑑 𝑖 𝑗)$ is proposed to explain observed flows between locations, as usually done in O-D matrix studies under a maximum entropy assumption, its predictive results need to be statistically tested against the full expected distribution of $𝑡 𝑖 𝑗$ , and not only against the (biased) statistic of observed or existing occupation numbers.

C. Fixed relative strength sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, 𝑇$

We consider now the case in which the only given constraint is the strength sequence $⃗ 𝑠$ . In this case, the constraints read

𝐶 𝑖 (𝑡 𝑖 𝑗) = \sum 𝑖 𝑡 𝑖 𝑗 = 𝑠 o u t 𝑖 𝐶 𝑗 (𝑡 𝑖 𝑗) = \sum 𝑗 𝑡 𝑖 𝑗 = 𝑠 i n 𝑗,

(18)

and Eq. becomes

𝛺 = \int 𝑒 - \sum 𝑖 𝛼 𝑖 𝑠 o u t 𝑖 𝑒 - \sum 𝑗 𝛽 𝑗 𝑠 i n 𝑗 (𝑁 \prod 𝑖 𝐷 [𝛼 𝑖] 𝐷 [𝛽 𝑖]) \times \sum \sum 𝑡 𝑖 𝑗 = 𝑇 𝑇! \prod 𝑖 𝑗 𝑡 𝑖 𝑗! \prod 𝑖 𝑗 (𝑒 𝛼 𝑖 + 𝛽 𝑗 + ℎ 𝑖 𝑗) 𝑡 𝑖 𝑗 = \int (𝑁 \prod 𝑖 𝐷 [𝛼 𝑖] 𝐷 [𝛽 𝑖] 𝑒 - 𝛼 𝑖 𝑠 o u t 𝑖 𝑒 - 𝛽 𝑖 𝑠 i n 𝑖) (\sum 𝑖 𝑗 𝑒 𝛼 𝑖 + 𝛽 𝑗 + ℎ 𝑖 𝑗) 𝑇 = \int 𝐷 [{𝛼 𝑖, 𝛽 𝑗}] e x p - \sum 𝑖 𝛼 𝑖 𝑠 o u t 𝑖 - \sum 𝑗 𝛽 𝑗 𝑠 i n 𝑗 + 𝑇 l n (\sum 𝑖 \neq 𝑗 𝑒 𝛼 𝑖 + 𝛽 𝑗 + ℎ 𝑖 𝑗) = \int 𝐷 [{𝛼 𝑖, 𝛽 𝑗}] e x p 𝑓 ({𝛼 𝑖}, {𝛽 𝑖}, {ℎ 𝑖 𝑗}) .

Now we need to solve the $2 𝑁$ saddle point equations,

𝜕 𝛼 𝑖 𝑓 | ℎ 𝑖 𝑗 = 0 \forall 𝑖, 𝑗 = 0 \Rightarrow 𝑠 o u t 𝑖 = 𝑇 \sum 𝑗 \neq 𝑖 𝑒 𝛼 𝑖 + 𝛽 𝑗 \sum 𝑒 𝛼 𝑖 + 𝛽 𝑗, 𝜕 𝛽 𝑗 𝑓 | ℎ 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑠 i n 𝑗 = 𝑇 \sum 𝑖 \neq 𝑗 𝑒 𝛼 𝑖 + 𝛽 𝑗 \sum 𝑒 𝛼 𝑖 + 𝛽 𝑗 .

(19)

And we apply again to get

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑇 𝑥 𝑖 𝑦 𝑗 \sum 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 = 𝑇 𝑝 𝑖 𝑗,

(20)

where we have identified $𝑥 𝑖 = 𝑒 𝛼 𝑖, 𝑦 𝑗 = 𝑒 𝛽 𝑗$ and $𝑝 𝑖 𝑗$ has the same multinomial structure and properties as in previous examples except for the fact that these probabilities are state dependent. In fact, the average number of multiedges between two nodes factorizes in an uncorrelated form.

For large $𝑁$ we are led to

𝑠 o u t 𝑖 = 𝑇 𝑥 𝑖 \sum 𝑗 𝑦 𝑗 \sum 𝑖 \neq 𝑗 𝑥 𝑖 \sum 𝑗 𝑦 𝑗 = 𝑇 𝑥 𝑖 𝑋 - 𝑥 𝑖 \approx 𝑇 𝑥 𝑖 𝑋, 𝑋 \equiv \sum 𝑖 𝑥 𝑖; 𝑌 \equiv \sum 𝑗 𝑦 𝑗; 𝑠 o u t 𝑖 \propto 𝑥 𝑖; 𝑠 i n 𝑗 \propto 𝑦 𝑗,

(21)

and we recover the weighted configuration model ,

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑠 o u t 𝑖 𝑠 i n 𝑗 𝑇 .

(22)

This result recovers the expression in and is in accordance with a maximum likelihood principle (contrary to the case of weighted networks as explained in ).

Let us notice that these results allow a straightforward extension to the canonical ensemble. By identifying from ,

𝑝 𝑠 o u t 𝑖 = 𝑥 𝑖 𝑋 𝑝 𝑠 i n 𝑗 = 𝑦 𝑗 𝑌,

(23)

one can work with a strength sequence which is no longer fixed but a collection of fluctuating integer random values with a multinomial structure, since the partial grouping of multinomial random variables $𝑡 𝑖 𝑗$ preserve their multinomial character and by construction $\sum 𝑗 𝑝 𝑠 i n 𝑗 = \sum 𝑖 𝑝 𝑠 o u t 𝑖 = 1$ . Therefore,

⟨ 𝑠 o u t 𝑖 ⟩ = 𝑇 𝑝 𝑠 o u t 𝑖, ⟨ 𝑠 i n 𝑗 ⟩ = 𝑇 𝑝 𝑠 i n 𝑗, 𝜎 𝑠 o u t 𝑖 = 𝑇 𝑝 𝑠 o u t 𝑖 (1 - 𝑝 𝑠 o u t 𝑖), 𝜎 𝑠 i n 𝑗 = 𝑇 𝑝 𝑠 i n 𝑗 (1 - 𝑝 𝑠 i n 𝑗) .

(24)

Concerning the entropy in this canonical ensemble scenario, making use of we further obtain a closed expression in terms of the constraints of the problem,

𝑆 B G ≃ - \sum 𝑖 𝑠 o u t 𝑖 l n 𝑠 o u t 𝑖 𝑇 - \sum 𝑗 𝑠 i n 𝑗 l n 𝑠 i n 𝑗 𝑇 + 𝑇 l n \sum 𝑖 \neq 𝑗 𝑠 o u t 𝑖 𝑠 i n 𝑗 𝑇 2 = - 𝑇 (\sum 𝑖 𝑝 𝑠 o u t 𝑖 l n 𝑝 𝑠 o u t 𝑖 + \sum 𝑗 𝑝 𝑠 i n 𝑗 l n 𝑝 𝑠 i n 𝑗), 𝑆 B G 𝑇 = - \sum 𝑖 𝑝 𝑠 o u t 𝑖 l n 𝑝 𝑠 o u t 𝑖 - \sum 𝑗 𝑝 𝑠 i n 𝑗 l n 𝑝 𝑠 i n 𝑗 = 𝑆 S H 𝑇 .

(25)

Let us remind that in this context $𝑝 𝑠 𝑖 = \sum 𝑗 𝑝 𝑖 𝑗$ represents the probability of a certain node $𝑖$ to accumulate $⟨ 𝑠 𝑖 ⟩$ (incoming or outgoing) events on average. We therefore recover a Shannon form for the entropy in both microcanonical and canonical ensembles, which scales with the total number of events $𝑇$ and is node (but also state) specific [from Eq. ].

Further additional constraints can be added leading to different models. An example would be to merge the last two examples of a network living in a metric space (or any network where we can give a cost to the edges) with fixed “accessibility” of each node (distance-weighted strength). We could also fix the average cost per trip $‾ 𝑐$ (global constraint) and the strength of each node (node local constraint). In this case one would obtain the popular doubly constrained gravity model in several forms , also Stouffers’ intervening opportunities model and even the newly proposed radiation model (by choosing an appropriate form for the cost function; see for extended discussion).

The case of the gravity law models deserves a closer look since an entropy maximization approach yields a form $⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑥 𝑖 𝑦 𝑗 𝑓 (𝛾, 𝑑 𝑖 𝑗)$ which is not equivalent (in general) to $⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑠 𝛼 𝑠 𝛽 𝑓 (𝛾, 𝑑 𝑖 𝑗)$ because the values of the multipliers depend on the particular spatial distribution of the considered nodes and their relative strengths. Hence, despite the success attained by these kinds of models to reproduce empirical data an entropy maximization approach could unify the different sets of exponents observed in each study (see ).

IV. MULTIEDGE NETWORK WITH GIVEN NONLINEAR CONSTRAINTS ON THE BINARY PROJECTION OF THE OCCUPATION NUMBERS $𝛩 (𝑡 𝑖 𝑗)$

Let us consider now more complex situations such as the case of nonlinear constraints on $𝑡 𝑖 𝑗$ . One of the most relevant objects to look at when dealing with complex networks is the degree distribution. It concerns only the existence of links between arbitrary nodes of the graph regardless the number of multiedges between them. In the framework presented in this paper it can be worked out as a function of the binary projection of the occupation numbers on the graph so, in general, such constraints can be expressed as

ˆ 𝐶 𝑞' = \sum 𝑖 𝑗 ̂ 𝑐 (𝑖 𝑗) 𝑞' 𝛩 (𝑡 𝑖 𝑗) .

(26)

The main technical difficulty in dealing with these types of constraints is that they do not allow the summation in a multinomial form of the terms in ${𝑡 𝑖 𝑗}$ inside the integral of the partition function for exactly fixed $𝑇$ . A workaround to perform the summation can be found, however: Instead of making the multinomial sum at once, we proceed in two steps. We first introduce a Kronecker delta in integral form inside the integral of the partition function with associated Lagrange multiplier $𝜃$ which fixes the total number of events. We secondly allow the sum inside the integral to cover all the available values of the phase space ( ${𝑡 𝑖 𝑗 \in [0, 𝑇] \forall 𝑖, 𝑗}$ . Finally, since the grouping of terms in the sum not fulfilling the constraints (including the constraint on the total expected number of events) will be penalized by the Kronecker deltas introduced earlier, we relax the limit on the sum over individual occupation number configurations from $𝑇 \to \infty$ what reminds the standard approach to the grandcanonical ensemble (in Appendix the calculation with finite $𝑇$ is also given).

Proceeding as explained we obtain a new version of Eq. :

𝛺 = \int 𝐷 [𝜃] 𝑒 - 𝜃 𝑇 𝑒 - \sum 𝑞' 𝜆 𝑞' ˆ 𝐶 𝑞' (\prod 𝑞' 𝐷 [𝜆 𝑞']) 𝑇! \prod {𝑖, 𝑗}' 𝑒 (𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞') 𝛩 (𝑡 𝑖 𝑗) 𝑒 (ℎ 𝑖 𝑗 + 𝜃) 𝑡 𝑖 𝑗 𝑡 𝑖 𝑗! = \int 𝐷 [𝜃] 𝑒 - 𝜃 𝑇 𝑒 - \sum 𝑞' 𝜆 𝑞' ˆ 𝐶 𝑞' (\prod 𝑞' 𝐷 [𝜆 𝑞']) 𝑇! [𝑒 (𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞') 𝑇 \sum 𝑡' 𝑖 𝑗 = 1 𝑒 (ℎ 𝑖 𝑗 + 𝜃) 𝑡' 𝑖 𝑗 𝑡' 𝑖 𝑗! + 1 0!] \sum \sum 𝑡 𝑖 𝑗 = 𝑇 - 𝑡' 𝑖 𝑗 \prod (𝑒 𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞') 𝛩 (𝑡 𝑖 𝑗) (𝑒 ℎ 𝑖 𝑗 + 𝜃) 𝑡 𝑖 𝑗 𝑡 𝑖 𝑗! = \int 𝐷 [𝜃] 𝑒 - 𝜃 𝑇 𝑒 - \sum 𝑞' 𝜆 𝑞' ˆ 𝐶 𝑞' (\prod 𝑞' 𝐷 [𝜆 𝑞']) 𝑇! e x p {\sum 𝑖 𝑗 l n {𝑒 𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 ℎ 𝑖 𝑗 + 𝜃 - 1) + 1}},

(27)

where we have introduced auxiliary fields $𝑏 𝑖 𝑗$ for the binary projections of the occupation numbers. Both the complete set of cumulants of the occupation numbers and their binary projections can now be recovered by differentiation. For the occupation numbers using and in the case of the binary projection,

⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ = 𝜕 𝑏 𝑖 𝑗 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖, 𝑗, 𝜎 2 𝛩 (𝑡 𝑖 𝑗) = 𝜕 𝑏 2 𝑖 𝑗 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖, 𝑗, 𝜎 2 𝛩 (𝑡 𝑖 𝑗), 𝛩 (𝑡 𝑘 𝑙) = 𝜕 𝑏 𝑖 𝑗 𝑏 𝑘 𝑙 l n 𝛺 ({𝑡 𝑖 𝑗}) | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖, 𝑗 .

(28)

From expression one can compute explicitly those values yielding

⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ \equiv ˆ 𝑝 𝑖 𝑗 = 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 𝜃 - 1) 𝛥 𝑖 𝑗, ⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑒 𝑒 𝜃 𝑒 𝜃 𝑒 𝑒 𝜃 - 1 ˆ 𝑝 𝑖 𝑗 \equiv ‾ 𝑡 + ˆ 𝑝 𝑖 𝑗, 𝜎 2 𝛩 (𝑡 𝑖 𝑗) = ˆ 𝑝 𝑖 𝑗 (1 - ˆ 𝑝 𝑖 𝑗); 𝜎 2 𝑡 𝑖 𝑗 = ⟨ 𝑡 𝑖 𝑗 ⟩ (1 + 𝑒 𝜃 - ⟨ 𝑡 𝑖 𝑗 ⟩), 𝜎 2 𝑡 𝑖 𝑗, 𝑡 𝑘 𝑙 = 𝜎 2 𝛩 (𝑡 𝑖 𝑗), 𝛩 (𝑡 𝑘 𝑙) = 0 if 𝑖 𝑗 \neq 𝑘 𝑙, 𝛥 𝑖 𝑗 \equiv 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 𝜃 - 1) + 1 .

(29)

For ease in notation, we perform the change $𝜌 \equiv 𝑒 𝜃$ and we identify $⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩$ with the probability $ˆ 𝑝 𝑖 𝑗$ of connection of nodes $𝑖$ and $𝑗$ and $‾ 𝑡 +$ with the graph-average occupation of existing links. In Appendix we prove that the resulting partition function is equal to the cumulant generating function of the outcome of $𝑁 (𝑁 - 1)$ independent zero-inflated Poisson processes (ZIP ) with individual associated probability,

𝑃 (𝑡 | ˆ 𝑝, 𝜌) = (1 - ˆ 𝑝) 1 - 𝛩 (𝑡) (ˆ 𝑝 𝑒 𝜌 - 1 𝜌 𝑡 𝑡!) 𝛩 (𝑡),

(30)

which in turn, regarding only the binary projection, corresponds to the outcome of independent Bernoulli processes with probabilities $ˆ 𝑝$ ,

𝑃 (𝛩 (𝑡) | 𝑝) = ˆ 𝑝 𝛩 (𝑡) (1 - ˆ 𝑝) 1 - 𝛩 (𝑡) .

(31)

From one sees that in this case the multiedge structure is completely determined by the binary constrained topology. Our coarse-grained description in terms of independent occupation numbers implies that for each state, two outcomes can be considered: Either the edge does not exist (obviously with 0 occupation) or it does exist, in which case the resulting (conditioned) statistics being Poisson with mean value $⟨ 𝑡 | 𝑡 \geq 1 ⟩ = 𝜌 𝑒 𝜌 𝑒 𝜌 - 1$ .

The constant relation of proportionality $⟨ 𝑡 𝑖 𝑗 ⟩ \propto ˆ 𝑝 𝑖 𝑗$ rapidly allows one to identify the graph-average occupation of the expected existing links $⟨ 𝐸 ⟩ = \sum 𝑖 𝑗 ˆ 𝑝 𝑖 𝑗$ ,

‾ 𝑡 + = 𝑇 ⟨ 𝐸 ⟩ = 𝑒 𝜌 𝜌 𝑒 𝜌 - 1 > 1; 𝜌 > 0,

(32)

which can be inverted leading to

𝜌 = 𝑊 (- 𝑒 - ‾ 𝑡 + ‾ 𝑡 +) + ‾ 𝑡 +,

(33)

where $𝑊 (𝑥)$ is the Lambert $𝑊$ function . Figure shows a plot of Eq. stressing the rapid asymptotical convergence $𝜌 \to ‾ 𝑡 +$ as $‾ 𝑡 + \to \infty$ (in fact, the approximation is clearly good as soon as $‾ 𝑡 + ≃ 5$ ).

FIG. 2

The result of Eq. is shown together with the equality line $𝜌 = ‾ 𝑡 +$ . One can see the rapid convergence: For $‾ 𝑡 + = 2.31$ we obtain $𝜌 ‾ 𝑡 + = 0.9$ and for $‾ 𝑡 + = 4.615$ one finds $𝜌 ‾ 𝑡 + = 0.99$ .

With and we can compute the relative fluctuations of both occupation numbers and binary links in the thermodynamic limit,

l i m 𝑇 \to \infty 𝜎 2 𝛩 (𝑡 𝑖 𝑗) ⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ 2 = l i m 𝑇 \to \infty (1 - ˆ 𝑝 𝑖 𝑗) ˆ 𝑝 𝑖 𝑗 = (1 - ˆ 𝑝 𝑖 𝑗) ˆ 𝑝 𝑖 𝑗, l i m 𝑇 \to \infty 𝜎 2 𝑡 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ 2 = l i m 𝑇 \to \infty (1 + 𝜌 - ‾ 𝑡 + ˆ 𝑝 𝑖 𝑗) ‾ 𝑡 + ˆ 𝑝 𝑖 𝑗 = 1 + ‾ 𝑡 + (1 - ˆ 𝑝 𝑖 𝑗) ‾ 𝑡 + ˆ 𝑝 𝑖 𝑗, 𝜎 2 𝑡 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ 2 \to 𝜎 2 𝛩 (𝑡 𝑖 𝑗) ⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ 2 as ˆ 𝑝 𝑖 𝑗 fixed, ‾ 𝑡 + = 𝑇 / \sum 𝑖 𝑗 ˆ 𝑝 𝑖 𝑗 \to \infty .

(34)

These expressions reflect the bimodal structure of the state statistics and explains the nonvanishing nature of the relative fluctuations: The variance of the occupation numbers has a maximum for $ˆ 𝑝 𝑖 𝑗 | m a x = 12 (1 + ‾ 𝑡 + - 1) \to 12$ , vanishes for the absence ( $ˆ 𝑝 \to 0$ ) of an edge and converges to Poisson statistics for edges that always exist ( $ˆ 𝑝 \to 1$ ). The existence of an edge is a binary event, hence the maximum variability corresponds to the draw situation ( $50 %$ chance). In such a case, approximately half of the times a graph is created the considered edge will have (on average) occupation $‾ 𝑡 +$ and the other half occupation 0, generating vast fluctuations on the overall statistics which are caused by the constrained binary structure of the graph.

Concerning the entropy, performing the steepest descent approximation on to first order as in we have

𝑆 B G ≃ - 𝑇 l n 𝜌 - \sum 𝑖 𝑗 ˆ 𝑝 𝑖 𝑗 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' + l n (𝑇!) + \sum 𝑖 𝑗 l n 𝛥 𝑖 𝑗 .

Using $𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝜌 - 1) = ˆ 𝑝 𝑖 𝑗 1 - ˆ 𝑝 𝑖 𝑗$ , Eq. , one is lead to

𝑆 B G = - {\sum 𝑖 𝑗 (ˆ 𝑝 𝑖 𝑗 l n ˆ 𝑝 𝑖 𝑗 + (1 - ˆ 𝑝 𝑖 𝑗) l n (1 - ˆ 𝑝 𝑖 𝑗))} - {𝑇 l n 𝜌 - l n 𝑇! - l n (𝑒 𝜌 - 1) ⟨ 𝐸 ⟩}, 𝑆 B G = 𝑆 b i n + 𝑆 d i s t .

(35)

This expression has two clear contributions: The first term $𝑆 b i n$ is the entropy corresponding to the binary constrained topology while the second one $𝑆 d i s t$ corresponds to the additional multiedge distinguishable structure. In other words, $𝑆 b i n$ counts all the possible ways to select $𝐸$ states out of a total number $𝑁 (𝑁 - 1)$ while $𝑆 d i s t$ refers to the possible ways to allocate the $𝑇$ events on these $𝐸$ surviving states.

The second term in can be explicitly evaluated in two limiting cases: The dense case which corresponds to the thermodynamic limit and the sparse for which the binary and weighted structure are equal $‾ 𝑡 + = 𝑇 / 𝐸 \to 1$ (hence $⟨ 𝐸 ⟩ = 𝐸$ is fixed).

Considering the sparse case, one has from $‾ 𝑡 + \to 1$ so $𝜌 \to 0$ and hence,

l i m 𝜌 \to 0, 𝑇 \to ⟨ 𝐸 ⟩ 𝑆 d i s t = l n 𝐸!,

(36)

which corresponds to the microcanonical counting of configurations coming from valid permutations of distinguishable multiedges over the fixed binary structure of $𝐸$ surviving occupied states.

The dense case corresponds to $𝑇 \to \infty$ , which implies from $𝜌 \to ‾ 𝑡 + = 𝑇 / ⟨ 𝐸 ⟩$ and $l n (𝑇!) ≃ 𝑇 l n 𝑇 - 𝑇$ ,

l i m 𝑇 \to \infty 𝑆 d i s t 𝑇 = l n ⟨ 𝐸 ⟩,

(37)

which has a Shannon form if we consider $𝑝 \equiv ⟨ 𝐸 ⟩ - 1$ which would be the probability associated with a multinomial process of sorting $𝑇$ events over the surviving $⟨ 𝐸 ⟩$ binary links with identical probability $𝑝$ . In this limit, the difference between sorting $⟨ 𝐸 ⟩$ independent Poisson processes with mean $𝑇 / ⟨ 𝐸 ⟩$ and sorting $⟨ 𝐸 ⟩$ Poisson processes excluding the zero-occupation events is negligible.

In the following, we present some examples to clarify the usage of this new methodology.

A. Fixed $𝐸$ , $𝑇$

We start by considering the most simple case where we only fix the total number of events $𝑇 = \sum 𝑖 𝑗 𝑡 𝑖 𝑗$ and the total number of existing binary links $𝐸 = \sum 𝑖 𝑗 𝛩 (𝑡 𝑖 𝑗)$ on the network, in analogy to the paradigmatic model of the Erdos-Renyi graph in binary networks .

In this case we introduce two Lagrange multipliers ( $𝜆, 𝜃$ ) to fix ( $𝐸, 𝑇$ ). Proceeding from Eq. , we readily obtain the saddle point equations,

𝐸 = \sum 𝑖 𝑗 ⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ = 𝑁 (𝑁 - 1) 𝜒 (𝑒 𝜌 - 1) 𝜒 (𝑒 𝜌 - 1) + 1, 𝑇 = \sum 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑁 (𝑁 - 1) 𝜒 𝑒 𝜌 𝜌 𝜒 (𝑒 𝜌 - 1) + 1,

(38)

where we identify $𝜒 = 𝑒 𝜆, 𝜌 = 𝑒 𝜃$ . We see that the average occupation numbers and edge existence probability are constant and their average values proportional as expected. Using we compute the relevant magnitudes,

⟨ 𝛩 (𝑡 𝑖 𝑗) ⟩ \equiv ˆ 𝑝 = 𝐸 𝑁 (𝑁 - 1) = C n t ⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑒 𝜌 𝜌 𝑒 𝜌 - 1 ˆ 𝑝 = ‾ 𝑡 + ˆ 𝑝 = 𝑇 𝐸 ˆ 𝑝 = 𝑇 𝑁 (𝑁 - 1) = C n t 𝜎 2 𝑡 = 𝑇 𝐸 ˆ 𝑝 (1 + 𝜌 - 𝑇 𝐸 ˆ 𝑝) 𝜎 2 𝛩 = ˆ 𝑝 (1 - ˆ 𝑝) .

(39)

We recover the binary structure of the well-known Erdös-Renyi graph as a result of the binary projection of a nontrivial multiedge structure, which on average values fulfills $⟨ 𝑡 ⟩ = 𝑇 𝐸 𝑝 = ‾ 𝑡 + 𝑝$ .

Despite the average occupation numbers over the ensemble being equal to the cases in Sec. , the underlying statistic is not. All the nodes (and states) in this case are statistically equivalent and their associated strengths and degrees (incoming and outgoing) are proportional on average (since they are fluctuating quantities not being fixed by the constraints):

⟨ 𝑠 ⟩ = ‾ 𝑠 = 𝑇 𝐸 ⟨ 𝑘 ⟩ = 𝑇 𝐸 ‾ 𝑘 = 𝑇 𝐸 𝐸 𝑁 = 𝑇 𝑁;   \sqrt 𝜎 2 𝑘 ⟨ 𝑘 ⟩ 2 =   \sqrt ‾ 𝑘 - 1 - \sum 𝑗 ˆ 𝑝 2 𝑖 𝑗 ‾ 𝑘 2; \sqrt 𝜎 2 𝑠 ⟨ 𝑠 ⟩ 2 =   \sqrt 1 ‾ 𝑡 + 𝑘 + \sum 𝑗 ˆ 𝑝 (1 - ˆ 𝑝) ‾ 𝑘 2 \to   \sqrt 𝜎 2 𝑘 ⟨ 𝑘 ⟩ 2 as 𝑇 \to \infty .

(40)

The entropy is readily computed from yielding,

𝑆 B G = 𝑆 𝐸 - 𝑅 + 𝑆 d i s t, 𝑆 𝐸 - 𝑅 = - 𝑁 (𝑁 - 1) (ˆ 𝑝 l n ˆ 𝑝 + (1 - ˆ 𝑝) l n (1 - ˆ 𝑝)), 𝑆 d i s t = l n (𝑇!) - 𝑇 l n (‾ 𝑡 + + 𝑊 (- ‾ 𝑡 + 𝑒 - ‾ 𝑡 +)) + 𝐸 l n (𝑒 ‾ 𝑡 + + 𝑊 (- ‾ 𝑡 + 𝑒 - ‾ 𝑡 +) - 1) .

(41)

Expression can also be computed using combinatorial arguments: Consider a process in which one selects $𝐸$ states out of $𝑁 (𝑁 - 1)$ and then populates each state with a single event chosen out of a set of $𝑇$ distinguishable entities, finally, the rest of the $𝑇 - 𝐸$ events are sorted in the $𝐸$ surviving states chosen in the first place. The counting of microstates reads

𝛤 (𝐸, 𝑇, 𝑁) = (𝑁 (𝑁 - 1) 𝐸) 𝑇! (𝑇 - 𝐸)! 𝐸 𝑇 - 𝐸,

(42)

and hence the microcanonical entropy is $𝑆 B G = l n 𝛤 (𝐸, 𝑇, 𝑁)$ . Here again one recovers the equivalence with in both the sparse and dense limits considered earlier.

B. Fixed degree sequence $⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}, 𝑇$

The next important case to consider is the one where the node binary connectivity of the graph is fixed. Such a situation is specially interesting as our framework permits one to understand binary networks as the projection (for instance, due to partial information or limited resolution) of a process generated by independent agents.

In this case the constraints read

ˆ 𝐶 𝑖 (𝑡 𝑖 𝑗) = \sum 𝑖 𝛩 (𝑡 𝑖 𝑗) = 𝑘 o u t 𝑖, ˆ 𝐶 𝑗 (𝑡 𝑖 𝑗) = \sum 𝑗 𝛩 (𝑡 𝑖 𝑗) = 𝑘 i n 𝑗 .

(43)

So we introduce two sets of Lagrange multipliers ( ${𝜆 𝑖, 𝛽 𝑖}$ ) and for expression we have

𝛺 = \int 𝐷 [{𝛼 𝑖, 𝛽 𝑗}] e x p {- \sum 𝑖 𝛼 𝑖 𝑘 o u t 𝑖 - \sum 𝑗 𝛽 𝑗 𝑘 i n 𝑗} \times e x p {- 𝜃 𝑇 + l n 𝑇! + \sum 𝑖 𝑗 l n [𝑒 𝛼 𝑖 + 𝛽 𝑗 + 𝑏 𝑖 𝑗 (𝑒 𝑒 ℎ 𝑖 𝑗 + 𝜃 - 1) + 1]} = \int 𝐷 [{𝛼 𝑖, 𝛽 𝑗}] e x p 𝑔 ({𝛼 𝑖, 𝛽 𝑗, ℎ 𝑖 𝑗, 𝑏 𝑖 𝑗}, 𝜃) .

And we solve the saddle point equations,

𝜕 𝜃 𝑔 | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑇 = \sum 𝑖, 𝑗 𝑥 𝑖 𝑦 𝑗 𝑒 𝜌 𝜌 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) + 1, 𝜕 𝛼 𝑖 𝑔 | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑘 o u t 𝑖 = \sum 𝑗 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) + 1, 𝜕 𝛽 𝑗 𝑔 | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑘 i n 𝑗 = \sum 𝑖 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) + 1,

(44)

where we have identified $𝑥 𝑖 \equiv 𝑒 𝛼 𝑖$ , $𝑦 𝑗 \equiv 𝑒 𝛽 𝑗$ , and $𝜌 \equiv 𝑒 𝜃$ . The occupation numbers and binary occupation probability, respectively, read

⟨ 𝑡 𝑖 𝑗 ⟩ = 𝑥 𝑖 𝑦 𝑗 𝑒 𝜌 𝜌 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) + 1, ˆ 𝑝 𝑖 𝑗 = 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) 𝑥 𝑖 𝑦 𝑗 (𝑒 𝜌 - 1) + 1 .

(45)

As expected we recover the proportionality of strengths and degrees for each node under this particular set of constraints.

⟨ 𝑡 𝑖 𝑗 ⟩ = ˆ 𝑝 𝑖 𝑗 𝑒 𝜌 𝜌 𝑒 𝜌 - 1 \Rightarrow ⟨ 𝑠 𝑖 ⟩ = ‾ 𝑡 𝑘 𝑖 = 𝑇 𝐸 𝑘 𝑖 .

(46)

Considering only the binary projection of the graph, one gets

ˆ 𝑝 𝑖 𝑗 = 𝜇 𝜅 𝑖 𝜆 𝑗 𝜇 𝜅 𝑖 𝜆 𝑗 + 1,

(47)

where we have identified $𝜅 𝑖 \equiv 𝑥 𝑖$ , $𝜆 𝑗 \equiv 𝑦 𝑗$ , and $𝜇 \equiv (𝑒 𝜌 - 1)$ .

The previous expression corresponds exactly with the expression for the so-called canonical ensemble of the random graph with any given degree distribution , where for fixed $𝑁$ , the $𝜇$ parameter controls the edge density . And since $𝑥 𝑖$ is a quantity related with node $𝑖$ (hence related with $𝑘 𝑖$ ), we obtain again the structural correlations of the configuration graph model .

V. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS DEPENDING ON BOTH THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$ AND THEIR BINARY PROJECTION $𝛩 (𝑡 𝑖 𝑗)$

We analyze for the sake of completeness the most general case where both types of considered constraints are fixed, and .

We introduce two sets of Lagrangian multipliers for the multiedge constraints ${𝜃 𝑞}$ and the binary ones ${𝜆 𝑞'}$ , plus an additional $𝜃$ corresponding to the constraint on the total number of events. The procedure then is analogous to the one in the previous section yielding finally from ,

𝛺 = \int 𝐷 [𝜃] 𝑒 - 𝜃 𝑇 𝑒 - \sum 𝑞' 𝜆 𝑞' ˆ 𝐶 𝑞' 𝑒 - \sum 𝑞 𝜃 𝑞 𝐶 𝑞 (\prod 𝑞 𝐷 [𝜃 𝑞]) \times (\prod 𝑞' 𝐷 [𝜆 𝑞']) 𝑇! e x p {\sum 𝑖 𝑗 l n \times {𝑒 𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 ℎ 𝑖 𝑗 + 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞 - 1) + 1}} .

(48)

Using and we obtain the statistic for the occupation numbers and their projections,

⟨ 𝑡 𝑖 𝑗 ⟩ = 1 𝛥 𝑖 𝑗 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' e x p (𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞) 𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞, 𝜎 𝑡 𝑖 𝑗 = ⟨ 𝑡 𝑖 𝑗 ⟩ (1 + 𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞 - ⟨ 𝑡 𝑖 𝑗 ⟩), ˆ 𝑝 𝑖 𝑗 = 1 𝛥 𝑖 𝑗 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞), 𝜎 𝛩 (𝑡 𝑖 𝑗) = ˆ 𝑝 𝑖 𝑗 (1 - ˆ 𝑝 𝑖 𝑗), 𝛥 𝑖 𝑗 = 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞 - 1) + 1 .

(49)

Assuming that the saddle point equations can be solved, which means that the imposed combination of binary and multiedge constraints is graphical, i.e., $𝑠 𝑖 \geq 𝑘 𝑖 \forall 𝑖 \in [1, 𝑁]$ for the case of fixed strength and degree sequence, for instance, then the probability to obtain a graph can still be written in terms of the Lagrange multipliers as a sum of independent ZIP processes with different parameters,

𝑃 (𝐓 | {ˆ 𝑝 𝑖 𝑗, 𝜇 𝑖 𝑗}) = \prod 𝑖, 𝑗 (1 - ˆ 𝑝 𝑖 𝑗) 1 - 𝛩 (𝑡 𝑖 𝑗) {ˆ 𝑝 𝑖 𝑗 𝑒 𝜇 𝑖 𝑗 - 1 𝜇 𝑡 𝑖 𝑗 𝑖 𝑗 𝑡 𝑖 𝑗!} 𝛩 (𝑡 𝑖 𝑗),

(50)

where $𝜇 𝑖 𝑗 = 𝑒 𝜃 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞$ is a quantity related to the average value of the occupation number of the given link, conditioned that this link exists,

\sum \infty 𝑡 = 1 𝑡 𝑖 𝑗 𝑃 (𝑡 𝑖 𝑗 | 𝜇 𝑖 𝑗, ˆ 𝑝 𝑖 𝑗, 𝑡 𝑖 𝑗 > 0) \sum \infty 𝑡 = 1 𝑃 (𝑡 𝑖 𝑗 | 𝜇 𝑖 𝑗, ˆ 𝑝 𝑖 𝑗, 𝑡 𝑖 𝑗 > 0) = 𝜇 𝑖 𝑗 𝑒 𝜇 𝑖 𝑗 𝑒 𝜇 𝑖 𝑗 - 1 = ⟨ 𝑡 + 𝑖 𝑗 ⟩ .

(51)

The relative fluctuations on occupation numbers do not vanish in the thermodynamic limit due to the strong constraints imposed by the binary structure. One can still express $𝜇 𝑖 𝑗 (⟨ 𝑡 + 𝑖 𝑗 ⟩)$ using ,

𝜇 𝑖 𝑗 (⟨ 𝑡 + 𝑖 𝑗 ⟩) = ⟨ 𝑡 + 𝑖 𝑗 ⟩ + 𝑊 (- ⟨ 𝑡 + 𝑖 𝑗 ⟩ 𝑒 - ⟨ 𝑡 + 𝑖 𝑗 ⟩) .

(52)

In the thermodynamic limit ( $𝑇 \to \infty$ which implies $⟨ 𝑡 + 𝑖 𝑗 ⟩ \to \infty$ ), expression converges to $𝜇 𝑖 𝑗 ≃ ⟨ 𝑡 + 𝑖 𝑗 ⟩$ .

Regarding the relation between expected occupation numbers and their binary projections, one has $⟨ 𝑡 𝑖 𝑗 ⟩ = ⟨ 𝑡 + 𝑖 𝑗 ⟩ ˆ 𝑝 𝑖 𝑗$ and the constant relation of proportionality is broken $⟨ 𝑡 𝑖 𝑗 ⟩ ˆ 𝑝 𝑖 𝑗 \neq C n t$ . This extends the well-known result that it is impossible to generate uncorrelated networks both at the level of strengths and degrees for multiedge networks or weighted networks .

Concerning the entropy, approximating by saddle point methods using we obtain the general expression that includes all the previous cases considered,

𝑆 B G = 𝑆 d i s t + 𝑆 b i n, 𝑆 d i s t = l n 𝑇! + \sum 𝑖 𝑗 ˆ 𝑝 𝑖 𝑗 l n (𝑒 𝜇 𝑖 𝑗 - 1) - \sum 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ l n 𝜇 𝑖 𝑗,

(53)

where $𝑆 b i n$ is still the binary contribution to the entropy [with the same form as in ].

The two limiting cases early considered can again be evaluated. The sparse case implies that $‾ 𝑡 + 𝑖 𝑗 = 𝑇 𝐸 \to 1 \forall 𝑖, 𝑗$ which means $𝜇 𝑖 𝑗 \to 0$ and $⟨ 𝑡 𝑖 𝑗 ⟩ \to ˆ 𝑝 𝑖 𝑗$ obtaining

l i m 𝑇 \to 𝐸 𝑆 d i s t = l n 𝐸!,

(54)

which is identical to the previous one (since it is equivalent to dropping the strength constraints).

For the thermodynamic limit (dense case), we have $𝑒 𝜇 𝑖 𝑗 - 1 \to 𝑒 𝜇 𝑖 𝑗 \to 𝑒 ⟨ 𝑡 + 𝑖 𝑗 ⟩$ and then,

l i m 𝑇 \to \infty 𝑆 d i s t 𝑇 = - \sum 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ 𝑇 l n ⟨ 𝑡 + 𝑖 𝑗 ⟩ 𝑇 = - \sum 𝑖 𝑗 ⟨ 𝑡 𝑖 𝑗 ⟩ 𝑇 l n ⟨ 𝑡 𝑖 𝑗 ⟩ 𝑇 ˆ 𝑝 𝑖 𝑗,

(55)

for which the previous expressions encountered are limiting cases. On one hand, if we relax the constraints on the occupation numbers, then $⟨ 𝑡 + 𝑖 𝑗 ⟩ = ‾ 𝑡 + = 𝑇 / 𝐸 \forall 𝑖, 𝑗$ and we recover expression . On the other hand, not fixing any binary related quantity implies that as $𝑇 \to \infty$ , $ˆ 𝑝 𝑖 𝑗 = 1 - 𝑒 - ⟨ 𝑡 𝑖 𝑗 ⟩ \to 1$ (fully connected topology) and we are lead to . Finally, not fixing any constraints, $⟨ 𝑡 + 𝑖 𝑗 ⟩ = ⟨ 𝑡 𝑖 𝑗 ⟩ = ‾ 𝑡$ , we recover .

For simplicity, the only example we report in this section corresponds to the very relevant case in which both strength and degree sequences are fixed, the rest of the cases being easily derivable from the general theory exposed.

A. Fixed relative strength sequence and fixed degree sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, ⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}$

We analyze here a situation where the constraints imposed are the strength and degree sequence [ and ]. The calculations are analogous to the previous section obtaining

𝜕 𝜃 ℎ | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑇 = \sum 𝑖 𝑗 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 𝑒 𝜃 𝑧 𝑖 𝑤 𝑗, 𝜕 𝛼 𝑖 ℎ | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑘 o u t 𝑖 = \sum 𝑗 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 (𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 - 1), 𝜕 𝛽 𝑗 ℎ | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑘 i n 𝑗 = \sum 𝑖 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 (𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 - 1), 𝜕 𝛾 𝑖 ℎ | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑠 o u t 𝑖 = \sum 𝑗 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 𝑧 𝑖 𝑤 𝑗 𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 𝑒 𝜃, 𝜕 𝜀 𝑗 ℎ | ℎ 𝑖 𝑗 = 𝑏 𝑖 𝑗 = 0 \forall 𝑖 𝑗 = 0 \Rightarrow 𝑠 i n 𝑗 = \sum 𝑖 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 𝑧 𝑖 𝑤 𝑗 𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 𝑒 𝜃, 𝛥 𝑖 𝑗 \equiv 𝑥 𝑖 𝑦 𝑗 {e x p (𝑒 𝜃 𝑧 𝑖 𝑤 𝑗) - 1} + 1,

where we have identified $𝑥 𝑖 \equiv 𝑒 𝜆 (𝑘 o u t) 𝑖, 𝑦 𝑗 \equiv 𝑒 𝜆 (𝑘 i n) 𝑗, 𝑧 𝑖 \equiv 𝑒 𝜃 (𝑠 o u t) 𝑖, 𝑤 𝑗 \equiv 𝑒 𝜃 (𝑠 i n) 𝑗$ corresponding to the $4 𝑁$ Lagrange multipliers introduced. We hence obtain the saddle point equations,

⟨ 𝑡 𝑖 𝑗 ⟩ = 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 e x p (𝑒 𝜃 𝑧 𝑖 𝑤 𝑗) 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃, 𝜎 𝑡 𝑖 𝑗 = ⟨ 𝑡 𝑖 𝑗 ⟩ (𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 + 1 - ⟨ 𝑡 𝑖 𝑗 ⟩), ˆ 𝑝 𝑖 𝑗 = 1 𝛥 𝑖 𝑗 𝑥 𝑖 𝑦 𝑗 (𝑒 𝑧 𝑖 𝑤 𝑗 𝑒 𝜃 - 1), 𝜎 𝛩 (𝑡 𝑖 𝑗) = ˆ 𝑝 𝑖 𝑗 (1 - ˆ 𝑝 𝑖 𝑗) .

(56)

Note that the expressions found in the previous cases are particular examples of this general problem and can be readily recovered by removing the appropriate constraints, i.e., making the Lagrange multipliers equal to zero, which in this case is equivalent to setting $𝑥 𝑖 = 𝑦 𝑗 = 1 \forall 𝑖, 𝑗$ or $𝑧 𝑖 = 𝑤 𝑗 = 1 \forall 𝑖, 𝑗$ or both.

We can revisit the case where only the strength sequence is fixed. Now, although the resulting statistics are Poisson and not multinomial, it can be proved that in the thermodynamic limit both descriptions are equivalent (see Appendix ). Additionally, in such a case one can obtain the statistics of the binary projection of the occupation numbers $ˆ 𝑝 𝑖 𝑗 = 1 - 𝑒 - ⟨ 𝑡 𝑖 𝑗 ⟩$ .

Unfortunately, the explicit form of the Lagrange multipliers for the degrees or the strengths cannot be solved, since the uncorrelated approximation is no longer valid,

𝑥 𝑖 𝑦 𝑗 (𝑒 𝑒 𝜃 𝑧 𝑖 𝑤 𝑗 - 1) ≃ 𝑒 𝑒 𝜃 𝑧 𝑖 𝑤 𝑗 𝑥 𝑖 𝑦 𝑗, if 𝑥 𝑖 𝑦 𝑗 𝑒 𝑒 𝜃 𝑧 𝑖 𝑤 𝑗 ≪ 1 \Rightarrow ˆ 𝑝 𝑖 𝑗 ≃ 𝑥 𝑖 𝑦 𝑗 𝑒 𝑒 𝜃 𝑧 𝑖 𝑤 𝑗 .

(57)

In this last expression the factorization of the connection probability in two node-dependent magnitudes is impossible, despite the approximation assumed. Hence one sees again that there is no way of generating uncorrelated networks at the level of degrees under the strict set of constraints considered.

VI. CONCLUSIONS

The present work deals with the statistical framework of multiedge networks, which in contrast to weighted networks, are entities formed by distinguishable events that can be grouped in edges weighted by integer numbers. It introduces flexible analytical null models that provide expectations for the main variables of a system represented in terms of a complex network (the entries of its adjacency matrix) thus complementing the existing studies on weighted networks. The decision upon which model to take depends on the kind of (physical) process generating the network at study.

We have started by properly defining the differences between weighted and multiedge networks based on the distinguishability or not of the elements forming a network. We have then properly set up a framework of multiedge networks in a statistical mechanics approach by defining ensembles of networks fulfilling some given sets of general constraints. Such a framework has been used to obtain analytical expressions considering cases with general constraints depending linearly on the number of events allocated to each edge as well as their binary projections, considering only the existence or not of a given edge. Some common interesting cases have been explicitly developed as examples such as fixed strength sequence, degree sequence, cost, strength sequence plus cost, total number of binary edges, and both strength and degree fixed sequences. Previous results found in the literature have been recovered, specifically the correlations of the configurational model (for binary constraints on the degree sequence) and the absence of correlation between occupation numbers and degree once the degree sequence is fixed among others. Our treatment uncovers explicit relations between the binary occupation probability of an edge and its expected occupation number, which allows one to fully characterize binary networks as projections of multiedge graph instances. These results permit an extensive treatment of the finite size effects present in this kind of network, since they are fully valid both in a dense scenario and intermediate cases.

Furthermore, we have explicitly derived general forms for the probability of obtaining a given graph with given constraints in terms of its adjacency matrix elements statistics, which can be understood as the canonical and grand-canonical description of the considered ensembles. Such expressions permit explicit entropy measures, useful for the analysis of similar systems for data inference . Additionally, the analytical character of the presented methodology allows for the explicit calculation of some network metrics using probabilistic examples such a clustering and strength-degree correlations, for instance. As a complement we have also introduced the main ideas which can lead to the efficient generation of multiedge graphs though its details are left for development in future work.

The applications of the theory developed here can be used in a wide variety of fields, mainly all those network studies applied to data generated by agent-based systems and in particular in the very active transportation research area, where networks are the result of complex, time-dependent, mobility patterns . It also opens the door to extensions of this kind of network to multiple layer scenarios .

Notes

¹ Please note that the framework we present permits self-edges if one desires; it is just a matter of extending the sums over available states over all values of $𝑖, 𝑗$ or excluding the term $𝑖 \neq 𝑗$ . It also is valid for both directed and undirected networks: One needs to perform such sums only over states for which $𝑗 < 𝑖$ in such cases.

² In the remainder of the paper, the $*$ signs will be omitted to simplify notation.

³ As done in for the binary case.

ACKNOWLEDGMENTS

This work has been partially supported by the Spanish DGICYT Grant No. FIS2009-13730-C02-01 and by the Generalitat de Catalunya 2009-SGR-00838. O.S. has been supported by the Generalitat de Catalunya through the FI Program.

Appendices

Appendix A: CUMULANT GENERATING FUNCTIONS

In this section we present the cumulant generating functions for the different models proposed (Bernoulli, multinomial, Poisson, and zero-inflated Poisson) and show that their close relationship to the expressions developed in the main text.

The cumulant generating function of the probability distribution $𝑃 (ℎ)$ of a variable $ℎ$ is defined as $𝐾 (ℎ, 𝑥) = l n 𝑀 (ℎ, 𝑥)$ , being $𝑀 (ℎ, 𝑥) = ⟨ 𝑒 ℎ 𝑥 ⟩$ its moment generating function and $𝑥$ and auxiliary field. Once $𝐾 (ℎ, 𝑥)$ is known, all central cumulants $𝜅 𝑘$ can be obtained by derivation, uniquely determining the distribution.

𝜅 𝑘 = 𝜕 𝑘 𝐾 (ℎ, 𝑥) | 𝑥 = 0 .

(A1)

Note that if we consider the joint distribution of two (or more) independent variables $ℎ 1, ℎ 2$ , being it a product of the individual distributions $𝑃 (ℎ 1), 𝑃 (ℎ 2)$ , then $𝑀 1, 2 (ℎ 1, ℎ 2, 𝑥 1, 𝑥 2) = 𝑀 1 (ℎ 1, 𝑥 1) 𝑀 2 (ℎ 2, 𝑥 2)$ and finally its joint cumulant generating function factorizes in the sum $𝐾 12 (ℎ 1, ℎ 2, 𝑥 1, 𝑥 2) = 𝐾 1 (ℎ 1, 𝑥 1) + 𝐾 2 (ℎ 2, 𝑥 2)$ .

Having introduced that, if we start at the microcanonical level [Eq. ] and identify

𝑃 ({𝑡 𝑖 𝑗}) = 𝛺 - 1 \prod 𝑞 𝛿 (𝐶 𝑞 - 𝐶 𝑞 ({𝑡 𝑖 𝑗})) = {\sum {𝑐} 𝑄 \prod 𝑞 𝛿 (𝐶 𝑞 - 𝐶 𝑞 ({𝑡 𝑖 𝑗}))} - 1 \times \prod 𝑞 𝛿 (𝐶 𝑞 - 𝐶 𝑞 ({𝑡 𝑖 𝑗})) \equiv 𝐶 - 1 \prod 𝑞 𝛿 (𝐶 𝑞 - 𝐶 𝑞 ({𝑡 𝑖 𝑗})), l n 𝛺 ({ℎ 𝑖 𝑗}) = l n {𝐶 \sum {𝑐} 𝑃 ({𝑡 𝑖 𝑗}) 𝑒 ℎ 𝑖 𝑗} = l n 𝐶 + 𝐾 ({𝑡 𝑖 𝑗}, {ℎ 𝑖 𝑗}),

(A2)

we clearly see that

𝜕 ℎ 𝑖 𝑗 l n 𝛺 ({𝑡 𝑖 𝑗}, {ℎ 𝑖 𝑗}) = 𝜕 ℎ 𝑖 𝑗 𝐾 ({𝑡 𝑖 𝑗}, {ℎ 𝑖 𝑗}),

(A3)

and the relation between both objects is apparent.

Starting at the level where only the number of events $𝑇$ is fixed (where all the distributions reduce to a multinomial form with associated probabilities ${𝑝 𝑖 𝑗}$ ), we take Eq. and perform the saddle point approximation to obtain

l n 𝛺 ≃ 𝑇 l n \sum 𝑖 𝑗 e x p (ℎ 𝑖 𝑗 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞) + 𝐹 ({𝜃 𝑞}, {𝐶 𝑞}) = 𝐻 ({ℎ 𝑖 𝑗}) + 𝐹 ({𝜃 𝑞}, {𝐶 𝑞}) .

(A4)

Despite having additional terms, with regards to differentiation with respect to $ℎ 𝑖 𝑗$ (needed to recover the moments of the distribution of $𝑡 𝑖 𝑗$ ), the form obtained is always of the type $𝐻 ({ℎ 𝑖 𝑗}) = 𝑇 l n \sum 𝜋 𝑖 𝑗 𝑒 ℎ 𝑖 𝑗$ (since $𝐹$ is a function not depending on the auxiliary fields ${ℎ 𝑖 𝑗}$ ) and hence the underlying statistic is multinomial because $\sum 𝑖 𝑗 𝜋 𝑖 𝑗 = 1$ .

Considering now only the most general case in Eq. we have again (adding external fields for the binary projection ${𝑏 𝑖 𝑗}$ ),

l n 𝛺 ≃ \sum 𝑖 𝑗 l n {𝑒 𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 ℎ 𝑖 𝑗 + \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞 - 1) + 1} + 𝐹 ({𝜃 𝑞}, {𝜆 𝑞'}, {𝐶 𝑞}, {𝐶 𝑞'}) .

(A5)

Dropping the constraints on the binary projection, i.e., $𝜆 𝑞' = 𝑏 𝑖 𝑗 = 0 \forall 𝑞, 𝑖 𝑗$ yields

𝐻 ({ℎ 𝑖 𝑗}) = \sum 𝑖, 𝑗 𝑒 ℎ 𝑖 𝑗 𝜇 𝑖 𝑗,

(A6)

where we identified $𝜇 𝑖 𝑗 = 𝑒 \sum 𝑞 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞$ . This expression (up to derivation with respect to $ℎ 𝑖 𝑗$ ) has the same form as a sum of independent Poisson cumulant generating functions ( $𝜇 𝑒 ℎ 𝑖 𝑗$ ).

Dropping the constraints on the multilink nature of the network $𝜃 𝑞 = ℎ 𝑖 𝑗 = 0 \forall 𝑞, 𝑖 𝑗$ (except the one on the total number of multilinks, $𝜃$ ), we have

𝐻 ({𝑞 𝑖 𝑗}) = \sum 𝑖, 𝑗 {l n (ˆ 𝑝 𝑖 𝑗 𝑒 𝑏 𝑖 𝑗 + (1 - ˆ 𝑝 𝑖 𝑗)) - l n (1 - ˆ 𝑝 𝑖 𝑗)},

(A7)

where we identified $ˆ 𝑝 𝑖 𝑗$ from Eq. . And we see that the prior expression is closely related to the cumulant generating function of independent Bernoulli processes (with respect again to derivation on ${𝑞 𝑖 𝑗}$ terms).

Finally, the most general case can be mapped to a mixed zero-inflated Poisson process as we shall prove: Imagine the outcome of a process in which we sort $𝑁 (𝑁 - 1)$ independent Bernoulli processes and from the result of it, if the outcome is positive, we sort a Poisson process on top of it (discarding the no-occurrence event). Since the processes are independent, we shall consider a single one of them and then write the overall probability of the events as the product of the different probabilities $𝑃 (𝑡)$ . The associated probability of the event just described is

𝑃 𝑏 - 𝑝 (𝑡) = (1 - 𝑝) 𝛩 (𝑡) (𝑝 𝑒 𝜇 - 1 𝜇 𝑡 𝑡!) 𝛩 (𝑡),

(A8)

which represents a probability measure over an integer quantity. We can compute the mean and variance of $𝑡$ yielding

⟨ 𝑡 ⟩ = 0 + 𝑝 𝑒 𝜇 - 1 \infty \sum 𝑡 = 1 𝑡 𝜇 𝑡 𝑡! = 𝑝 𝜇 𝑒 𝜇 𝑒 𝜇 - 1, 𝜎 2 𝑡 = ⟨ 𝑡 ⟩ (1 + 𝜇 - ⟨ 𝑡 ⟩), ⟨ 𝛩 (𝑡) ⟩ = 𝑝; 𝜎 2 𝛩 (𝑡) = 𝑝 (1 - 𝑝) .

(A9)

The obtained expressions need to be compared with , which allows one to identify $𝜇 = 𝑒 𝜃 \prod 𝑞 𝑒 𝜃 𝑞 𝑐 (𝑖 𝑗) 𝑞$ and $𝑝 = ˆ 𝑝 𝑖 𝑗$ . Moreover, concerning the cumulant generating function, one finds

l n ⟨ 𝑒 ℎ 𝑡 ⟩ = l n \infty \sum 𝑡 = 0 𝑃 𝑏 - 𝑝 (𝑡) 𝑒 ℎ 𝑡 = l n {1 - 𝑝 + 𝑝 𝑒 𝜇 - 1 (𝑒 𝜇 𝑒 ℎ - 1)} = l n (1 + 𝑝) + l n {1 + 𝑝 1 - 𝑝 (𝑒 𝜇 𝑒 ℎ - 1) 𝑒 𝜇 - 1},

(A10)

which is identical to the argument in the sum of Eq. (except for a linear constant) and captures the more general case considered.

Appendix B: ENSEMBLE EQUIVALENCE AND GRAPH GENERATION

Throughout this paper we have uncovered the mathematical expressions allowing one to generate networks under different ensembles using a probabilistic framework over $𝑡 𝑖 𝑗$ . Explicitly they can be summarized:

(1) Canonical ensemble (linear constraints on $𝑡 𝑖 𝑗$ ):

$𝑃 (𝐓 | 𝑇, {𝜃 𝑞}) = 𝑇! \prod 𝑖 𝑗 𝑡 𝑖 𝑗! \prod 𝑖 𝑗 𝑝 𝑡 𝑖 𝑗 𝑖 𝑗 .$

(B1)
(2) Grandcanonical ensemble [linear constraints on $𝛩 (𝑡 𝑖 𝑗)$ and/or on $𝑡 𝑖 𝑗$ ]:

$𝑃 (𝐓 | {ˆ 𝑝 𝑖 𝑗, 𝜇 𝑖 𝑗}) = \prod 𝑖, 𝑗 (1 - ˆ 𝑝 𝑖 𝑗) 1 - 𝛩 (𝑡 𝑖 𝑗) {ˆ 𝑝 𝑖 𝑗 𝑒 𝜇 𝑖 𝑗 - 1 𝜇 𝑡 𝑖 𝑗 𝑖 𝑗 𝑡 𝑖 𝑗!} 𝛩 (𝑡 𝑖 𝑗) .$

(B2)
(3) Microcanonical ensemble. This ensemble can be used by generating sequences of ${𝑡 𝑖 𝑗}$ using the two above expressions and discarding those not corresponding exactly with the imposed constraints.

We have shown already that the relative fluctuations of the linear constraints on the occupation numbers $𝑡 𝑖 𝑗$ vanish in the thermodynamic limit and that the binary depending constraints are nonvanishing in this limit. We finally prove here that the grand-canonical ensemble considering only linear constraints on occupation numbers is strictly equivalent to the canonical and the microcanonical in the thermodynamic limit.

To do so, we make use of the properties of the multinomial distribution to recover it under the cases where no constraints or only linear constraints on $𝑡 𝑖 𝑗$ are imposed. In this case $𝑝 𝑖 𝑗 = 1 - 𝑒 - 𝜇 𝑖 𝑗$ and reduce to the product of Poisson distributions with different mean parameters $𝜇 𝑖 𝑗$ .

The outcome of a process in which we sort different independent Poisson variables of parameters ${𝜇 𝑖 𝑗}$ can be equivalently expressed as a product of a multinomial process of $⟨ 𝑇 ⟩ = \sum 𝜇 𝑖 𝑗$ multinomial trials with associated probabilities ${𝑝 𝑖 𝑗 = 𝜇 𝑖 𝑗 ⟨ 𝑇 ⟩}$ . Hence we have that

𝑃 (𝐓) = Mult ({𝑝 𝑖 𝑗} | 𝑇) Pois (𝜆 𝑇 = \sum 𝜇 𝑖 𝑗) .

(B3)

And since the resulting occupation numbers statistics derive to a Poisson distribution, which has vanishing fluctuations on the thermodynamic limit ( $𝜎 2 𝑡 𝑖 𝑗 / 𝜇 2 𝑖 𝑗 = 𝜇 - 1 𝑖 𝑗 \propto 𝑇 - 1$ , then the equivalence between the two presented calculations is completely proved.

Appendix C: FINITE T IN EQ.

Let us assume that instead of having $𝑇 \to \infty$ in we consider $𝑇$ finite. Therefore,

𝛺 = \int 𝑒 - 𝜃 𝑇 𝑒 - \sum 𝑞' 𝜆 𝑞' 𝐶 𝑞' 𝐷 [𝜃] (\prod 𝑞' 𝐷 [𝜆 𝑞']) \times 𝑇! [𝑒 𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆' 𝑞 ̂ 𝑐 (𝑖 𝑗) 𝑞' 𝑇 \sum 𝑡' 𝑖 𝑗 = 1 𝑒 (ℎ 𝑖 𝑗 + 𝜃) 𝑡' 𝑖 𝑗 𝑡' 𝑖 𝑗! + 1 0!], \times \sum \sum 𝑡 𝑖 𝑗 = 𝑇 - 𝑡' 𝑖 𝑗 \prod (𝑒 (𝑏 𝑖 𝑗 + \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞') 𝛩 (𝑡 𝑖 𝑗) + ℎ 𝑖 𝑗 + 𝜃) 𝑡 𝑖 𝑗 𝑡 𝑖 𝑗!,

(C1)

and we now perform the finite sum inside the integral,

𝑇 \sum 1 𝑧 𝑡 𝑡! = 𝛤 (𝑇 + 1, 𝑧) 𝑒 𝑧 𝛤 (𝑇 + 1) .

(C2)

And for the occupation numbers and edge probability obtain

⟨ 𝑡 𝑖 𝑗 ⟩ = 1 𝛥 𝑖 𝑗 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' 𝑒 𝜃 (𝑒 𝑒 𝜃 𝛤 (𝑇 + 1, 𝑒 𝜃) 𝛤 (𝑇 + 1) - 𝑒 𝑇 𝜃 𝛤 (𝑇 + 1)), ˆ 𝑝 𝑖 𝑗 = 1 𝛥 𝑖 𝑗 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖 𝑗) 𝑞' (𝑒 𝑒 𝜃 𝛤 (𝑇 + 1, 𝑒 𝜃) 𝛤 (𝑇 + 1) - 1), 𝛥 𝑖 𝑗 = 𝑒 \sum 𝑞' 𝜆 𝑞' ̂ 𝑐 (𝑖, 𝑗) 𝑞' (𝑒 𝑒 𝜃 - 1) + 1,

(C3)

which converge extremely quickly to the obtained results for $𝑇 \to \infty$ (see Fig. ).

FIG. 3

$\sum 𝑇 1 𝜌 𝑡 𝑡! 𝑒 𝜌 = 𝛤 (𝑇 + 1, 𝜌 (‾ 𝑡 +)) 𝛤 (𝑇 + 1)$ for different values of $‾ 𝑡 +$ and $𝑇$ , we observe the extreme rapid convergence to unity as $𝑇$ grows (note that $𝑇 \geq ‾ 𝑡 + = 𝑇 / 𝐸$ ).

Appendix D: ADDITIONAL EXAMPLES. FIXED BINNED DISTRIBUTION OF COSTS ${𝑐 𝑛, 𝑁 𝑐 𝑛}$

We here report an additional example which may be of interest for studies on transportation origin-destination matrices, where some forms of trip-cost distribution have been discussed .

Starting from Sec. (), it is a matter of considering the additional term,

e x p (\sum 𝑛 𝜅 𝑛 𝜉 𝑛 (𝑑 𝑖 𝑗)),

on the equations, where $𝜅 𝑛$ are additional Lagrange multipliers satisfying that

𝑁 𝑛 = \sum 𝑖 𝑗 𝜅 𝑛 𝜉 𝑛 (𝑑 𝑖 𝑗) 𝑡 𝑖 𝑗,

where $𝑁 𝑛$ is the number of trips whose distance is in the interval $[𝑑 𝑛 - 1, 𝑑 𝑛)$ and $𝜉 𝑛$ is the indicator function of such an event. The size of the bins needs to be chosen in an appropriate manner as to give consistency to the distribution obtained. Such an example is particularly important for its importance to assess whether an observed O-D is caused by a particular tendency of the agents that create it to move or conversely, the space where the move shapes the form of the obtained O-D.

The expressions of $𝑝 𝑖 𝑗$ in this case read

𝑝 𝑖 𝑗 = 𝑒 \sum 𝑛 𝜅 𝑛 𝜉 (𝑑 𝑖 𝑗) \sum 𝑖 𝑗 𝑒 \sum 𝑛 𝜅 𝑛 𝜉 (𝑑 𝑖 𝑗) .

And the prior considerations are also valid. Note also that keeping the multinomial framework, all the quantities considered are intensive, while the expected strength and average occupation numbers remain extensive variables. Additionally, in the thermodynamic limit these quantities have vanishing relative fluctuations.

Statistical mechanics of multiedge networks

Abstract

Article Text

I. INTRODUCTION

II. ENSEMBLE APPROACH TO MULTIEDGE NETWORKS: NOMENCLATURE AND DEFINITIONS

III. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS ON THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$

A. No constraints

B. Fixed average event cost $‾ 𝑐 = 𝐶 / 𝑇$

C. Fixed relative strength sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, 𝑇$

IV. MULTIEDGE NETWORK WITH GIVEN NONLINEAR CONSTRAINTS ON THE BINARY PROJECTION OF THE OCCUPATION NUMBERS $𝛩 (𝑡 𝑖 𝑗)$

A. Fixed $𝐸$ , $𝑇$

B. Fixed degree sequence $⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}, 𝑇$

V. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS DEPENDING ON BOTH THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$ AND THEIR BINARY PROJECTION $𝛩 (𝑡 𝑖 𝑗)$

A. Fixed relative strength sequence and fixed degree sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, ⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}$

VI. CONCLUSIONS

Notes

ACKNOWLEDGMENTS

Appendices

Appendix A: CUMULANT GENERATING FUNCTIONS

Appendix B: ENSEMBLE EQUIVALENCE AND GRAPH GENERATION

Appendix C: FINITE T IN EQ.

Appendix D: ADDITIONAL EXAMPLES. FIXED BINNED DISTRIBUTION OF COSTS ${𝑐 𝑛, 𝑁 𝑐 𝑛}$

References (40)

Statistical mechanics of multiedge networks

Abstract

Article Text

I. INTRODUCTION

II. ENSEMBLE APPROACH TO MULTIEDGE NETWORKS: NOMENCLATURE AND DEFINITIONS

III. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS ON THE OCCUPATION NUMBERS 𝑡𝑖⁢𝑗

A. No constraints

B. Fixed average event cost ‾𝑐=𝐶/𝑇

C. Fixed relative strength sequence ⃗𝑠={(𝑠out,𝑠in)𝑖},𝑇

IV. MULTIEDGE NETWORK WITH GIVEN NONLINEAR CONSTRAINTS ON THE BINARY PROJECTION OF THE OCCUPATION NUMBERS 𝛩⁡(𝑡𝑖⁢𝑗)

A. Fixed 𝐸,𝑇

B. Fixed degree sequence ⃗𝑘={(𝑘out,𝑘in)𝑖},𝑇

V. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS DEPENDING ON BOTH THE OCCUPATION NUMBERS 𝑡𝑖⁢𝑗 AND THEIR BINARY PROJECTION 𝛩⁡(𝑡𝑖⁢𝑗)

A. Fixed relative strength sequence and fixed degree sequence ⃗𝑠={(𝑠out,𝑠in)𝑖},⃗𝑘={(𝑘out,𝑘in)𝑖}

VI. CONCLUSIONS

Notes

ACKNOWLEDGMENTS

Appendices

Appendix A: CUMULANT GENERATING FUNCTIONS

Appendix B: ENSEMBLE EQUIVALENCE AND GRAPH GENERATION

Appendix C: FINITE T IN EQ. (27)

Appendix D: ADDITIONAL EXAMPLES. FIXED BINNED DISTRIBUTION OF COSTS {𝑐𝑛,𝑁𝑐𝑛}

References (40)

III. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS ON THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$

B. Fixed average event cost $‾ 𝑐 = 𝐶 / 𝑇$

C. Fixed relative strength sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, 𝑇$

IV. MULTIEDGE NETWORK WITH GIVEN NONLINEAR CONSTRAINTS ON THE BINARY PROJECTION OF THE OCCUPATION NUMBERS $𝛩 (𝑡 𝑖 𝑗)$

A. Fixed $𝐸$ , $𝑇$

B. Fixed degree sequence $⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}, 𝑇$

V. MULTIEDGE NETWORK WITH GIVEN LINEAR CONSTRAINTS DEPENDING ON BOTH THE OCCUPATION NUMBERS $𝑡 𝑖 𝑗$ AND THEIR BINARY PROJECTION $𝛩 (𝑡 𝑖 𝑗)$

A. Fixed relative strength sequence and fixed degree sequence $⃗ 𝑠 = {(𝑠 o u t, 𝑠 i n) 𝑖}, ⃗ 𝑘 = {(𝑘 o u t, 𝑘 i n) 𝑖}$

Appendix C: FINITE T IN EQ.

Appendix D: ADDITIONAL EXAMPLES. FIXED BINNED DISTRIBUTION OF COSTS ${𝑐 𝑛, 𝑁 𝑐 𝑛}$