A polynomial eigenvalue approach for multiplex networks

We explore the block nature of the matrix representation of multiplex networks, introducing a new formalism to deal with its spectral properties as a function of the inter-layer coupling parameter. This approach allows us to derive interesting results based on an interpretation of the traditional eigenvalue problem. More specifically, we reduce the dimensionality of our matrices but increase the power of the characteristic polynomial, i.e, a polynomial eigenvalue problem. Such an approach may sound counterintuitive at first glance, but it allows us to relate the quadratic problem for a 2-Layer multiplex system with the spectra of the aggregated network and to derive bounds for the spectra, among many other interesting analytical insights. Furthermore, it also permits to directly obtain analytical and numerical insights on the eigenvalue behavior as a function of the coupling between layers. Our study includes the supra-adjacency, supra-Laplacian, and the probability transition matrices, which enable us to put our results under the perspective of structural phases in multiplex networks. We believe that this formalism and the results reported will make it possible to derive new results for multiplex networks in the future.


I. INTRODUCTION
Complex network theory has become one of the main tools for the analysis of complex systems, allowing the representation of a wide range of systems composed by interacting discrete elements [1]. However, real-life systems are also organized in layers, which represent different channels of interaction. In order to incorporate these characteristics, one should work with multilayer networks, which allows for a proper representation of multiplex and interconnected systems [2][3][4]. The introduction of this extra level of complexity also imposes new challenges on the analysis of its structural and dynamical properties. Furthermore, a key element on the analysis of networks is their spectral properties [5]. In fact, they play an important role in explaining the connection between structure and dynamics. For instance, in epidemic spreading the critical point below which the infection prevalence is null is predicted to be the inverse of the leading eigenvalue of the adjacency matrix, in both, single [6] and multiplex networks [7]. Additionally, its nature also seems to be connected to those properties [7][8][9]. Although the literature about the spectra of singlelayer networks is well developed [5], the theory of spectral properties of multiplex networks is still in its infancy. This motivates us to propose a different formalism aimed at filling this gap.
In this paper, we will consider the matrix representation of multiplex networks, constraining ourselves to finite matrices. First of all, we are interested in weighting differently inter and intra-layer edges. This implies that those matrices will also be a function of the inter-layer coupling parameter, here called p. Consequently, the associated eigensystem will be a function of that same parameter. Additionally, the matrix approach is especially interesting in this context since it allows us to directly use linear algebra and spectral graph results already available.
When varying the coupling parameter, a multiplex system might present different structural phases, which are characterized in terms of eigenvalue crossings and eigengaps and are intuitively defined as: (i) decoupled phase, for small values of p, where the layers are virtually decoupled and act by themselves, with a negligible interaction between layers, (ii) multiplex/multilayer phase, where the system is coupled and the intra-layer edges play an important role and (iii) the aggregate network phase, where the system behaves as the superposition of all layers. It is clear that a good understanding of the eigenvalues' behavior might be useful since we could move our system into different structural regimes, aiming at different goals such as improved robustness, better performance regarding diffusion or spreading, among many other possible applications.
The different structural phases above are related through the interlacing properties of quotient graphs [10,11]. More specifically, in [10] the authors showed that the spectra of different scales of a multiplex network (aggregated network, the network of layers and individual layers) characterize the three phases. In practical terms, the interlacing provides us bounds for the spectra [10,11], but also emphasizes that the different scales are intrinsically connected. Indeed, it is impossible to tune the leading eigenvalue of the network of layers without also increasing the leading eigenvalue of the whole multi-plex. Furthermore, in [12] the authors characterized multiple topological scales using the supra-Laplacian matrix. More specifically, they analyzed eigengaps to characterize them.
Following a different approach, in [13], the author evaluated the normalized Laplacian matrix, which, in fact, shares the same set of eigenvalues of the probability transition matrix (see Section IV C), and proposed a similar classification. However, it is worth mentioning that in [13], a different nomenclature was used and a fourth phase was defined. Namely, the proposed structural regimes were: (i) bipartite phase, (ii) decoupled phase, (iii) indistinguishable, where the author argues that the system is topologically and dynamically indistinguishable [13], and (iv) a mixed phase, (called BD in [13] -bipartite and decoupled phases) where the layers are structurally and dynamically distinguishable. Although here we do not make a distinction between the regimes (iii) and (iv) and consider both as multiplex regimes, we acknowledge the differences pointed out in [13]. It is also noteworthy that [13] considered structural correlations for the analysis, which is a key ingredient for the reported results. On the other hand, here we focus on uncorrelated networks and a multiplex structure.
The paper is organized as follows. In Section II, we present the polynomial eigenvalue formalism, giving its general definitions and properties in Section II A. Next, we formalize the 2-Layer problem into a quadratic eigenvalue problem in Section III, analytically exploring its behavior as a function of the coupling parameter. In Section III A, obtaining some bounds, in Section III B, and discussing the simplified symmetric problem in Section III C. We present our main applications in Section IV, where we explore the supra-Laplacian matrix, in Section IV A, the supra-adjacency matrix in Section IV B and the probability transition matrix in Section IV C. To round off this paper, we discuss the physical consequences of our findings, summarize our main results and perspectives in Section V.

II. POLYNOMIAL EIGENVALUE PROBLEM
In this section, we formally define the polynomial eigenvalue problem and present some of its fundamental properties. The aim of this section is to generically define our main mathematical object, establishing its basic properties. Thus, this will allow us to properly study the matrices associated with 2-Layer multiplex networks, which will naturally appear as a consequence of a simple manipulation of a linear system describing the network. From this simple approach, we expect to provide a different perspective on the spectral properties of multiplex networks.

A. General definition and properties
A matrix polynomial of order l is a matrix-valued function of a complex variable of the form [14] where M 0 , M 1 , . . . , M l are n × n matrices and they are said coefficient matrices. If M l = I the identity matrix the polynomial R is said to be monic. The eigenvalues of R are the solution to the characteristic equation Right and left eigenvectors are defined by where x and y are the right and left eigenvectors associated to the eigenvalue λ. It reduces to the standard eigenvalue problem when R(l) is a monic matrix polynomial with l = 1 and M 0 = −A, for any matrix A.
A generalization of the Jordan form theory to a general matrix polynomial is possible and is briefly presented in the Appendix A.
A particular case of general interest is the Quadratic Eigenvalue Problem (QEP), which is directly related to the 2-layer case of our interest in this work. A quadratic matrix polynomial can be written as [14][15][16] Q(λ) = Aλ 2 + Bλ + C.
Besides, it is worth mentioning that, without loss of generality, we assume in the following that the eigenvectors are unitary. Note, however, that if the right eigenvector x is normalized, then the left eigenvector y is not. In the following, we assume that x is unitary to simplify the equations.
Furthermore, a special class of problems is obtained if A, B and C are Hermitian, called Hyperbolic Quadratic Eigenvalue Problem (HQEP) [16]. Unfortunately, most of the problems in our context does not fall in this class.

III. A GENERAL 2-LAYER CASE: A BLOCK MATRICIAL PROBLEM
The general form of any matrix associated to a multiplex network composed by two layers (supra-adjacency, supra-Laplacian and transition matrices for example) can be written as a block matrix. The resulting eigenvalue problem is the following where M 12 = M T 21 . Interpreting it as a system of equations and isolating v 1 on the second row, we have Finally, inserting it in the first row we have This expression defines a QEP, whose coefficient matrices are which poses a restriction on the inter-layer coupling matrix M 12 , i.e., it must be invertible. Additionally, in our context, exchanging M 11 and M 22 does not change the system, neither the solutions. Note that this operation is equivalent to relabeling the layers. However, if the polynomial of the first is Q(λ), then, for the second it is Q(λ) T . In this way we found a relation between the right and left eigenvectors and these two possible configurations of our system. Formally, such observation implies x = v 2 and y = v 1 . As usual, we consider coupling matrices that are functions of a coupling parameter, p, i.e., M ij = M ij (p), for i = j. In fact, throughout this paper we explore how the spectral properties of our network evolve as we change such coupling parameter. As a constraint, we should mention that we only consider finite matrices. Furthermore, note that B in equation 9 is intimately related to the aggregated and the loop-less aggregated networks of the original multiplex network (for more, see [10] or Section 2.3.2 of [17]). More specifically, if the coupling matrices M 12 and M 21 are the identity matrix (or proportional to this matrix), thus B = − (M 11 + M 22 ), which is proportional to the loop-less aggregated network. Besides, note that a network of layers in the twolayer multiplex is a simple line graph with two nodes.
A. Spectral analysis as a function of p So far, we have defined our main mathematical objects, making as less constraints as possible. Now we restrict ourselves to diagonal coupling matrices -i.e., multiplex networks -and assume a linear function of the parameter p > 0 [27], M 12 = pD, where D is a diagonal invertible matrix (such constraint will be relaxed later). Then, defining the scalar equation that describes each eigenvalue as the product of Q(λ), by its left and right eigenvectors we have where a(y T , x) = y T Ax, b(y T , x) = y T Bx and c(y T , x) = y T Cx. The solution of this equation is given by where ∆(y T , x) = b(y T , x) 2 − 4a(y T , x)c(y T , x). Note that for each pair of right and left eigenvectors we have two possible solutions, but just one of them is an eigenvalue of Q(λ). Additionally, differentiating equation 12 by p we obtain information on how the eigenvalues change as p changes. Formally we have Note that the eigenvalues and eigenvectors are also a function of p. For continuity, two different eigenvalues may cross each other when varying p. Observe that for non-crossing points the relations dy T dp Q(λ)x = 0 and y T Q(λ) dx dp = 0 holds, since the derivatives are bounded for non-crossing points. However, on the crossings we have two eigenvectors associated to the same eigenvalue, which imply two solutions for the derivatives. Then, isolating the derivative of λ we have dλ dp = Such relation can be applied to drive a system through different regimes. For instance, considering the adjacency matrix, one can use this equation in order to chose an edge or set of edges to be removed (or weighted) in order to optimally reduce or increase the leading eigenvalue and consequently the critical point of spreading processes, such as epidemic spreading. Obviously the matrix under study depends on the process. Another application is to design a numerical method to follow the correct eigenvalues as a function of p in a problem that might present eigenvalues crossings [7,18].

B. Bounds
Aiming to find bounds to equation 13 we study the scalar polynomial defined by x T Q(λ)x = 0, where x is an eigenvector (left or right), which guarantees that the polynomial is equal to zero. In order to simplify the problem we multiply Q(λ) by D, obtaining a monic matrix polynomial , then we must bound the terms b(x T , x) and ∆(x T , x), which allow us to bound both solutions. Those terms can be bounded by the numerical range of the matrices to which they are related. The numerical range is formally defined for any matrix X as F (X) = {x T Xx : x ∈ C and x T x = 1}. Additionally, σ(X) ⊆ F (X), where σ(X) is in the set of eigenvalues of X. Moreover, if X is an Hermitian matrix x T Xx is the Rayleigh quotient of X, which implies λ 1 (X) ≤ x T Xx ≤ λ N (X). Finally, to bound a non-Hermitian matrix we use the relation of the spectral norm and the numerical range, given as 1 2 |||X||| 2 ≤ r(X) ≤ |||X||| 2 , where r(X) is its numerical radius, defined as r(X) = max First, consider the term b(x T , x), which is bounded by however, in many cases, B is an Hermitian matrix, allowing us to improve this bound to More precisely, observe that B is often related to the aggregated network, connecting both scales.
x and defining the matrix ∆ = B 2 − 4C, we can focus on the problem x T ∆x instead of the initial definition of ∆(x T , x), since both have the same bounds. Besides, since in most of the problems on networks we are dealing with symmetric matrices (undirected networks), we might also impose that ∆(x T , x) ≥ 0 because we already know that the spectra is real in this case. consequently, we have Observe that those bounds can be further improved when applied to the analysis of particular matrices (supraadjacency, supra-Laplacian, and probability transition) since their particularities also impose constraints on the solutions and could be explored to improve the bounds.

C. Comments on symmetric problems: HQEP
As previously mentioned, if A, B and C are Hermitian, we have a special class of problems called Hyperbolic Quadratic Eigenvalue Problem (HQEP) [16]. The HQEP has interesting properties, for instance, if x is a right eigenvector associated with the eigenvalue λ, then it is also a left eigenvector of the same eigenvalue [16]. In order to take advantage of those properties one can interpret the original problem as an HQEP plus asymmetric perturbation. Thus, the matrix polynomial defined by the matrix coefficients in Equation 9 is not symmetric in most cases. However, a class of problems that arise naturally is defined by M 12 = p I and in this case the matrices A and B are Hermitian. Observe that C still might be asymmetric. However, we can use the Toeplitz decomposition [19] in order to analyze a simplified problem. Such decomposition states that any square matrix can be uniquely written as the sum of an Hermitian (X = X * ) and a skew Hermitian matrix (X = −X * ) as . This allows us to In this way we can re-write our QEP into two parts, one composed by Hermitian matrices, which is a HQEP, and a skew Hermitian matrix, that can be interpreted as a perturbation. The natural consequence from the perturbation theory is that the matrix pC of the HQEP is perturbed by 1 2 (M 11 M 22 − M 22 M 11 ) and such matrix norm goes to zero as the layers are more similar. From the Bauer and Fike theorem [19] we can write a quality function for the approximation of the perturbed matrix C as whereλ is the eigenvalue of C = C H + C S , C H = UΛU −1 and κ(·) is the condition number with respect to the matrix norm |||·|||. Considering the spectral norm is near 1, small perturbations imply small changes on the eigenvalues. On the other hand, large values of κ(U) suggest a poor approximation. Observe that such analysis concerns only the matrix C and not the whole QEP, however, it can be an estimate of the quality of the approximation and show that the general solution interpolates between a HQEP and a general QEP.
In addition to the HQEP properties, the perturbation analysis also emphasizes an important multiplex property. We must note that the more similar the layers are, the closer to zero the norm 1 2 (M 11 M 22 − M 22 M 11 ) is. Moreover, we also have another criteria which is based on the commutativity of the matrices M 11 and M 22 . Moreover, observe the role of correlations in this approximation. If both layers are identical, they are obviously correlated and the problem is symmetric.

D. Limits for sparse inter-layer coupling: singular
So far we have assumed a node-aligned multiplex, i.e., a multiplex network in which each node has a counterpart on every layer [2], fulfilling the invertibility of M 12 , which is necessary to formally define the problem, however, we can use the limit of D ii = → 0 to obtain an approximation of the sparse coupling. Observe that equation 5 can be analyzed in two different steps, first calculating the limit of decoupled edges and secondly the rest of the system. The first limit is analyzed as follows. From 5 the absent edges are factorized as whereD(i) = D −1 . Multiplying equation 21 by p and using the following limit we havẽ where the term of order p D vanishes in the limit of → 0. Observe thatD = lim →0 D −1 = I if both layers are decoupled and the polynomial equation can be factorized as (M 11 − λ I) (M 22 − λ I) = 0, whose solutions are the union of the solution of the standard eigenvalue problem of each layer. An important observation is that the number of nodes that are not connected to the other layer is also the number of eigenvalues that do not change as a function of p.
Equation 23 presents the solution for nodes that do not have any counterpart on the other layer. In order to calculate the remaining solutions we have to redefine the original problem in terms of the Moore -Penrose pseudoinverse, denoted by X † , for a matrix X. Denoting byD = p −1 D † we haveD jj = p −1 D −1 jj if D jj = 0 andD jj = 0 otherwise. Note that the zeros ofD jj are ones inD jj . For the sake of simplicity, in the following we assume that M 12 is invertible, however, the strategy mentioned above can be applied if it is not the case. From the computational point of view, we can reduce the cost to calculate the whole spectra as a function of a closed range of p by separating it into two components, where a subset is constant and the remaining subset varies.

IV. APPLICATIONS
We next apply our main formalism to study the supra-Laplacian and the supra-adjacency matrices. For the sake of completeness, let us explicitly define the supraadjacency matrix in terms of its block matrices (adjacency matrix of the individual layers). Formally, the supra-adjacency matrix is defined as where we weight differently the intra and inter-layer edges. The definition of the supra-Laplacian matrix is whereD is a diagonal matrix whose elements are D ii = A ij and the Laplacian matrices of the individual layers are denoted as L a and L b . Those matrices are related to many dynamical processes. The supra-Laplacian is used to describe diffusion and synchronization of coupled oscillators, while the supra-adjacency matrix is intimately related to epidemic and information spreading. It is also noteworthy that many structural metrics are also directly extracted from the spectral properties of those matrices. For instance, the communicability, which can be easily written as a matrix function, or more specifically, as the exponential of the adjacency or supra-adjacency matrix.
Here we are going to focus on the spectral properties of these matrices, and their behavior as a function of the coupling parameter p under different conditions.
In addition to the supra-adjacency and supra-Laplacian matrix, we also analyze the probability transition matrices in Section IV C, which can be used to describe classical random walks on networks. The analysis of such a matrix is left to the last section since it is mainly numerical. Note that the probability transition matrix has a well bounded spectra where [20]. This characteristic imposes an extra challenge on the derivation of the bounds. Although we could not improve those bounds, we report an interesting spectral behavior found numerically.

A. Supra-Laplacian matrix
The simplest supra-Laplacian matrix can be built considering a diagonal coupling matrix M 12 = −p I, where each node has a counterpart on the other layer and the coupling is homogeneous. This implies that the QEP is defined with the following coefficient matrices It is noteworthy that the aggregated network, L + = L a + L b , appears naturally under this formalism. This is interesting since it is physically understandable. On the other hand, the term L a L b , in the definition of C, is of not so direct interpretation. This system presents a structural transition, which can be directly derived from our formalism. This derivation is presented in Section IV A 1. Additionally, we can also obtain bounds for the spectra using the ideas discussed in Section III B. Those improved bounds are derived in Section IV A 2, where we use the particular properties of a Laplacian matrix to improve our previous results. In Section IV A 3, we evaluate the spectra of the supra-Laplacian matrix as a function of p and also compare our previous results with sparse and heterogeneous couplings. Specifically, on the heterogeneous case we consider a coupling matrix The analysis of such QEP is not trivial, since the matrices are not symmetric, however we can explore it numerically and compare with the homogeneous case, M 12 = −p I.

Structural transitions
Firstly, we discuss the structural transition presented in [21] on the Laplacian matrix. Here we calculate the exact transition points using the QEP formulation. We can easily derive such transition points using our formalism. It is noteworthy that those transition points were also calculated in [18] using two different methods: eigenvalue sensitivity analysis and a ShurâĂŹs complement approach. Both derivations are quite complicated, contrasting with our approach, where the solutions are given using simple arguments. Note, however, that our approach presents a different expression if compared to the method presented in [18], but both expressions yield the same final result. We do not prove the equivalence mathematically, but verified their equivalence numerically.
To begin with, it is well known that λ = 2p is an eigenvalue of the supra-Laplacian and the crossing points are a consequence of this eigenvalue crossing the bounded part of the supra-Laplacian spectra, producing the so-called structural transitions. In this way, from our definition of QEP, we have that which has two possible solutions: (i) det (L a + L b ) = 0, which is always true, since the sum of two Laplacian matrices is also the Laplacian of the aggregated network and also has determinant equal to zero and (ii) the crossing points or eigenvalues of multiplicity larger than one. Since it is also an eigenvalue problem in terms of p, we have that the crossing points are expressed as There are N possible values of p that solve Equation 29, each one representing one crossing. The first crossing is trivial, at p = 0, the second is the one called structural transition in [21], which is relevant for some dynamical processes [22]. As said before, this expression is different from the previous one presented in the literature, however both give the same result.

Bounds
Here, exploiting the ideas presented in Section III B, we improve the former bounds using specific Laplacian properties, such as its semi-positiveness. Thus, the QEP of the supra Laplacian can be bounded considering the individual bounds of B, which is a semi-positive definite Hermitian matrix, leading to Besides, the discriminant function is also bounded by where the upper bound can be defined as a function of the spectral properties of (L a − L b ) 2 . On the other hand, regarding the lower bound, it can be improved by realizing that the matrix ∆ = (L a − L b ) 2 − 2p (L a + L b ) + 4p 2 I, defined on Section III B, is semi-positive definite for undirected networks, ∆ 0. In this way, [28]. From these properties, we can establish the lower bound as 4p 2 . Formally, The previous bounds imply that in the asymptotic analysis formalism we have ∆(x T , x) ∈ Θ(p 2 ). Moreover, observe that the lower and the upper bounds converge to each other as the layers become similar. On the extreme case of identical layers we have ∆(x T , x) = 4p 2 . Finally, combining the formerly obtained bounds we have, and Interestingly, these bounds can be analyzed in terms of their asymptotic behavior (approximation), where for a sufficiently large value of p they can be approximated to Observe that from the asymptotic point of view we have λ − (x) ∈ Θ(1) and λ + (x) ∈ Θ(p).
As an example, in Figure 1 we present the evaluation of the eigenvalues as a function of the coupling parameter p of a multiplex network composed by two Erdös Renyi layers with n = 10 3 nodes. The first layer has an average degree k = 12, while the second has k = 16.

Spectral properties as a function of the coupling p
In this section, we focus on the behavior of the eigenvalues as a function of the coupling parameter, λ i (p). First of all, we apply the concepts of Section III A regarding the derivative of Q(λ). Consider the simplest case, where D = I. In such a case, we have a monic polynomial matrix, where B depends on the aggregated network, which is semi-positive definite. Besides, C is a matrix that contains the product of both layers and accounts for similarities between them. In this way, Equation 16 can be expressed as whereb(y T , x) = y T (L a + L b )x and y T x = cos(θ) is the cosine of the angle between left and right eigenvectors of our QEP. Observe that part of the spectra has dλ dp → 0, while the other part has dλ dp → 2 as p increases, which can be proved as follows. Firstly, suppose that λ is constant as a function of p, then dλ dp → 0 because the denominator grows as a function of p and the numerator is bounded, as supposed. Secondly, suppose that λ grows with p r , where r < 1. In this case dλ dp → 0, by the same arguments as before, since the linear function of the denominator dominates it. However, if r = 1 we have dλ dp → 2, since both, the numerator and the denominator, grow linearly. Finally, with p r , where r > 1, both the numerator and denominator are dominated by p r , which imply that for large p the derivative dλ dp → 1, which is also a contradiction, since it was supposed to be a linear function of p. In this way, we conclude that the derivatives of λ, for large values of p, cannot grow faster than linearly and their growth will be one of two values, 0 or 2. These results are in agreement with the previously obtained bounds. Additionally, as an example, in Figure 1 we also observe such a behavior.
Although for the simplified case we have two possible solutions at large p, observe that the above arguments fail for the case of general coupling matrix. From Equation 16 and the definition of the Laplacian QEP we conclude that only the denominator of Equation 16 changes for a different choice of D, since the terms that have dependencies on D vanish in the partial derivatives of the numerator. The denominator follows the general form y T 2λD −1 − 2p I −L a D −1 − D −1 L b x. In this way, different coupling weights can change the behavior of each eigenvalue differently for large p. For instance, if D = I the spectral distribution for large p is bimodal, however, if D = p n n i i diag(1, 2, ..., n) [29], this behavior changes completely and the eigenvalues change with different rates, presenting a "continuous" bulk. This argument is valid for infinity size networks, since for finite size networks for a large p gaps between eigenvalues may appear due to different rates of growth (as a function of p). An example of this is shown in Figure 2. Furthermore, we also found an empirical function that seems to bound the spectra as a function of p in this experiment. The lower bound is trivial, since it is a semi-positive definite matrix. The upper bound can be obtained correcting p = max{diag(D)}p, hence From Figure 2 we observe that such bound is not as close to the largest eigenvalue as the homogeneous case. In addition to a non-homogeneous coupling matrix, the last case studied is the sparse coupling. The analytical part of this study was presented in Section III D. As predicted, each uncoupled node implies a pair of eigenvalues that does not depend on p. Due to the nature of the Laplacian matrix, where just one eigenvalue varies with p, while the other remains bounded, the set of bounded eigenvalues increases by one. For example, if we haveñ uncoupled nodes, the bounded part have n +ñ eigenvalues, while the "unbounded" part have n −ñ eigenvalues. Note that the upper bound for the bounded part is not 1 2 λ max (L a + L b ) anymore. However, the general upper bound for D = I seems to be also an upper bound for the sparse problem, as we numerically verified. The figures for these are not shown since they are visually similar to Figure 1.

B. Supra-adjacency matrix
Similarly to the supra-Laplacian case, here we also begin with the simplest case, i.e., the diagonal homogeneous coupling, and increase the level of complexity considering heterogeneous inter-layer weights and sparsity. Thus, in the simplest case we have M 12 = p I, therefore, the QEP, Equation 5, is defined by the following coefficient matrices Note that, in a way similar to the Laplacian, B is also defined in terms of the aggregated network. On the other hand, the physical interpretation of C is still difficult due to the product A a A b . In Section IV B 1 we improve the bounds proposed in Section III B. Then, in Section IV B 2, we evaluate the spectral properties of the supra-adjacency matrix as a function of p in three different contexts: (i) diagonal homogeneous coupling, (ii) diagonal heterogeneous coupling and (iii) sparse diagonal homogeneous coupling. Note that, in order to analyze the heterogeneous coupling we must consider the general QEP with the following coef-

Bounds
Similarly to the analysis performed for the supra-Laplacian, here we also extend the ideas presented in Section III B to the supra-adjacency matrix. First of all, regarding the diagonal heterogeneous coupling case, D = I, we can also find bounds for the spectral distribution of the adjacency matrix. Beginning with B, we can bound it based on its eigenvalues as Interestingly, those are the eigenvalues of the aggregated network, which have a clear physical meaning. Similarly, for the discriminant we have Finally, combining those bounds we can bound both solutions by which asymptotically converge (as an approximation) to (45) In other words, the spectral density of the adjacency matrix is bimodal and part of the eigenvalues grows linearly with p, while the other part decreases at the same rate.

Spectral properties as a function of the coupling p
In its general form, the first derivative is given as where x and y T are the right and left eigenvectors associated with λ. Firstly, focusing on the case where D = I and using a similar approach as that applied to the Laplacian case, we can suppose that λ is a constant function of p or a function of p r with r < 1, however, it would give us dλ dp ∼ p, which is a contradiction. Next, we can suppose that it is a linear function of p, which implies dλ dp → ±1, depending on the sign of the linear coefficient. Finally, supposing that it is a function of p r with r > 1 we obtain that dλ dp → 0, since the denominator grows faster than the numerator, which again is a contradiction. In this way, based on such analysis we infer that the first derivative of λ can assume only dλ dp → ±1.
Secondly, for the general case observe that both, the numerator and the denominator of Equation 46 vary as a function of D. Additionally, D weights the product of the components of the eigenvectors, which allows the derivatives to assume more values, even a "continuous bulk" instead of the bimodal distribution of the diagonal homogeneous case, similarly to the case discussed for the supra-Laplacian matrix. Here we also use the coupling matrix D = p n n i i diag(1, 2, ..., n) . We show the spectral evolution as a function of p for the non-homogeneous case in Figure 4. Similarly to the Laplacian case, there are evidences that the bounds can be corrected using p = max{diag(D)}p, hence Here we can also obtain a similar conclusion as for the supra-Laplacian case. From Figure 4 we observe that the corrected bounds are not as close to the homogeneous case. Finally, we evaluate the sparse coupling case, whose analytical study was presented in Section III D. As predicted, each uncoupled node implies a pair of eigenvalues that does not depend on p. In this way, if we haveñ uncoupled nodes, the central part of the spectra will have 2ñ eigenvalues that do not change as a function of p. Next, n−ñ grows linearly with p, while the other n−ñ eigenvalues with −p. This is illustrated in Figure 5. Note that in this figure the horizontal lines bounding the central part of the spectra are not calculated, but numerically obtained and are only shown to serve as a reference.

C. Probability transition matrix
In this section, we evaluate the probability transition matrix, mainly focusing on its spectral properties as a function of the coupling parameter p. Due to the probabilistic nature of this matrix, we were not able to improve its bounds. Therefore, we mainly report numerical results.
Formally, the probability transition matrix is defined whereD X ii = k X i + p, and X = {A, B} represents the label of each layer. It is known that this matrix models the classical random walk, where the walker choses a neighbor based on the weights of its surrounding edges.
It is important to mention that in [23] the authors studied random walks on top of multiplex networks and analyzed them in terms of the normalized supra-Laplacian matrix. This matrix is defined as where L is the supra-Laplacian and P is the probability transition matrix. Note that the normalized supra-Laplacian matrix is intimately related to the probability transition matrix. In fact, their spectra are trivially related. Furthermore, we can also relate the spectra of the normalized Laplacian as follows where S =D − 1 2 AD − 1 2 has the same set of eigenvalues as P and if v is an eigenvector of S, thenD −1 v is an eigenvector of P associated with the same eigenvalue [20]. Note, however, that S is symmetric [20]. In the context of random walks in multiplex networks, in [13] the author used the normalized supra-Laplacian matrix. Here, in this section, we will study P, defined in Equation 49.
Next, following our formalism, from Equation 49, we can define our QEP in its monic form as Note that such quadratic polynomial present some similarities with the one for the supra adjacency matrix, however the probability transition matrix is not symmetric and the matricesD X presents a dependency on p. This fact, associated with the natural bound for stochastic matrices, make the derivation of the spectral bounds more complicated than the previous cases. Here we focus on the spectral properties of the probability transition matrix as a function of the coupling parameter p.

Spectral properties as a function of the coupling p
For the sake of completeness, let us study the spectral properties of the transition matrix as a function of the coupling strength p. This exercise is much more of an example than a practical application since we already know that the spectra are bounded on stochastic matrices, which does not allow unbounded grow. In Figure 6 we present the spectra as a function of the coupling parameter p. The first observation is that, aside from being bounded, the growth rate of the eigenvalues is quite different from what was observed for the Laplacian and adjacency cases.
Firstly, lets proceed with the analysis of Equation 16 aiming for an approximation, which qualitatively describes the λ(p). First of all, the partial derivative of B can be expressed as where the term ∂ ∂p Next, expanding the partial derivative of C we have All the expressions obtained so far are quite complicated to be analyzed in its exact form. Thus, we will proceed with an asymptotic analysis, aiming for a hypothesis of a possible formula that qualitatively describes the behavior of λ i as a function of p. In other words, we propose a formula that fits the expected asymptotic behavior, but we also expect it to work for smaller values of p. We must remark that this analysis is an approximation and, in order to verify its validity we perform numerical fittings and evaluate the obtained errors.
From the previously mentioned perspectives, the asymptotic behavior are ∂C Next, in Figure 7 Note that we considered a single value of k x , without considering products between different constants. Besides, just one term is considered. We remark that the main goal of this exercise is to have insights on the qualitative behavior of more complicated functions, such as Equation 57. In other words, the performed approximations are not expected to quantitatively predict those terms, but qualitative represent and "catch" the main behavior of those functions.
Firstly, recall that the spectra on stochastic matrices is bounded, which consequently restricts its derivatives. In other words, λ i ∈ O(1). However, for the sake of the argument, let us suppose that λ i = c 1 p r +O(p r−1 ), hence dλ i dp = c 1 rp (r−1) + O(p r−2 ), where r is an integer. Thus, comparing with Equation 16, we have that can be rewritten as which simplifies to which implies that dλ dp → 0 and r ≤ 0 since on the lefthand side we have a function in O(max{2r − 1, −1}), while, on the right-hand side we have a function in O(max{r − 2, −1}). This simple analysis suggests that r ≤ 0, for consistency. Note that we are not inferring anything regarding its "velocity" (how fast it goes to zero). Such arguments reinforce that λ i ∈ O(1), as previously mentioned. However, there are a huge class of functions that satisfies such restriction. In order to satisfy the so far established restrictions, let us suppose that which is a function that satisfies our previous analysis. Thus, it also implies that which yields dλ i dp ∈ O(p −2 ). Next, from Equation 16 we Note that it allows a set of possible solutions and, among them, it allows our initial supposition, Equation 61, i.e., Next, we proceed with a numerical experiment, extracting some eigenvalues presented in Figure 6 we perform a fitting aiming to obtain the same curve. We chose 5 eigenvalues: (i) the leading eigenvalue, λ 1 = 1, just as a reference and to emphasize that our proposed equation also works for that case, (ii) λ 3 , the first eigenvalue on the bulk (note that there can be a crossing between λ 3 and λ 2 , which would change the index of the eigenvalue -here we are not going to enter into details of this possible crossing behavior and, in order to avoid that, we   error of these two curves is expected to be larger than the previous cases. It is important to remark that, as previously mentioned, there can be a crossing between λ 3 and λ 2 , but here we are looking at the main global behavior and such a change would not be a big source of error. In this way, we are showing that there is a set of parameters that approximates the spectra using the proposed equations. The proposed mentioned experiment does not serve as a proof, but it does serve as an evidence of such, or a similar, behavior. Following the proposed pipeline, firstly, in order to obtain the fittings, we used the nonlinear least squares method, the Levenberg-Marquardt algorithm [24][25][26] and the least absolute residual (LAR) robust regression. Additionally, all the initial conditions were set to one. In Figure 8 we show the obtained fittings and the numerically obtained eigenvalues. Complementary, in Table I we present the fitted parameters. Interestingly, we observe that the proposed approximation fits really well the observed curves, which can be objectively measured by means of the Sum of Squares Due to Error (SSE), whose values are also reported in Table I. Thus, the behavior of λ i assumed in Equation 61 seems to be a very good guess. Besides, we also observe that there seems to be a symmetry on the obtained parameters for λ 3 and λ N , which are close. The only exception is k 0 , since both have the same modulus, but with a different sign, as expected. The last important observation also regards the parameter |k 0 |. Note that such a parameter is very close to one on all the fittings, suggesting some underlying property of our formulation.
Finally, for the sake of completeness and for comparison reasons, we numerically evaluate the spectra of the probability transition matrix for the sparse and heterogeneous coupling cases. Regarding the sparsity, in Figure 9, we present a similar experiment as done for the supra-Laplacian and supra-adjacency cases. Similarly to those experiments, here we also observe a group of eigenvalues that do not change as a function of p. Moreover, we also verified that forn = 100 decoupled nodes, we have 2n = 200 eigenvalues that remain constant, validat-ing the insights we obtained in Section III D. Although the behavior observed is similar to the previously studied matrices, for the probability transition matrix we observe a slightly different behavior for intermediate values of p (here 1 < p < 10), where the intermediate eigenvalues change, forming the "central bulk".
Furthermore, we remark that for the heterogeneous coupling (D = p n n i i diag (1, 2, ..., n) ), if compared with the supra-Laplacian and supra-adjacency, a completely different behavior emerged. In the probability transition matrix case, the spectra seem to be always bimodal. This effect is shown in Figure 10, where, for a large enough value of p, the eigenvalues tend to a constant. It is also noteworthy that the rate at which this phenomenon takes place is much slower than the rate of the homogeneous case, shown in Figure 6.

V. DISCUSSION AND CONCLUSIONS
From the developed theory, we applied and analyzed three different matrices: (i) the supra-Laplacian, (ii) the supra-adjacency and (iii) the probability transition matrix. In all these cases we have considered three different coupling schemes: (a) diagonal homogeneous coupling, D = p I, (b) diagonal homogeneous sparse coupling and (c) diagonal heterogeneous coupling, D = p n n i i diag(1, 2, ..., n) . Regarding the supra-Laplacian and the supra-adjacency matrices, on the first scenario, (a), we were able to extract some analytical results regarding the derivatives of the eigenvalues, which suggested a different behavior for the other two cases, (b) and (c). On the other hand, regarding the probability transition matrix, due to its stochastic nature, we were not able to go further with the analytical analysis. However, we followed an asymptotic analysis, proposing a function that describes the eigenvalues behavior. This function was validated with numerical fittings of the original spectra. Although it is just an approximation, it also helps us understand the nature of the phenomena behind this structure. Furthermore, we also reported the differences between the spectral distributions for large p, where we can have bimodal, multi-modal or even a continuous bulk for the adjacency and Laplacian cases, just changing the coupling matrices. On the sparse case, this analysis was analytically supported, while the other cases were explored numerically.
Our analysis pointed out some important features about multiplex systems. As a general observation, as we increase p we will find (roughly) three different structural phases, which might take place at different points for each structure and matrix. Thus, the structural phases of a multiplex network can be defined as: (i) decoupled phase, for small values of p, where the layers are virtually decoupled and act by themselves, with a neglectable interaction, (ii) multiplex/multilayer phase, where the system is coupled and the intra-layer edges play an important role and (iii) a network of layers phase, where the structure of the network of layers plays the major role. Note that, from the perturbation theory point of view, in the decoupled phase the eigenvalues are basically the union of the eigenvalues of the individual layers plus some perturbations. On the other extreme, in the network of layers phase, the intra-layer edges might be understood as the perturbation, since p 1. In this case, we can interpret the system as a set of n virtually disconnected small networks, whose structure is given by the network of layers (considering a multiplex case, where each node has a counterpart on the other layers). Finally, the most interesting scenario is the multiplex phase, where the inter and intra-layer topologies play a fundamental role on the dynamics.
Throughout our analysis, we were able to verify these regimes in the different matrices we evaluated. Although all of them showed this behavior, the differences between those matrices are also evident. Considering the supra-Laplacian and supra-adjacency matrices with diagonal homogeneous coupling, we observe that, for a large enough p, the spectral distribution is bi-modal, while for the sparse case we have three bulk's, where the central one results from the nodes that do not have an inter-layer edge. Finally, on the diagonal heterogeneous coupling, we observe a completely different behavior, where the eigenvalues are distributed into a single bulk. We remark that, although we were not able to analytically quantify this last phenomenon, our analysis suggested such a behavior.
Furthermore, comparing those results with the ones obtained using the probability transition matrix, we observed a completely different behavior. On the diagonal, homogeneous or heterogeneous cases, the spectra seem to be bi-modal for a sufficiently large p. Note that for the homogeneous case this convergence to the bulks is much faster than the heterogeneous case. This is an interesting phenomenon since it contrasts with the supra-adjacency and supra-Laplacian cases, where the heterogeneous coupling implies a "continuous" bulk. It is noteworthy that our predictions for a central bulk for the uncoupled nodes are also fulfilled for the probability transition matrix.
The analysis performed here emphasize the importance of a proper study of the structural phases in different contexts. The sparse and heterogeneous cases might change completely the spectra (depending on p). Obviously, the analysis should also take into account the correct matrix since the structural changes are different from case to case. In other words, different matrices present their phases in different intervals (values of p). In this context, dynamical insights can also be useful for a better understanding.
Since all the analyzed matrices are also related to dynamical processes, the results reported here will directly impact on these processes too. Note that in [22] the authors analyzed diffusion processes in multiplex networks and found the so-called superdiffusion. This process is described by the supra-Laplacian matrix and it is intrinsically connected to the so-called structural transition of this matrix as pointed in [22] and latter discussed in [18], where the authors found the exact structural transition point. Furthermore, in [7], while studying epidemic spreading in multiplex networks, the authors verified this structural behavior in the analysis of the supra-adjacency matrix. Besides, it was also shown that it is intimately related to spreading processes and the layer-localization phenomena [7]. Thus, in the mentioned cases, the dynamical regimes can be understood as a consequence of the structural changes.
In summary, we have proposed a new mathematical formalism for the analysis of spectral properties in multiplex networks using the polynomial eigenvalue problem. This approach, we reduces the dimensionality of our matrices (coefficient matrices) at the cost of a higher order of the characteristic polynomial. This technique might seem counterintuitive at a first glance, but it reveals an underlying relationship between the eigenvalues of the matrices associated to multiplex structures. In contrast to single layer networks, multiplex networks are defined as matrix functions since we are interested in weighting inter and intra-layer edges differently. Thus, they depend on the coupling parameter, here denoted by p. Therefore, it is of utmost importance to derive spectral bounds as a function of p and to obtain insights on their asymptotic behavior.
We hope this work motivates the community to study in further details the structural behavior of multiplex and multilayer systems and their dynamical consequences. Besides, other matrices and processes might also be studied and evaluated. In this case, we believe that our formalism might also be helpful. Finally, we also hope to motivate studies on the analysis of eigenvectors, which is still lacking in the literature and are of great importance as they were shown to play a major role in dynamical processes [7,8].
Note that the Jordan matrix, J, the diagonal blocks are the Jordan blocks. Besides, if all the eigenvalues are simple, J is a diagonal matrix, where J ii = λ i . Additionally, observe that X ∈ R n×nl , where X = [x 1 , ..., x nl ] is composed by the right eigenvectors, while Y ∈ R nl×n , where Y = [y 1 , ..., y nl ] T is composed by its left eigenvectors. Finally, J ∈ R nl×nl .