Belief-propagation algorithm and the Ising model on networks with arbitrary distributions of motifs

# Belief-propagation algorithm and the Ising model on networks with arbitrary distributions of motifs

S. Yoon Departamento de Física da Universidade de Aveiro, I3N, 3810-193 Aveiro, Portugal    A. V. Goltsev Departamento de Física da Universidade de Aveiro, I3N, 3810-193 Aveiro, Portugal A. F. Ioffe Physico-Technical Institute, 194021 St. Petersburg, Russia    S. N. Dorogovtsev Departamento de Física da Universidade de Aveiro, I3N, 3810-193 Aveiro, Portugal A. F. Ioffe Physico-Technical Institute, 194021 St. Petersburg, Russia    J. F. F. Mendes Departamento de Física da Universidade de Aveiro, I3N, 3810-193 Aveiro, Portugal
July 14, 2019
###### Abstract

We generalize the belief-propagation algorithm to sparse random networks with arbitrary distributions of motifs (triangles, loops, etc.). Each vertex in these networks belongs to a given set of motifs (generalization of the configuration model). These networks can be treated as sparse uncorrelated hypergraphs in which hyperedges represent motifs. Here a hypergraph is a generalization of a graph, where a hyperedge can connect any number of vertices. These uncorrelated hypergraphs are tree-like (hypertrees), which crucially simplifies the problem and allows us to apply the belief-propagation algorithm to these loopy networks with arbitrary motifs. As natural examples, we consider motifs in the form of finite loops and cliques. We apply the belief-propagation algorithm to the ferromagnetic Ising model with pairwise interactions on the resulting random networks and obtain an exact solution of this model. We find an exact critical temperature of the ferromagnetic phase transition and demonstrate that with increasing the clustering coefficient and the loop size, the critical temperature increases compared to ordinary tree-like complex networks. However, weak clustering does not change the critical behavior qualitatively. Our solution also gives the birth point of the giant connected component in these loopy networks.

###### pacs:
05.10.-a, 05.40.-a, 05.50.+q, 87.18.Sn

## I Introduction

The belief-propagation algorithm is an effective numerical method for solving inference problems on sparse graphs. It was originally proposed by J. Pearl Pearl:p88 () for tree-like graphs. The belief-propagation algorithm was applied to study diverse systems in computer science, physics, and biology. Among its numerous applications are computer vision problems, decoding of high performance turbo codes and many others, see Frey:98 (); McEliece:mmc98 (). Empirically, it was found that this algorithm works surprisingly good even for graphs with loops. J. S. Yedidia et al. Yedidia:yfw01 () found that the belief-propagation algorithm actually coincides with the minimization of the Bethe free energy. This discovery renewed interest in the Bethe-Peierls approximation and related methods Pretti:pp03 (); Mooij:mk05 (); Hartmann:hw05 (). In statistical physics the belief-propagation algorithm is equivalent to the so-called cavity method proposed by Mézard, Parisi, and Virasoro mp2003 (). The belief-propagation algorithm is applicable to models with both discrete and continuous variables mooij04 (); ohkubo05 (); Dorogovtsev:dgm08 (). Recently the belief-propagation algorithm was applied to diverse problems in random networks: spread of disease karrer10 (), the calculation of the size of the giant component sk10 (), counting large loops in directed random uncorrelated networks bianconi08 (), the graph bi-partitioning problem sulc10 (), analysis of the contributions of individual nodes and groups of nodes to the network structure ras10 (), and networks of spiking neurons smd09 ().

These investigations showed that the belief-propagation algorithm is exact for a model on a graph with locally tree-like structure in the infinite size limit. However, only an approximate solution was obtained by use of this method for networks with short loops. Real technological, social and biological networks have numerous short and large loops and other complex subgraphs or motifs which lead to essentially non-tree-like neighborhoods ab01a (); dm01c (); Dorogovtsev:2010 (); Newman:n03a (); Milo2002 (); Milo2004 (); Sporn2004 (); Alon2007 (). That is why it is important to develop a method which takes into account finite loops and more complex motifs. Different ways were proposed recently to compute loop corrections to the Bethe approximation by use of the belief-propagation algorithm mr05 (); ps05 (). These methods, however, are exact only if a graph contains a single loop. Loop expansions were proposed in Refs. Chertkov:cc06a (); Chertkov:cc06b (). However, they were not applied to complex network yet.

Recently, Newman and Miller newman:n09 (); m09 (); kn10b () independently introduced a model of random graphs with arbitrary distributions of subgraphs or motifs. In this natural generalization of the configuration model, each vertex belongs to a given set of motifs [e.g., vertex is a member of motifs of type , motifs of type , and so on], and apart from this constraint, the network is uniformly random. In the original configuration model, the sequence of vertex degrees is fixed, where . In this generalization, the sequence of generalized vertex degrees is given. For example, motifs can be triangles, loops, chains, cliques (fully connected subgraphs), single edges that do not enter other motifs, and, in general, arbitrary finite clusters. The resulting networks can be treated as uncorrelated hypergraphs in which hyperedges represent motifs. In graph theory, a hypergraph is a generalization of a graph, where a hyperedge can connect any number of vertices Berge73 (). Because of the complex motifs, the large sparse networks under consideration can have loops, clustering, and correlations, but the underlying hypergraphs are locally tree-like (hypertrees) and uncorrelated similarly to the original sparse configuration model. Our approach is based on the hypertree-like structure of these highly structured sparse networks. To demonstrate our approach, we apply the generalized belief-propagation algorithm to the ferromagnetic Ising model on the sparse networks in which motifs are finite loops or cliques. We obtain an exact solution of the ferromagnetic Ising model and the birth point of the giant connected component in this kind of highly structured networks. Note that the Ising model on these networks is a more complex problem than the percolation problem because we must account for spin interactions between spins both inside and between motifs. We find an exact critical temperature of the ferromagnetic phase transition and demonstrate that finite loops increase the critical temperature in comparison to ordinary tree-like networks.

## Ii Ensemble of random networks with motifs

Let us introduce a statistical ensemble of random networks with a given distribution of motifs in which each vertex belongs to a given set of motifs and apart of this constraint, the networks are uniformly random. These networks can be treated as uncorrelated hypergraphs, in which motifs play a role of hyperedges, so the number of hyperdegrees are equal to the number of specific motifs attached to a vertex. In principal, one can choose any subgraph with an arbitrary number of vertices as a motif, see Fig. 1. In the present paper, for simplicity, we only consider simple motifs such as single edges, finite loops, and cliques.

Let us first note how one can describe a statistical ensemble of random networks with the simplest motifs, namely simple edges (see, for example, Refs. Dorogovtsev:dgm08 (); Bianconi09 () and references therein). We define the probability that an edge between vertices and is present () or absent (),

 p2(aij)=⟨Q2⟩N−1δ(aij−1)+(1−⟨Q2⟩N−1)δ(aij) (1)

where are entries of the symmetrical adjacency matrix, is the mean number of edges attached to a vertex. It is well known that the probability of the realization of the Erdős-Rényi graph with a given adjacency matrix , is the product

 G2({aij})=N−1∏i=1N∏j=i+1p2(aij). (2)

The degree distribution is the Poisson distribution. In the configuration model, the probability of the realization of a given graph with a given sequence of degrees, , , ,…, , is

 G2({aij})=1AN∏i=1δ(Q2(i)−N∑j=1aij)∏i

The delta-function fixes the number of edges attached to vertex . is a normalization constant. The distribution function of degrees is determined by the sequence ,

 P2(Q2)=1N∑iδ(Q2−Q2(i)). (4)

The second simplest motif, the triangle, plays the role of a hyperedge that interconnects a triple of vertices. Let us introduce the probability that a hyperedge (triangle) among vertices , , and is present or absent, i.e., or , respectively:

 p3(aijk)=pδ(aijk−1)+(1−p)δ(aijk), (5)

where

 p=2⟨Q3⟩(N−1)(N−2) (6)

is the probability that vertices and form a triangle. is the mean number of triangles attached to a randomly chosen vertex. are entries of the adjacency matrix of the hypergraph. This matrix is symmetrical with respect to permutations of the indices and ,   .

For example, one can introduce the ensemble of the Erdős-Rényi hypergraphs. Given that the matrix elements are independent and uncorrelated random parameters, the probability of realization of a graph with a given adjacency matrix is the product of probabilities over different triples of vertices:

 G3({aijk})=N−2∏i=1N−1∏j=i+1N∏k=j+1p3(aijk). (7)

This function describes the ensemble of the Erdős-Rényi hypergraphs with the Poisson degree distributions of the number of triangles attached to vertices,

 P3(Q3)=(⟨Q3⟩)Q3Q3!e−⟨Q3⟩. (8)

In the configuration model, the probability of the realization of a given graph with a sequence of the number of triangles, , , ,…,, attached to vertices is defined by an equation,

 G3({aijk})=1AN∏i=1δ(Q3(i)−12∑j,kaijk)∏i

Here is given by Eq. (5). The delta-function fixes the number of triangles attached to vertex . is a normalization constant. The distribution function of triangles is

 P3(Q3)=1N∑iδ(Q3−Q3(i)). (10)

One can further generalize equations Eqs. (5) and (6) and introduce the probability that vertices and form a hyperedge of size 4 (a clique or loop of size 4). In the case of loops, one can arrange these vertices in order of increasing vertex index. Then, for a given sequence of squares or cliques, , attached to vertices , one can introduce the probability of the realization of a given graph with a given sequence similar to Eq. (9), and so on. The network ensemble of the configuration model with a given sequences of edges , triangles , squares and other motifs is described by a product of the corresponding probabilities,

 G({aij},{aijk},{aijkl},...)=G2({aij})G3({aijk})... (11)

The average of a quantity over the network ensemble is

 ⟨K⟩en=∫K({aij},{aijk},{aijkl}...)× G({aij},{aijk},{aijkl},...)∏i

where we integrate over all possible edges and hyperedges.

One can prove that the probability that different motifs have a common edge, i.e., they are overlapping, tends to zero in the infinite size limit . For example, the total number of edges that overlap with triangles equals

 ⟨12∑i,j,kaijaijk⟩en=N⟨Q2⟩⟨Q3⟩N−1. (13)

This number of double connections is finite in the limit . Because the total number of edges is of the order of , the fraction of double connections tends to zero in the thermodynamic limit. Thus, one can neglect overlapping among different motifs. Using the standard methods in complex network theory Bianconi:bc03 (); Bianconi:bm05a (); Newman:n03b (), one can show that in the configuration model the number of finite loops formed by hyperedges is finite in the thermodynamic limit. Therefore, sparse random uncorrelated hypergraphs have a hypertree-like structure.

Uncorrelated random hypergraphs are described by a distribution function which is the probability that a randomly chosen vertex has edges, triangles, squares, and other motifs attached to this vertex. In the case of uncorrelated motifs, we have

 P(Q2,Q3,Q4,…)=∞∏ℓ=2Pℓ(Qℓ). (14)

Here, is the distribution function of loops of size ,

 Pℓ(Qℓ)=1N∑iδ(Qℓ−Qℓ(i)). (15)

In this kind of network, the number of nearest neighboring vertices of vertex is equal to

 Q(i)=Q2(i)+2Q3(i)+2Q4(i)+…. (16)

Therefore, the mean degree is

 ⟨Q⟩=1NN∑i=1Q(i)=∑Q2,Q3,…QP(Q2,Q3,Q4,…). (17)

The local clustering coefficient of node with degree , Eq. (16), is determined by the number of triangles as follows,

 C(i)=2Q3(i)Q(i)(Q(i)−1). (18)

Therefore, is maximum if a network consists only of triangles. In this case, we have and therefore

 C(i)=1Q(i)−1. (19)

On the other hand, if a network only consists of -cliques with a given (fully connected subgraphs of size ), then the local clustering coefficient of vertex with attached cliques equals

 C(i)=ℓ−2Q(i)−1, (20)

where is degree of vertex . If other motifs (edges, squares, and larger finite loops) are present in the network, then according to Eqs. (18) and (20) the local clustering coefficient is smaller than . In Refs. Serrano2006a (); Serrano2006b (); Serrano2006c (), it was shown that properties of networks with ”weak clustering”, , may differ from properties of networks with ”strong clustering” when decreases slower than . The latter networks are beyond the scope of the present article.

## Iii Belief-propagation algorithm

Let us consider the Ising model with pairwise interactions between nearest neighboring spins on a sparse random network with arbitrary distributions of motifs. The Hamiltonian of the model is

 E= −∑iHiSi−∑i

Here we sum the energies of spin clusters corresponding to motifs in the network. The second sum corresponds to simple edges. The third sum corresponds to triangles. The forth sum corresponds to motifs of size 4, and so on. is a local magnetic field at vertex . The coupling determines the energy of pairwise interaction between spins and . In general, one can introduce multi-spin interactions between spins in hyperedges, for example, , , and so on. However, this kind of interaction is beyond the scope of our article. Generally, the local fields and the couplings can be random quantities.

In order to solve this model, we generalize the belief-propagation algorithm. At first, we note the belief-propagation algorithm in application for the Ising model on tree-like uncorrelated complex networks without short loops (for more details, see Ref. Dorogovtsev:dgm08 ()). In this case, there are only edges determined by the adjacency matrix and the clustering coefficient is zero in the thermodynamic limit. Consider spin . A nearest neighboring spin sends a so-called message to spin that we denote as . This message is normalized,

 ∑Si=±1μji(Si)=1. (22)

Within the belief-propagation algorithm, the probability that spin is in a state is determined by the normalized product of incoming messages to spin and the probabilistic factor due to a local field :

 pi(Si)=1AeβHiSi∏jμji(Si), (23)

where is a normalization constant and is the reciprocal temperature. Now we can calculate the mean moment of spin ,

 ⟨Si⟩=∑Si=±1Sipi(Si). (24)

A message can be written in a general form,

 μji(Si)=exp(βhjiSi)/[2cosh(βhji)], (25)

Using this representation, we rewrite Eq. (24) in a physically clear form,

 ⟨Si⟩=tanh(βHi+β∑jajihji). (26)

This equation shows that a parameter plays the role of an effective field produced by spin at site . Messages obey a self-consistent equation that is called update rule,

 B∑Sj=±1e−βE(j,i)∏n≠iμnj(Sj)=μji(Si), (27)

where the index numerates nearest neighbors of vertex and is a normalization constant. The diagram representation of this equation is shown in Fig. 2 . The probabilistic factor takes account of the interaction energy of spins and and a local field ,

 E(j,i)=−HjSj−ajiJjiSiSj. (28)

Thus, a message from to is determined by the coupling between these spins and messages received by from its neighbors except . Multiplying Eq. (27) by and summing over , we obtain a self-consistent equation for the effective fields ,

 tanh(βhji)=tanh(βajiJji)tanh[β(Hj+∑n(≠i)anjhnj)]. (29)

In the thermodynamic limit, this equation is exact for a tree-like graph. The Bethe-Peierls approach and the Baxter recurrence method give exactly the same result Dorogovtsev:dgm08 (); bbook82 (). In numerical calculations, Eqs. (27) and (29) may be solved by use of iterations. In a general case, an analytical solution of these equations is unknown. In the case of all-to-all interactions with random couplings (the Sherrington-Kirkpatrick model), this set of equations is reduced to well-known TAP equations Thouless:tap77 () that are exact in the thermodynamic limit. These equations also give an exact solution of the ferromagnetic Ising model with uniform couplings on a random uncorrelated graph with zero clustering coefficient dgm02 (); lvvz02 (); Dorogovtsev:dgm08 ().

We now generalize the belief-propagation algorithm to the case of networks with given distributions of motifs described in Sec. II. In a network that consists only of edges, a message goes along an edge from the spin at one edge end to the spin at the other edge end. For the networks with motifs, it is natural to introduce a message that is sent by a motif attached to a vertex. This message goes to a spin to which this motif (hyperedge) is attached. Different motifs may be attached to a vertex, see Fig. 1. Let denote a cluster of size attached to vertex that we will call the destination vertex. Vertices together with the destination vertex form this motif. We introduce a message to spin from a motif . This message is normalized and can be written as follows,

 μMℓ(i)(Si)=exp(βhMℓ(i)Si)/[2cosh(βhMℓ(i))]. (30)

As above, the probability that spin is in a state is determined by the normalized product of incoming messages from motifs attached to spin and the probabilistic factor ,

 pi(Si)=1AeβHiSi∏{Mℓ(i)}μMℓ(i)(Si). (31)

Using Eq. (24), we find that the mean moment is determined by effective fields acting on spin from motifs , , attached to , see Fig. 1,

 ⟨Si⟩=tanh(βHi+β∑{Mℓ(i)}hMℓ(i)). (32)

Here the sum is taken over motifs attached to .

Let us find an update rule for a message from a given motif to vertex . We introduce the following update rule:

 B∑{Sj=±1}e−βE(Mℓ(i))∏j∏{Mn(j)≠Mℓ(i)}μMn(j)(Sj)=μMℓ(i)(Si). (33)

This rule shows that the message is equal to the product of incoming messages from motifs attached to all spins in the motif except spin and the motif itself (see Fig. 2. In Eq. (33) we average over all spin configurations of the spins , and in the motif. is a normalization constant. is an energy of the interaction between spins in the motif . This update rule also takes account of local fields acting on all spins except the destination spin . We note that Eq. (33) is valid for arbitrary motifs even with a complex internal structure. The only assumption is that a sparse network under consideration, in terms of hypergraphs, has a hypertree-like structure. Messages can be found numerically by use of iterations of Eq. (33) that start from an initial distribution of the messages.

For the sake of simplicity, let motifs be finite loops of size . Then

 E(Mℓ(i))=−ℓ−1∑n=1HjnSjn−ℓ−1∑n=0Jjnjn+1SjnSjn+1, (34)

where . If motifs are -cliques, then the energy takes account of interactions between all spins in these cliques,

 E(Mℓ(i))=−∑j(≠i)HjSj−12∑i,j∈Mℓ(i)Ji,jSiSj. (35)

Multiplying Eq. (33) by and summing over all spin configurations, we obtain an equation for the effective field ,

 tanh(βhMℓ(i))=⟨Si⟩Mℓ(i). (36)

The function on the right hand side is

 ⟨Si⟩Mℓ(i)=1Z(Mℓ(i))∑{Si,Sj1,...=±1}Sie−β˜E(Mℓ(i)). (37)

This function has a meaning of the mean moment of spin in the cluster with an energy

 ˜E(Mℓ(i))=−ℓ−1∑n=1Ht(jn)Sjn−ℓ∑n=0Jjnjn+1SjnSjn+1. (38)

This energy takes into account both the interaction between spins in this motif and total effective fields acting on these spins. In turn, the field acting on is the sum of a local field and effective fields produced by incoming messages from motifs attached to vertex except the motif , see Fig. 3,

 Ht(j)=Hj+∑{Mm(j)(≠Mℓ(i))}hMm(j). (39)

is the partition function of the cluster ,

 Z[Mℓ(i)]=∑{Si,Sj1,...=±1}e−β˜E(Mℓ(i)). (40)

Thus, the function is a function of the total fields acting on all spins in the motif except spin . For a given network with motifs, it is necessary to solve Eq. (36) with respect to the effective fields . One then can calculate the mean magnetic moments from Eq. (32). Note that the only condition we used to derive the equations above was hypertree-like structure of the networks. The absence of correlations in these hypertree-like networks was not needed. Equations (33) and  (36) are our main result that actually generalizes the Bethe-Peierls approach and the Baxter recurrence method.

Let us study the ferromagnetic Ising model with uniform couplings, , at zero magnetic field, i.e., , on a uncorrelated hypergraph which consists of edges and finite loops. In a random network, effective fields are also random variable and fluctuate from vertex to vertex. For a given , we introduce the distribution function of fields ,

 Ψℓ(h)=1AN∑i=1∑{Mℓ(i)}δ(h−hMℓ(i)). (41)

Here, we sum over vertices and attached motifs of size . is the normalization constant. We assume that in the thermodynamic limit, , the average over vertices in the network is equivalent to the average over the statistical network ensemble, Eq. (12). Using Eq. (36), we obtain an equation for the distribution function ,

 Ψℓ(h)=∫δ(h−Ttanh−1[⟨S⟩Mℓ])ℓ−1∏j=1Φℓ(Ht(j))dHt(j). (42)

Here is the function defined by Eq. (37). is a total field, Eq. (39), acting on vertex in the motif , see Fig. 3. The integration is over fields with the distribution function . Using Eq. (39), we can find this function,

 Φℓ(H)=∑Q2,Q3,...P(Q2,Q3,...)Qℓ⟨Qℓ⟩× ∫δ(H−∞∑m(≠ℓ)Qm∑α=1hMm(α)−Qℓ−1∑α=1hMℓ(α))× ∞∏m(≠ℓ)(Qm∏α=1Ψm(hMm(α))dhMm(α))Qℓ−1∏α=1Ψℓ(hMℓ(α))dhMℓ(α).

Here is given by Eq. (14). This is the probability that a destination vertex has edges, triangles, and so on. In turn, incoming messages to the destination spin from attached loops of size are numerated by the index , . Only for incoming messages from loops of size do we have , because we should not account for the motif . The integration is over incoming messages to all vertices in the motif (hyperedge) except the destination vertex , see Fig. 3. Note that we have an additional factor , because there is the probability in the configuration model that if we arrive at a destination vertex along a hyperedge , then this vertex has outgoing hyperedges . Equations  (42) and  (LABEL:phi) represent a set of self-consistent equations for the functions , . Using Eq. (32), one then can find a distribution function of spontaneous magnetic moments in the network. Equations  (42) and  (LABEL:phi) are also valid for networks with cliques.

It should be noted that our approach is similar to the generalized belief-propagation algorithm proposed in Ref. Yedidia:yfw01 (). It was shown that the belief-propagation algorithm is equivalent to the cluster variation method (CVM) introduced by R. Kikuchi (see Ref. Pelizzola05 () for more details). Note that in the model of complex networks in the present work, two motifs can only overlap each other by a single node. They have no joint links. Motifs here play the role of “hyperedges” in contrast to the approach in Ref. Yedidia:yfw01 () where clusters are considered as “super-nodes”. Finally, we assumed that networks under consideration have a hypertree-like structure in the infinite size limit.

## Iv Critical temperature and critical behavior

In the paramagnetic phase at zero magnetic field, there is no spontaneous magnetization and the effective fields are equal to zero, i.e., . Therefore, we obtain

 Ψℓ(h)=δ(h). (44)

for any . This is the only solution of Eqs. (42) and  (LABEL:phi) in the paramagnetic phase. Below a critical temperature , this solution becomes unstable and a non-trivial solution for the function emerges. In this case, a mean value of the effective field ,

 ⟨hMℓ⟩T=∫hΨℓ(h)dh, (45)

becomes non-zero. Here denotes the thermodynamic average. In the case of a continuous phase transition, tends to zero if the temperature tends to from below. Let use this fact. From Eq. (36) we obtain

 ∫tanh(βh)Ψℓ(h)dh=∫⟨S⟩Mℓℓ−1∏j=1Φℓ(Ht(j))dHt(j). (46)

We expand the functions and on the left and right hand sides of this equation in small and , respectively:

 tanh(βh)=βh+O(h3), ⟨Si⟩Mℓ=ℓ−1∑j=1χℓ(ij)Ht(j)+O(H3t(j)), (47)

where we define

 χℓ(ij)≡∂⟨Si⟩Mℓ∂Ht(j)∣∣Ht(j)=0=β⟨SiSj⟩Mℓ∣∣Ht(j)=0. (48)

is a non-local susceptibility in a spin cluster at zero magnetic field. Note that . Assuming that at the first moment of the function , i.e., , is much larger than higher moments, i.e., , in the leading order we obtain a linear equation,

 ⟨hMℓ⟩T=Fℓ(T)⟨Qℓ⟩[⟨Qℓ(Qℓ−1)⟩⟨hMℓ⟩T +∞∑m(≠ℓ)⟨QℓQm⟩⟨hMm⟩T], (49)

where defines an average over the distribution function given by Eq. (14). The function is defined as follows,

 Fℓ(T)≡Tℓ−1∑j=1χℓ(ij)=Tχℓ−1, (50)

is the total zero-field magnetic susceptibility of the Ising model on a ring of size . Simple calculations give

 F2(T)=t,ℓ=2, Fℓ(T)=2t(1−tℓ−1)(1−t)(1+tℓ),ℓ≥3, (51)

where . In the paramagnetic phase, the set of linear equations (49) for parameters has only a trivial solution, . A non-trivial solution appears at a temperature at which

 detˆM=0, (52)

where the matrix is defined as follows:

 Mmn=⟨QmQn⟩,m≠n, Mmm=⟨Qm(Qm−1)⟩−⟨Qm⟩Fm(T), (53)

for . Equation (52) determines the critical point of the continuous phase transition. Below , spontaneous effective fields appear, i.e., . Therefore, there is a non-zero spontaneous magnetic moment. If all motifs are loops of equal length , then the critical temperature is determined by the equation

 BℓFℓ(T)=1 (54)

where is the average branching coefficient for these hyperedges (-loops in this network),

 Bℓ=⟨Qℓ(Qℓ−1)⟩⟨Qℓ⟩. (55)

In particular, if there are only edges, i.e., , then Eq. (54) gives

 Tc=2J/ln(B2+1B2−1). (56)

This result was found for uncorrelated random complex networks with arbitrary degree distributions dgm02 (); lvvz02 (); Dorogovtsev:dgm08 (). In the limit , we find that

 Tc(ℓ=2)/J≅Bℓ+o(1), (57) Tc(ℓ≥3)/J≅2Bℓ+1+o(1), (58)

Figure 4 shows the dependence of the critical temperature on the mean number of nearest neighbors in the networks with the Poissonian distribution of -loop motifs. In this case the average branching coefficient is . Notice also that at we have according to Eqs. (16) and (17). Comparison between Eqs. (57) and (58) shows that at a given mean degree for the Poisson distribution of hyperdegrees the critical temperature is higher than only by 1. This shift is the only effect of finite loops.

If there are only two motifs, for example, loops of size and , then Eq. (52) leads to the equation

 [Bℓ−1Fℓ(T)][Bℓ′−1Fℓ′(T)]=⟨QℓQℓ′⟩2⟨Qℓ⟩⟨Qℓ′⟩. (59)

From Eqs. (51) and (54), it follows that if an uncorrelated random hypergraph has a divergent second moment for any , then becomes infinite in the infinite size limit. This result was obtained for ordinary uncorrelated random complex network in Refs. dgm02 (); lvvz02 (); Dorogovtsev:dgm08 ().

Another example of complex motifs are -cliques. In this kind of complex networks the clustering coefficient is larger than in networks with loops of size (compare between Eqs. (19) and (20)). Calculating the susceptibility of a cluster of spins on this subgraph, we can find the function in Eq. (50) at ,

 Fℓ(T)=(ℓ−1)[1−4Iℓ−2(T)Iℓ(T)]. (60)

Here we introduced the function

 Iℓ(T)=ℓ∑n=0Cℓnexp[12βJ(2n−ℓ)2], (61)

where is the binomial coefficient. If a network only consists of uncorrelated -cliques then the critical temperature of the ferromagnetic Ising model is given by Eq. (54) with the function Eq. (60). In the limit we obtain an asymptotic result,

 Tc≈(ℓ−1)Bℓ+(ℓ−2)(1+o(1)). (62)

Let us consider the Poisson distribution. In this case we have and the mean degree is . Critical temperatures of networks with cliques and loops are compared in Fig. 4. Comparison between Eqs. (62), (57) and (58) shows that at a given mean degree the critical temperature in Eq. (62) is larger than the critical temperature Eq. (57) only by . This shift is the only effect of clustering. It gives a relative correction of order , because . This result shows that the internal structure of motifs may change the critical temperature of the Ising model although the percolation threshold may be the same (see below).

Thus in the considered random networks the average branching coefficient is crucially important. If , clustering and finite loops lead to a relatively small shift of the critical temperature in comparison to networks without loops but with the same average branching coefficient. However, in networks with a small branching coefficient of motifs, influence of clustering and finite loops is stronger. Indeed, at the critical temperature, Eq. (56), exceeds zero, , if the average branching coefficient (for the Poisson distribution this corresponds to the mean degree ). This is due to the fact that only in this case there is a giant connected component [see Eq. (67)]. In networks with finite loops of size , the percolation threshold from Eq. (67) is . Therefore, the phase transition appears at the average branching coefficient which is much smaller than 1 if (for the Poisson distribution this corresponds to the mean degree ).

Does clustering influence the critical behavior in the network with motifs? In order to answer this question we consider networks that consist of cliques of a given size with the mean number of cliques attached to a node. It is convenient to assume that the distribution function of cliques has the following asymptotic behavior: . In Eq. (LABEL:phi), we use the ”effective medium” approximation introduced in Ref. dgm02 ():

 Qℓ−1∑α=1hMℓ(α)≈(Qℓ−1)⟨hMℓ⟩T (63)

where is the mean effective field acting on a node from an attached -clique. As was shown in Ref. dgm02 (), this approximation takes into account the most dangerous highly connected nodes and gives exact critical behavior for random uncorrelated tree-like networks (see also Refs. lvvz02 (); Dorogovtsev:dgm08 ()). As a result, we obtain the distribution function of the total field acting on a node from attached -cliques as a function of ,

 Φℓ(H)≈∑QℓPℓ(Qℓ)Qℓ⟨Qℓ⟩δ(H−(Qℓ−1)⟨hMℓ⟩T). (64)

Substituting Eqs. (42) and (LABEL:phi) into Eq. (45), we obtain a self consistent equation for the mean effective field ,

 ⟨hMℓ⟩T=G(⟨hMℓ⟩T), (65)

where is the right hand side of Eq. (45). Critical behavior of the model is determined by analytical behavior of the function at small dgm02 (). Analysis of analytical properties of the function at zero magnetic field shows that if the forth moment of degree distribution is finite, i.e., (scale-free networks with the degree exponent ), then where and are certain functions of temperature . Solving Eq. (65), we find that the spontaneous magnetization and the mean effective field behave as below with the standard mean-field exponent . If diverges but the second moment is finite (scale-free networks with ) then . In this case, the critical exponent becomes dependent on the asymptotic behavior of the distribution function , namely, . Finally, if diverges but the mean degree is finite (scale-free networks with ), then in the thermodynamic limit the critical temperature tends to infinity, i.e., at any finite temperature the Ising model is in the ordered phase. This is the critical behavior that was found for the ferromagnetic Ising model on random uncorrelated complex networks dgm02 (); lvvz02 (); Dorogovtsev:dgm08 ().

Our approach is explicit for sparse random networks that have local hypertree-like structure. Two motifs can only overlap each other by a single node. If two clusters have common edges, one can combine them and introduce a new motif. The internal structure of motifs does not influence the critical behavior of the ferromagnetic Ising model. However this structure may be essential for models with frustrations (spin glasses). Relaxing the hypertree-like assumption may change crucially the critical behavior of the Ising model. Furthermore, as we have showed above, asymptotic behavior of distribution function of motifs is significant for the critical behavior of the Ising model and other models of statistical physics on complex networks Dorogovtsev:dgm08 (). In our analysis of critical behavior we also assumed that there are no correlations in the distribution of motifs over a network. The role of degree-degree correlations in critical behavior for percolation on complex networks was demonstrated in Ref. gdm08 (). Weak clustering in the sense of Sec. II does not influence the critical behavior of the Ising model. The influence of strong clustering demands further investigation.

## V Percolation threshold

Equation (52) permits us to find the percolation threshold in these loopy networks. Below the percolation threshold, a network consists of finite clusters and there is no giant component. In this case, there is no phase transition in the Ising model. This spin system is in a paramagnetic state at any . At a point of the birth of a giant component, there is a giant connected cluster and the critical temperature is . Above the percolation threshold, the critical temperature is non-zero, . Therefore, in a general case for an arbitrary distribution function