Spectral properties of the Laplacian of multiplex networks

One of the more challenging tasks in the understanding of dynamical properties of models on top of complex networks is to capture the precise role of multiplex topologies. In a recent paper, Gomez et al. [Phys. Rev. Lett. 101, 028701 (2013)] proposed a framework for the study of diffusion processes in such networks. Here, we extend the previous framework to deal with general configurations in several layers of networks, and analyze the behavior of the spectrum of the Laplacian of the full multiplex. We derive an interesting decoupling of the problem that allow us to unravel the role played by the interconnections of the multiplex in the dynamical processes on top of them. Capitalizing on this decoupling we perform an asymptotic analysis that allow us to derive analytical expressions for the full spectrum of eigenvalues. This spectrum is used to gain insight into physical phenomena on top of multiplex, specifically, diffusion processes and synchronizability.


I. INTRODUCTION
One of the major lines of research in complex networks has focused in the comprehension of the relationship between network topologies and the behavior of processes occurring on them. The general approach to model a network is not general enough to ascertain the true interdependence that arises in some complex systems. Examples of such a systems are: user relationships in different social networks [1], transportation systems [2], or the learning organization in the brain [3]. Particularly interesting are those topologies of interconnected networks called multiplex, see Fig. 1, where each object is univocally represented in each independent layer and so the interconnectivity pattern among layers become one-toone [4][5][6][7][8][9][10], allowing the simultaneous study of different interconnected patterns between the same objects.
The behavior of any linearized dynamical process on a complex system is related to the Laplacian matrix of the underlying network and particularly to its second smallest eigenvalue. This eigenvalue, also called algebraic connectivity λ 2 , turns out to be essential to understand, for example, the time required to synchronize phase oscillators [11], or to converge to the maximum entropy state in a diffusion process [12]. Moreover, the largest eigenvalue of the Laplacian matrix plays a determinant role in the assessment of the stability of the synchronization manifold in networks of coupled oscillators [13,14]. In this article, we rely on the particularities of the multiplex networks to model its Laplacian matrix in terms of a decomposition between intra-and interlayer structure. This decomposition allows us to characterize the spectrum of the Laplacian, using perturbation theory, and hence the behavior of several dynamic processes. In particular, we are able to assess the diffusion time scales in any multiplex structure, and we can also infer the optimal value of the synchronization ratio in terms of the master stability function.
The paper is structured as follows: in Sec. II, we present the structural decomposition of the Laplacian of the multiplex (from now on supra-Laplacian) into intralayer and interlayer networks contributions; in Sec. III, we analyze the role of the interlayer network; Sec. IV is devoted to the perturbative analysis of the eigenvectors of the supra-Laplacian; in Sec. V, we expose the implications of the findings in terms of the physics of dynamical processes in multiplex networks, and finally we state the conclusions.

II. MULTIPLEX SUPRA-LAPLACIAN MATRIX
Let us consider a multiplex network consisting of M layers and N nodes per layer. The intralayer connectivity of layer α is expressed as an adjacency or strengths matrix W (α) ∈ R N ×N whose corresponding Laplacian is is the diagonal matrix of the nodes' intralayer strengths. In multiplex networks it is supposed that the interlayer connectivity is identical for all nodes [12], thus we may define the interlayer network W I ∈ R M×M whose components represent the strength of the connection between every pair of layers, and the associated interlayer Laplacian is L I = S I − W I . For the sake of simplicity, we will assume that the interlayer and intralayer networks are undirected, globally connected and without self-loops.
The supra-Laplacian L of the whole multiplex [12] may be separated in two contributions, where L L stands for the supra-Laplacian of the independent layers and L I for the interlayer supra-Laplacian. The first one is just the direct sum of the intralayer Laplacians, Figure 1: Sketch of a multiplex networks of four layers. Note that the same nodes are represented at each layer, although with different connectivity in each one of them. The interlayer connections corresponds to the links between the different categorical layers each node has. In the upper panel we observe the elements that will end in the direct sum in Eq.2. In the middle panel we show the structure corresponding to Eq.3. Finally, in the lower panel, we show the structure that will correspond to the network of layers and its laplacian L I .
while the interlayer supra-Laplacian may be expressed as the Kronecker (or tensorial) product of the interlayer Laplacian and the N × N identity matrix I, The decomposition of the supra-Laplacian given in Eq. (1) is fundamental for the discovery of several spectral properties of the multiplex, which will be uncovered in the following sections.

III. THE ROLE OF THE INTERLAYER NETWORK
It is well-known that the spectrum of the Kronecker product of two matrices is formed by the products of the eigenvalues of the individual matrices, and the associated eigenvectors are obtained by the Kronecker products of the eigenvectors. We can apply this property to Eq. (3), thus the eigenvalues of L I are equal to the eigenvalues of the interlayer network Laplacian L I , since the identity matrix I only has eigenvalues equal to one. Moreover, the multiplicity of each eigenvalue of L I is just N times its multiplicity in L I . Let x I ∈ R M be any of the eigenvectors of L I , with eigenvalue λ. The vector x ≡ x I ⊗ 1, being 1 ≡ (1, . . . , 1) T ∈ R N , is an eigenvector of L I with eigenvalue λ, as explained above, but it is also an eigenvector of L L with zero eigenvalue: Therefore, i.e. x is an eigenvector of the full supra-Laplacian. Consequently, denoting by Λ(L I ) = {λ I 1 = 0, λ I 2 , . . . , λ I M } the spectrum of the interlayer Laplacian, with the eigenvalues sorted in ascending order, λ I 1 ≤ λ I 2 ≤ · · · ≤ λ I M , we have shown that, Namely, all eigenvalues of the interlayer Laplacian are also eigenvalues of the supra-Laplacian. In general, these eigenvalues cannot be easily calculated and numeric computation is needed. In the following subsections we analyze several particular cases in which they can be easily computed.

A. Solvable case: small number of layers
When the number of layers is small it is possible to derive analytical expressions for the eigenvalues of L I . For example, if there are two layers, the interlayer network is characterized by, with spectrum Λ( where D x is a parameter that controls the relative strength between the interlayer and the intralayer contributions, and then the non-zero eigenvalues λ I 2 and λ I 3 are given by where D aβ (a, β = 1, 2, 3) denotes the weight of the links between the layers a and β. With four and five layers this kind of analytical expressions also exist (e.g. based on Cardano's method to solve cubic and quartic equations), but they are too involved to be useful. And for six or more layers, no general algebraic expressions exist due to the limitations imposed by Galois theorem.

B. Solvable case: Uniform interlayer weights
Another solvable case is the one in which the interlayer network is fully connected (without self-loops) and with all the weights equal to the same value, If the interlayer network is uniform but not fully connected, all the eigenvalues are proportional to the common weight D x , but the coefficients depend on the precise chosen topology. Many of them are also solvable, e.g. regular graphs, cycles, paths, etc.

IV. WEAK AND STRONG INTERLAYER INTERACTIONS: PERTURBATIVE ANALYSIS
In the previous section we have shown that the eigenvalues of the interlayer network are also eigenvalues of the full supra-Laplacian. Here we show, using perturbation theory, that these eigenvalues also impose important restrictions on the shape of the whole spectrum of the supra-Laplacian. These characteristics are specially useful to understand physical phenomena on multiplex networks.
To account for the relative strength between the interand intralayer connections, we denote by D x the quotient between their maximum values, thus we can write, where and the largest components of L L andL I are of the same order of magnitude. From now on we will use hats above objects where a factor D x has been extracted. According to Eqs. (6) and (10), the spectrum of the supra-Laplacian L contains the eigenvalues whereλ I α are the eigenvalues ofL I . On the other hand, the eigenvalues of the intralayer supra-Laplacian L L in Eq. (2) are the union of the eigenvalues of the intralayer networks,

A. Weak interlayer networks
For small values of the interlayer network, D x ≪ 1, the smallest eigenvalues are those given in Eq. (12), while the largest ones are perturbations of the non-zero intralayer eigenvalues in Eq. (13). For the calculation of the shape of these perturbations, let us select any of the eigenvectors x (γ) ∈ R N of the γ-th intralayer Laplacian, The corresponding eigenvector of the intralayer supra- where e γ ∈ R M denotes the canonical vector with a unity in the γ-th component and zeros elsewhere. By substituting the first order perturbation, and Eq. (10) into the eigenvalue equation, we recover Eq. (15) at 0-th order, and, at O(D x ). Multiplying both sides by the left with the transpose of v (γ) , and taking into account the symmetry of the supra-Laplacians and Eq. (15), we get wherel I αβ denotes the (α, β)-component of the interlayer LaplacianL I . Hence, the first order perturbation of the γ-th intralayer Laplacian eigenvalues at D x ≪ 1 are given by wherel I γγ =ŝ I γ represents the strength of node γ in the interlayer networkŴ I . This means that, for small values of D x , all the non-zero eigenvalues of the intralayer networks are shifted by a valueŝ I γ D x , which only depends on the layers. For the zero eigenvalues no perturbation analysis is needed since we know the exact values, given in Eq. (12).

B. Strong interlayer networks
The analysis in the limit of strong interlayer interactions, D x ≫ 1, is analogous to the weak interactions case. Defining ǫ = 1/D x , we can write the supra-Laplacian as, and consider the behavior at ǫ ≪ 1. The eigenvectors and eigenvalues to be perturbed are those from the interlayer supra-LaplacianL. Calling x I andλ I any eigenvector and eigenvalue pair ofL I , and using Eq. (3), we realize that, for any vector u ∈ R N , Therefore, we apply the following perturbation, leads, at O(ǫ), to Its α block row is equal to Now we multiply by x I α , sum over α β α and use Eq. (24) to obtain Among the sorted eigenvalues Λ(L I ) = {λ I 1 , . . . ,λ I M } it is convenient to studyλ I 1 = 0 separately from the rest of the non-zero ones. In this case, we know that x I = 1, and Eq. (32) reduces to This means that, for D x ≫ 1, the associated eigenvalues λ of the supra-Laplacian L, are the eigenvalues of the Laplacian L AV of the average network This result was introduced in [12] only for the particular case of the two-layer multiplex, M = 2. For the non-zero eigenvaluesλ I = 0 therefore the corresponding eigenvalues of the supra-Laplacian diverge linearly with D x , with shifts given by the eigenvalues of Eq. (32). One particular shift isλ ′ = 0 when u = 1, thus recovering the previously found exact eigenvalues given in Eq. (12).

C. Global structure of the supra-Laplacian spectrum
The global picture of the supra-Laplacian spectrum, which can be deduced from the previous subsections, may be summarized in the following list: This structure is completely general provided the interlayer and intralayer networks are undirected, connected and without self-loops. The changes produced by breaking the last two conditions are not important, but for directed networks complex eigenvalues may appear, thus modifying significantly the properties of the supra-Laplacian spectrum. Figure 2 illustrates the accuracy of the approximation for weak and strong interlayer networks. The figure shows a subset of the eigenvalues, together with the proposed approximation, of a multiplex network of tree layers. In each layer we have a different toy network of five nodes.
In the following section we will present the exploitation of these results in the context of physical processes running on top of general multiplex networks.

A. On the timescale of diffusion dynamics
Equivalently to [12], we consider a particular diffusion dynamics where nodes are linearly coupled between them. Under the context of a multiplex network, as commented in the introduction, this coupling is differentiated between intralayer and interlayer. We consider that the interlayer coupling constant is the same for all nodes between two layers. Under this setting, it is possible to model the evolution of the states of node x iα (node i at layer α) with the following differential equation: i,j is the weight of the connectivity between nodes x iα and x jα and ℓ αβ is the interlayer coupling constant.
The diffusion equation defined by Eq. (37) can be represented in matrix form by: where L L and L I = D xL I represent the intralayer and interlayer supra-Laplacians and D x represents the quotient defined in Sec. IV. The solution of this equation in terms of normal modes is given by φ i (t) = φ(0)e −λit , where λ i are the eigenvalues of the supra-Laplacian matrix L. We can see that the eigenvalue that governs the convergence of the process is given by the second eigenvalue. Concluding that the time scale for the diffusion process becomes τ ∝ λ −1 2 . For small couplings between the different layers (i.e. D x ≪ 1) λ 2 (L) will correspond to the second eigenvalue of the interlayer network, L I . In that case, τ ∝ (λ I 2 D x ) −1 (see Sec. IV A).
For large interlayer coupling (i.e. D x ≫ 1), λ 2 (L) can be approximated by the second eigenvalue of the average network defined in Eq. (35). Thus for large D x values τ ∝ (λ 2 (W AV )) −1 (see Sec. IV B). Figure 3 shows the wellness of the proposed approximation for small and large D x coefficients for a multiplex of four layers. See Supplemental Material for more examples.

B. On the stability of the synchronization manifold of coupled phase oscillators
Now, we focus on another ubiquitous physical process, the synchronization of phase oscillators in networks [14]. Let us assume that the phase oscillators are embedded in the nodes of a multiplex network structure, and that the phases (states) are different at different layers, meaning that we have a multi-phase oscillator. Extending previous results on the synchronizability (or strictly speaking, the stability of the synchronization manifold) to multiplex and taking advantage of the Master Stability Function [13,15], we reduce the problem of assessing synchronizability to that of computing the eigenratio R = λ N /λ 2 , in the multiplex network, where λ N is the largest eigenvalue. In particular, we analyze the asymp-  totic behavior of R for the cases of weak interlayer coupling values (D x ≪ 1) and strong interlayer coupling values (D x ≫ 1).
For weak interlayer coupling values, as Sec. IV A states, the second eigenvalue of the supra-Laplacian becomes the second eigenvalue of the network of layers, i.e. λ 2 (L) ≈ λ 2 (L I ). The largest eigenvalue in this case corresponds to the largest eigenvalue of the Laplacians of the different layers, i.e. λ max (L) = max α (λ max (L (α) ) +ŝ I α D x ). Thus, the eigenratio can be approximated as, For strong interlayer coupling the second eigenvalue can be approximated by the second eigenvalue of the average network, that is L AV (see Sec. IV B). The largest eigenvalue, in this case, is a function of the largest eigenvalue of the network of layers, λ max (L I ), with an additional offset. This offset depends on the individual layers, L (α) , weighted by the entries of the eigenvector x ′ I corresponding to the maximum eigenvalue of the network of layers. Thus, the largest eigenvalue is where the laplacian L WA max is given by Thus, the eigenratio for large interlayer coupling can be approximated as, To illustrate the proposed approximation, we computed the eigenratio and the proposed approximation for a multiplex of four layers. Figure 4 shows the obtained results. See Supplemental material for additional plots with different multiplex topologies.

VI. CONCLUSIONS
We have presented the asymptotic analysis of the spectrum of the Laplacian of multiplex networks. We found analytical expressions that allow us to infer the behavior of dynamical processes, such as diffusion or synchronization, on top of multiplex networks. The findings reveal physical implications of the multiplex structure that have no counterpart in monoplex (one layer) networks. For example, in diffusive processes we find a super-diffusive behavior where the time scales of diffusion associated to the multiplex are shorter than in any particular individual layer network. In the analysis of the synchronizabity ratio R, we have found the existence of an optimal value of the coupling D x for which the synchronization of the full structure is the most stable. The applicability of the mathematical findings on the features of the Laplacian of multiplex networks sure go beyond this particular cases, and eventually can help to understand any process whose behavior is linked to the spectral properties of the Laplacian or other akin matrices.  Plot of the eigen-ratio and the proposed approximation for a multiplex of 3 layers. Each layers contains a network with 4 communities, each community corresponds to an Erdös-Rényi network with edge probability 0.5, and the inter-community edge probability is 0,1. The communities between different layers strongly overlap. Comparison between the second eigenvalue of the different laplacians for a multiplex of 4 layers. Each layers contains a network with 4 communities, each community corresponds to an Erdös-Rényi network with edge probability 0.5, and the inter-community edge probability is 0,05. The communities between different layers strongly overlap. Comparison between the second eigenvalue of the different laplacians for a multiplex of 4 layers. Two of the layers contain a network with 4 communities, each community corresponds to an Erdös-Rényi network with edge probability 0.5, and the inter-community edge probability is 0,05. The two communities strongly overlap. The other two layers contain an Erdös-Rényi network of 200 nodes with edge probability of 0.5.  Comparison between the second eigenvalue of the different laplacians for a multiplex of 4 layers. Each layers contains a scale-free network of 200 nodes generated using the Barábasi-Albert model. The first two layers have been generated attaching, at each step, the new node to a 3 existing nodes, the third and fourth layers attaching the new node to a 7 existing nodes.