Unique superdiffusion induced by directionality in multiplex networks

The multilayer network framework has served to describe and uncover a number of novel and unforeseen physical behaviors and regimes in interacting complex systems. However, the majority of existing studies are built on undirected multilayer networks while most complex systems in nature exhibit directed interactions. Here, we propose a framework to analyze diffusive dynamics on multilayer networks consisting of at least one directed layer. We rigorously demonstrate that directionality in multilayer networks can fundamentally change the behavior of diffusive dynamics: from monotonic (in undirected systems) to non-monotonic diffusion with respect to the interlayer coupling strength. Moreover, for certain multilayer network configurations, the directionality can induce a unique superdiffusion regime for intermediate values of the interlayer coupling, wherein the diffusion is even faster than that corresponding to the theoretical limit for undirected systems, i.e., the diffusion in the integrated network obtained from the aggregation of each layer. We theoretically and numerically show that the existence of superdiffusion is fully determined by the directionality of each layer and the topological overlap between layers. We further provide a formulation of multilayer networks displaying superdiffusion. Our results highlight the significance of incorporating the interacting directionality in multilevel networked systems and provide a framework to analyze dynamical processes on interconnected complex systems with directionality.

nonsymmetrical graph matrices, which mostly features complex and non-orthogonal spectrum. Consequently, spectrum related analysis on undirected multilayer networks might be proven as invalid or inaccurate. Few recent studies on directed multilayer networks reveal new and unexpected physical behaviors [16,19], as well as increased sensitivity to structural changes [20].
The contribution of this paper is to provide a theoretical framework for diffusive dynamics on multilayer networks with directed layers. We show that although certain undirected multilayer networks can result in enhanced diffusive (faster than the slowest layer) or superdiffusive (faster than any of the layers) regimes for large values of the interlayer coupling, the diffusive rate of the overall system never exceeds the diffusion on the equivalent aggregated system (the direct sum of each layer). We also show that certain directed multilayer networks can also show enhanced diffusive and superdiffusive regimes for large values of interlayer coupling. Remarkably, directed multilayer networks, differently from their undirected counterparts, can exhibit a unique superdiffusive behaviour for intermediate values of coupling, where diffusion dynamics are even faster than the corresponding aggregated system. This unique superdiffusion regime only emerges when certain structural conditions are satisfied. We analytically and numerically uncover conditions driving the unique superdiffusion in directed multilayer networks. In addition, we propose a model to configure directed multilayer networks that achieves the unique superdiffusion with tunable magnitude.
The paper is organized as follows. Diffusive dynamics on multilayer networks with direction is described in the Model section 2. The section 3 of Results presents the conditions for diffusion regimes below the diffusivity in the corresponding aggregated system in section 3.1 and for the superdiffusion regime above the diffusivity in the aggregated system in section 3.2, followed by a model for the construction of a multilayer network with superdiffusion in section 3.3. Section 4 concludes the work.

Model
Multilayer networks consist of nodes interacting both within the same layer and across different layers via intralinks and interlinks which encodes multiple types of interactions. In this study, we focus on a particular type of multilayer networks called multiplex networks, characterized by layers consisting of the same set of nodes, but possibly different connectivity (layer topology); and layers interacting with each other only via counterpart node. The topological structure of the multiplex is described by the socalled supra-adjacency matrix A = (a ij ), which is a block diagonal matrix. Each block is a N × N matrix, where N is the number of nodes per layer. Each diagonal block corresponds to the adjacency matrix of a layer (intra-layer connectivity, and therefore an entry a ij = 1 in such a block corresponds to an interlayer link), while the offdiagonal blocks, given the definition of a multiplex network, are just N × N identity matrices, representing the inter-layer links between replica nodes across layers. Note that a directed multiplex, is a multiplex where at least one of the diagonal blocks is not symmetric.
Let Q be the corresponding Laplacian matrix, defined as Q = D − A where D = diag (d i ) and d i = N i=1 a ij denoting the out-going degree of a node i. We consider without loss of generality a multiplex network consisting of two directed layers and interconnected via links of weight p ≥ 0, the corresponding Laplacian matrix Q can be written in a block form as where Q i = D i − A i , i = 1, 2, denotes the Laplacian matrix for the directed graph in each layer and pI encodes the weighted interconnections across layers. The above defined Laplacian matrix for multilayer networks with directed layers is non-symmetrical, which might lead to complex rather than real spectrum. As the Laplacian satisfies a zero row sum, i.e., Qu = 0 with u denoting the all one vector, value 0 is an eigenvalue with right eigenvector u and left eigenvector denoted as y. For a nonsymmetric matrix, left and right eigenvectors are in general not the same.
To analyze the effect of network directionality on dynamical processes, we consider the standard diffusive dynamics on multilayer networks with directed layers and determine the convergence rate to a fixed point solution. Let x i (t) denotes the state of a node i at time t. Changes of the state x i (t) with respect to time is governed bẏ The diffusive dynamics on a multilayer network can be written in a matrix form as where () T denotes transposition, and x(t) is the vector encoding the state of all nodes at time t. The solution of the above governing equation follows The fixed point x * T at whichẋ * T = 0 is calculated employing the conservation law x(0) T u = x(∞) T u, where all one vector u is the right eigenvector corresponding to eigenvalue 0 of the Laplacian matrix Q. Let y denote the left eigenvector of Q corresponding to eigenvalue 0. Combining the conservation law x(0) T u = x * T u and y T Q = 0 yields the solution of fixed point We show in Appendix A a convergence to the fixed point in Eq. (3) is guaranteed for a connected directed graph. The convergence rate to the fixed point x * is characterized by the second smallest real part [16], denoted as Reλ 2 (Q), of the eigenspectrum of the Laplacian matrix.

General analysis for diffusive behavior
To analyze the convergence behavior of diffusive dynamics on multiplex networks with directed layers, we study the spectrum of the governing Laplacian matrix and compare with the integrated system. The characteristic polynomial for the Laplacian matrix Q in Eq. (1) can be calculated by, applying the Schur complement theorem on the block Laplacian matrix, The eigen-pair of eigenvalue 2p and eigenvector (u, −u) holds for directed multilayer networks, which is previously found in undirected multilayered networks [11,21]. The joint effect of the Laplacian matrices Q 1 and Q 2 for each layer, encoded in the matrix ∆Q, determines the deviation of the convergence rate of the whole system from that of the integrated multilayer. Depending on relations between the diffusive rate of the system as a whole and that of the integrated system, we refer two complementary regimes, superdiffusive and non-superdiffusive, as follows: (i) if Reλ 2 (Q) of the overall system is greater than λ 2 ((Q 1 + Q 2 )/2) of the integrated system, we refer to as superdiffusive, (ii) otherwise as non-superdiffusive.
We mention that the superdiffusive regime defined in this work slightly differs from the definition in literature [14,16], where superdiffusion implies a faster diffusion on the overall system than the diffusion on the fastest layer. Here, we mainly focus on diffusive behavior comparing the overall system and the correspondent integrated system, given the later serves as the asymptotic diffusion as intercoupling strength goes to infinity and, most importantly, serves as a theoretical upper limit for diffusion on undirected multiplex networks (as shown in the following Eq. (8)). Under the special case of Q 1 = Q 2 , solving Eq. (4) yields Reλ 2 (Q) = min (λ 2 ((Q 1 + Q 2 )/2), 2p) and therefore only the non-superdiffusive regime is exhibited. General cases of Q 1 = Q 2 result in possibilities of both non-superdiffusive and superdiffusive regimes. Given the rich dynamical behavior of diffusion in directed multiplex, it is of paramount importance to identify the key structural properties underpinning each of the observed regimes. In this study, we mainly focus on identifying the underlying conditions driving the non-superdiffusive and superdiffusive regimes employing the principal submatrix approach and the Cauchy interlacing theorem [22] (or the Poincare separation theorem). In addition, we employ the eigenvalue property of a normal matrix, which satisfies XX T = X T X. For a normal matrix, the real parts of eigenvalues are the eigenvalues of the Hermitian part ReX = (X + X T )/2 of the matrix X.
The Laplacian matrix is similar to the following matrix Q in a transformed basis which can be simplified as The integrated Laplacian Q 1 +Q 2 2 appears as a principal submatrix of the transformed Laplacian Q of the whole system. In addition, eigenvalues of the matrix Q are the same with eigenvalues of Q due to the similarity relation between Q and Q. Applying the Cauchy interlacing theorem on the transformed matrix Q and the integrated system as a principal submatrix, undirected multilayer networks always satisfy the following interlacing relation Hence, there never occurs a superdiffusion in undirected multilayer networks, regardless of the strength of intercoupling. Eq. (8) serves as a theoretical limit and implies that although multilayer structure of undirected layers enhances the less diffusive layer, the diffusive rate is upper bounded by the integrated system.

Conditions for the non-superdiffusive regime
For multilayer networks with directed layers, the conditions for non-superdiffusive and superdiffusive regimes are translated into conditions when the Cauchy interlacing theorem between a non-symmetric matrix and its principle submatrix is satisfied. Starting from the special case of coupled identical layers, it follows which implies a non-superdiffusive regime for directed multilayer networks with identical directed layers for p ≥ λ 2 ((Q 1 + Q 2 )/2) /2. When the transformed Laplacian Q and all the principal matrices including the integrated matrix (Q 2 + Q 1 )/2 are normal, the interlacing is also satisfied. Because for a principally normal matrix, all eigenvalues of matrices Q and (Q 2 + Q 1 )/2 lie on the same complex line with an angle θ in the complex plane [23,24]. For a principally normal matrix, the partial order of generalized interlacing on a complex plain is obtained: As the real part of eigenvalue reads |λ 2 (Q) |e iθ and θ is the same for all eigenvalues, it implies an interlacing between real parts of eigenvalues and accordingly a nonsuperdiffusive regime for multilayer networks with principally normal Laplacian matrices.
When the requirement of normality of all principal submatrices is relaxed to the normality of a single principal submatrix, generalized interlacing of lexicographic order can be applied to Q and (Q 2 + Q 1 )/2. For a normal Laplacian matrix and a normal principal submatrix Q 1 +Q 2 2 , we have that where ≤ θ represents the lexicographic order [25] characterized by the positive cone H := {a+bi : a > 0 or a = 0 and b > 0}. Eq. (11) means that e −iθ λ 2 −e −iθ λ 2 (Q) lies in a positive cone, which implies λ 2 ≥ Reλ 2 (Q). Therefore, directed multilayer networks with normal Laplacian matrix and normal principal submatrix exhibit only the non-superdiffusive regime.

Conditions for the superdiffusive regime
Diffusion on undirected multilayer never exceeds the diffusion of integrated part. However, certain multilayer networks with directed layers break the constraint and exhibit a unique superdiffusion. In this subsection, we analyze conditions for the occurrence of such superdiffusion. For the generalized interlacing theorem to hold, it requires both the whole system and the integrated system to be normal or simultaneously close to normal. When the condition is relaxed such that the integrated system is normal or close to normal, while the whole system deviates from normality, it might occur with superdiffusion. We provide analytical evidences by firstly relating the real part of eigenvalues and eigenvalues of the Hermitian part of the original matrix, which follows where Re Q denotes the Hermitian part of Q. Because N j=1 Reλ j = N j=1 λ j Re Q and λ 1 Re Q = 0 for a normal subgraph, we have that Reλ 2 (Q) ≥ λ 2 Re Q . The equality is achieved if the matrix is normal. Therefore, if the integrated system is normal, it implies λ 2 . Secondly, we analyze the superdiffusion condition by relating eigenvalues of the Hermitian part of the integrated system and Hermitian part of the whole system. On one hand, observing that the Hermitian part of the integrated system coincides with the principal submatrix of the Hermitian part of the whole system, we have that λ 2 Re Q ≤ λ 2 Re Q 1 +Q 2 2 due to the Cauchy interlacing theorem for Hermitian matrices. On the other hand, A superdiffusion is thus established if (i) the integrated system is or close to normal and (ii) the structure difference of coupled layers ReQ 2 − ReQ 1 is negligibly small. The equality of λ 2 Re Q = λ 2 Re Q 1 +Q 2 2 is reached for a multilayer network consisting of a directed graph coupled with its reversed graph, Q 2 = Q T 1 . In this case, it follows ReQ 2 = ReQ 1 and there exists a superdiffusion of Reλ 2 (Q) > λ 2 . It was reported that the superdiffusion occurs when a fully connected network coupled with a network with a certain level of directionality, quantified by a metric called network directionality index [16]. Here, we prove in Appendix B an open problem of the normative property of the Network directionality index (NDI). In addition, we show that the condition of sufficient directionality as quantified by NDI is in line with the derived normality conditions in this study (as shown in Figure 4). The derived condition for superdiffusion in this work, i.e., a close to normal and a deviation from normal (upper left panel in Fig. 4) for the integrated and overall multiplex, respectively, is equivalent to the condition of a sufficient level of NDI. The nonnormality condition for the non-superdiffusive regime (lower right panel in Fig. 4) also agrees with the NDI condition. However, the nornormality condition is widely applicable for multiplex networks with any number of directed layers, compared to a single directed layer dealt by NDI condition.

Formulation of multilayer networks with superdiffusion
To interpret the normality conditions for superdiffusion, we propose a model to construct multiplayer networks exhibiting superdiffusion with a tunable level of magnitude. When the coupled layers have identical Hermitian part of Laplacian, the minimization of the normality level of the integrated structure is translated into with proof in Appendix C. Achieving superdiffusion by minimizing Eq. (14) is in general dauntingly difficult, even for undirected graphs of symmetrical Laplacian matrices. The simultaneous diagonalization, which always leads to the commutativity of matrices Q 1 and Q T 2 , i.e., Q 1 Q T 2 = Q T 2 Q 1 , was posted as an open problem by Hiriart-Urruty [27] more than a decade ago. Extensive research is performed due to the close association with quadratically constrained quadratic programming and certain advances are made for real symmetrical matrices [28]. However, results on nonsymmetrical matrices are extremely scarce. Nonetheless, we propose a construction model to generate a directed multilayer network by constructing one matrix (e.g., Q 2 ) from a given matrix (e.g., Q 1 ), such that superdiffusion occurs.
We decompose the construction model into two parts: (i) constructing a graph G 2 from the replication of graph G 1 but after reversing the direction of directed links; (ii) preserving the out-degree of each node in graph G 1 such that graph Laplacian Q 2 has the same diagonal elements as Q 1 . To this end, we propose to reverse links in a cycle fashion, i.e., simultaneously reversing the direction of all links involved in a cycle. In addition, the number of links reversed in a cycle fashion is closely associated with the magnitude Reλ 2 (Q) − λ 2 of superdiffusion. For two directed graphs G 1 and G 2 , we define a variable q as direction overlap, the fraction of links having the same direction and calculated by q = |E(G 1 )∩E(G 2 )| |E(G 1 )∪E(G 2 )| or in a matrix form as The denominator normalizes the direction overlap, in which q = 0 corresponds to Q 2 = Q T 1 and q = 1 corresponds to Q 1 = Q 2 . The magnitude of superdiffusion can thus be tuned by reversing different fraction of links, parameterized by 1 − p. Figure 3 shows diffusive behaviors on constructed multiplex networks consisting of (a) a circulant directed graph and the duplicate digraph with 1 − q reversed links and (b) a real-world genetic and protein interactions network of the Epstein-Barr virus [29,30] and the duplicate digraph with 1 − q reversed links. The coupled digraph and its duplicate without link reversing (q = 1) displays a non-superdiffusive behavior. For 0 ≤ q < 1, the constructed multiplex built upon both synthetic and real-world directed graph shows superdiffusion regimes. In addition, a higher fraction of reversed links (a smaller q) results in a higher magnitude of superdiffusion.

Conclusions
In this paper, we study the diffusive dynamics on multiplex networks consisting of at least one directed layers. Though multilayer structure of undirected layers enhances diffusion, it is restricted to a theoretical limit by the aggregation of each layer. We show that ubiquitous property of directionality in each layer could break this limit and achieve a unique superdiffusion regime, wherein diffusion is even faster than the corresponding aggregated system. We analytically and numerically uncover that the unique superdiffusion is driven by the directionality and the underlying structure of layers rather than how strongly layers are coupled. In particular, diffusive behaviors are associated with the non-normality level of the integrated and the overall system: when both the integrated and the overall system are normal, it is assured that the multiplex exhibits only a non-superdiffusive regime; when the integrated system is normal or close to normal while the whole system is deviated from normal, it exhibits the unique superdiffusive regime. Additionally, we provide a model to construct multiplex networks, achieving superdiffusion with a tunable level of magnitude. We show that directionality induces a unique superdiffusion, a regime not existed in the counterpart of undirected multilayer, and suggest the key role of network directionality in analyzing diffusive dynamics on real-world systems which oftentimes are structurally directed and multilayered.

Interlayer coupling
Interlayer coupling Interlayer coupling Smallest nonzero eigenvalue Smallest nonzero eigenvalue Smallest nonzero eigenvalue a) b) c) ). Panel (b) shows a multiplex network with coupled layers in a reversed direction, exhibiting a superdiffusion (Reλ 2 (Q) > λ 2 ) after the transition from a linear growth of 2p. Panel (c) shows a multiplex with coupled layers in the same direction, exhibiting a non-superdiffusive regime (Reλ 2 (Q) = λ 2 ).

Appendix A. Convergence to the fixed point
The convergence to the fixed point is characterized by ∆x(t) where λ min is the minimum eigenvalue of the Laplacian Q. For a connected directed graph, λ min equals to zero which means a guaranteed convergence to the steady state x * calculated by Eq. (3).
Appendix B. Proof that network directionality index is normalized.
In the appendix, we prove that the metric quantifying the level of directionality is normalized to 1. Rewrite the definition [16] of Network Directionality Index (NDI) as or equivalently, Interlayer coupling Interlayer coupling  : The normality condition is in line with the condition of sufficient directionality for the occurrence of superdiffusion. The size of each circle is proportional to the value of network directionality index, NDI. The simulation is performed on a two-layered network of N = 200 nodes wherein one layer consists of a fully connected, undirected graph and the other layer consists of a directed graph whose NDI is modified progressively. The intercoupling weight is p = 200. A high level of NDI is translated to a low deviation from normality of the integrated network and a high deviation of normality of the multiplex network as a whole.
where ∆d ij = |d i→j − d j→i |.
Theorem 1. For a directed cycle network consisting of N nodes and N directed links in a clockwise or counterclockwise direction, the NDI can be simplified as where N j=1 |d i→j − d j→i |/N is the average path asymmetry and N j=1 (d i→j + d j→i ) /2N is the average path length starting from an arbitrary node i and arriving at each of the remaining N − 1 nodes.
Proof. For a directed circulant graph of N nodes and N directed links in a clockwise direction, the shortest directed path from node i to node j and back to i forms a closed cycle of length N , regardless of choices of pairs of nodes (i = j). Thus, we have that Applying the symmetry of ∆d ij = ∆d ji into the numerator of (B.2), we have that For a directed cycle consisting of N nodes and N directed links in a clockwise direction, the sum N j=1 ∆d ij is the same for each node i and we arrive at Theorem 2. Given a directed cycle consisting of N nodes and N directed links in a clockwise direction, adding a single or multiple directed links between any pair of nodes i and j generates shortcuts either in forward (i → j) or backward (j → i) direction between that pair of nodes and thus forms k (k ≥ 1) sub-cycles with each consisting of c k (c k ≤ N ) nodes and c k links. The NDI after adding directed links is smaller than or equal to the NDI of the original directed cycle.
Proof. Adding a directed link connecting node j to i forms a sub-cycle k consisting of c k nodes and c k directed links in a clockwise direction, as shown in Figure B2. For a pair of nodes i and j in the sub-cycle k, the difference of path asymmetry ∆d ij before and after adding directed links reads It has been shown [16] that a directed cycle of N nodes has NDI < 1. Subtracting the same number j N − (m) j both in the numerator and denominator of (B.2) does not increase the NDI of the original cycle consisting of N nodes.