Exotic Critical Behavior of Weak Multiplex Percolation
G. J. Baxter, R. A. da Costa, S. N. Dorogovtsev, J. F. F. Mendes
EExotic Critical Behavior of Weak Multiplex Percolation
G. J. Baxter, R. A. da Costa, S. N. Dorogovtsev, and J. F. F. Mendes Department of Physics, University of Aveiro & I3N,Campus Universit´ario de Santiago, 3810-193 Aveiro, Portugal (Dated: June 2, 2020)We describe the critical behavior of weak multiplex percolation, a generalization of percolation tomultiplex or interdependent networks. A node can determine its active or inactive status simply byreferencing neighboring nodes. This is not the case for the more commonly studied generalization ofpercolation to multiplex networks, the mutually connected clusters, which requires an interconnect-ing path within each layer between any two vertices in the giant mutually connected component.We study the emergence of a giant connected component of active nodes under the weak percola-tion rule, finding several non-typical phenomena. In two layers, the giant component emerges witha continuos phase transition, but with quadratic growth above the critical threshold. In three ormore layers, a discontinuous hybrid transition occurs, similar to that found in the giant mutuallyconnected component. In networks with asymptotically powerlaw degree distributions, defined bythe decay exponent γ , the discontinuity vanishes but at γ = 1 . γ = 1 + 1 / ( M −
1) in M layers. I. INTRODUCTION
Complex systems with interdependent sub-systemsmay be modeled as a multiplex (or colored) network,where links of different types (colors) represent connec-tions within different sub-systems, while nodes havingmore than one type of connection encompass interdepen-dencies between subsystems [1]. Equivalently, one mayuse a multiplex network, with a layer for each sub-systemand links between nodes in different layers representinginterdependencies [1, 2]. To study the resilience of suchsystems, one typically considers a generalization of perco-lation. The concept of connected cluster in a single-layernetwork generalizes to mutually connected clusters, de-fined as a set of nodes each pair of which is connected byat least one path in all of the layers in which they partic-ipate [2–4]. The interdependency between layers (colors)leads to increased fragility of the system, and under ran-dom damage, the giant mutually connected componentcollapses discontinuously [2] in a hybrid phase transitionof the k -core type [3, 5]. The collapse occurs due to longrange cascading failures in the system [3].This percolation process applies a global condition toidentify surviving nodes: a path of every color must ex-ist between every pair of nodes in a mutually connectedcluster in order for the members of the cluster to sur-vive. One may identify the mutually connected clustersby a global pruning process, iteratively removing clustersin each layer that do not have a counterpart cluster ineach other layer, until a stable equilibrium is reached.Alternatively, one may identify the connected clusters ineach layer, and remove any non-overlapping parts, thenrepeating the process until no more nodes are pruned.This percolation process has been extensively studied,with works considering effects of partial interdependence[6] overlapping edges [7–9], multiple dependencies [10]correlations [11] among many others [12–15]. As in manynetwork processes, heavy-tailed degree distributions havea strong effect on the phase transition, with the point a b FIG. 1. Illustration of the difference between weak percola-tion and mutually connected clusters. In (a) all nodes haveconnections of both solid and dashed types, and belong to thesame weak percolating cluster. This graph contains no mutu-ally connected clusters. (b) In order for a set of nodes to forma mutually connected cluster, there must be a path of bothkinds between every pair of vertices. The four vertices showsform a mutually connected cluster and a weak percolatingcluster. at which the giant mutually connected cluster emerges p c → − p is the fraction of nodes initiallydamaged) as the powerlaw exponent γ →
2. Approach-ing the same point, the size of the discontinuity decaysrapidly as S c ∼ − / ( γ − / ( γ −
2) [3].An alternative definition for percolation on multiplexnetworks was proposed in Ref. [16], in which the survivalof a node is established by a strictly local rule: if it has atleast one connection to another surviving node in everylayer in which it participates. In this rule the survival ofan agent depends only on its immediate neighborhood.Due to the less restrictive definition, we refer to thispercolation process as weak multiplex percolation. De-spite its purely local (and hence more realistic) process,it nevertheless may undergo the same discontinuous hy-brid transition as found in the ‘strong’ rule describedabove. The phase diagram in networks with rapidly de-caying degree distributions was delineated in Ref. [16].This process was further explored in Ref. [17], and therelationship with the stronger rule elaborated in Ref.[4].In ( M ≥ a r X i v : . [ c ond - m a t . d i s - nn ] M a y (1 − − ... − γ , the discontinuous transitiondisappears at γ = 1 . / ( M −
1) in M layres) in contrast to the limit γ = 2 found for example in ordinary percolation [19, 20]and the mutually connected cluster [3].The remainder of this paper is organized in the follow-ing way. In Section II we define the problem and givethe general self-consistency equations which allow for acomplete solution. We consider the continuous transitionwhich occurs in two layer networks in Section III. In Sec-tion IV we consider the discontinuous hybrid transitionwhich occurs in three or more layers. Finally discussionand conclusions are given in Section V. II. PROBLEM AND GENERAL ANALYSIS ba ∞ ∞ ≥ ≥ S = ∑ = ∞ ∞ ≥ ∑ ∞ ≥ = ∞ ∞ ≥ ∑ ∞ ≥ FIG. 2. Diagrammatic representation of self-consistencyequations for M = 2 layers, a and b . (a) In a tree-like net-work, a node belongs to the giant weak-percolation cluster(giant component) if it has at least one connection via anedge of type a to an infinite subtree satisfying the property Z a (represented by a solid edge leading to an infinity sym-bol), and one of type b leading to a satisfying Z b (dashededge leading to an infinity symbol), see Eq. (4). (b) Therecursive relations obeyed by the probabilities Z a (left) and Z b (right), see Eq. (3). For Z a an edge in layer a leads toa node with at least one edge of type b satisfying Z b , andmay also have connections satisfying Z a , and similarly for Z b , exchanging the labels. Let us consider a large sparse random multiplex net-work, consisting of N nodes connected in one or moreof M layers (each having its own unique type of edge).Note that a node does not necessarily participate in alllayers. The analysis which follows is not impeded by the presence of degree correlations between layers, so we con-sider a generalised configuration model network definedby its joint degree distribution P ( q , q , ..., q M ). A nodemay be considered active if it retains at least one connec-tion to other active nodes in each of the layers in whichit participates. A weak percolating cluster is then a setof such active nodes which are connected to each other(each member is connected to at least one other memberin at least one layer).In the large size limit N → ∞ we can use the locallytree-like property of the network to write self-consistencyequations for the probability that a randomly selectednode is active. This is equal to the the relative size ofthe giant weak-percolation cluster (we will from now onuse “giant component” as a shorthand for this cluster) insuch networks. The size of giant component in an infinitesparse network is given by S = (cid:88) q ,q ,...,q M P ( q , q , ..., q M ) M (cid:89) α =1 [1 − (1 − Z α ) q α ] (1)where the probabilities Z , Z , ..., Z M are given by Z α = (cid:88) q ,q ,...,q M q α P ( q , q , ..., q M ) (cid:104) q α (cid:105) (cid:89) β (cid:54) = α [1 − (1 − Z β ) q β ](2)for α = 1 , , ..., M . They represent the probability that,upon following an edge of type α , we encounter a vertexwith at least one edge of type β satisfying Z β for all layers β (cid:54) = α . These equations are illustrated diagrammaticallyfor two layers in Fig. 2.One may then obtain the size of the giant componentby solving, Eqs. (2) and substituting the solution intoEq. (1). In a two layer network with rapidly decayingdegree distribution (such as an Erd˝os-R´enyi network), asconnectivity increases, the giant component appears con-tinuously with a second-order phase transition, which dif-fers from the standard percolation transition as the giantcomponent grows quadratically above the critical thresh-old. For three or more layers, on the other hand, one findsthat the giant component appears with a discontinuoushybrid transition, similar to that seen in k -core percola-tion. The size of the giant component S jumps from zeroto a finite value at a critical threshold, and then growsas the square root of the distance above the threshold. III. TWO LAYERS - CONTINUOUSTRANSITION
Let us first consider the case of two layers. In this casethe giant percolating cluster emerges with a continuoustransition, but with different characteristics than the or-dinary percolation transition.Equation (2) becomes Z a = (cid:88) q a ,q b q a P ( q a , q b ) (cid:104) q a (cid:105) − (1 − Z b ) q b ] ,Z b = (cid:88) q a ,q b q b P ( q a , q b ) (cid:104) q b (cid:105) [1 − (1 − Z a ) q a ] (3)and the expression for the size of the giant weak perco-lation cluster (which corresponds to the (1 − k -core) is S = (cid:88) q a ,q b P ( q a , q b )[1 − (1 − Z a ) q a ][1 − (1 − Z b ) q b ] . (4) A. Rapidly decaying degree distributions
If the moments (cid:104) q a (cid:105) , (cid:104) q b (cid:105) , (cid:104) q a q b (cid:105) , (cid:104) q a q b (cid:105) , and (cid:104) q a q b (cid:105) are finite, which is the case, for example, for Erd˝os-R´enyinetwork layers, we may expand Eqs. (3) for small Z a and Z b , finding Z a ∼ = 1 (cid:104) q a (cid:105) (cid:2) (cid:104) q a q b (cid:105) Z b − (cid:104) q a q b ( q b − (cid:105) Z b (cid:3) ,Z b ∼ = 1 (cid:104) q b (cid:105) (cid:2) (cid:104) q a q b (cid:105) Z a − (cid:104) q a q b ( q a − (cid:105) Z a (cid:3) (5)which indicates that the continuous transition occurswhen (cid:104) q a (cid:105)(cid:104) q b (cid:105) = (cid:104) q a q b (cid:105) . (6)Near the transition point, we can write S ∼ = (cid:104) q a q b (cid:105) Z a Z b , (7)and using Eq. (5) we find the size of the giant clusternear the critical point: S ∼ = 4 (cid:104) q a (cid:105)(cid:104) q b (cid:105) (cid:0) (cid:104) q a q b (cid:105) − (cid:104) q a (cid:105)(cid:104) q b (cid:105) (cid:1) (cid:104) q a q b (cid:105) Q a Q b . (8)where for compactness we have defined Q a ≡ (cid:2) (cid:104) q a q b ( q a − (cid:105)(cid:104) q b (cid:105) + (cid:104) q a q b ( q b − (cid:105)(cid:104) q a q b (cid:105) (cid:3) , (9) Q b ≡ (cid:2) (cid:104) q a q b ( q b − (cid:105)(cid:104) q a (cid:105) + (cid:104) q a q b ( q a − (cid:105)(cid:104) q a q b (cid:105) (cid:3) . (10)The giant weakly percolating component grows as thesquare of the distance from the critical point, i.e. β = 2,as opposed to the usual percolation transition which has β = 1. To illustrate this, let us consider the simplifiedcase in which a node’s degrees in each layer are indepen-dent, P ( q a , q b ) = P a ( q a ) P b ( q b ). Then (cid:104) q a q b (cid:105) = (cid:104) q a (cid:105)(cid:104) q b (cid:105) ,and the condition of Eq. (6) may be reduced to (assum-ing still that (cid:104) q a (cid:105) , (cid:104) q b (cid:105) < ∞ ) (cid:104) q a (cid:105)(cid:104) q b (cid:105) = 1 . (11) Then near the transition, S = Z a Z b ∼ = 4 (cid:104) q a (cid:105) ( (cid:104) q a (cid:105)(cid:104) q b (cid:105) − [( (cid:104) q b (cid:105)(cid:104) q a (cid:105) − (cid:104) q a (cid:105) + (cid:104) q a (cid:105) − (cid:104) q a (cid:105) ] , (12)where we have used that (cid:104) q b (cid:105) = 1 / (cid:104) q a (cid:105) at the threshold.Further, if the network is symmetric ( (cid:104) q a (cid:105) = (cid:104) q b (cid:105) ≡ (cid:104) q (cid:105) ),then the giant component exists for (cid:104) q (cid:105) >
1. Assuming,moreover, (cid:104) q (cid:105) is finite, we arrive at the following relativesize of the giant component near the transition S ∼ = 4 ( (cid:104) q (cid:105) − ( (cid:104) q (cid:105) − . (13)Here ( (cid:104) q (cid:105) −
1) plays the role of a control parameter, andwe see immediately that the growth is quadratic. In thesymmetric Erd˝os–R´enyi situation, S coincides with thesquare of the relative size of the giant connected compo-nent in an individual layer. B. Heavy-tailed degree distributions
For strongly heterogeneous degree distributions thecondition that the leading moments are finite may notbe met. In this case we may use generating func-tions to study the asymptotics of the solutions. When P ( q a , q b ) = P a ( q a ) P b ( q b ), we may rewrite Eqs. (3) usinggenerating functions (see Appendix B) as Z a = 1 − G b (1 − Z b ) ,Z b = 1 − G a (1 − Z a ) . (14)while the size of the giant percolating cluster is simply S = Z a Z b . (15)For concreteness, we will consider uncorrelated power-law tailed degree distributions of the form P ( q ) = Aq − γ (16)for each layer. Note that, as we will see, the exponentfound above for rapidly decaying degree distributions, β = 2, applies for γ > β = 1 applies only for γ > P ( q a ) = P ( q b ) ≡ P ( q ). When 2 < γ <
3, using the asymptoticbehaviour of the generating function, see Appendix B, Z ∼ = (cid:104) q (cid:105) Z − A Γ(1 − γ ) Z γ − , (17)so the size of the giant component equals S = Z ∼ = (cid:104) (cid:104) q (cid:105) − A Γ(1 − γ ) (cid:105) / ( γ − . (18)In Fig. 3 we compare this theoretical calculation of S with simulated networks containing N = 10 and 10 . . . . . . . . . h q i . . . . . . S Poisson γ =2 . γ =2 . γ =2 . γ =2 . γ =3 γ =4 FIG. 3. Relative size of the largest weak-percolating cluster S as a function of mean degree (cid:104) q (cid:105) in symmetric two layermultiplex networks, M = 2, with each layer having a pow-erlaw degree distribution P ( q ) ∼ Aq − γ generated using thestatic model, containing N = 10 and N = 10 nodes (resultsare virtually identical), see Appendix C. Results for Poissondegree distributions are shown for comparison. Black solidcurves in each case are numerical solutions of Eqs. (1) and(2) . nodes, above the threshold ( c ≡ (cid:104) q (cid:105) = 1), showing per-fect agreement. Each layer is an independently generatedconfiguration model network, generated according to thestatic model degree distribution, which is asymptoticallypowerlaw [21, 22]. See Appendix C for more details. Thisfigure illustrates the very slow growth for values of γ closeto 2, and the approach to the quadratic growth at γ ≥ < γ <
2, the mean degree diverges and onemust proceed with caution, as the conditions requiredfor the self consistency equations to be exact may nothold. Nevertheless we find, after extensive comparisonagainst numerical simulations, that our equations giveaccurate results. In this region we can no longer use themean degree as a control parameter. Instead, we mayconsider applying random damage to the network. Ifedges are retained with probability p and removed withprobability 1 − p , the tail of the degree distribution retainsthe same powerlaw exponent γ , with a reduced coefficient A p = A p γ − , as we show in Appendix A. In this casewe find Z ∼ = − A p γ − Γ(1 − γ ) Z γ − , (19)so, the size of the giant component is S = Z ∼ = [ − A Γ(1 − γ )] / (2 − γ ) p γ − / (2 − γ ) . (20)For site percolation, we have Z ∼ = p [ − A Γ(1 − γ )] Z γ − , (21) and so S = pZ ∼ = [ − A Γ(1 − γ )] / (2 − γ ) p (4 − γ ) / (2 − γ ) . (22)Thus, for both edge and vertex removal, the giant com-ponent appears immediately from p c = 0.We now consider the cases in which the exponents foreach layer are different, P a ( q a ) = Aq − γ a a and P b ( q b ) = Aq − γ b b . In general, the critical behaviour is determinedby the smaller of the two degree distribution exponents.Without loss of generality, let us assume γ a > γ b . Resultsfor the opposite case can be obtained by simply exchang-ing the subscripts a and b . If both exponents are greaterthan 3, we have the behavior described in the previousSection.We first consider the case γ a >
3, 2 < γ b <
3. Then γ a − >
2, so the leading terms in the expansion of Z b are linear and quadratic. We may neglect the quadraticterm, so we have Z a ∼ = (cid:104) q b (cid:105) Z b − A b Γ(1 − γ b ) Z γ b − b ,Z b ∼ = (cid:104) q a (cid:105) Z a . (23)The solution is Z a ∼ = (cid:20) (cid:104) q a (cid:105)(cid:104) q b (cid:105) − A b Γ(1 − γ b ) (cid:21) / ( γ b − (cid:104) q a (cid:105) − ( γ b − / ( γ b − , (24)so S = (cid:20) (cid:104) q a (cid:105)(cid:104) q b (cid:105) − A b Γ(1 − γ b ) (cid:21) / ( γ b − (cid:104) q a (cid:105) − γ b / ( γ b − . (25)Note that, somewhat counterintuitively, the fatter-taileddegree distribution (smaller exponent) determines the be-havior, in contrast to what would be the case for moretraditional percolation problems.If both exponents are less than three, 2 < γ a , γ b < Z a ∼ = (cid:104) q b (cid:105) Z b − A b Γ(1 − γ b ) Z γ b − b ,Z b ∼ = (cid:104) q a (cid:105) Z a − A a Γ(1 − γ a ) Z γ a − a . (26)Substituting the second line into the first,( (cid:104) q a (cid:105)(cid:104) q b (cid:105) − Z a ∼ = (cid:104) q b (cid:105) A a Γ(1 − γ a ) Z γ a − a + (cid:104) q a (cid:105) γ b − A b Γ(1 − γ b ) Z γ b − a . (27)If γ a > γ b , then the first term on the right-hand side ofEq. (27) should be neglected, and we obtain Z a ∼ = (cid:20) (cid:104) q a (cid:105)(cid:104) q b (cid:105) − (cid:104) q a (cid:105) γ b − A b Γ(1 − γ b ) (cid:21) / ( γ b − , (28)thus S ∼ = (cid:104) q a (cid:105) − γ b / ( γ b − (cid:20) (cid:104) q a (cid:105)(cid:104) q b (cid:105) − A b Γ(1 − γ b ) (cid:21) / ( γ b − . (29)Now let us consider 1 < γ b < γ a >
2. We obtain Z a ∼ = − A b Γ(1 − γ b ) Z γ − b ,Z b ∼ = (cid:104) q a (cid:105) Z a . (30)Note that, given the first equation, we can neglect higherorder terms in the second equation, so it applies both for2 < γ a < γ a >
3, and we treat both these casestogether.As before, we apply random damage, using the re-tention probability p as the control parameter. Thedegree distribution for layer b is then asymptotically A b p γ b − q − γ b , where A b, ≡ A b ( p =1). The solution toEqs. (30) is Z a ∼ =[ − A b Γ(1 − γ b )] / ( γ b − (cid:104) q a (cid:105) − ( γ b − / (2 − γ b ) ∝ p − ( γ b − / (2 − γ b ) , (31)so that S = (cid:104) q a (cid:105) Z a ∼ = [ − A b, Γ(1 − γ b )] / ( γ b − (cid:104) q a (cid:105) − γ b / (2 − γ b ) p − ( γ b − / (2 − γ b ) . (32)Note that this expression contains only (cid:104) q a (cid:105) of layer a .Again, the behaviour is determined by the fatter taileddegree distribution.Finally, when both exponents are small, 1 < γ a , γ b < Z a and Z b have the form: Z a ∼ = − A b Γ(1 − γ b ) Z γ b − b ,Z b ∼ = − A a Γ(1 − γ a ) Z γ a − a . (33)So S = Z a Z b ∼ = (cid:104) − A b Γ(1 − γ b ) (cid:105) γ a / [1 − ( γ b − γ a − × (cid:104) − A a Γ(1 − γ a ) (cid:105) γ b / [1 − ( γ b − γ a − ∼ = (cid:104) − A a Γ(1 − γ a ) (cid:105) γ b / [1 − ( γ b − γ a − × (cid:104) − A b Γ(1 − γ b ) (cid:105) γ a / [1 − ( γ b − γ a − × p ( γ a − γ b / [1 − ( γ b − γ a − a p ( γ b − γ a / [1 − ( γ b − γ a − b . (34)Here we assumed that the fractions of retained edges inlayers a and b are p a and p b respectively.Finally, for completeness, we may make the same cal-culation for the case of vertex removal. Vertices survivewith probability p . A factor of p is added to Eq. (33)[compare Eq. (21)]. Solving for Z a and Z b then substi-tuting into S = pZ a Z b we find S ∼ = (cid:104) − A a Γ(1 − γ a ) (cid:105) γ b / [1 − ( γ b − γ a − × (cid:104) − A b Γ(1 − γ b ) (cid:105) γ a / [1 − ( γ b − γ a − × p ( γ a γ b ) / [1 − ( γ b − γ a − . (35) IV. HIGHER NUMBER OF LAYERS
For more than two layers, the giant weakly percolatingcomponent typically appears with a discontinuous hybridtransition [4, 16], of the same type observed in the mu-tually connected cluster [3] and in k -core percolation [5].For M ≥ Z α = M (cid:89) β (cid:54) = α [1 − G β (1 − Z β )] (36)for α = 1 , , ..., M , and S = M (cid:89) β =1 [1 − G β (1 − Z β )] = (cid:16) M (cid:89) β =1 Z β (cid:17) / ( M − . (37)In the symmetric case, in which every layer is arandom network with the same degree distribution, P ( q a , q b , q c , ... ) = P ( q a ) P ( q b ) P ( q c ) ... , we have Z = [1 − G (1 − Z )] M − (38)and S = Z M/ ( M − . (39)In this situation, for Poisson degree distributions with c ≡ (cid:104) q (cid:105) , Eq. (38) leads to the equation: Z = [1 − e − cZ ] M − . (40)This equation is practically identical to the one obtainedin Ref. [23] for the relative size S ∗ of the giant mutuallyconnected component in M -layer multiplex Erd˝os-R´enyinetworks: S ∗ = [1 − e − cS ∗ ] M . (41)Comparing Eqs. (39) and (40) with Eq. (41), we obtainthe following relation between these two problems: S ( c, M ) = S ∗ M/ ( M − ( c, M − ,c c ( M ) = c ∗ c ( M − , (42)where M ≥
3, and c c ( M ) and c ∗ c ( M ) are the critical valueof the average degree for the giant mutually connectedcomponent and weak percolation, respectively, for M -layer multiplex Erd˝os-R´enyi networks. Thus the weakpercolation problem on a M -layer multiplex Erd˝os-R´enyinetwork is equivalent to the problem of giant mutuallyconnected component in the corresponding M − M , the asymptotics of these quantities arethe same in both problems: c c ∼ = ln M + ln ln M + 1 + ln ln M ln M ,S c ∼ = 1 − M + ln ln M ln M . (43) h q i . . . . . . S Poisson γ =4 γ =3 γ =2 . γ =2 . FIG. 4. Relative size of the giant component S for threeidentically powerlaw distributed layers ( M = 3) as a functionof mean degree (cid:104) q (cid:105) for various values of γ greater than 2.For details of the simulations see Appendix C. Black curvesshow analytic results (from numerical solution of Eqs. (1) and(2)), symbols show measurements averaged over 100 syntheticnetworks of N = 10 nodes (circles) and N = 10 nodes(squares). A. Effect of heavy-tailed degree distributions
The discontinuous hybrid transition always maintainsthe same square-root scaling above the transition, how-ever the size of the discontinuity and the location of thecritical point for such transitions may be strongly affectedby the degree distribution [3]. For orientation, we againbegin with the symmetric case, P ( q ) = Aq − γ . As can beseen in Fig. 4, the location and size of the discontinuitydepends strongly on the powerlaw exponent γ . Simula-tion results were obtained using independently generatedconfiguration model networks following the static modeldegree distribution for each layer, just as in Fig. 3, seeAppendix C. As γ approaches 2, the mean degree di-verges and finite size effects become particularly promi-nent, as evidenced by the divergence between theoreticaland numerical results, and between numerical results fornetworks of different sizes at γ = 2 . γ = 2, we must againproceed with caution. As shown in Figs. 5 and 6, weverify all results against numerical simulations, and findthat the numerical measurements converge to the analyt-ical results as the system size increases. We again intro-duce random damage, and use the undamaged fractionof edges or vertices p as a control parameter. The hybridtransition disappears at γ = 1 + 1 / ( M − γ = 2.To find the size and location of the jump, we modify . . . . . . p . . . . . . S γ =1 . γ =1 . γ =2 . γ =2 . γ =3 . γ =4 . γ =6 . γ =10 . FIG. 5. Relative size of the giant component S for threeidentically powerlaw distributed layers ( P ( q ) = Aq − γ , q ≥ p for various values of γ greater than 1 .
5. Black curves show analytic results, symbolsshow measurements for synthetic networks of N = 10 nodes(circles, 100 realisations) and N = 10 nodes (squares, onerealisation), see Appendix C. .
00 0 .
05 0 .
10 0 .
15 0 .
20 0 . p . . . . . . S γ =1 . γ =1 . γ =1 . γ =1 . γ =1 . FIG. 6. Relative size of the giant component S for threeidentically powerlaw distributed layers ( P ( q ) = Aq − γ , q ≥
4) as a function of undamaged fraction p for various valuesof γ less than or equal to 1 .
5. Black curves show analyticresults, with the finite threshold for γ = 1 . N = 10 nodes (circles) and N = 10 nodes (squares), see Appendix C. Eq. (37) for symmetrical layers, obtaining S = [1 − G (1 − Z )] M (44)where Z obeys Eq. (38), so that S = Z M/ ( M − . (45)We look for an expansion of G (1 − Z ) for small Z ,and keep the two leading orders, see Eqs. (B6)-(B9) inAppendix B. When γ < γ − < Z is then(for Z (cid:28) Z ∼ = (cid:8) − A Γ(1 − γ ) Z γ − − BZ (cid:9) M − ≡ Ψ( Z ) (46)where the coefficient B , given by Eq. (B10), depends onthe specific form of the degree distribution.A hybrid transition occurs when the line Z is tangentto ψ ( Z ), which occurs when (cid:18) Ψ Z (cid:19) (cid:48) = 1 Z (cid:20) Ψ (cid:48) − Ψ Z (cid:21) = 0 . (47)Assuming that Z (cid:54) = 0 this gives us the value of Z abovethe discontinuity: Z c = (cid:26) − A Γ(1 − γ ) [( M − γ − − B ( M − (cid:27) / (2 − γ ) . (48)We see that Z c (and hence S c ) tends to zero at γ =1 + 1 / ( M − γ = 1 + 1 / ( M −
1) + δ , we cansee that Z c near δ = 0 behaves as Z c ∼ = (cid:20) − A ( M − B ( M −
2) Γ (cid:18) − M − (cid:19) δ (cid:21) ( M − / ( M − . (49)To account for the damage applied to the network, asvertex removal, we use the original degree distribution,but modify Eqs. (44) and (38) as follows S = p [1 − G (1 − Z )] M (50)and Z = p [1 − G (1 − Z )] M − , (51)so that S = p (cid:20) Zp (cid:21) M/ ( M − . (52)The self-consistency equation for Z is now (for Z (cid:28) Z ∼ = p (cid:8) − A Γ(1 − γ ) Z γ − − BZ (cid:9) M − ≡ Ψ( Z ) . (53)Following the same procedure, applying the condition forthe hybrid transition point Eq. (47) gives us again Eq.(49).Substituting Eq. (49) back into Eq. (53) gives p c ∼ = (cid:20) − A Γ (cid:18) − M − (cid:19)(cid:21) − ( M − × (cid:20) ( M − B ( M − δ (cid:21) − δ ( M − / ( M − . (54) Taking the limit δ → p c = (cid:20) − A Γ (cid:18) − M − (cid:19)(cid:21) − ( M − (55)which depends on the degree distribution only throughthe amplitude A . We see in Fig. 5 that, although thesize of the discontinuity tends to zero, the critical pointremains finite. For example, for M = 3, for the distribu-tion used in Figs. 5 and 6, p c → A Γ( − . ≈ . ... as γ → . + . (56)This point is marked with a black circle in Fig. 6.Finally, using Eq. (52) gives, for site removal: S c ∼ = (cid:20) − A Γ (cid:18) − M − (cid:19)(cid:21) M − / ( M − × (cid:20) ( M − B ( M − δ (cid:21) M/ ( M − . (57)If, instead, we wish to consider edge removal, we usethe amplitude A p = A p γ − as given by Eq. (A3). ThenEq. (49) becomes Z c ∼ = (cid:20) − A ( M − B ( M −
2) Γ (cid:18) − M − (cid:19) δ (cid:21) ( M − / ( M − × p / ( M − c . (58)Combining Eq. (58) with Eq. (46), gives the same criticalpoint, Eq. (54).Finally, using Eq. (45) gives, for edge removal: S c ∼ = (cid:20) ( M − B ( M − δ (cid:21) M/ ( M − . (59)For 1 < γ < / ( M −
1) the transition is continu-ous, and the critical point is always zero, with extremelyslow growth of the giant component as shown in Fig.6. Finite size effects are even more significant close to γ = 1 + 1 / ( M − S as thesize of the network increases.Keeping only the leading order in Eq. (46) we have Z = [ − A p Γ(1 − γ )] M − Z ( M − γ − , (60)so Z = [ − A p Γ(1 − γ )] ( M − / [1 − ( M − γ − ∝ p ( γ − M − / [1 − ( M − γ − , (61)where we recalled that A p = A p γ − . The exponent( γ − M − / [1 − ( M − γ − γ < M − . This gives immediately S = [ − A p Γ(1 − γ )] M/ [1 − ( M − γ − =[ − A Γ(1 − γ )] M/ [1 − ( M − γ − p ( γ − M/ [1 − ( M − γ − . (62)We see that S grows as a power of p , so p c = 0 in thisregion. The exponent diverges as we approach γ = 1 + M − from below, so S grows extremely slowly in thislimit. This is illustrated in Fig. 6. Specifically for M = 3, S ∼ p γ − / [1 − γ − . (63)This gives S ∼ p at γ = 1 . Z = p [ − A Γ(1 − γ )] M − Z ( M − γ − , (64)so Z = [ − A Γ(1 − γ )] ( M − / [1 − ( M − γ − p / [1 − ( M − γ − . (65)Hence S = pZ M/ ( M − =[ − A Γ(1 − γ )] M/ [1 − ( M − γ − × p [2 M − − ( M − ( γ − / { ( M − − ( M − γ − } . (66)The interval in which the hybrid transition is absent,1 < γ < / ( M − M → ∞ . Thus only rather fat-tailed degree distributionscan maintain this singularity.For the non-symmetric multiplex networks, Eqs. (36)and (37), let us consider the case of M = 3, 1 <γ a , γ b γ c <
2. Then at small Z we have Z a ∼ Z γ b − b Z γ c − c ,Z b ∼ Z γ a − a Z γ c − c ,Z c ∼ Z γ a − a Z γ b − b (67)From this system of equations we get Z a ∼ Z ( γ a − γ b − γ a γ b ( γ c − /γ c a . (68)The discontinuity is absent (the transition, i.e., a singu-larity, in this case is at zero—hyper-resilience) if the ex-ponent of the right-hand side of Eq. (68) is smaller than1. This leads to the following condition for the absenceof the discontinuity:2( γ a − γ b − γ c − γ a − γ b −
1) + ( γ b − γ c −
1) + ( γ c − γ a − < . (69) V. CONCLUSIONS
Multi-level networks have received significant attentionin recent years. The structure and resilience of such net-works has generally been studied by generalizing the con-cepte of connected clusters to mutually connected clus-ters, in which there must exist a path connecting everypair of vertices in all the layers in which they participate.A node is active if it belongs to such a mutually con-nected cluster. This yields an exotic percolation phasetransition that is discontinuous yet retains some featuresof a second-order transition. The same type of transitionhas been found in k -core percolation. However, identi-fying the mutually connected clusters requires a globalview of the multiplex network: vertices must belong tothe same connected cluster in each layer. An alterna-tive definition of multiplex percolation was introduced inRef. [16]. Under this definition, a node is active if itmaintains connections to active nodes in each of the lay-ers to which it belongs. Thus the state of a node canbe determined by examining the state of its neighbors.Despite this simplicity, in this paper we have shown thatthis ”weak multiplex percolation” exhibits a complex setof critical phenomena.When the network consists of two layers, we en-counter a continuous second-order transition, as in or-dinary percolation. We have shown, however that forrapidly decaying degree distributions, the giant compo-nent grows quadratically rather than linearly above thecritical threshold. When the degree distributions of thelayers are heavy-tailed, such as powerlaw distributed, wefind that the giant component grows nonlinearly abovethe critical point, with an exponent that depends onthe powerlaw decay exponent of the degree distribution.When this exponent is different in the two layers, thecritical behavior depends on the smaller of the two, ex-cept when both are smaller than 2, in which case thegrowth of the giant component above the critical point isdetermined by both powerlaw exponents.In networks with three or more layers, the giant weakmultiplex-percolation component emerges with the samediscontinuous hybrid transition found in the mutuallyconnected component. This transition consists of a dis-continuity, with a square root singularity above the crit-ical threshold. The weak percolation problem on an M -layer multiplex Erd˝os-R´enyi network is equivalent to theproblem of the giant mutually connected component inthe corresponding ( M − γ = 2 the threshold and discontinuity are still finite.Exponents smaller than 2 are not usually investigated,as the diverging mean degree means the usual locallytree-like assumptions for configuration model networksdo not strictly hold. However, by carefully comparingwith large scale simulations, we show that our equationsgive meaningful and accurate results in this regime. Weshow that the discontinuity and critical point don’t be-come zero until γ = 1+1 / ( M − M layer network.This differs sharply from the mutual connected compo-nent, for which p c = 0 at γ = 2 in two layers [3], andnormal percolation, where p c = 0 at γ = 3. In the rangeof powerlaw exponents γ >
2, for a large number of lay-ers, the weak percolation behavior becomes essentiallythe same as that of the mutually connected component.The weak multiplex percolation process has the advan-tage of being locally decidable and thus corresponds to adifferent type of process than mutually connected com-ponents, being defined by physical bonding to neighborsin all layers vs connectivity to all members of a clusterwithin each layer. While sharing many of the same crit-ical phenomena, the two processes differ significantly innetworks with heavy-tailed degree distributions. Choos-ing the appropriate model for a given system is thereforeimportant for making correct predictions about its re-silience and critical behavior.
ACKNOWLEDGMENTS
This work was developed within the scope of theproject i3N, UIDB/50025/2020 and UIDP/50025/2020,financed by national funds through the FCT/MEC. Thiswork was also supported by National Funds throughFCT, I. P. Project No. IF/00726/2015 and Project No.EXPL/FIS-NAN/1275/2013. R. A. d. C. acknowledgesthe FCT Grants No. SFRH/BPD/123077/2016 and No.CEECIND/04697/2017.
Appendix A: Modification of degree distribution byremoval of edges or vertices
Let p be the fraction of undeleted edges (or vertices).For the original, undamaged, network, p = 1, we have thetail of the degree distribution A k − γ and its first moment (cid:104) k (cid:105) , where the amplitude A p =1 ≡ A . The first momentfor p < (cid:104) k (cid:105) p = p (cid:104) k (cid:105) . (A1)The damage acts on the scale-free degree distributionin the following way. The low-degree part of the dis-tribution increases. In particular, (additional) verticesof degree 0 and especially importantly ones of degree 1emerge. The high-degree asymptotics stays ∝ k − γ , butits amplitude decreases.For establishing the relation between A p and A for p (cid:28)
1, note the following. Under edge removal, for large k , the number of (surviving) vertices with degrees q >k in the damaged network should be equal to the thenumber of vertices with degrees q > k/p in the original network. After integration, this gives A p k − γ ∼ = A p γ − k − γ . (A2)So A p ∼ = A p γ − . (A3)For the case of vertex removal, only a fraction p of ver-tices survive, and these then keep each edge with proba-bility p . This gives an extra factor of p :˜ A p ∼ = A p γ . (A4) Appendix B: Generating functions
The generating function of the degree distribution P ( q )is defined as G ( x ) = (cid:88) q P ( q ) x q . (B1)For a Poisson degree distribution with mean degree (cid:104) q (cid:105) ≡ c , G ( x ) = e − c (1 − x ) , (B2) G (cid:48) ( x ) = ce − c (1 − x ) . (B3)For a scale-free degree distribution P ( q ) = Aq − γ , witha minimum degree q , we can write: G (1 − x ) = (cid:88) q Aq − γ (1 − x ) q ≈ (cid:90) ∞ q Aq − γ e − qx dq . (B4)Let y = qx , then G (1 − x ) ≈ Ax γ − (cid:90) ∞ xq y − γ e − y dy . (B5)Integrating by parts twice gives G (1 − x ) ∼ = 1 − (cid:104) q (cid:105) x + A Γ(1 − γ ) x γ − + O ( x ) . (B6)The term in order γ − γ . Keepingonly the leading two terms in x (after the constant), wehave that:(i) if γ >
3, then G (1 − x ) ∼ = 1 − (cid:104) q (cid:105) x + 12 (cid:104) q ( q − (cid:105) x , (B7)(ii) if 2 < γ <
3, then G (1 − x ) ∼ = 1 − (cid:104) q (cid:105) x + A Γ(1 − γ ) x γ − , (B8)0(iii) if 1 < γ < G (1 − x ) ∼ = 1 + A Γ(1 − γ ) x γ − + Bx (B9)where the coefficient B of the linear term is no longerequal to the mean degree, which diverges, but insteaddepends on the specific form of the distribution, B = − (cid:40)(cid:88) q q [ P ( q ) − Aq − γ ] − ζ ( γ − (cid:41) . (B10)Note that Γ( z ) < z ∈ ( − ,
0) while Γ( z ) > z ∈ ( − , − Appendix C: Numerical simulations
In this Appendix we describe the numerical proceduresused in our simulations. To calculate each data pointwe use a configuration model method to generate mul-tiple realizations of networks with the same degree dis-tributions. In the illustrative examples of Figs. 3 and 4we use the same distribution in all the layers withoutdegree-degree correlations or correlations between layers.In each layer, and in each realization, we set the degreeof the nodes independently at random according to thefollowing distribution: P SM ( q )= [ (cid:104) q (cid:105) ( γ − γ − ( γ − γ − Γ ( q +1 − γ, (cid:104) q (cid:105) [ γ − / [ γ − q + 1) (C1) ∼ = [ (cid:104) q (cid:105) ( γ − γ − ( γ − γ − q − γ , (C2)where Γ( s ) = (cid:82) ∞ t s − e − t dt is the gamma function, andΓ( s, x ) = (cid:82) ∞ x t s − e − t dt is the incomplete gamma func-tion. The total degree must be even, so, if the sum of alldegrees is odd, we add 1 to the degree of a random node.Each edge is shared by two nodes, so the degree can beseen as the number of ‘half-edges’ belonging to a node.Then, the configuration model inserts edges by joininguniformly at random pairs of ‘half-edges’. This config-uration model imposes no restrictions on the emergenceof self-loops and multiple edges in the network, which,for γ ≤
3, is necessary in order for the degree-degreedistribution to remain uncorrelated [24].Equation (C1) is the exact degree distribution gener-ated by the static model for infinite N [21]. The moti-vation for using this distribution is that, in the small q region, it contains deviations to the asymptotic form ofEq. (C2), which are more realistic than a pure power-lawdistribution.We cannot, however, use the distribution of Eq. (C1)to generate networks with exponent 1 < γ ≤
2, becausethe mean degree (cid:104) q (cid:105) diverges. Notice that our main re-sults are obtained in terms of the amplitude A and theexponent γ of the asymptotics of P ( q ) ∼ = Aq − γ , and, although Eq. (C1) cannot describe them, there are stilldistributions with 1 < γ ≤ A .Figures 5 and 6 show results of simulations for valuesof γ ∈ (1 , γ , we generatethe node’s degrees from a pure power-law distributionwith a minimum degree q : P PL ( q ) = (cid:40) q < q , (cid:16) ζ ( γ ) − (cid:80) q − q =1 q − γ (cid:17) − q − γ q ≥ q , (C3)where ζ ( s ) = (cid:80) ∞ n =1 n − γ is the Riemann zeta function.Notice that, unlike for Figs. 3 and 4, for Figs. 5 and 6we cannot use (cid:104) q (cid:105) as control parameter because of itsdivergence; instead we first generate networks using thedistribution of Eq. (C3), which depends only on γ , andlater apply damage by removing a fraction p of edges atrandom.For our simulations, we choose the minimum degree q = 4 in Eq. (C3) to ensure that the transition is wellobserved at a value of p smaller than 1 for all γ > q sufficiently large.In the range 1 < γ ≤ (cid:104) q (cid:105) leadsto a dramatic increase of the amount of CPU time andmemory required by simulations. Additionally, in thisextreme range of γ , there is another issue in simulationsof (necessarily) finite systems, namely, a single node ac-cumulates a large fraction of all the edges.To elucidate this point, let us consider the effects oftruncating the degree distribution at some cutoff degree ∼ N α when γ <
2. In this case, the average degreeis ∼ N α (2 − γ ) , and the total degree of the network is E = N (cid:104) q (cid:105) ∼ N α (2 − γ ) . The expected number of self-loops of a node of degree q is proportional to q /E , i.e., q times the probability that a ‘half-edge belonging to thatnode is picked uniformly at random out of E possibilities.In particular, for the highest-degree node present in thesystem, with degree q max ∼ N α , the number of self-loopsis ∼ q /E ∼ N αγ − . For the amount of self-loops tobe a vanishingly small fraction of all the edges, the ratiobetween the number of self-loops of the highest-degreenode and the total degree E , which is ∼ N α ( γ − − ,must go to zero as N → ∞ . Then, by using an exponentof the truncation cutoff α < / ( γ − γ ≤
2, while for γ > α that can be used in thewhole range of γ < α = 1, in the simulations ofFigs. 5 and 6 we generated the degrees from distributionstruncated at N , i.e., P ( q < N ) = P PL ( q ) / (cid:80) q 2, which avoids the explosion of self-loops, also requires much lower amounts of CPU time and memory, allowing us to explore larger system sizes N . [1] S.-W. Son, G. Bizhani, C. Christensen, P. Grassberger,and M. Paczuski, Percolation theory on interdependentnetworks based on epidemic spreading, EPL , 16006(2012).[2] S. V. Buldyrev, R. Parshani, R. Paul, H. E. Stanley, andS. Havlin, Catastrophic cascade of failures in interdepen-dent networks, Nature , 1025 (2010).[3] G. J. Baxter, S. N. Dorogovtsev, A. V. Goltsev, andJ. F. F. Mendes, Avalanche collapse of interdependentnetworks, Phys. Rev. Lett. , 248701 (2012).[4] G. J. Baxter, D. Cellai, S. N. Dorogovtsev, A. V. Goltsev,and J. F. Mendes, A unified approach to percolation pro-cesses on multiplex networks, in Interconnected Networks (Springer, 2016) pp. 101–123.[5] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, k -core organization of complex networks, Phys. Rev. Lett. , 040601 (2006).[6] G. Dong, J. Gao, L. Tian, R. Du, and Y. He, Percola-tion of partially interdependent networks under targetedattack, Physical Review E , 016112 (2012).[7] B. Min, S. Lee, K.-M. Lee, and K.-I. Goh, Link overlap,viability, and mutual percolation in multiplex networks,Chaos, Solitons & Fractals , 49 (2015).[8] G. J. Baxter, G. Bianconi, R. A. da Costa, S. N. Doro-govtsev, and J. F. F. Mendes, Correlated edge overlapsin multiplex networks, Phys. Rev. E , 012303 (2016).[9] D. Cellai, S. N. Dorogovtsev, and G. Bianconi, Messagepassing theory for percolation models on multiplex net-works with link overlap, Phys. Rev. E , 032301 (2016).[10] J. Shao, S. V. Buldyrev, S. Havlin, and H. E. Stan-ley, Cascade of failures in coupled network systems withmultiple support-dependence relations, Phys. Rev. E ,036116 (2011).[11] Y. Hu, D. Zhou, R. Zhang, Z. Han, C. Rozenblat, andS. Havlin, Percolation of interdependent networks withintersimilarity, Phys. Rev. E , 052805 (2013).[12] G. Bianconi, Multilayer Networks: Structure and Func-tion (Oxford University Press, Oxford, 2018). [13] M. Kivel¨a, A. Arenas, M. Barthelemy, J. P. Gleeson,Y. Moreno, and M. A. Porter, Multilayer networks, J.Complex Networks , 203 (2014).[14] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Ge-nio, J. G´omez-Garde˜nes, M. Romance, I. Sendi˜na-Nadal,Z. Wang, and M. Zanin, The structure and dynamics ofmultilayer networks, Phys. Rep. , 1 (2014).[15] E. Cozzo, G. F. De Arruda, F. A. Rodrigues, andY. Moreno, Multiplex Networks: Basic Formalism andStructural Properties (Springer, Berlin, 2018).[16] G. J. Baxter, S. N. Dorogovtsev, J. F. F. Mendes, andD. Cellai, Weak percolation on multiplex networks, Phys.Rev. E , 042801 (2014).[17] B. Min and K.-I. Goh, Multiple resource demands andviability in multiplex networks, Physical Review E ,040802 (2014).[18] N. Azimi-Tafreshi, J. G´omez-Gardenes, and S. N. Doro-govtsev, k -core percolation on multiplex networks, Phys.Rev. E , 032816 (2014).[19] R. Cohen, D. Ben-Avraham, and S. Havlin, Percolationcritical exponents in scale-free networks, Physical ReviewE , 036113 (2002).[20] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes,Critical phenomena in complex networks, Rev. Mod.Phys. , 1275 (2008).[21] M. Catanzaro and R. Pastor-Satorras, Analytic solutionof a static scale-free network model, The European Phys-ical Journal B-Condensed Matter and Complex Systems , 241 (2005).[22] K.-I. Goh, B. Kahng, and D. Kim, Universal behavior ofload distribution in scale-free networks, Phys. Rev. Lett. , 278701 (2001).[23] J. Gao, S. V. Buldyrev, S. Havlin, and H. E. Stanley,Robustness of a network of networks, Phys. Rev. Lett. , 195701 (2011).[24] M. Catanzaro, M. Bogun´a, and R. Pastor-Satorras, Gen-eration of uncorrelated random scale-free networks, Phys.Rev. E71