Bond percolation in coloured and multiplex networks

Percolation in complex networks is a process that mimics network degradation and a tool that reveals peculiarities of the network structure. During the course of percolation, the emergent properties of networks undergo non-trivial transformations, which include a phase transition in the connectivity, and in some special cases, multiple phase transitions. Such global transformations are caused by only subtle changes in the degree distribution, which locally describe the network. Here we establish a generic analytic theory that describes how structure and sizes of all connected components in the network are affected by simple and colour-dependent bond percolations. This theory predicts locations of the phase transitions, existence of wide critical regimes that do not vanish in the thermodynamic limit, and a phenomenon of colour switching in small components. These results may be used to design percolation-like processes, optimise network response to percolation, and detect subtle signals preceding network collapse.

One of the richest tools for exploring properties of complex networks is probing these networks by randomly removing links.This process is called percolation, and on many occasions percolation has shaped our understanding of physical and socio-economic phenomena [1].Naturally, percolation is related to network resilience.This connection was exploited, for instance, in studies on communication, transportation, and supply networks [2][3][4][5].Percolation has defined the modern view on disease epidemics [2,[6][7][8][9], as well as on other spreading processes [10][11][12].Percolation is instrumental in material science, where the gelation [13][14][15][16] and jamming [17,18] have been both connected to this process.
It is known that during percolation, complex networks undergo a series of non-trivial transformations that include splitting of connected components and criticality in the global connectivity [5,[19][20][21][22].Remarkably, even the simplest models that define a network solely by its degree distribution, the configuration models, do also feature these phenomena.In many studies, this observation justified the use of the configuration model as the null-model, which one would compare against any real network, in order to reveal presence of trends that are, or are not, explainable by the degree distribution alone.
In the edge-coloured configuration model, the edges are categorised into different groups or colours.This distinction leads to a more realistic representation of complex networks.Indeed, most of real networks do feature different types of interactions: chemical bonds, communication channels, social interactions between infected individuals, among other examples, all require a partition of the edges into discrete categories [9,[23][24][25].Furthermore, a whole new world of possibilities unfolds when the number of colours is virtually unlimited: different edge colours may encode all combinations of one interaction occurring between different types of nodes, as was proposed in the assortative mixing model [26], or one type of interaction with discrete strengths [27].In modular networks, interactions within each community can be represented with a unique colour.In contrast to unicoloured networks [28], less is known about how percolation happens in coloured networks.For instance, presence of multiple phase transitions has been reported [20,21,24], but no complete theory could explain their nature nor predict locations where these phase transitions occur.It has been observed that negative correlation in the degree distribution leads to multiple phase transitions [24].This paper demonstrates that in order to explain the percolation in coloured networks, one has to study not only the giant component but also the rest of the system.In fact, the tail asymptote for the sizes of connected components already covers sufficient amount of information.By building upon this result, we establish universal criteria of phase transitions for two percolation processes: simple bond percolation, which removes every edge with equal probability, and the colour-dependent bond percolation, in which edges of different colours have different probability to be removed.This theory detects multiple phase transitions, if such occur, and in the case of colour-dependant percolation, the theory yields a manifold containing all critical points.Surprisingly, the theory predicts existence of wide critical windows that do not vanish in infinite systems and colour switching that occurs in small connected components during percolation.

I. RESULTS
A. The N -colour asymptotic theory Newman et al. observed that in a unicoloured network defined by the degree distribution, the size distribution of connected components can be found numerically [19].His ideas have been later formalised into an analytical theory in our previous work [29] (see SI 1 for a brief summary).The current work takes this analytical concept further by considering a network with an arbitrary number of colours.Suppose every edge has one of N colours, so that a randomly chosen node bears k 1 edges of colour one, k 2 edges of colour two, and so on.A state of a node is parametrised by a vector of colour numbers k = (k 1 , . . ., k N ).Let u(k) be the probability that a randomly chosen node has state k.In this paper, we often refer not to u(k) directly, but to its parameters µ 0 , M , and T i that are defined via expectations of u(k): , i, j, l = 1, . . ., N.
The probabilistic interpretation of µ 0 , M , T i and the concept of the expected value is explained in the Methods.In our previous paper [30], the formal expression for the size distribution of connected components w(n) was derived by applying Joyal's theory of species, where operations f * g and f * k denote, respectively, the N -dimensional convolution product and the convolution power [30], whereas D is the auxiliary function defined in SI 2. Although the computations of equation (S7) require large complexity of O(n N log n), it is possible to write the tail asymptote for arbitrary number of colours N : Here coefficients C 1 > 0, C 2 ≥ 0 are defined in terms of µ 0 , M and T i .The Methods provide the expressions for these coefficients, and SI 3 contains the derivation of equation (5).Coefficient C 2 is especially important as it defines how fast the size distribution decreases: C 2 0 implies exponentially rapid decrease; C 2 = 0, triggers a scale-free behaviour, and close-to-zero values of C 2 are associated with transient scale-free behaviour that is eventually overrun by an exponentially fast decrease.The universality of − 3  2 constant has been noticed by Newman et al. in unicoloured networks [19].Remarkably, the exponent − 3 2 also holds universally for arbitrary number of colours N and an arbitrary multi-degree distribution that features finite mixed moments In the configuration models, a connected component of size n has a tree-like structure, and therefore, contains n − 1 edges.Let v i denote the numbers of i−coloured edges that are present in a component of size n 1.By normalising with the total number of edges f i := vi n−1 , one obtains the respective colour fractions.The probability of a particular configuration of colour fractions P[f 1 , f 2 , . . ., f N ] is the Gaussian function with independent of n mean m and covariance matrix 1  n Σ that vanishes when the component size tends to infinity.The expressions for m and Σ are given in terms of M and T i in the Methods, and the necessary derivations can be found in SI 7.

B. Criticality and simple bond percolation
There exists a class of degree distributions for which C 2 = 0 and therefore the asymptote ( 5) is scale-free.We call these degree distributions critical.As shown in SI 4, the critical degree distribution can be identified by testing with a following criterion: Before testing with criterion (6), it is useful to check first if det[M −I] = 0, which is only a necessary condition for criticality but not a sufficient one.Equation ( 6) generalises many known relations, for instance, Molloy-Reed criterion, equation 33 in [19]; the criterion for the multiplex giant component, equation 11 in [31], or the message passing criticality, equation 27 in Ref. [32].Now, let us investigate a simple process that continuously changes the degree distribution in such a way that the latter traverses though one of the critical points.Consider a process on the edge-coloured network that thinners the network by randomly removing edges with probability 1 − p, or equivalently, by keeping edges with probability p.Such removal of edges alters the degree distribution, so that: µ 0 = pµ 0 , M = pM and (T i ) j,l = p 2 (T i ) j,l + p(1 − p)M j,i δ j,l .As shown in SI 5 by plugging M into the criticality criterion (6) one obtains the following p-dependant criterion: the edge-coloured network features the critical behaviour at p = p c ∈ (0, 1] if there is vector v for which, An alternative way of looking at equation ( 7) is reformulating this criterion as an eigenvalue problem: the edge-coloured network features the critical behaviour at p c = λ −1 if the following conditions hold true: v is non-negative when normalised, i.e. v |v| ≥ 0. These three conditions provide a generalisation to a percolation criterion that has been introduced elsewhere [24] (equation 7).Indeed, suppose M i,j > 0, then according to the Perron-Frobenius theorem, M has exactly one positive eigenvector and this eigenvector corresponds to the largest eigenvalue so that the third condition of the above-mentioned list is automatically satisfied.In which case, one can write a simplified criterion: The essential difference between condition (8) and the generalisation (7) lies in the fact that equation (8) predicts at most one critical point p c , whereas equation ( 7) recovers multiple critical points, if such exist.

C. The hierarchy of critical points and secondary phase transitions
To this end, our notion of connectivity was based on a multicoloured path, i.e. the path that consists of any colours form index set {1, . . ., N }.By excluding some of the colours, one may define stronger notions of connectivities.Any non-empty subset S ⊂ {1, . . ., N } gives rise to a valid definition of a path, and all possibilities collected together comprise a power set P{1, . . ., N } \ {∅}, i.e. the set of all subsets.This power set features a natural hierarchy as defined by inclusion relation "⊂".For example, let N = 3, then if two nodes are connected with a path having colours from {1, 2}, this nodes are also connected by a weaker notion of the path that may include colours from {1, 2, 3}, since {1, 2} ⊂ {1, 2, 3}.Every connectivity notion S = {i 1 , i 2 , . . ., i k } ⊂ {1, . . ., N } is associated with a distinct notion of criticality, and the inclusion of these subsets imposes a certain ordering of critical points p S , for instance, in the above example: p {1,2,3} ≤ p {1,2} .Criterion (7) can be adjusted to detect these S-criticalities: let x S is the indicator vector for colour subset S and (λ, v) is an eigenpair of then the model features S-criticality at p S = λ −1 if v |v| ≥ 0. In this way, one has means of identifying all critical points p S and all colour subsets that correspond to them.
Although, recovering all p S requires solving 2 N −1 eigenvalue problems as defined by equation ( 9), much can be said on how p S are distributed when M is a diagonal-dominant matrix.In this case, the critical points are grouped into batches, and the locations of this batches are given by the Gershgorin circle theorem: each critical point p S belongs to at least one of the following intervals:

D. Colour dependant percolation
Instead of single percolation parameter p, consider a vector p composed of p i , the probabilities that i-coloured edge is not removed.Colour-dependant percolation alters the degree distribution parameters in the following fashion: µ 0 = diag{p}µ 0 , M = diag{p}M , T i = diag{p}T i diag{p} + diag{p}diag{1 − p}diag{M 1,j , . . ., M N,j }.The derivations of this transformations can be found in SI 6.By plugging M into equation ( 6) one obtains the criterion for colour-dependent percolation: edge-coloured network features the critical behaviour at 0 < p < 1 if and only if The latter criterion can be viewed as a parameter equation for a surface placed in an N -dimensional space, and, unlike in the case of simple percolation, one cannot reduce criterion (11) to an eigenvalue problem.

E. Description of the giant component
The node-size of the giant component g node is the probability that a randomly sampled node belongs to the giant component.This probability can be written in terms of expected values of the degree distribution, where s i are obtained by solving and e i are the standard basis vectors.In a similar fashion, the edge-size of the giant component, the probability that a randomly sampled edge of colour i is a part of the giant component, is given by g i = 1 − s 2 i .The weight-average size of finite connected components is given by: where } and X(s) is a matrix function with the following elements , i, j = 1, . . ., N.
When subjected to the bond percolation, g node and w avg become functions of percolation parameter p: and In this way, the original degree distribution provides enough information to describe how the size of the giant component and the average size of finite components evolve during the percolation progress.Equation ( 14) has been also derived in refs.[24,31], and (15) generalises similar relation for uncoloured networks [16].The derivations for equations ( 12)-( 15) are given in SI 8.

F. Stochastic simulations of percolating coloured networks
In this section, we compare the asymptotic theory and stochastic simulations for a few peculiar examples.Consider a configuration model with three colours that is defined by where Poiss k, λ := e −λ λ k k! is the Poisson mass density function and C ensures the appropriate normalisation.From equations ( 2),(3) one obtains matrices µ 0 , M , T 1 , T 2 and T 3 , which define the asymptotical properties of the network, and from equation ( 5) one obtains the asymptote of the size distribution: Figure 1a compares this asymptote to the size distribution obtained from the simulations.This asymptote appears to be almost scale-free.This is in accordance with the close-to-zero value of C 2 = 0.00043: the network is very close to the critical point, however, one cannot assess if the network is just below or just above the phase transition by simply looking at the asymptote.Since Perron-Frobenius conditions are satisfied for M , we can safely apply the simplified test (8): the network is below its critical point since ρ(M ) < 1. Simple bond percolation with occupancy probabilities p = 0.9 and p = 0.8 gives a progressively faster decreasing size distribution.The colour fractions in the non-giant components settle down on an uneven proportion in components of large size n 1.According to equation ( 25), a connected component contains, in average, 19% edges of colour 1, 74% edges of colour 2, and 7% edges of colour 3. The covariance matrix of this distribution, as defined in equation ( 26), features negative correlation and vanishes when n is large: Figure 1b compares the theory with simulated colour fractions.At first glance, this tendency may seem to be contra-intuitive: the larger component is, the less random structure it attains.The intuition behind this phenomena lies in the fact that since the size of the component is conditioned to be a large number n, the model finds the optimal balance between two contradictory forces: selecting as many nodes of high degree as possible versus selecting nodes that are most abundant according to the degree distribution, which are in this case the nodes with small degree and specific configuration of colours.The variance of the spread present in the scattered data in Fig. 1b is well explained by the theory.This spread is an outcome of the fact that finite components in infinite systems feature fluctuations in their internal structure.
In the next example, we consider degree distribution When α = 0, each node has solely links of one colour, so that the whole system is a composition of three unicoloured networks.This fact results in a somewhat exotic situation when multiple giant components can stably coexist, which is an attractive phenomenon from the perspective of many applied disciplines [33].When α > 0, some nodes have links with many colours and the whole network contains, if any, one giant component.Nevertheless, as Fig. 2a shows, the network is modular in both cases, which points towards the utility of coloured edges for modelling networks with community structures.In the case α = 0, the dependence of C 2 on the percolation probability p, as given by the solid line in Fig. 2d, reveals that such criticalities occur three times at p c = 1 5 , 2 5 , and 2 3 .Fig. 2e shows that the weight-average component size is singular at these critical points.The values of p c can be easily deduced from matrix M , which is diagonal when α = 0: In this case, the diagonal elements of M are also its eigenvalues, and all three eigenvectors are positive when normalised.So that, according to the criterion (7), all three eigenvalues are associated with valid phase transitions.Since all off-diagonal elements are zeros, one can invoke equation (10) to show that all secondary phase transitions coincide with the primary critical points, the complete list is given in SI 9. Note, that when α = 0, matrix M contains zero elements, and therefore the simplified phase transition criterion (8) is not applicable.
Setting α = 0.1 in degree distribution (19) perturbs M so that it becomes a full matrix with small off-diagonal elements, Although the total size of the giant components, as shown in Fig. 2e and Fig. 2f differs only a little, the eigenvalue decomposition of M shows that there is only one positive eigenvector, and therefore, one phase transition.This observation is supported by the fact that C 2 (p) = 0 only once, as shown by the dashed line in Fig. 2d, and that the weight-average component size is singular only at p ≈ 0.21, as indicated in Fig. 2f.One can yet observe that C 2 (p) has two local minima at the locations where the case of α = 0 features phase transitions.
Since M is strongly diagonal dominant for α = 0.1, the secondary critical points appear in groups.The full list, as obtained from criterion (9), is given in Fig. 2b.Phase transitions associated with three colours feature the hierarchy as indicated by the partially ordered set in Fig. 2c, one may also think of Fig. 2b as a linear sorting of this partial order.Probably the most surprising implication of the equation describing the internal structure of connected components is not that the colour fractions settle on an uneven ratio, as depicted in Fig. 1b, but that this ratio peculiarly evolves as a function of p.For instance, in the case of α = 0.1, the colour fractions f 1 , f 2 , and f 3 feature a switching behaviour.Non of the switching points coincides with the phase transitions, but as Fig. 3a reveals, they are rather equidistant from the critical points.This trend may be exploited to device early warning strategies that, similar to those devised in [34], detect proximity of a phase transition in empirical networks.In contrast to large components, in which the colour fractions converge to precise values as shown in Fig. 1b, the structure of small components is not deterministic but features fluctuations.The variance of these fluctuations is predicted by the theory and does not vanish in large networks.Figure 3b compares the simulated data for components of size n = 50 against the theoretical standard deviation.In colour-dependent percolation, one investigates the properties of the percolated network as a function of probability vector p = (p 1 , p 2 , p 3 ).All configurations of p amount to the volume of a unit cube.Figure 4 presents the regions of this parameter space where the network becomes critical, that is C 2 = 0, or close-to-critical, C 2 is small.This configurations are recovered by numerically solving equation ( 11) that parametrises the corresponding manifolds.One can see that what appeared as single curves in Fig. 2d,e for simple percolation, is now a surface placed in the unit cube, see Fig. 4. When p 1 = p 2 = p 3 , which corresponds to a diagonal of this unit cube as indicated by the yellow line, the colour-dependant and simple percolations are equivalent.Fig. 4b shows that in the case α = 0.1, the critical points form a box-shaped surface, whereas α = 0 changes this surface into a more complex shape, see Fig. 4a.Note, in Fig. 4a the yellow line intersects this surface three times, which corresponds to the three phase-transition points that are also observed in Fig. 2e.
Although there is a conceptual difference between the cases presented in Fig. 4a,b, one would expect that the actual networks should be close in some sense.Indeed, Fig. 4b is obtained by perturbing the network represented in Fig. 4a with a small parameter α.This similarity can be highlighted if we depart from our strict definition of criticality, C 2 = 0, in favour of a weaker criterion: C 2 < 0.02.All configurations of p that satisfy C 2 < 0.02 form together the three-dimensional manifolds presented in Fig. 4c,d: one can see that these manifolds do bear resemblance.Now, the time evolution the colour-dependant percolation can be represented with a path p(t), t ∈ [0, 1] that starts at p(0) = (1, 1, 1) , which corresponds to the intact network, and ends at p(1) = (0, 0, 0) -a completely disintegrated one.Furthermore, the question of how an individual network responds to percolation is about how this path relates to the above-described geometric structures, and it turns out that colour-dependant percolation can lead to completely new behaviours that are not observed in simple percolation.In order to demonstrate this fact, we devise a path p(t) in such a way that C 2 (p) FIG. 5. Colour-dependant percolation that features a wide critical window.Networks corresponding to the degree distribution defined by equation (19) with: a,b,c, α = 0 and d,e,f, α = 0.1.a,d, optimal percolation paths p(t) together with the reference surfaces on which C2 = 0.02.b,e, the C2-profiles corresponding to the optimal paths.c,f, comparisons of the theoretical weighted-average component size profiles against the corresponding profiles extracted from stochastically generated networks with 10 5 nodes.
is minimal on it.This is achieved by performing minimisation  ,d, whereas Figs 5b,e depict the corresponding profiles of C 2 (t).One immediately notices that in these examples the profiles of C 2 (t) vanish not at discrete points, as was the case with simple percolation, but on intervals of continuum.Everywhere on this intervals the network is critical: that is the sizes of connected components feature the scale-free behaviour and thus the weight-average component size is infinite.As Fig. 5c,f shows, this theoretical considerations are also supported by numerically generated networks.Generated networks comprise of finite number of nodes, and therefore they cannot feature infinite average component size.However, this quantity features macroscopic fluctuations within the critical window and diverges to infinity with growing system size.Remarkably, the width of this critical window does not shrink to zero in infinite systems.This phenomenon constitutes a new type of phase transitions in coloured networks: the phase transitions with wide critical windows.Within this critical window the average size of a connected component features macroscopic fluctuations.

II. CONCLUSIONS
In this paper, edge colours provide an abstraction of an additional layer of information that a network can be equipped with, but also an abstraction of various network structures.The colours may represent an affiliation to communities, multiplexity, different types of interactions, assortative/disassortative relationships, and other aspects prevalent to complex networks.We have shown that colour dependencies may amplify failures in connectivity or make the collapse of the network less sudden.This is why this theory provides the foundation for answering the key question of complex networks science: "How we can we economically design robust multilayer networks of infrastructures or financial networks with trustworthy links?" [35].In many ways a modeller exploit the abovedescribed geometric interpretation of colour-dependant percolation in sake of network design and control.By following similar motivations to the ones in Ref. [4], one might aim to optimise the percolation path so that minimal/maximal number of edges is removed before percolation reaches the phase transition, or for a fixed path, one might optimise the network robustness so that the phase transition manifold avoids the path as much as possible.Another objective, reducing the sharpness of the phase transition [36,37], can be achieved by reducing the angle between the path and the surface at the intersection point.In fact, as we have demonstrated by an example, one may even construct a path that do not immediately intersect the manifold but stays inside for a long time, and therefore, keeps the network in the critical window.The latter observation shifts the paradigm of the critical window itself, as it demonstrates that even in infinite systems, the critical window may feature a non-vanishing width.

Expectations of the multi-degree distribution
In a network with N colours, the degree distribution u(k), k = (k 1 , . . ., k N ) is the probability that a randomly selected node bears k i edges of colour i.The expected value of f (k) given degree distribution u(k), is defined by the following sum: It is convenient to gather all expected values for k i into the expected vector, If instead of selecting a node at random, one selects a node at the end of a randomly chosen, icoloured edge, the corresponding degree distribution is called i-biased and given by (k + e i ) u(k+ei) E[ki] .The expected vectors of the i-biased degree distributions are given by the columns of the following N × N matrix, and the covariance matrices of i-biased degree distributions are given by

Coefficients of the asymptote
Consider an axillary function where z is a vector of dimension Let z * = arg min f (z), z * > 0, |z * | ≤ 1, is the point where f (z) reaches its minimum value, then the size distribution (S7) features the following asymptote at the tail: where the coefficients are given by Here T i e j , and R i,j = µ 0 σ −1 z * adj(M − I)T i e j are square matrices of size N × N , and 24) is the Hessian matrix of f (z) having size N − 1 by N − 1.The standard basis vectors are denoted by (e i ) j = δ i,j .The derivation of these expressions is given in SI 3.

Colour-fraction configuration of a connected component of size n
Since the total number of edges is |v| = n − 1, one writes v 1 = 1 − N i=2 v i , and the configurations for the rest of edge types v 2 , v 3 , . . .v N obey the following law of large numbers: Let f i = vi n−1 denote colour fractions in components of size n.The configurations of the colourfraction vector, (f 1 , f 2 , . . ., f N ), are normally distributed with mean m = z * and variance matrix where a = The derivation of equations ( 25)-( 26) is given in SI 7.

(S6)
In a similar fashion to how Eq. (S1) was derived, Ref. [30] applies Joyal's theory of species to write the exact expression for the size distribution of weakly connected components: where the , is defined as The sum in Eq. (S8) runs over all partitions of the N -dimensional vector n into two non-negative terms, j and k, det * (D) refers to the determinant computed with the multiplication replaced by the convolution, and matrix D has the following elements 1) , i, j = 1, . . ., N. (S9)

SI6. DERIVATIONS FOR THE SIZE DISTRIBUTION ASYMPTOTE
This section derives the asymptote for the equation (S7).Note that in Eq. (S7), the arguments of the convolved functions and the convolution powers are both taken from the same vector k.In order to circumvent this issue, let us decouple the argument from the powers.We thus define To introduce the idea of the asymptotic analysis, we first start with a simplified version of (S7).
A. Asymptote for the product of convolution powers Let φ(ω) denotes the characteristic function of u(k), and φ i (ω) denote the characteristic functions of u i (k).Then, the characteristic function of G(k, k ) is given by the product however, as follows from the central limit theorem, φ ki i (ω) approaches a limit function when k i is large: where and µ i and T i are the expected values and covariance matrices as defined by (S5) and (S6).Applying similar procedure to the equation (S11), we obtain: where Inverse flourier transform of Eq. ( S13) is easy to obtain since ψ ∞ (ω, k) itself is the characteristic function of the multivariate Gaussian function: Consider the following matrix and the transformation it induces z = 1 n Ak.Since , which replaces the functions appearing in the exponent of (S14) with and where K = (M − I)A −1 , and are independent of n.The Taylor expansion of f (z) around the point where this function reaches its minimum value, gives: where the N − 1-by-N − 1 matrix H z * is the Hessian of f (z) at z = z * : S18) and R(z, z * ) is the expansion residue.Note that the gradient ∇f (z * ) = 0 since z * is a minimum of f (z), and by plugging these expansion into (S14) we obtain:

Now, by multiplying
and rewriting the sum of exponents as a product of exponential functions isolates a multivariate Gaussian function in the expression of F n (z): Let us now turn back to the asymptotic analysis of the equation (S10), which becomes written as The latter summation is performed over so that, on one hand, the sum form (S20) can be approximated with the integral On another hand, the first fraction in Eq. (S19) is a properly normalised multivariate Gaussian function with mean z * and shrinking variance 1 n H −1 z .This Gaussian function approaches the Dirac's delta δ(z − z * ) in the large n limit: its variance vanishes as O( 1 n ) while the mean remains constant.By combining these two observations together we obtain: where Note that by the definition, the residue vanishes at the expansion point: R(z * , z * ) = 0. Finally, by moving 2π and 1 n outside the determinants, one obtains where

B. Complete size distribution asymptote
The characteristic function for (S9) is given by: so that the characteristic function of the full expression appearing under the sum in (S7) is given by where The limiting functions for φ i (ω) are given in (S12), so that one can write which feature the following partial derivatives By replacing φ i (ω) and i N ki ∂ ∂ωj φ i (ω) with their limiting functions (S24),(S25) we obtain a chain of transformations: Determinant det[D ] can be now rewritten in the matrix form: where Let us expand determinant det[D ] into the sum over the set S N of all permutations of {1, 2, . . ., N }: Since the gradient of ψ ∞ (ω, k) is given by one can express the iωψ ∞ (ω, k) from the latter equation as: which is the characteristic function for Since series (S27) is the sum of powers of iωψ ∞ (ω, k), this series is the characteristic function for an identical expression in which iωψ ∞ (ω, k) is substituted by (S28): The latter expression is the product of F n (z) and a polynomial in an indeterminate y = 1 n : Instead of this polynomial, it is convenient to consider its collapsed form that we obtain by observing that the series from Eqs. (S27) and (S29) coincide under the substitution iω → x n (z), therefore performing the same substitution in Eq. (S26) leads to: where p(n) is independent of z * : and By taking into the account the asymptotical behaviour of F n (z) as given by (S23), the complete expression for the asymptote of (S7) reads: where and

SI7. NECESSARY AND SUFFICIENT CONDITIONS FOR THE SCALE-FREE BEHAVIOUR OF w(n)
In this section we show that C 2 = 0 if and only if According to the definition (S15), σ z is a linear combination of covariance matrices, and therefore σ −1 z is a positive definite matrix that features the Cholesky factorisation: Suppose C 2 = 0.Then, according to the definition of C 2 , f (z * ) = 0 and Recall that K = (M − I)A −1 and since det L = 0 we have So that v ∈ ker(M − I).

According to the definition of
We will now demonstrate that v > 0. Assume that z * / ∈ Ω ∞ , then for any z ∈ Ω ∞ , f (z) > 0, and therefore for arbitrary α > 0 and non-singular function γ(z) the following product vanish in the limit of large n, lim n→∞ n α Ωn γ(z)e − 1 2 nf (z) dz → 0, which contradicts our assumption that the asymptote is scale-free.One therefore concludes that z * ∈ Ω ∞ , so that the set of inequalities (S21) holds for z * .This inequalities imply: v 1 = 1 − |z * | > 0 and v i = z i−1 > 0, for i = 2, 3, . . ., N , and therefore v |v| = v > 0. This completes the proof of the forward implication.
Consider the reverse implication: suppose v ∈ ker(M −I), and v |v| > 0. Without loss of generality, assume |v| = 1.We will demonstrate that vector z * = (v 2 , v 3 , . . ., v N ) gives minimum to f (z), and f (z * ) = 0. On the one hand, we have a chain of transformations that is reverse to (S35): On the other hand, since f (z) ≥ 0 for z ≥ 0, z * must also be a minimum of f (z).Furthermore, this minimum belongs to Ω ∞ .Indeed, The latter inequalities imply that z * ∈ Ω ∞ , which guarantees validity of asymptote (S32), and since C 2 = f (z * ) = 0, this asymptote is a scale-free one.
and plugging these into Eqs.(S5)-(S6) allows us to express µ, M , T i as functions of p: By plugging M into Eq.(S33) one obtains the criterion for criticality: edge-coloured network features a critical behaviour at percolation probability vector 0 < p < 1 if and only if: In order to recover the manifold containing all critical points p numerically, one may follow this practical procedure: 1. for all the points from a discretised unit hypercube test if det(diag{p}M − P ) = 0; 2. for those points that pass the first test, find the eigenpair (v, 0) of diag{p}M − P and check if v |v| > 0.

SI10. DERIVATION OF COLOUR FRACTIONS IN CONNECTED COMPONENTS OF SIZE n
For the same reason as in the conventional configuration model, every connected component of finite size n has a tree-like structure, and therefore this component contains n − 1 edges.The generalised model operates with N types of edges and it is interesting to see how these n − 1 edges are partitioned among N edge types.It turns out that z * , as defined in (S16), plays an essential role in defining this partition.Let v i , i = 1, . . ., N denotes the number of edges of i th type in a component of size n.Since the total number of edges is |v| and the probabilities of configurations for the rest of edge types v 2 , v 3 , . . .v N follow form F n (z * ) as defined in Eq. (S19).In fact, since we are interested in the conditional probability given component size is n, it is enough to consider the first fraction appearing in (S19).So that the following law of large numbers holds: Evidently, the multivariate stochastic variable ( v2 n−1 , v3 n−1 , . . ., v N n−1 ) is normally distributed with mean z * and variance 1 n H z * , and the whole vector v n−1 is normally distributed with mean and covariance matrix where a =

SI11. DERIVATIONS FOR THE GIANT COMPONENT AND WEIGHT-AVERAGE COMPONENT SIZES
Let random variable n is the size of the component containing a randomly selected node.The generating function for n, W (x) = E[x n ], satisfies the system of N equations: These equations constitute a generalisation of the corresponding system introduced for uncoloured networks by Newman et al. [19].In Eq. (S45), W i (x) generate probabilities that a randomly selected node is connected to a component of size m on either of its sides.Since the network may contain infinitely large components, n and m are improper random variables.We formalise this facts by writing: Plugging these definitions into (S45) yields the following equations: and Here, e i are the standard basis vectors, and vector power is evaluated element-wisely: s ki i .The quantities g node and s have straightforward interpretations: g node , or the node size of the ginat component, is the probability that a randomly sampled node belongs to the giant component; s i are the probabilities that a randomly sampled i-coloured edge is not connected to a giant component on at least one side.So that the edge size of the giant component is written as the probability that a randomly chosen edge of colour i belongs to the giant component: i , i = 1, . . ., N.
By weighting these numbers with the total fractions of coloured edges c i = E[ki] / N i=1 E[ki], one obtains the vector of colour fractions in the giant component, v * = (v * 1 , v * 2 , . . ., v * N ), v * i = gici g c , which is the giant-component analog of (S43).Since the giant component is infinite by definition, v * is not a stochastic variable as is the case with colour fractions in finite components (S43).
We will now give a quantitive estimate for the size of average finite connected component.Formally speaking, we wish to extract some properties of random variable n , which is the size of randomly chosen component.It is important to note the subtle difference between n and n : they differ in the method of sampling.For the former, we first chose a node, and then take the size of the component the node belongs to, for the latter -we chose the component itself.
Given the tools at hand, it is not easy to derive the expression for average component size E[n ].With little efforts, however, we can derive one for the ratio In polymer literature, this ratio is commonly referred to as the weight-average size [38] and we will conveniently borrow this terminology.The weight-average size of finite connected components w avg is found by first expressing the derivatives of W i (x) and W (x) from equation (S47).Let y = (y 1 , y 2 , . . ., y N ) , y i = d dx W i (x)| x=1 , then See also Ref [16] for the derivation of the wight average size of connected components in the case of unicoloured networks.Finally, we show how expressions (S49) and (S48) change as a result of simple bond percolation.In this case, the degree distribution becomes dependant on percolation probability p as defined in Eq. (S36), which induces the following transformation of the giant component size: (S50)

FIG. 1 .
FIG. 1. Simple bond percolation in a sub-critical coloured network.a. Simulated size distributions of connected components (scatter plots) are compared with the theoretical asymptotes (solid lines).Inset: dependance of the exponent coefficient C2 on the percolation probability p. b.Fraction of coloured edges fi that belong to components of size n.Scatter plots: simulated networks (1-red, 2-blue, 3-yellow).Dashed lines: theoretical mean values.Solid lines: theoretical two-sigma standard deviation.

FIG. 2 .
FIG. 2.A coloured network with multiple phase transitions.a, Samples of network topologies corresponding to the degree distribution given by equation (19) with α = 0 and α = 0.1.b, The list of all secondary phase transitions and the colour subset indicators that correspond to them.c, The partial order of all secondary phase transitions in a network with three colours.An arrow points from an earlier-occurring phase transition to a later-occurring one.d, Dependance of the exponent coefficient C2 on the percolation parameter p for two cases of the degree distribution: α = 0 (solid line), and α = 0.1 (dashed line).e,f, Left axis: The total size of all giant components as a function of percolation parameter p as obtained from the theory (red solid line) and simulations (black dashed line).Right axis: The weight-average size of finite connected components as obtained form the theory (black solid line) and simulations (grey dashed line).These results are obtained for degree distributions with: e, α = 0, and f, α = 0.1.The vertical guide lines provide theoretical vales for phase transition points.

FIG. 3 .
FIG. 3. Colour switching in finite components.a, Theoretical expected faction of colours in finite components as a function of percolation probability p features a switching behaviour in the coloured network with α = 0.1.b, Comparison of colour fractions f1, f2, and f3 in components of size n = 50: a generated network (scatter plot) versus theoretical two-sigma standard deviation (shaded areas).The vertical guide lines provide theoretical vales for phase transition points.

FIG. 4 .
FIG. 4. Colour-dependant percolation in a network with three colours.a,b, Surfaces containing critical configurations of vector p in the network model with: a, α = 0 and b, α = 0.1.c,d, Manifolds containing p-configurations for which C2 < 0.02 are computed for: c, α = 0 and d. α = 0.1.The redblue colouring is used for better visual distinction between the surface's sides in a, and b, and manifold's interior/exterior in c, and d.

we have z 1
= 1 and therefore we may write z = 1 z , where z is a vector of dimension N − 1.Let us now set k = k and introduce the transformation of variables k = nA −1 1 z R i,j = µ 0 σ −1 z * adj(M − I)T i e j , are square matrices of size N ×N , and (e i ) j = δ i,j are the standard basis vectors.One can derive the coefficients of the expansion (S30) by applying differential operator 1 k! ∂ k ∂y k | y=0 to (S31).For k = 0, 1 this procedure yields, a 0 = det(M − I) det[Q] = 0, and a 1 = det(M − I)tr(adj(Q)R) > 0.