Unique superdiffusion induced by directionality in multiplex networks

PAPER • OPEN ACCESS Unique superdiffusion induced by directionality in multiplex networks To cite this article: Xiangrong Wang et al 2021 New J. Phys. 23 013016 View the article online for updates and enhancements. This content was downloaded from IP address 155.210.59.210 on 04/03/2021 at 12:21 New J. Phys. 23 (2021) 013016 https://doi.org/10.1088/1367-2630/abdb71 OPEN ACCESS RECEIVED 3 November 2020 REVISED 7 January 2021 ACCEPTED FOR PUBLICATION 13 January 2021 PUBLISHED 29 January 2021 Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. PAPER Unique superdiffusion induced by directionality in multiplex networks Xiangrong Wang1,2, Alejandro Tejedor3,4,YiWang 1,2and Yamir Moreno5,6,7,∗ 1Institute of Future Networks, Southern University of Science and Technology, Shenzhen, People’s Republic of China 2Research Center of Networks and Communications, Peng Cheng Laboratory, Shenzhen, People’s Republic of China 3Department of Science and Engineering, Sorbonne University Abu Dhabi, Abu Dhabi, United Arab Emirates 4Department of Civil and Environmental Engineering, University of California Irvine, Irvine, CA, United States of America 5Institute for Biocomputation and Physics of Complex Systems, University of Zaragoza, Zaragoza, Spain 6Department of Theoretical Physics, University of Zaragoza, Zaragoza, Spain 7ISI Foundation, Turin, Italy ∗Author to whom any correspondence should be addressed. E-mail: yamir.moren[email protected]om Keywords: networks, multiplex, directionality, complex systems, spectral properties Abstract 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. 1. Introduction Multilayer networks provide a proper mathematical representation of complex systems with different types of interactions, e.g. social interactions among individuals across different social platforms such as Facebook and Twitter, and enable the understanding of dynamics acting on or evolving within those systems [1,2]. A number of dynamical processes are studied on the framework of multilayer networks: for example, innovation/information diffusion [3], the control of formation and task allocation of swarms [4–6], synchronization of oscillators [7], the understanding of functional brain connectivity [8], and many others. More specifically, recent studies based on a multilayer network representation have revealed unforeseen dynamical regimes and behaviors, such as new spreading regimes when more than one diseases in spreading in a given population [9–12], and enhanced stability of synchronized states [13]. Moreover, conversely from building up a multilayer structure, gradually dismantling a multilayer network triggers an abrupt transition [14], suggesting a non-additive effect of a multilayer structure from the integration of its layers. In terms of diffusion dynamics, multilayer network studies revealed that for any multilayer configuration, the diffusive behavior of the overall system (multilayer network) could be tuned to be faster than that of the slowest layer © 2021 The Author(s). Published by IOP Publishing Ltd on behalf of the Institute of Physics and Deutsche Physikalische Gesellschaft New J. Phys. 23 (2021) 013016 XWanget al for significant values of interlayer diffusivity (in comparison to intralayer diffusivity). More surprisingly, but only for certain multilayer network configurations, the overall system can exhibit a superdiffusive behavior for large values of interlayer diffusivity, where the diffusive dynamics on the multilayer network is faster than in any of the individual layers, when those are considered in independently [15–17]. Most of the studies made in the field of multilayer networks adopt undirected structures, i.e. the interactions represented by both interlayer and intralayer links are assumed to be undirected. However, most real world systems are inherently directed. Examples are the World Wide Web, in which hyperlinks runinonedirectionfromoneWebpagetoanother;foodwebs,inwhichenergyflowsfrompreyto predator; phone call networks; metabolic networks; citation networks; social networks of followers and followee; gene regulation networks, to name a few. Additionally, some dynamical processes are governed by physical laws and inherently evolve in a directional manner despite of the underlying undirected connectivity topology, like current/information flow from high voltage to low voltage, virus transmissions from infected to susceptible individuals. Moreover, multi-leveled unidirectional interactions are commonly found in interacting systems. For example, a gene regulation network aims to properly characterize both protein-DNA and protein-protein directed interactions between genes and proteins [18]. Therein, multilayer networks with directionality provide a realistic representation for such systems consisting of multilevel and directed interactions. Multilayer networks where at least one of the layers is directed give rise to new theoretical challenges in analyzing their dynamical behaviors. As diffusion dynamics on undirected multilayer networks is fully characterized by the spectrum of network matrices, like supra-Laplacian [15], network directionality leads to 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 behavior 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.Thesection3of 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.Section4concludes the work. 2. 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 so-called supra-adjacency matrix A=aij, which is a block diagonal matrix. Each block is an N×Nmatrix, where Nis the number of nodes per layer. Each diagonal block corresponds to the adjacency matrix of a layer (intra-layer connectivity, and therefore an entry aij =1 in such a block corresponds to an interlayer link), while the off-diagonal blocks, given the definition of a multiplex network, are just N×Nidentity 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 Qbe the corresponding Laplacian matrix, defined as Q=D−Awhere D=diag (di)and di=N i=1aij 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 2 New J. Phys. 23 (2021) 013016 XWanget al Laplacian matrix Qcan be written in a block form as Q=Q1+pI −pI −pI Q2+pI(1) where Qi=Di−Ai,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 =0withudenoting the all one vector, value 0 is an eigenvalue with right eigenvector uand 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 xi(t) denotes the state of a node iat time t. Changes of the state xi(t)withrespectto time is governed by ˙ xi(t)=−dixi(t)+iaijxi(t). The diffusive dynamics on a multilayer network can be writteninamatrixformas dx(t)T dt=−x(t)TQ(2) where ()Tdenotes transposition, and x(t) is the vector encoding the state of all nodes at time t. The solution of the above governing equation follows x(t)T=x(0)Te−Qt.Thefixedpointx∗Tat which ˙ x∗T=0is calculated employing the conservation law x(0)Tu=x(∞)Tu,whereallonevectoruis the right eigenvector corresponding to eigenvalue 0 of the Laplacian matrix Q.Letydenote the left eigenvector of Q corresponding to eigenvalue 0. Combining the conservation law x(0)Tu=x∗Tuand yTQ=0yieldsthe solution of fixed point x∗=x(0)Tu yTuy.(3) We show in appendix Aa convergence to the fixed point in equation (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. 2.1. 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 Qin equation (1) can be calculated by, applying the Schur complement theorem on the block Laplacian matrix, det (λI−Q)=det λI−Q1+Q2 2(λ−2p)I−Q1+Q2 2−ΔQ(4) where 4ΔQ=(Q1−Q2)2−2(Q1Q2−Q2Q1). The eigen-pair of eigenvalue 2pand eigenvector (u,−u) holds for directed multilayer networks, which is previously found in undirected multilayered networks [12,21]. The joint effect of the Laplacian matrices Q1and Q2for 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((Q1+Q2)/2) of the integrated system, we refer to as superdiffusive, (ii) otherwise as non-superdiffusive. Re λ2(Q)⎧ ⎨ ⎩ >λ 2((Q1+Q2)/2) Superdiffusive ⩽λ2((Q1+Q2)/2) Non-superdiffusive (5) We mention that the superdiffusive regime defined in this work slightly differs from the definition in literature [15,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 equation (8)]. Under the special case of Q1=Q2, solving equation (4)yieldsReλ2(Q)=min λ2((Q1+Q2)/2), 2p and therefore only the non-superdiffusive regime is exhibited. General cases of Q1=Q2result in possibilities of both non-superdiffusive and superdiffusive regimes. Figure 1exemplifies three typologies 3 New J. Phys. 23 (2021) 013016 XWanget al Figure 1. Superdiffusive and non-superdiffusive regimes on multiplex networks with at least one directed layers. Panel (a) shows a multiplex network consisting of a directed layer and an undirected layer displaying a non-superdiffusive regime (Re λ2(Q)<λ 2Q1+Q2 2). Panel (b) shows a multiplex network with coupled layers in a reversed direction, exhibiting a superdiffusion (Re λ2(Q)>λ 2Q1+Q2 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Q1+Q2 2). including (a) coupled directed and undirected layers, (b) coupled directed layers with the opposite direction, and (c) coupled directed layers with the same direction. Figure 1(a) and (c) show only non-superdiffusive regimes and figure 1(b) exhibits superdiffusion. Figures 2(a) and (c) show the non-superdiffusive and superdiffusive behavior on real-world multiplex networks, respectively. 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 XXT=XTX. For a normal matrix, the real parts of eigenvalues are the eigenvalues of the Hermitian part Re X=(X+XT)/2ofthematrixX. The Laplacian matrix is similar to the following matrix  Qin a transformed basis  Q=1 √2II −II QI−I II 1 √2(6) which can be simplified as  Q=⎡ ⎢ ⎣ Q1+Q2 2 Q2−Q1 2 Q2−Q1 2 Q1+Q2 2+2pI ⎤ ⎥ ⎦.(7) The integrated Laplacian Q1+Q2 2appears as a principal submatrix of the transformed Laplacian  Qof the whole system. In addition, eigenvalues of the matrix  Qare the same with eigenvalues of Qdue to the similarity relation between Qand  Q. Applying the Cauchy interlacing theorem on the transformed matrix  Qand the integrated system Q1+Q2 2as a principal submatrix, undirected multilayer networks always satisfy the following interlacing relation λ2(Q)⩽λ2Q1+Q2 2.(8) Hence, there never occurs a superdiffusion in undirected multilayer networks, regardless of the strength of intercoupling. Equation (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. 4 New J. Phys. 23 (2021) 013016 XWanget al Figure 2. Non-superdiffusive and superdiffusive real-world multiplex networks. Panel (a) shows a non-superdiffusive regime on the Vickers–Chan multiplex social networks [31] with directed layers. Panel (c) shows a superdiffusion on the multilayer network where one layer is from Vickers–Chan multiplex and the other layer is the reverse of the first layer. Panels (b) and (d) show the distribution of all eigenvalues of the Laplacian Qfor the overall multiplex and that of the integrated system (Q1+Q2)/2 corresponding to networks in panels (a) and (c), respectively. 3. Results 3.1. 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 Re λ2(Q)=λ2Q1+Q2 2(9) which implies a non-superdiffusive regime for directed multilayer networks with identical directed layers for p⩾λ2(Q1+Q2)/2/2. When the transformed Laplacian  Qand all the principal matrices including the integrated matrix (Q2+Q1)/2 are normal, the interlacing is also satisfied. Because for a principally normal matrix, all eigenvalues of matrices Qand (Q2+Q1)/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: |λ2(Q)|⩽λ2Q1+Q2 2.(10) As the real part of eigenvalue reads |λ2(Q)|eiθand θis the same for all eigenvalues, it implies an interlacing between real parts of eigenvalues and accordingly a non-superdiffusive 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  Qand (Q2+Q1)/2. For a normal Laplacian matrix and a normal principal submatrix Q1+Q2 2,wehavethat λ2(Q)⩽θλ2Q1+Q2 2(11) where ⩽θrepresents the lexicographic order [25] characterized by the positive cone H:={a+bi : a>0ora=0andb>0}.Equation(11)meansthate −iθλ2Q1+Q2 2−e−iθλ2(Q)lies in a positive cone, 5 New J. Phys. 23 (2021) 013016 XWanget al which implies λ2Q1+Q2 2⩾Re λ2(Q). Therefore, directed multilayer networks with normal Laplacian matrix and normal principal submatrix Q1+Q2 2exhibit only the non-superdiffusive regime. 3.2. 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 N  j=N−k Re λj= N  j=N−k x∗ jRe  Qxj⩽ N  j=N−k λjRe  Q(12) where Re  Qdenotes the Hermitian part of matrix  Q.BecauseN j=1Re λ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.Theequality, λ2Q1+Q2 2=λ2Re Q1+Q2 2, is achieved if the integrated system is normal. 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 Q1+Q2 2due to the Cauchy interlacing theorem for Hermitian matrices. On the other hand, λ2Re  Q⩾λ2Re Q1+Q2 2+λ2Re Q2−Q1 2due to Ky Fan majorization [26]theorem λ(A+B)≺λ(A)+λ(B). Therefore, if the difference Re (Q2−Q1)2is bounded, it leads to λ2Re  Q≈λ2Re Q1+Q2 2. A superdiffusion is thus established Re λ2(Q)λ2Q1+Q2 2(13) if (i) the integrated system is or close to normal and (ii) the structure difference of coupled layers Re Q2−Re Q1is negligibly small. The equality of λ2Re  Q=λ2Re Q1+Q2 2is reached for a multilayer network consisting of a directed graph coupled with its reversed graph, Q2=QT 1. In this case, it follows Re Q2=Re Q1and there exists a superdiffusion of Re λ2(Q)>λ 2Q1+Q2 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 Ban 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 3). The derived condition for superdiffusion in this work, i.e. a close to normal and a deviation from normal (upper left panel in figure 3) 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 figure 3)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. Though the results are demonstrated for two-layer multiplex networks, numerical results on multiplex networks with a larger number of layers show analogous superdiffusive and non-superdiffusive behaviors, also depending on the nonnormality and NDI conditions of the structure of each layer. However, in those cases it is far more challenging to provide a complete mathematical proof due to the difficulty in isolating the Laplacian of each layer from the block Laplacian. 3.3. 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 6 New J. Phys. 23 (2021) 013016 XWanget al Figure 3. 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. Figure 4. The difference Re λ2(Q)−λ2Q1+Q2 2as a function of the interlayer coupling strength pfor (a) a multilayer network consisting of directed circulant graph of 100 nodes and average degree 4 and the constructed layer by the proposed model; and (b) a multilayer network consisting of a genetic and protein interactions network of the Epstein–Barr virus network and the constructed layer. translated into min Q1QT 2−QT 2Q12(14) with proof in appendix C. Achieving superdiffusion by minimizing equation (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 Q1and QT 2,i.e.Q1QT 2=QT 2Q1, 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. Q2) from a given matrix (e.g. Q1), such that superdiffusion occurs. We decompose the construction model into two parts: (i) constructing a graph G2from the replication of graph G1but after reversing the direction of directed links; (ii) preserving the out-degree of each node in graph G1such that graph Laplacian Q2has the same diagonal elements as Q1. To ensure the second smallest real part of the Laplacian eigenvalue to be nonzero, we assume a single requirement of graph G1to be a strongly connected directed graph. 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Q1+Q2 2of superdiffusion. For two directed graphs G1and G2,wedefineavariableqas direction overlap, the fraction 7 New J. Phys. 23 (2021) 013016 XWanget al of links having the same direction and calculated by q=|E(G1)∩E(G2)| |E(G1)∪E(G2)|or in a matrix form as q=i,j(A1)ij(A2)ij i,j(A1)ij +i,j(A2)ij −i,j(A1)ij(A2)ij .(15) The denominator normalizes the direction overlap, in which q=0 corresponds to Q2=QT 1and q=1 corresponds to Q1=Q2. The magnitude of superdiffusion can thus be tuned by reversing different fraction of links, parameterized by 1 −p. Figure 4shows diffusive behaviors on constructed multiplex networks consisting of (a) a circulant directed graph and the duplicate digraph with 1 −qreversed links and (b) a real-world genetic and protein interactions network of the Epstein–Barr virus [29,30] and the duplicate digraph with 1 −qreversed 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)resultsina higher magnitude of superdiffusion. 4. 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. Acknowledgments XW acknowledges the project 62003156 supported by NSFC and project ‘PCL Future Greater-Bay Area Network Facilities for Large-scale Experiments and Applications (LZC0019)’. AT acknowledges partial support from NSF Grants EAR-181190 and the UK Research and Innovation Global Challenges Research Fund Living Deltas Hub Grant NES0089261. YM acknowledges partial support from the Government of Arag´ on, Spain through a grant to the group FENOL (E36-20R), by MINECO and FEDER funds (Grant FIS2017-87519-P) and by Intesa Sanpaolo Innovation Center. The funders had no role in study design, data collection and analysis, or preparation of the manuscript. Appendix A. Convergence to the fixed point The convergence to the fixed point is characterized by Δx(t)T=(x(t)−x∗)T.Therateofconvergenceis regulated by ˙ Δx(t)T=d(x(t)−x∗)T dt=−(Δx(t)T+x∗T)Q=−Δx(t)TQ.ThenormΔx(t)Tof the difference vector is bounded by Δx(t)T=Δx(0)Te−Qt⩽Δx(0)Te−Qt.(A.1) If one mainly concerns with the convergence behavior for the decrease of Δx(t), the following holds for a bounded operator Qin a Hilbert space lim t→∞ log e−Qt t=−λmin (A.2) 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 equation (3). 8