{'fraction_non_alphanumeric': 0.03023332237923102, 'fraction_numerical': 0.009924416694051922, 'mean_word_length': 4.547028800583303, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Large deviation theory of percolation on multiplex networks

Recently increasing attention has been addressed to the fluctuations observed in percolation defined in single and multiplex networks. These fluctuations are extremely important to characterize the robustness of real finite networks but cannot be captured by the traditionally adopted mean-field theory of percolation. Here we propose a theoretical framework and a message passing algorithm that is able to fully capture the large deviation of percolation in interdependent multiplex networks with a locally tree-like structure. This framework is here applied to study the robustness of single instance multiplex networks and compared to the results obtained using extensive simulations of the initial damage. For simplicity the method is here developed for interdependent multiplex networks without link overlap, however it can be generalized to treat multiplex networks with link overlap.

Multiplex networks are often characterized by interdependencies [3] existing between nodes of different layers implying that the failure of one node in one layer necessarily causes the failure of the interdependent nodes in the other layers independently of the rest of the network. Interestingly interdependent multiplex networks are much more fragile than single networks [3][4][5][6]. In particular the percolation threshold p = p c of multiplex networks is larger than the percolation threshold of their single layers taken in isolation. Moreover if the layers are not formed by identical networks, the percolation transition is discontinuous and hybrid and characterized by large avalanches of failures going back and forth among different layers. These large avalanches of failure events are the most notable effect of interdependencies and they are a strong signal of the increased fragility of multiplex networks.
The order parameter of the percolation defined on interdependent multiplex networks is given by the Mutually Connected Giant Component (MCGC) [3,5] which generalizes the giant component and is formed by the set of nodes such that each pair of nodes is connected by a path in at least one layer. In infinite interdependent multiplex networks the percolation transition can be studied using a suitable mean-field approach to percolation [3][4][5][6]. In fact in multiplex network as well as in single networks, percolation theory is self-averaging and in the infinite network limit the mean-field approach gives exact results. In finite interdependent multiplex networks when the initial damage configuration is known, the percolation transition on interdependent multiplex networks can be described by a message passing algorithm [4,6,21,22] that is able to predict which nodes are in the MCGC after the initial damage provided the multiplex network is locally tree-like. Interestingly this message passing algorithm has a definition that depends whether the multiplex network displays or not link overlap, i.e. it depends wheather there are pairs of nodes connected in more than one layer [10][11][12][13]22]. When the specific nature of the initial damage configuration is not known the message passing algorithm can be averaged over the distribution of the initial damage. This mean-field approach provides the probability that each node is in the MCGC when the initial damage configuration is drawn from the known distribution (for instance when each node is damaged with probability f = 1 − p). However this approach is not able to capture the fluctuations that can be observed in the size of the MCGC for different realizations of the initial damage.
Recently, there has been a surge of interest in characterizing the risk of dramatic failure of events as a consequence of a random damage [23][24][25][26][27][28][29]. In fact real networks are typically finite and often sufficiently far from the thermodynamic limit. Therefore the outcome of a initial damage configuration drawn from a given distribution can have large fluctuations not predicted by the mean-field approach to percolation. Characterizing these fluctuations is of fundamental importance for applications as the average response to perturbation captured by the mean-field theory of percolation can be dramatically misleading [28,29]. Moreover characterizing the node sets that are responsible for safeguarding the cohesiveness of a network [28] or for more efficiently disrupting a network (optimal percolation) [30][31][32][33] have a number of applications. These range from the design principles for robust infrastructures to the control of epidemic processes on networks. From the theoretical perspective characterizing these fluctuations is very interesting as well as it sheds new light on the nature of the percolation transition [23][24][25][26][27].
These phenomena are naturally addressed using the theory of large deviations [34,35]. On single networks the characterization of the large deviation properties of random graphs has been first tackled using the properties of the Ising model [36] and the Potts model [37,38]. While the study of the Potts model sheads light on the large deviations of the number of components [37,38] the properties of the Ising model can characterize the large deviations of the giant component in random Erdös and Rényi graphs [36] with constant average degree. However the analytical method developed in Ref. [36] is not able to predict the distribution of the giant component in a real networks or in a network model different from a Erdös and Rényi graph. In order to address this problem a sofisticated numerical technique based on a Markov Chain Monte Carlo method has been shown [39] to be able to uncover the full distribution of the giant component down to rare events with probabilities 10 −100 . This method has subsequently been also used for numerically evaluating the distribution of the diameter [40] and of the largest biconnected component of networks [41].
Recently [24] a large deviation theory of percolation has been proposed to capture the fluctuations observed in percolation of single finite networks. This approach is based on message-passing algorithm and specifically on Belief Propagation [42] and can be applied to single instances of networks provided that they are locally tree-like. This work reveals that when one considers aggravating initial damage configurations, the percolation transition can become discontinuous also in percolation defined on single networks.
Here we generalize the large deviation theory of percolation to multiplex networks. For simplicity we consider only multiplex networks without link overlap. However it is possible to generalize the method to multiplex networks with link overlap. We show that the percolation transition remains discontinuous for both aggravating and buffering configuration of the initial damage. Interestingly we can measure the fluctuations in the size of the MCGC and we observe that these fluctuations can remain very significant up to the percolation transition in the typical scenario. Finally we use the large deviation theory to theoretically predict the convex envelop of the rate function I(R) of the distribution π(R) the size R of the MCGC finding very good agreement with extensive numerical simulations.
The paper is structured as follows. In Sec. 2 we present the mean-field theory of percolation of interdependent duplex networks without link overlap. In Sec. 3 we present the large deviation theory of percolation. In Sec. 4 we show how the Belief Propagation approach can be used to fully characterize the large deviation theory of percolation. In Sec. 5 we compare the numerical results obtained with the Belief Propagation approach with extensive numerical simulations. Finally in Sec. 6 we give the conclusions.
2. The mean-field theory of percolation in interdependent multiplex networks

Multiplex networks and mutually connected giant component
We consider an interdependent duplex network G = (G [1] , G [2] ) formed by two networks [6,20]. Each network α is defined by the adjacency matrix a [α] . For each node i we consider two replica nodes (i, 1) and (i, 2) indicating the identity of node i in layer 1 and in layer 2 respectively. Each pair of replica nodes is connected by an interlink that has a different valence with respect to the links present in each layer. In fact in interdependent multiplex networks interlinks indicate interdependencies whose role in determining the robustness properties of the duplex network is explained in the following. Duplex networks can be classified depending on the presence or the absence of link overlap. The total overlap between the two layers of a duplex network is given by the number of pairs of nodes that are connected in both layers, i.e.
ij a [2] ij . (1) We therefore call duplex networks without link overlap the duplex networks in which the total overlap vanishes and we call duplex networks with link overlap the duplex networks in which the total overlap is greater than zero. We monitor the robustness of the considered duplex networks by studying the size of the Mutually Connected Giant Component (MCGC) [3] after an initial damage of a fraction f = 1 − p of the nodes. Due to the presence of interdependencies, replica nodes cannot be in the MCGC if their corresponding replica node is not included in it. Therefore the MCGC can be evaluated by propagating the initial failure back and forth among the different layers of the duplex network [3,6]. Specifically the algorithm that define the MCGC prescribes to first evaluate the giant component in each individual layer and subsequently to damage each replica node whose corresponding interdependent replica node is not in the giant component. This algorithm is then iterated until the damage does not propagate any more among different layers. Alternatively the MCGC can be defined as the giant component formed by pairs of nodes connected to each other at least by a path in each layer [5].

Message passing algorithm on single network and single realization of damage
On a locally tree-like interdependent duplex network without link overlap percolation can be described by the following message passing equations. We assume that initially each node of the duplex network is damaged randomly with probability f = 1 − p, and we associated a variable x i = 1 to a node i that is not initially damaged and a variable x i = 0 to a node i that is initially damaged. Additionally we assign to each node i a variable ρ i indicating if the node belongs to the MCGC (ρ i = 1) or not (ρ i = 0). The indicator function ρ i depends on the configuration of the initial damage x = {x i |i ∈ {1, 2, . . . , N }} and the messages σ α i→j that each node i send to a neighbour node j in layer α. In particular a node i belongs to the MCGC if it is not initially damaged and receives at least one positive message in each layer,i.e.
The messages σ α i→j ∈ {0, 1} take the value σ α i→j = 1 if the node i is not initially damaged and if in each layer node i receives at least a positive message from a neighbour nodes different from node j, i.e.
where here and in the following β indicates the layer β ∈ {1, 2} with β = α. The size R of the MCGC is given by the number of nodes that belong to it, i.e.
We note that different initial damages can yield MCGC of different sizes. The goal of this paper is to uncover the fluctuations that can be observed in the sizes of the MCGC resulting from comparable configurations of the initial damage.
We stress here that Eqs. (2), (3) are only valid in absence of link overlap. In fact when the multiplex network display link overlap these equations only define the Directed Mutually Connected Component (DMCGC) which can be considered as the outcome of an epidemic spreading dynamics in which each node can become infected only if is in contact with at least an infected individual in each layer [10,11]. In presence of link overalp the MCGC is instead captured by a message passing algorithm in which the messages are not scalars like σ α i→j but they are instead M -dimensional vectors n i→j of elements n α i→j ∈ {0, 1}. Here n α i→j indicates wherether (n α i→j = 1) or not (n α i→j = 0) node i connects node j to the MCGC through layer α as long as node j is assumed to be in the MCGC [11,12]. These message elements are set of one (n α i→j = 1) if the following two conditions are met: (i) node i is in the MCGC; (ii) node i is connected in layer α to at least one node of the MCGC different from node j.
Additionally node i is in the MCGC if for any layer of the multiplex network node i is either connected to node j or connected to the MCGC through nodes different from j. Therefore the message passing equations for n α i→j read where δ(c, d) indicates the Kronecker delta and v i→j is defined as Therefore the message elements n α i→j are correlated and this property of these message passing equations considerably complicates the algorithm. Finally in presence of link overlap the indicator function σ i determining if a node is in the MCGC satisfies Naturally a close inspection to the message passing equations in presence of link overlap show that these equations reduce to the Eqs. (2) and (3) in absence of link overlap. Moreover we note that in absence of link overlap the equations for the MCGC and for the DMCGC coincide. We note that the message passing Eqs. (5), (6) and (7) can be also written in terms of the indicator functions σ n i→j ∈ {0, 1} that indicates if node i sends to node j the message n (σ n i→j = 1) or not (σ n i→j = 0). For a detailed account of the message passing theory for percolation of interdependent multiplex networks see Ref. [6,11].
Here and in the following we will consider exclusively multiplex network without link overlap. In the following we will indicate with σ the set of all the messages and with σ i the set of all the messages starting or ending at node i, i.e.
Additionally we will indicate with x the configuration of the initial damage, i.e.

Random realization of the damage and typical behaviour
Let us consider initial damage configurations x obtained by damaging each node with probability 1 − p, i.e. each configuration x is drawn from a distributioñ In the mean-field theory of percolation the expected sizeR of the MCGC given bŷ is obtained by averaging the original message passing algorithm over the distributioñ P (x). On a locally tree-like multiplex network without link overlap this procedure reduces to considering a message passing algorithm where the messagesσ α i→j correspond to the average of the messages σ α i→j over the distributionP (x), i.e.
Similarly the probabilityρ i that node i is in the MCGC can be obtained starting from the messagesσ α i→j aŝ and corresponds to the average of the indicator function ρ i over the distributionP (x), Finally the expected size of the MCGCR can be calculated directly aŝ Therefore by monitoringR versus p it is possible to characterize the typical response of the duplex network to random damage of the nodes

Random network ensemble and random realization of the initial damage
In the previous paragraph we have considered the message passing algorithm that allows to study percolation the mean-field theory of percolation on single instances of duplex networks when the initial damage configuration is unknown but drawn from the distributionP (x) given by Eq. (10). Here we consider the case in which also the duplex network is unknown but drawn from ensemble of duplex networks formed by two layers each one formed by an uncorrelated network with degree distribution P (k). In this case we cannot any more consider a message passing algorithm, rather we should consider the mean-field equation [3][4][5][6] for percolation of interdependent duplex networks without link overlap obtained by averaging Eqs. (12) and (14) over the considered duplex network ensemble. We indicate with R the probability that a node is in the MCGC. This quantity is given by the average ofρ i over the considered ensemble of duplex networks. Similarly we indicate with S the probability that by following a link in one layer we reach a node in the MCGC. This quantity is given by the average ofσ i→j over the duplex network ensemble. Using this notation the mean field equations for percolation of independent networks without link overlap [3,4,6] read where G 0 (x) and G 1 (x) are the generating functions These equations yield a discontinuous phase transition at the percolation threshold p = p c determined [3][4][5][6] by the equations (17) together with the equation 3. Large deviation theory of percolation

The large deviation approach to percolation
While in infinite networks the percolation process is self-averaging, i.e. the typical behaviour characterizes a set of measure one of random instances of the initial damage, in finite networks deviations from the typical behaviour can be observed. It is therefore necessary to establish the large deviation theory of percolation. Here our aim is to generalize the framework proposed in Ref. [24] for characterizing the large deviation in single networks in order to treat the large deviation of percolation in interdependent multiplex networks. The large deviation theory of percolation of multiplex networks aims at characterizing the probability distribution π(R) of size R of the MCGC when the initial damage configuration x is chosen with probabilityP (x), i.e.
where with δ(c, d) we indicate the Kronecker delta. For large network sizes N 1 and given value of p the probability π(R) has a scaling with the network size N determined by the rate function I(R) ≥ 0. In particular we observe the scaling [34] π(R) ∼ e −N I(R) .
This expression implies that for N → ∞ the percolation is self-averaging and π(R) is non-zero only for s R =R for which I(R) takes its minimum value I(R).Therefore in the infinite network limit all realizations of the initial damage yield almost surely the same size of the MCGC. However the rate function I(R) captures the fluctuations in the size of the MCGC that can be observed in finite multiplex networks. In order to find I(R) we adopt a canonical approach and we introduce the partition function Z = Z(ω) given by Using the definition of π(R) given by Eq. (20) it can be easily shown that Z is the generating function of π(R) as Z can be written as The corresponding free-energy F and the free energy density f can be calculated as The Legendre-Fenchel transform of the rate function I(R) [34] can be expressed in terms of the free energy density as ωf (ω) and we have Additionally as long ωf (ω) is differentiable, the Legendre-Fenchel transformÎ(R) of ωf (ω) given bŷ fully determines the rate function I(R), i.e. I(R) =Î(R). However when I(R) is non-convex ωf (ω) is not differentiable and the Legendre-Fenchel transform of ωf (ω) given byÎ(R) only provides the convex envelop of the rate function I(R) [34].
Therefore the knowledge of the biased distributionπ(R) which can be sampled with the Markov Chain Monte Carlo method is in principle sufficient to reconstruct the full distribution π(R). However this approach is entirely numerical. In the following we will consider the first approach and we will evaluate the partition function Z using the Belief Propagation algorithm.

The Gibbs measure over messages
In order to calculate the partition function Z we follow the approach proposed in Ref. [24] and we consider a Gibbs measure P (σ) over the set σ of all messages. The probability distribution P (σ) weights the configuration of the messages σ according to the probability of the corresponding initial damage configurations. Moreover we introduce a Lagrangian multiplier ω conjugated to the size of the MCGC R that is able to tune the relative weight of the configurations of the messages corresponding to different sizes of the MCGC. Therefore the Gibbs measure P (σ) is given by where the function χ(σ, x) enforces the message passing Eqs. (3) , i.e.
The partition function Z in Eq. (29) clearly reduces to Z defined in Eq. (23). In fact we have Therefore characterizing the partition function Z and the free energy F (given by Eq. (24)) corresponding to the Gibbs measure defined in Eq. (29) allows us to directly calculate the Legendre-Fenchel transform of the rate function I(R). We note that the Gibbs measure P (σ) can be interpreted as a canonical ensemble, where ω plays the role of the inverse temperature and R plays the role of the energy. Since R is the sum of the node variable ρ i and ρ i can only take two values (ρ i = 0the node does not belong to the MCGC or ρ i = 1 the node does belong to the MCGC) the Gibbs measure can be interpreted as a a statistical mechanics problem of a two level system. It follows that in this case we can investigate the properties of the Gibbs measure for values of ω that can be also negative.
For ω < 0, the Gibbs measure weights more the buffering configurations of the initial damage resulting in a MCGC larger than the typical one. On the contrary for ω > 0 the Gibbs measure weights more the aggravating configurations of the initial damage resulting in a MCGC smaller than the typical one. For ω = 0 we recover the typical scenario.
The Gibbs measure P (σ) given by Eq. (29) can also be expressed as where the set of constraints ψ i (σ i , ω) for i = 1, 2, . . . , N defined over all the messages σ i starting or ending to node i read Here δ(m, n) indicates the Kronecker delta andρ i is given bỹ Finally, using Eq. (31) it can be easily shown that the partition function Z can be also written as In the following sections we will characterize the large deviation properties of percolation on multiplex networks by calculating the Gibbs measure using Belief Propagation (BP).

The Belief Propagation equations
On a locally tree-like duplex network without link overlap the Gibbs distribution P (σ) can be expressed explicitly using the Belief Propagation (BP) algorithm [42] by finding the messagesP α i→j (σ α i→j , σ α j→i ) that each node i sends to the generic neighbour node j in layer α. These messages satisfy the following recursive BP equations (see Appendix A for their explicit expression) i→j are normalization constants enforcing the normalization condition In the Bethe approximation, valid on locally tree-like networks the probability distribution P (σ) is given by [42] where P i (σ i ) and P α ij (σ α i→j , σ α j→i ) indicate the marginal distribution of nodes and links and are given by [42] with C i and C

[α]
ij indicating normalization constants (see Appendix B for their explicit expression).

Free energy
The free energy F given by Eq. (24) can be found by minimizing the Gibbs free energy F Gibbs given by where ψ(σ) indicates the set of constraints In fact it can be easily shown that the Gibbs free energy F Gibbs is minimal when calculated over the probability distribution P (σ) given by Eq. (31) and that its minimum value is On a locally tree-like duplex network the Gibbs measure P (σ) reduces to Eq. (35), and it can be easily see that the free energy can be expressed as Given the explicit expression of C [α] ij and of C i in terms of the messages (see Appendix B), the free energy F can be calculated easily when the BP equations have been solved.

Energy and Specific Heat
The energy and the specific heat corresponding to the Gibbs measure P (σ) have also a very important interpretation in terms of the underlying percolation process. The energy is given by the the average size of the MCGC R. In fact we have The mean-field percolation transition correspond to the phase transition from a nonpercolating phase with R = 0 to a percolating phase R > 0 at p = p c and ω = 0. However the present approach allows to consider the full line of critical points (p , ω ) at which the transition is observed when the large deviations are considered. Following a statistical mechanics definition, we can also define the specific heat C as The specific heat has the immediate interpretation in terms of the variance in the size of the MCGC, i.e.
Both R and C/ω 2 can be derived from the message passing algorithm. Indeed we have where indicating the probability that node i is in the giant component is given by where the explicit expression of z i and C i is given in Appendix B. The quantity C/ω 2 can be also interpreted as the expected fraction of nodes that given two random realizations of the initial damage are found in the MCGC in one realization but not in the other. This quantity generalizes the measure proposed in Ref. [23] to characterize the fluctuations in percolation in single networks.

The typical scenario
It is instructive to see that the proposed large deviation theory of percolation reduces to the mean-field theory of percolation in the typical scenario obtained by putting ω = 0.
In this case we obtain that the BP equations have solution witĥ and the BP equations for ω = 0 reduces to Eqs.(12) when we put Similarly it is easy to show that r i reduces toρ i given by Eq. (14).

Numerical results
Here we consider the results obtained by running the BP algorithm over a Poisson duplex network with average degree z = 6 and network size N = 100 (see Figure 1). The typical scenario observed for ω = 0 gives a discontinuous transition of the average size R of the MCGC at p = p c = 0.409235 . . . with R = R c = 0.209405 . . . for the considered duplex network as predicted by the mean-field theory of percolation. The fluctuations of the size of the MCGC measured by C/ω 2 for ω → 0 have also a jump at p = p c . Therefore as we approach the critical percolation threshold from above (i.e. for p → p + c ) we observe significant fluctuations in the size of the MCGC. However these fluctuations are maximal only for larger values of p.
Interestingly the BP results allow us to investigate also how the average size of the MCGC R and its fluctuations C/ω 2 change if we deviate from the typical scenario, i.e. for ω = 0. For ω < 0 when we consider buffering configuration of the initial damage we observe a transition for lower values of p, for ω > 0 we observe a transition for larger value of p. In both cases the transition remains discontinuous in the investigated range of values of ω. Additionally we note that for buffering configurations, as ω decreases the discontinuous jump in the size of the MCGC R appears to reach a constant value. This suggests that the discontinuity might be preserved even beyond the observed range of values of ω. The fluctuations in the size of the MCGC C/ω 2 observed at the transition point increase for aggravating configuration of the damage. Moreover for aggravating configuration of the damage these fluctuations can achieve their maximum at the transition point itself.
In Figure 2 we compare the results obtained for the duplex networks to the results that can be obtained by studying the large deviations of percolation on single networks as discussed in Ref. [24]. In the duplex network we observe a discontinuous jump of C/ω 2 also in the typical scenario ω = 0 while in single networks C/ω 2 is continuous at the transition point. Therefore the observed discontinuity of C/ω 2 at p = p c and ω = 0 is a purely multiplex network phenomenon not observed in percolation of single networks. In fact the presence of significant fluctuations of the percolation order parameter at the percolation transition can be only observed if the order parameter has a discontinuous jump and in this case only for p → p + c when R → R c > 0. Finally we note that in the multiplex case the discontinuity of C/ω 2 extend also for negative value of ω corresponding to buffering configurations of the initial damage while in single networks we observe a continuous behaviour of C/ω 2 .
In order to test the validity of the BP algorithm, we have compared the rate function I(R) measured starting from 10 6 initial damage configurations performed on the same duplex network instance on which the BP algorithm is run, withÎ(R) provided by the BP algorithm finding very good results (see Figure 3). We notice that while for large values of p the rate function I(R) is convex and therefore I(R) =Î(R), as we approach the percolation transition the rate function I(R) becomes non-convex andÎ(R) only provides the convex envelop of I(R).

Conclusions
In this paper we have characterized the large deviation of percolation of interdependent multiplex networks without link overlap using a Belief Propagation algorithm. In the typical scenario, well captured by the mean-field theory of percolation we observe a discontinuous percolation transition. Our analysis reveals that when we depart from the typical scenario the percolation transition can occur for values of p significantly distant from the percolation threshold in the typical scenario while remaining discontinuous both for buffering configurations of the random damage and aggravating ones. Interestingly our study show a significant difference of interdependent percolation with respect to percolation in single layers. In fact when we consider percolation of isolated networks the fluctuations in the size of the giant component go to zero at the critical point. However in interdependent percolation we can observe significant fluctuations of the size of the MCGC when we approach the critical point from above, i.e. for p → p + c . Finally we have compared the rate function I(R) of observing a MCGC of size R with the predictedÎ(R) obtained using the proposed canonical BP algorithm. As observed in the single layer scenario the rate function I(R) becomes bimodal for small values of p, and therefore the proposed canonical BP algorithm can only capture its convex envelop, i.e. we haveÎ(R) = I(R).  The proposed method can be extended in different ways by considering percolation and directed percolation in multiplex network with link overlap and by developing a microcanonical approach in which instead of introducing the Lagrangian multiplier ω fixing the size of the MCGC in average, we consider the corresponding hard constraint. When for arbitrary degree k In this appendix we give the explict expression of C α ij , C i and z i in terms of the messages that solve the BP equations. In particular by normalizing the marginals in Eqs. (36) we obtain C α ij =P α i→j (0, 0)P α j→i (0, 0) +P α i→j (0, 1)P α j→i (1, 0) +P α i→j (1, 0)P α j→i (0, 1) +P α i→j (1, 1)P α j→i (1, 1),