Driving interconnected networks to supercriticality

Networks in the real world do not exist as isolated entities, but they are often part of more complicated structures composed of many interconnected network layers. Recent studies have shown that such mutual dependence makes real networked systems potentially exposed to atypical structural and dynamical behaviors, and thus there is a urgent necessity to better understand the mechanisms at the basis of these anomalies. Previous research has mainly focused on the emergence of atypical properties in relation with the moments of the intra- and inter-layer degree distributions. In this paper, we show that an additional ingredient plays a fundamental role for the possible scenario that an interconnected network can face: the correlation between intra- and inter-layer degrees. For sufficiently high amounts of correlation, an interconnected network can be tuned, by varying the moments of the intra- and inter-layer degree distributions, in distinct topological and dynamical regimes. When instead the correlation between intra- and inter-layer degrees is lower than a critical value, the system enters in a supercricritical regime where dynamical and topological phases are not longer distinguishable.

The traditional study of networks as isolated entities has been recently overcome by a more realistic approach that accounts for interactions between networks [1].Networks in the real world are in fact often, if not always, mutually connected: social networks (e.g., Facebook, Twitter) are coupled because they share the same actors [2]; transportation networks are composed of different layers (e.g., buses, airplanes) with common nodes standing for the same geographic locations [3]; the functioning of communication and power grid systems depend one on the other [1].As the properties of an isolated network are not trivially deducible from those of the individual vertices that are part of it, at the same strength decomposing an interconnected network and studying each component in isolation does not allow to understand the whole system and its dynamics.Indeed, interconnected networks often exhibit properties that strongly differ form those typical of isolated networks: for example, structural [1] and dynamical [4] transitions, that are usually continuous in isolated networks [5], may become discontinuous in interconnected networks.So far, theoretical studies have pointed out that phase transitions in random interconnected networks become abrupt only if the strength (or density) of the interconnections is sufficiently large when compared to the first moment of the strength (or degree) distribution of the whole network [4,[6][7][8][9].In this paper, we show that the situation is not so simple, but depending on the combination of basic structural features -strength of interconnections, first two moments of the degree distribution of the entire network, and correlation between intra-and inter-layer degrees -an interconnected network explores a complicated scenario where its critical properties can drastically mutate.
In the following, we focus our attention on the case of * Electronic address: filiradi@indiana.edu.
two arbitrarily interconnected network layers, and study the spectral properties of its associated normalized laplacian [10].The choice of this matrix is not arbitrary.The normalized laplacian of a network is in fact an object of fundamental importance for the understanding of its structural and dynamical properties, sometimes even more than the adjacency matrix.For example, the spectrum of the normalized laplacian is used in spectral graph clustering to determine the internal organization of a graph [11,12], and many useful measures, such as graph energy [13], graph conductance and resistance [14], and the Randić index [15], are quantifiable in terms of the eigenvalues of the normalized laplacian.Also, the normalized laplacian fully describes the behavior of one of the most important dynamical processes studied in science, from biology [16] to computer science [17], from ecology [18] to finance [19] and physics [20]: random walk dynamics.The normalized laplacian is representative for both the classical [10] and the quantum [21] versions of random walk in networks, and more in general for every time reversible Markov chain [22].
Let us consider a symmetric and weighted interconnected network formed by two interconnected layers of identical size N whose adjacency matrix G is written in the block form (1) A, B and C are N × N square matrices.A and B are symmetric matrices, while C is not necessarily symmetric.According to Eq. ( 1), the adjacency matrix G is a convex linear combination of the intra-and inter-layer adjacency matrices G in and G out , and the parameter p ∈ [0, 1] serves to continuously tune the entire network from a bipartite graph (p = 0) to a perfectly decoupled networked system (p = 1).The normalized graph laplacian associated with the adjacency matrix of Eq. ( 1) is As the weight of the intra-layer layer connections p varies, the system explores different regimes clearly visible from the spectrum of L. A For p = 0.3, only the largest eigenvalue ν2N is well separated from the rest of the spectrum.The system is said to be in the "bipartite" regime (B-phase).B For p = 0.5, there are no detached eigenvalues.The system is in the "indistinguishable" regime (I-phase).C For p = 0.7, the second smallest eigenvalue ν2 is well separated for the rest of the spectrum.The system is in the "decoupled" regime (D-phase).D In the B-phase, the components of largest eigenvector w2N corresponding to a single network layer (blue squares) are linearly correlated with the square root of the respective node degrees √ k, while those of the second smallest eigenvector w2 (red circles) are not.E In the I-phase, there is no correlation between the components of w2N (w2) and √ k.F In the D-phase, the components of second smallest eigenvector w2 are linearly correlated with √ k, while the components of w2N are not.G The transition points p − c 0.42 and p + c 0.58 (delimiting the vertical gray area) between the different regimes correspond to the points in which our prediction [Eq.( 5), black dashed line] enters the continuous band of the spectrum [Eq.( 4), horizontal gray area].Red circles refer to ν2, while blue squares refer to ν2N .H The sums of the components of w2 (red circles) and w2N (blue squares) corresponding to nodes of a single network layer can be used as a order parameter to monitor the transition between the different regimes (u is the vector whose components are all equal to one).I As noted above, in the D-phase (B-phase), the components of w2 (w2N ) are perfectly correlated with the square root of the node degrees.

defined as
In the definition of L, 1 is the identity matrix, and D −1/2 is a diagonal matrix whose diagonal elements are equal to the inverse of the square root of the node strengths [10].Since G is symmetric, all the eigenvalues of L are real numbers in the range [0, 2].For simplicity, we sort them in ascending order such that 0 and we indicate with w i the eigenvector associated to ν i .The smallest eigenpair (ν 1 , w 1 ) is generally said to be "trivial" because it depends only on the strength sequence of the graph.We always have ν 1 = 0 and w T 1 = 1/ 2N s √ s 1 , . . ., √ s 2N , with s first moment of the strength distribution, and s i strength of the node i.When translated into the language of random walk dynamics, w 1 tells us that, independently of the initial conditions and the topology of the network, the stationary probability to find the random walker on a given node is linearly proportional to its strength.The other two external eigenvalues ν 2 and ν 2N , and their associated eigenvectors w 2 and w 2N , are much more meaningful from the structural and dynamical points of view.ν 2 is always strictly larger than zero in graphs composed of a single connected component.The associated eigenvector w 2 is typically used in graph clustering to determine the bisection corresponding to the minimum of the normalized cut of the graph [11].In terms of random walk dynamics, such eigenvector is also very meaningful because the signs of its components define the so-called almost invariant sets of the graph [12,23].In essence, a random walker spends most its time moving among nodes whose correspondent components in w 2 have the same sign, and less frequently jumps between nodes whose correspondent components in w 2 have different signs.The largest eigenvalue ν 2N = 2 only if the graph is bipartite, while ν 2N < 2 in all other cases, and the components of w 2N identify the bipartite components of the graph.A random walker typically jumps between pairs of nodes that correspond to components of w 2N with different signs.
To understand the spectral properties of L, we restrict our attention to the case of random network models generated according to the following procedure.Intraand inter-layer degrees of both network layers A and B are respectively given by the same identical intra-and inter-layer degree sequences The numbers appearing in the two degree sequences are also identical except for their order of appearance.The intra-and inter-layer degree distributions are thus identical P (k in ) = P (k out ) = P (k), while the alignment between the two sequences regulates the joint probability distribution P (k in , k out ) of both network layers.In the construction of the network, internal and external connections are randomly placed with the only constraints of preserving the a priori provided degree sequences and generating network layers composed of a single connected component, so that ν 2 > 0 and ν 2N < 2 for every p ∈ (0, 1).It is worth to briefly discuss the meaning, limitations and justifications of the previous assumptions.The fact that intraand inter-layer connections are randomly placed allows us to say that no substructure is present in layers A and B. Although this is not a good assumption for real networks that instead often exhibit modular structure at single layer level [24], it is a necessary assumption to concentrate our attention only on the effect of the layer structure on the spectrum of L. The fact that the intra-and inter-layer degree distributions are identical, and both network layers have the same degree distribution and joint probability function P (k in , k out ) are also non realistic assumptions.In real world interconnected networks, in fact, one should expect the various distributions to be different.Although our general solution requires a much weaker assumption (see SI), here we decided to impose these constraints for illustrative purposes.The scenario is in fact already so much rich that having an additional symmetry in the system makes easier the explanation of the various behaviors that the system can exhibits.
If we think in terms of random walk dynamics, our network models are symmetric in two ways.First, the stationary probability to find the random walker in one of the two layers is equal to 1/2.Second, they introduce a symmetry around the point p = 1/2 for the various regimes in which the system can be tuned by varying p. Intuitively, we expect the following: (i) For p < 1/2, the network should be in a "bipartite" regime (B-phase) with the random walker moving more often from one layer to the other, and less frequently diffusing in the same layer; (ii) For p > 1/2, the network should be in a "decoupled" regime (D-phase) with the random walker more likely diffusing between nodes of the same layer, and less likely jumping between layers; (iii) For p 1/2, the network should be in a "indistinguishable" regime (I-phase), with the random walker moving between and among layers with equal probability.Please note that, in the I-phase, the system is topologically and dynamically identical to a single-layer network, thus in a potentially safe state.On contrary, when layers are structurally and dynamically distinguishable (i.e., Band D-phases), then the network behaves as an interconnected network, and, as such, is exposed to potentially catastrophic failures or anomalous dynamical behavior.Intuitively, one should expect that the network is driven in a continuous way among the different phases as p is varied.In reality, the scenario is much more complicated and appealing.Let us first focus on the exactly solvable case in which the intra-and inter-layer degree sequences of both layers are ordered in the same way (k in i = k out i = k i , for all nodes i).Please note that in this case the node strengths equal the node degrees, and thus they do not change with p In addition to the trivial eigenpair, L has always another eigenpair (ν * , w * ) given by where c i = c, for i = 1, . . ., N , and c i = −c, for i = N + 1, . . ., 2N , and c = 1/ 2N k with k first moment of the degree distribution (see SI). Eq. ( 3) tells us that the eigenvector w * is able to perfectly reflect the layer structure of the system.Such eigenvector is asso-  ], the results of numerical computations (symbols) are in very good agreement with our analytic predictions (dashed lines).The critical point p + c is determined as the point for which the dashed lines enter the spectral band (gray area).For ξ = ξmin, our predictions fail to correctly describe how ν2 changes as a function of p. B As in Fig. 1, we monitor the transition between the D-and the I-phase, by looking at the order parameter u T w2/ N/2.As ξ decreases, the change of the order parameter becomes smoother, and the value of the critical point gets closer to p = 0.5.For ξ ∈ [ k 2 , k 2 ], the order parameter goes to zero if p < p + c .For ξ = ξmin instead, the order parameter is larger than zero even for p = 0.5.C The same behavior described for the order parameter is visible if we monitor the transition through the linear correlation coefficient between the eigenvector components and the r.h.s of Eq. ( 6).D Phase diagram representing the expected behavior of the system in the thermodynamic limit.The red line stands for the critical value p + c , while the blue line represents p − c [Eq. ( 9)].These critical values are plotted as functions of the quantity k 2 − ξ, ranging between zero and k 2 − ξmin.As a reference, the dashed line represents the point at which ξ = k 2 .Although the curves corresponding to p + c and p − c do not intersect, for values of the correlation term ξ sufficiently low, their analytic continuations suggest the disappearance of the I-phase, and the appearance of a supercritical regime where both the bipartite and the decoupled phases are present (i.e., BD-phase).
ciated to an eigenvalue that spans the entire range [0, 2], and therefore must intersect all the other eigenvalues of L. At each point of intersection, two orthogonal eigenvectors correspond to the same eigenvalue, and thus, in physical terms, each point of intersection corresponds to a degenerate energy level.Such phenomenology is typical of discontinuous phase transitions [25].Please note that this happens in networks of any size, and for any degree distribution.The changes between the various phases that we described above are of the same nature: as p is increased, the system moves discontinuously from the B-phase to the I-phase, and from the I-phase to the D-phase.The B-phase can be identified as the regime in which ν * corresponds to the largest eigenvalue of L. Similarly, the D-phase can be determined as the regime in which ν * represents the second smallest eigenvalue of L. Finally in the I-phase, ν * is inside the spectrum of L. The values of p for which phase transitions occur can be approximately estimated with the following argument.Numerical computations of the spectrum of L show that, whereas the eigenvalue ν * is rapidly changing as a function of p, all other eigenvalues vary much more slowly (see upper panels of Fig. 1).The graph obtained for p = 1/2 can be thus be used as a representative for this almost constant behavior.We can then approximate, thanks to the prediction of Chung et al. [26], the expected spectral radius of the normalized laplacian of a random graph with prescribed degree distribution as Eq. ( 4) generally provides very good estimates for the spectral radius of L for random networks with poissonian and power-law degree distributions.By comparing ν * to the r.h.s. of Eq. ( 4), we finally find that the two transition points are In Eq. ( 5), the solution with the minus sign corresponds to critical point of the transition between the Band Iphases, while the solution with the plus sign corresponds to threshold between the Iand D-phases.As expected, such predictions are in very good agreement with the results of numerical experiments (see Fig. 1).
Next, we relax the former assumption about perfectly aligned degree sequences, and we let k out correspond to a slightly modified permutation of the sequence k in .Unfortunately, our theory works only for Poisson degree distributions, but the results we obtain allows us to understand also the qualitative behavior of networks with power-law degree distributions.If degrees are random variates extracted from a Poisson distribution, we are able to numerically show that the i-th component of the second smallest (only in the D-phase) or the largest eigenvector (only in the B-phase) w * obeys the relation where c i = c, for i = 1, . . ., N , and c i = −c, for i = N + 1, . . ., 2N , and c is a proper normalization constant.This essentially means that the i-th component of the eigenvector w * is proportional to the difference between the internal and external strengths divided by the square root of the total strength of node i (Fig. S1).
The validity of Eq. ( 6) is testified by the very high linear correlation coefficients that one can measure between the components of the eigenvectors (numerically computed), and the r.h.s. of Eq. ( 6).Using Eq. ( 6) as ansatz for the solution of the eigenvalue problem, we are able to determine that the eigenvalue of L associated to the eigenvector w * of Eq. ( 6) is given by and thus by the ratio between a quadratic and a linear function of p.The function appearing at the denominator depends only on the second moment k 2 of the degree distribution.The function appearing at the numerator, instead, depends on the second moment of the degree distribution, and the correlation between intra-and interlayer degrees By comparing Eq. ( 7) with the expected spectral radius of Eq. ( 4), we finally arrive to the prediction of the critical points When the term inside the square root of Eq. ( 9) is larger than zero, we have an intersection between ν * and the rest of spectrum of L. This indicates the presence of an abrupt transition between the corresponding phases.Strictly speaking, these eigenvalues do not intersect in finite size random networks, as predicted by the Wignervon Neumann noncrossing rule [27,28], however, the abrupt nature of the transition appears clear already for networks of moderately large size.When the term in the square root of Eq. ( 9) is equal to zero, the non trivial eigenvalue ν * touches the rest of spectrum of L tangentially at p c .In this case, the corresponding transition becomes continuous.The non existence of a real solution in Eq. ( 9) reveals instead the absence of a crossing between ν * and the rest of the spectrum even in infinite size networks, and thus is the sign of a dramatic change in the behavior of the system.
Eqs. (7) and (9) show that the the entire scenario is characterized by two main factors.On one end, the degree distribution, with its first two moments, play a fundamental role for the determination of the various phases.
On the other hand, the degree distribution alone is not able to explain the behavior of the system.At parity of moments in fact, the features of the system can be still drastically changed by the correlation term defined in Eq. ( 8).To better understand how the various possibilities are regulated by the interplay between the correlation term ξ and the moments of the degree distribution, we vary the correlation term by regulating the alignment between the intra-and inter-layer degree sequences k in and k out .If the order of the sequences k in and k out is identical, ξ = k 2 , and Eqs. ( 7) and ( 9) correctly reduce to Eqs. ( 3) and ( 5), respectively.If from this initial alignment we randomly shuffle a selected portion of pairs of entries in the inter-layer degree sequence k out , we decrease the value of ξ.In particular, if we randomly mix all entries, then k out corresponds to a random permutation of k in and the correlation term reads ξ = k 2 .Finally, if the entries of k in are ordered in ascending (descending) order, while those of k out are sorted in descending (ascending) order, the correlation term reaches the minimum value ξ min [29].When ξ is in the range [ k 2 , k 2 ], the term in the square root of Eq. ( 9) is strictly larger than zero: all regimes exist, and the corresponding phase transitions are discontinuous (see Fig. 2).At the same time, the phase diagram suggests that a further decrement of ξ leads to a drastic mutation in the features of the system.For lower values of ξ, the analytic continuation of our predictions tell us: (i) phase transitions change their nature and become continuous; (ii) the Iphase is not longer present, and leaves space to an hybrid BD-phase where the system is simultaneously in both the bipartite and decoupled regimes.The presence of this hybrid phase appears evident by looking at the spectral properties of the normalized laplacian even of finite size systems: for ξ = ξ min and p = 0.5, the order parameters associated to the Band D-phases are in fact simultaneously larger than zero (see Figs. 2 and S2).
The previous effect is amplified when we consider networks whose degrees obey a power-law degree distribution P (k) ∼ k −γ , with γ < 3. Differently from the case of poissonian networks in fact, the correlation term ξ has a much larger range of variability.Although the ansatz of Eq. ( 6) breaks down for scale-free networks, and thus the entire analytic approach is not working, Eq. ( 9) still tells us what we should expect to see: unless the correlation term ξ is very close to its maximal value k 2 , the term in the square root is negative; this means that the BD-phase is reached already for little variations from the perfect alignment between the intra-and inter-layer sequences k in and k out .This is indeed verified in our numerical experiments as shown in Figs. 3 and S3.For p = 0.5, the presence of the hybrid BD-phase is testified by the fact that the eigenvalues ν 2 and ν 2N are simultaneously well separated from the other eigenvalues of L, and the associated eigenvectors w 2 and w 2N have components that are very well correlated with the r.h.s of Eq. ( 6).The mixing between the two degree Spectral analysis of the normalized laplacian L for an interconnected network with degree sequence identical to the one analyzed in Fig. 1.In this case, however, intra-and inter-layer degree sequences are aligned in such a way that the correlation term is ξ 1/2( k 2 + k 2 ).A For p = 0.5, the system is in the BDphase, where both the eigenvalues ν2 and ν2N are well separated from the other eigenvalues of L. B The components of both eigenvectors w2 and w2N are correlated with the r.h.s. of Eq. ( 6).C The eigenvalues ν2 (red circles) and ν2N (blue squares) touch the continuous band tangentially.D The order parameter clearly shows the presence of supercritical regime, where both the B-and D-phases are simultaneously present, for a wide range of values of p.Only for very low (large) values of p, the system is in the pure B-phase (D-phase).
sequences causes a mismatch between intra-and interlayer hubs which do not correspond to the same nodes as in the case of perfectly aligned degree sequences.From the structural point of view, the total disappearance of the so-called indistinguishable phase can be interpreted as the absence of a buffer of structural safety.Scalefree networks seem therefore to be constantly at risks of potentially catastrophic failures.There are, however, also dramatic dynamical consequences.While for high values of ξ there is a neat distinction between the Band D-phases and a random walker moves in the network performing either intra-or inter-layer diffusion, in the BD-phase instead both motions happen simultaneously (Fig. 4).Such ability makes the walker faster in the exploration of the entire system (Fig. S4).
To summarize, the spectral properties of the normalized laplacian of random interconnected networks show that these systems face a scenario that is much richer than the one valid for isolated networks.Depending on the combination of basic structural properties -degree distribu- tion, degree correlation, and strength difference between intra-and inter-layer connections -topological and dynamical phase transitions associated to the spectrum of the normalized laplacian can be either discontinuous or continuous, and different regimes can disappear or even coexist.This allows us, somehow, to divide the space of structural parameters in regions where the "classical" theory of isolated network holds, and those in which interdependent networks behave in anomalous manner.As a final remark, we would like to stress an appealing analogy between our findings and the typical thermodynamic behavior of any substance, e.g., water, carbon dioxide and methane.In normal conditions, a substance changes its phase from liquid to gas in a discontinuous manner.However, above its critical temperature and pressure, the substance ceases to exhibit distinct liquid-or gas-like state, and it becomes a so-called supercritical fluid [30].Since their properties can be fine-tuned, supercritical fluids have many industrial and scientific applications.We believe that the possibility to drive interconnected networks to supercritical regimes could also be used in a fruitful way to better design and control real world networked systems.

Normalized Laplacian
Consider now the normalized laplacian matrix derived from the adjacency matrix of Eq. (S1), given by In order to solve in a simpler way the eigenvalue problem it is better to consider the matrix where 1 is the identity matrix, and D −1 G is the so-called random walk matrix.Since For our purposes, it is even more convenient to rewrite the eigenvalue problem as follows Please note that previous equation has one trivial solution given by given by 1 − ν = 0 and |v = |1 , thus any other eigenvector with different eigenvalue must be orthogonal to |1 .
If we write the eigenvector |v = |v A , v B as the composition of two vectors |v A and |v B with N entries each corresponding to the nodes of layers A and B, respectively, the previous equation is equivalent to the coupled equations If we multiply these equations for 1|, we have In-and out-degrees sequences are randomly mixed so that the correlation term reads ξ = k 2 .The description of the various panels is identical to the one of Fig. 1 of the main text.In panels D, E and F, the components of the eigenvectors w2 and w2N are plotted as functions of the components of the vectors appearing on the r.h.s. of Eq. (S6).In Panel G, the dashed line stands for Eq.(S9).
If we take their difference, we have Please note that the graphical closure of the graph imposes that When the internal strength of the both layers is comparable, i.e., then numerical computations of the spectrum of L, such as those presented in Fig. S1, show that, the second smallest eigenvector (largest eigenvector) in the decoupled regime (bipartite regime) of the normalized laplacian is approximately given by where n A and n B are normalization constants.This means that Since |v A , v B is orthogonal to |1, 1 , we can write So in conclusion, the eigenvalue corresponding to this eigenvector is given by Please note that where m in 2,A is the second moment of the in-strength distribution of layer A, m out 2,A is the second moment of the out-strength distribution of layer A, and ξ A is the correlation between the in-and out-strength of nodes in layer A.

Similarly, one can write ∆s
where m in 1,A and m out 1,A are the first moments of the in-and out-strength distribution, respectively.The same relations are also valid for layer B. Eq. (S8) can be thus rewritten as .
If we finally assume that both network layers have identical in-and out-strength distributions and identical correlations between in-and out-strengths, the former expression further simplifies to Please note that if we consider the adjacency matrix of Eq. 1 of the main text, we have to perform the following substitutions: We thus recover Eq. ( 7) of the main text.

Critical points
The critical points of the transitions are calculated as the values of p for which Eq. ( 7) intersects the spectral band of Eq. ( 4) predicted by Chung et al.For p = 1/2, this leads to the equation and finally Such quadratic equation has discriminant given by and thus solution given by The solutions of this equations closer to p = 1/2 are not not interesting for us because ν * has already entered the spectrum at point.The effective critical points are thus given by Eq. ( 9) of the main text.

Positively/neutrally correlated network layers
A simple strategy to explore this regime of correlation is the following.We sort both the intra-layer and inter-layer degree sequences in ascending (descending) order, but we then shuffle, with probability q, each entry of the sequence with another randomly chosen entry.For q = 0, the two sequences are identical and the correlation term reads ξ = k 2 .In this case, we say that network layers are positively correlated.For q = 1, the out-degree sequence is effectively a random permutation of k in , and the correlation term reads ξ = k 2 .We thus say that network layers are neutrally correlated.For general values of q, the correlation term obeys The discriminant of Eq. (S10) becomes , and thus The value q c of q for which ∆ = 0 is therefore the solution of a quadratic equation in q.This is given by from which we finally find Please note that , thus the only solution that is potentially smaller or equal to one For poissonian graphs, the former expression shows that q c > 1.Thus for any value of q ∈ [0, 1] ∆ > 0. For scale-free graphs with degree exponent γ < 3 instead, 1, and q c ≤ 1.

Negatively/neutrally correlated network layers
To explore this regime of correlation, we sort both the intra-layer in ascending (descending) order and inter-layer degree sequences in descending (ascending) order, but we then shuffle, with probability t, each entry of the sequence with another randomly chosen entry.For t = 0, the correlation term reads ξ = ξ min and the layers are negatively correlated.For t = 1, the out-degree sequence is effectively a random permutation of k in , and the correlation term reads ξ = k 2 .Network layers are neutrally correlated.For general values of t, the correlation term obeys Thus the critical value of t is given by The solution with the plus sign is always negative, thus As explained in the text, the continuous parameters q and t are used to control the alignment between the entries of the two degree sequences sequences, so that the correlation term reads ξ = (1 − q) k 2 + q k 2 or ξ = (1 − t)ξmin + t k 2 .A For any value of q, the largest eigenvalue ν2N intersects the continuous band (gray area).The corresponding transition, in the thermodynamic limit, is expected to be discontinuous.The results of numerical computations (symbols) are in very good agreement with our analytic predictions (dashed lines).For t < 1, our predictions fail to correctly describe how ν2N changes as a function of the weight of the intra-layer connections.B Same as for panel A, but for the second smallest eigenvalue ν2.C Sum of the components of w2N corresponding to nodes of a single layer.D Same as for panel C, but for the second smallest eigenvector w2.

Positively correlated layers
This is a special case in which we can determine an additional eigenpair of L.
This means that the vector |1, −1 is eigenvector of K −1 G with eigenvalue (2p − 1).Thus the vector K 1/2 |1, −1 is A When the weight of the intra-layer connections p is sufficiently low, the system is in a pure bipartite regime.B After that, the system enters in a hybrid phase, where both eigenvalues ν2 and ν2N are simultaneously outside the continuous band.
C For large values of p, ν2N gets absorbed in the continuous band, and only ν2 is detached from the rest of the spectrum.D In the pure bipartite regime, the components w2N (blue squares) of a single layer are correlated with the r.h.s. of Eq. (S6), while those of w2 (red circles) are not.E In the hybrid regime, the components of both eigenvectors w2 and w2N are correlated with the r.h.s. of Eq. (S6).F In the pure decoupled regime, only the components of w2 are correlated with the r.h.s. of Eq. (S6).G The eigenvalues ν2 (red circles) and ν2N (blue squares) touch the continuous band tangentially.This is an indication of continuous transitions between the various regimes.H The sum of the eigenvector components of a single network layer clearly show the coexistence of the two dynamical phases.I The same is also visible if we monitor the state of the system with the correlation coefficient between the components of the eigenvectors and the r.h.s. of Eq. (S6).
eigenvector of L with eigenvalue 2(1 − p).Please note that this eigenvector should be normalized such that and thus the normalization reads n = 1 2 1|k .
In conclusion L has always associated the eigenpair

Autocorrelation function of random walk dynamics
Indicate with |π the vector whose i-th component represents the stationary distribution to find the random walker on node i, i.e., |π i = |s i 1|s .
Suppose now the walker starts its walk at position j.The probability to find the walker at a given position after t steps is given contained in the vector where we indicated with |j the vector whose components are all equal to zero except for its j-th component that equals one, and with D −1 G t the t-th power of the transition matrix D −1 G.
Suppose we monitor the position of the random walker at time τ throught the variable x(τ ).We define x(τ ) = +1 if the walker is in one of the nodes of layer A at stage τ , and x(τ ) = −1, otherwise.The average value of x(τ |j) given that the initial position is on node j is thus given by where the vector |1, −1 is a vector whose first N components are all equal to 1, and its second N components are equal to −1.Since, at stationarity, the proability to start at position j is equal to |π j , we can write

7 Figure 1 :
Figure1: Abrupt transition in interconnected networks.Spectral analysis of the normalized laplacian L of an interconnected network composed of two network layers of size N = 1024.The degree sequence is composed of integers extracted from a power-law probability distribution with exponent γ = 2.5 defined on the support [32, N ].Intra-and inter-layer degree sequences are perfectly aligned so that the correlation term reads ξ = k in ,k out k in k out P (k in , k out ) = k 2 .As the weight of the intra-layer layer connections p varies, the system explores different regimes clearly visible from the spectrum of L. A For p = 0.3, only the largest eigenvalue ν2N is well separated from the rest of the spectrum.The system is said to be in the "bipartite" regime (B-phase).B For p = 0.5, there are no detached eigenvalues.The system is in the "indistinguishable" regime (I-phase).C For p = 0.7, the second smallest eigenvalue ν2 is well separated for the rest of the spectrum.The system is in the "decoupled" regime (D-phase).D In the B-phase, the components of largest eigenvector w2N corresponding to a single network layer (blue squares) are linearly correlated with the square root of the respective node degrees √ k, while those of the second smallest eigenvector w2 (red circles) are not.E In the I-phase, there is no correlation between the components of w2N (w2) and √ k.F In the D-phase, the components of second smallest eigenvector w2 are linearly correlated with √ k, while the components of w2N are not.G The transition points p −

Figure 2 :
Figure 2: Phase diagram.We consider interconnected networks composed of two layers of size N = 8192 with intraand inter-layer degrees obeying a Poisson distribution with average k = 32.A Second smallest eigenvalue ν2 as a function of the weight of the intra-layer connections p.For values of the correlation term ξ in the range [ k 2 , k 2 ], the results of numerical computations (symbols) are in very good agreement with our analytic predictions (dashed lines).The critical point p + c is determined as the point for which the dashed lines enter the spectral band (gray area).For ξ = ξmin, our predictions fail to correctly describe how ν2 changes as a function of p. B As in Fig.1, we monitor the transition between the D-and the I-phase, by looking at the order parameter u T w2/ N/2.As ξ decreases, the change of the order parameter becomes smoother, and the value of the critical point gets closer to p = 0.5.For ξ ∈ [ k 2 , k 2 ], the order parameter goes to zero if p < p + c .For ξ = ξmin instead, the order parameter is larger than zero even for p = 0.5.C The same behavior described for the order parameter is visible if we monitor the transition through the linear correlation coefficient between the eigenvector components and the r.h.s of Eq. (6).D Phase diagram representing the expected behavior of the system in the thermodynamic limit.The red line stands for the critical value p + c , while the blue line represents p − c [Eq. (9)].These critical values are plotted as functions of the quantity k 2 − ξ, ranging between zero and k 2 − ξmin.As a reference, the dashed line represents the point at which ξ = k 2 .Although the curves corresponding to p + c and p − c do not intersect, for values of the correlation term ξ sufficiently low, their analytic continuations suggest the disappearance of the I-phase, and the appearance of a supercritical regime where both the bipartite and the decoupled phases are present (i.e., BD-phase).

[ 5 BFigure 3 :
Figure3: Supercritical regime.Spectral analysis of the normalized laplacian L for an interconnected network with degree sequence identical to the one analyzed in Fig.1.In this case, however, intra-and inter-layer degree sequences are aligned in such a way that the correlation term is ξ 1/2( k 2 + k 2 ).A For p = 0.5, the system is in the BDphase, where both the eigenvalues ν2 and ν2N are well separated from the other eigenvalues of L. B The components of both eigenvectors w2 and w2N are correlated with the r.h.s. of Eq. (6).C The eigenvalues ν2 (red circles) and ν2N (blue squares) touch the continuous band tangentially.D The order parameter clearly shows the presence of supercritical regime, where both the B-and D-phases are simultaneously present, for a wide range of values of p.Only for very low (large) values of p, the system is in the pure B-phase (D-phase).

ξ 5 -Figure 4 :
Figure4: Random walk dynamics.We consider two interconnected networks composed of N = 8192 nodes.Degrees are extracted from a power-law distribution with exponent γ = 2.5 defined on the support [32, N ].We monitor the trajectory of a random walker moving on these systems for different values of the correlation term ξ and different values of the weight of the intra-layer connections p.A Autocorrelation function R(τ ) = x(t)x(t + τ ) as a function of the time delay τ .In the definition of R(τ ), x(t) = +1 if the random walker is on layer A at time t, and x(t) = −1, otherwise.B Spectral density S(ω) of the autocorrelation functions of panel A. Intra-layer movements are characterized by the presence of a peak of S(ω) at ω = 0, while inter-layer jumps are emphasized by the presence of peaks at ω = ±1.
) D A and D B are two diagonal matrices such that the i-th diagonal elements are (D A ) i,i = (|s A ) i (i.e., the i-th component of the vector |s A ) and (D B ) i,i = (|s B ) i (i.e., the i-th component of the vector |s B ).
L and 1 − D −1 G are related by a similarity transformation, they share the same spectrum.The right eigenvector |w of L is related to the left eigenvector of |u of 1 − D −1 G by |w = D −1/2 |u , while the right eigenvector |w of L is related to the right eigenvector of |v of 1 − D −1 G by |w = D 1/2 |v , and consequently the left eigenvector |u of 1 − D −1 G and the right eigenvector of |v of 1 − D −1 G are related by |v = D −1 |u .

[[[ 7 Figure S1 :
Figure S1: Spectral analysis of the normalized laplacian L of an interconnected network composed of two network layers of size N = 1024 with intra-and inter-layer degrees obeying a Poisson distribution with average degree k = 32.In-and out-degrees sequences are randomly mixed so that the correlation term reads ξ = k 2 .The description of the various panels is identical to the one of Fig.1of the main text.In panels D, E and F, the components of the eigenvectors w2 and w2N are plotted as functions of the components of the vectors appearing on the r.h.s. of Eq. (S6).In Panel G, the dashed line stands for Eq.(S9).

2 CFigure S2 :
FigureS2: We consider interconnected networks composed of two network layers of size N = 8192 with intra-and inter-layer degrees obeying a Poisson distribution with average degree k = 32.As explained in the text, the continuous parameters q and t are used to control the alignment between the entries of the two degree sequences sequences, so that the correlation term reads ξ = (1 − q) k 2 + q k 2 or ξ = (1 − t)ξmin + t k 2 .A For any value of q, the largest eigenvalue ν2N intersects the continuous band (gray area).The corresponding transition, in the thermodynamic limit, is expected to be discontinuous.The results of numerical computations (symbols) are in very good agreement with our analytic predictions (dashed lines).For t < 1, our predictions fail to correctly describe how ν2N changes as a function of the weight of the intra-layer connections.B Same as for panel A, but for the second smallest eigenvalue ν2.C Sum of the components of w2N corresponding to nodes of a single layer.D Same as for panel C, but for the second smallest eigenvector w2.

[[[ 9 FFigure S3 :
Figure S3: Spectral analysis of the normalized laplacian L the scale-free interconnected network analyzed in Fig.3of the main text.A When the weight of the intra-layer connections p is sufficiently low, the system is in a pure bipartite regime.B After that, the system enters in a hybrid phase, where both eigenvalues ν2 and ν2N are simultaneously outside the continuous band.C For large values of p, ν2N gets absorbed in the continuous band, and only ν2 is detached from the rest of the spectrum.D In the pure bipartite regime, the components w2N (blue squares) of a single layer are correlated with the r.h.s. of Eq. (S6), while those of w2 (red circles) are not.E In the hybrid regime, the components of both eigenvectors w2 and w2N are correlated with the r.h.s. of Eq. (S6).F In the pure decoupled regime, only the components of w2 are correlated with the r.h.s. of Eq. (S6).G The eigenvalues ν2 (red circles) and ν2N (blue squares) touch the continuous band tangentially.This is an indication of continuous transitions between the various regimes.H The sum of the eigenvector components of a single network layer clearly show the coexistence of the two dynamical phases.I The same is also visible if we monitor the state of the system with the correlation coefficient between the components of the eigenvectors and the r.h.s. of Eq. (S6).

2 8 2 9 2 10 2 11 2 12 2 13 2 14 2 15 2 16 sizeξ = k 2 ξ = k 2 ξ
Figure S4: Cover time Tcover of a random walker moving in an interconnected network composed of two layers of N nodes.The degrees are extracted from a power-law distribution with exponent γ = 2.5 defined on the support [5, N ].Cover time is divided by the sum of the strengths of all nodes in the network.The weight p of intra-layer connections is set equal to 0.5.Each point represents the average value of the normalized cover time obtained in at least 100 realizations of the network models.The dashed lines stand for fits with the function Tcover = a + b log N .We measure b = 0.20 (red), b = 0.19 (blue) and b = 0.15 (green).