{'fraction_non_alphanumeric': 0.03291344281504713, 'fraction_numerical': 0.010925195154314074, 'mean_word_length': 4.566084788029925, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 18, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Diffusion Dynamics and Optimal Coupling in Multiplex Networks with Directed Layers

Multiplex networks have been intensively studied during the last few years as they offer a more realistic representation of many interdependent and multilevel complex networked systems. However, even if most real networks have some degree of directionality, the vast majority of the existent literature deals with multiplex networks where all layers are undirected. Here, we study the dynamics of diffusion processes acting on coupled multilayer networks where at least one layer consists of a directed graph; we call these directed multiplex networks . We reveal a new and unexpected signature of diffusion dynamics on directed multiplex networks, namely, that different from their undirected counterparts, they can exhibit a nonmonotonic rate of convergence to steady state as a function of the degree of coupling, resulting in a faster diffusion at an intermediate degree of coupling than when the two layers are fully coupled. We use synthetic multiplex examples and real-world topologies to illustrate the characteristics of the underlying dynamics that give rise to a regime in which an optimal coupling exists. We further provide analytical and numerical evidence that this new phenomenon is solely a property of directed multiplex, where at least one of the layers exhibits sufficient directionality quantified by a normalized metric of asymmetry in directional path lengths. Given the ubiquity of both directed and multilayer networks in nature, our results have important implications for studying the dynamics of multilevel complex systems. DOI:


I. INTRODUCTION
Studying the structural and process-based dynamics of complex systems operating in natural, biological, climatic, and social environments using network-based theories has been a popular area of research over the past decades, often revealing dynamics and emergent behavior not accessible otherwise [1][2][3][4][5][6]. Lately, the realization that many systems consist of processes that operate on many levels of connectivity and dynamics has spurred research on the so-called multilayer networks, i.e., a class of multidimensional networks which includes networks of networks, interdependent networks, multiplex networks, etc. [7,8]. Multiplex networks are coupled multilayer networks where each layer consists of the same set of nodes but possibly different topologies and layers interact with each other only via counterpart nodes in each layer [7][8][9]. Multiplex networks have been shown useful for the study of diverse processes including social [10][11][12], transportation [13], and biochemical networks [14,15], among others. Recent studies have shown that the coupling of the layers in a multiplex network can result in emergent structural [16] and dynamical behavior such as enhanced diffusion (superdiffusion) [17], increased resilience to random failure [18], and emergence of critical points in the dynamics of coupled spreading processes [10,12]. These richer dynamics arise as a direct consequence of the emergence of more paths between every pair of nodes brought about by layer switching via an interlayer link.
Most of the theory for multiplex networks has been developed when all layers consist of undirected networks [9]. However, more often than not, real social, biological, and natural networks are structurally directed. Additionally, even if the underlying topology is not directed, the functional and dynamical connectivity of undirected networks is often directed due to gradients or the directionality in the flow of mass or information. Examples include geophysical processes on tributary river networks [19][20][21][22] and river delta channel networks [23][24][25][26][27][28], ecohydrology [29,30], food webs [31][32][33], gene regulation networks [31,34], and social dynamics [35], to name a few.
In this paper, we study diffusion processes on multiplex networks wherein the connectivity of at least one of the layers forms a directed graph; we term them directed multiplex networks. We document for the first time a nonmonotonic increase in the rate of convergence to the steady state as a function of the degree of coupling between layers, which translates to the emergence of a prime regime for intermediate values of coupling where transport processes are enhanced and diffusion is faster than when layers are fully coupled. Within the prime regime, an optimal coupling exists for which the convergence time to the steady state of the diffusion process is minimal (fastest transport). We demonstrate that this behavior is inherent to directed multiplexes by providing a general analytical proof of the monotonically increasing nature of the rate of convergence of the diffusion process as a function of interlayer coupling for the case of undirected multiplexes. From a phenomenological point of view, the new phenomenon unveiled can be traced back to the anisotropy introduced by the directionality of the edges (of at least one of the layers) and the emergence of directional paths at the network scale, which is the main factor that induces much richer dynamics than in the undirected counterpart. The combination of two factors-(i) a drastically different behavior of the expected system dynamics of directed vs undirected multiplexes for a range of between-layer levels of coupling and (ii) the fact that intermediate coupling scenarios are most likely to dominate the interaction between many intralayer processes (e.g., multiplatform transportation systems)-make the results of this study fundamental to understanding the emergent system response in many natural and engineered systems where directionality and multiprocess connectivity hold.
The rest of the paper is organized as follows. We first describe diffusion processes in both undirected and directed networks. For the former case, we first show analytically that it is not possible to observe a regime in which the diffusion is faster than when the layers of the multilayer system are fully coupled. Next, we describe the new phenomenology numerically for a few synthetic network configurations. We prove analytically the existence of the nonmonotonic behavior in the rate of convergence to steady state and the existence of the prime regime. We show that the emergence of the prime regime and the optimal coupling requires sufficient network-scale directionality in at least one of the layers of the directed multiplex, as quantified by the proposed network directionality index (NDI). Finally, we illustrate the validity and potential of our results by analyzing two real-world topologies. We round off the paper by discussing the phenomenological implications of our results and providing concluding remarks.

II. DIFFUSIONLIKE DYNAMICS ON NETWORKS
Let xðtÞ represent the N × 1 vector of concentration associated with the N nodes of a network at time t (throughout, vectors are thought of as column vectors). The diffusion dynamics on an undirected single-layer network (monoplex) can be described by where the vector _ xðtÞ is the temporal derivative of xðtÞ, and D is a scalar that represents the diffusion constant. The L ¼ ½l ij N i;j¼1 matrix is the Laplacian of the network, which is defined as L ¼ S − W. Here, the N × N matrix W is the weighted adjacency matrix, whose entries w ij ¼ w ji ∈ R þ represent the strength of the connectivity between nodes i and j, and S is the N × N diagonal matrix with diagonal entries s ii ¼ P N j¼1 w ij . For an unweighted network where the entries of W are binary (w ij ¼ 1 if there exists a link between nodes i and j, and 0 otherwise), the matrix S is the degree matrix and its diagonal entries s ii correspond to the number of links connected to node i, and the Laplacian is known as the combinatorial Laplacian (see Ref. [36] and references therein).
The negative Laplacian −L (or more generally, −DL, when the diffusion constant D ≠ 1) can be interpreted as the transition-rate matrix of a continuous-time Markov chain (CTMC); its entries l ij represent the rate at which the transition from node i to node j takes place [36,37]. The dynamics of the corresponding continuous-time random walk (CTRW) is governed by where ðÞ T denotes transposition, and pðtÞ is a vector whose ith component is the probability that the CTRW visits node i at time t. Note that for undirected networks, L is symmetric, and, therefore, Eqs. (1) and (2) are equivalent for D ¼ 1. Thus, the stationary distribution ½ _ pðtÞ T ¼ 0 of the CTRW is identical to the stationary solution of the diffusion process p s T ¼ 1 N ð1; 1; 1; …; 1Þ, which is the eigenvector of L that corresponds to the zero eigenvalue. The characteristic timescale of convergence to that solution is inversely proportional to the smallest nonzero eigenvalue λ 2 of L assuming that the graph is connected; i.e., there exists a single component and, therefore, λ 2 > 0 [38,39].
For directed networks, the weighted adjacency matrix W is not symmetric; i.e., the strength w ij of the directed link starting at node i and ending at node j may differ from w ji . Consequently, the Laplacian is not symmetric either. Two different Laplacians can be naturally defined by extension of the undirected definition: (i) out-Laplacian where entries of the diagonal out-and in-strength matrices are defined as s out ii ¼ P N j¼1 w ij and s in jj ¼ P N i¼1 w ij , respectively. Both Laplacians (out and in) have 1 (column vector with all entries equal to 1) as the eigenvector corresponding to the zero eigenvalue, making it tempting to assume that _ xðtÞ ¼ −D · L out xðtÞ as the equation governing the diffusion process in directed networks. However, due to the asymmetry of L out , 1 T is not in general a left eigenvector of L out for the zero eigenvalue, and, therefore, continuity (conservation of mass) is not guaranteed in the process ( d dt P x . On the other hand, it is meaningful to interpret −L out as the transition-rate matrix of a CTMC in a directed network, where the dynamics of the random walk are described by In this case, continuity of the process (conservation of probabilities or mass) is assured since 1 is a right eigenvector corresponding to the zero eigenvalue of −L out . Furthermore, if the directed network is strongly connected (i.e., there exists a directed path between every pair of nodes), a unique stationary distribution p T s exists (unique left eigenvector of the −L out corresponding to eigenvalue zero), which is a proper probability vector (the Perron-Frobenius theorem guarantees that all the vector entries are positive given that the network is strongly connected), The spectrum of L out is in general complex, and the convergence towards the stationary distribution is exponential (asymptotic) with rate Reðλ 2 Þ where λ 2 is the eigenvalue with the smallest nonzero real part [36,38].

III. UNDIRECTED MULTIPLEX
Gómez et al. [17] generalized the study of the characteristic timescale of diffusive processes to (undirected) multilayer networks for different degrees of coupling across two different layers. The three main assumptions in Ref. [17] were (1) the same set of nodes, albeit with different connectivity, forms the networks at each layer, (2) the connectivity in each layer (intralayer connectivity) consists of undirected networks forming a single-connected component, and (3) the interlayer connectivity consists of undirected links between counterpart nodes in the different layers, i.e., multiplex. The authors showed that the characteristic time of convergence of the diffusion process τ is inversely proportional to the smallest nonzero eigenvalue Λ 2 of the Laplacian of the multiplex called supra-Laplacian L (i.e., τ ∼ 1 Λ 2 ). Let L be the supra-Laplacian of a two-layer multiplex, each layer consisting of a singlecomponent undirected network [17]. L can be written as where L 1 , L 2 are the Laplacian matrices of the respective layers, and D X is the interlayer diffusion coefficient, which controls the degree of coupling between the layers. We are interested in the spectrum of the supra-Laplacian, i.e., the set of eigenvalues Λ i of L, and more specifically, in its second smallest eigenvalue Λ 2 , as this dictates the rate of convergence to the steady-state dynamics. Here, we prove the monotonic increase of Λ 2 as a function of D X . To this end, we use a transformation of L in a proper basis, which also allows us to obtain the asymptotic behavior of Λ 2 as D X → 0 and as D X → ∞. These asymptotic behaviors were previously derived in Gómez et al. [17], but we obtain them again for the sake of completeness and to capitalize on the fact that the procedure used here is applicable to both directed and undirected multiplexes. LetL be the transformation of L in a new basis, where L þ ¼ 1 2 ðL 1 þ L 2 Þ and L − ¼ 1 2 ðL 2 − L 1 Þ. Notice that L andL are similar matrices, and, therefore, they have the same eigenvalues. The spectrum of L is obtained as the set of solutions to the characteristic equation jL − ΛIj ¼ 0.
Using the Schur complement [40], we can rewrite this determinant as follows: From Eq. (7), we can now deduce the asymptotic behavior of the spectrum ofL and, more specifically, its second smallest eigenvalue Λ 2 as D X → 0 and as D X → ∞. These two cases are explained below.
(i) Case D X → 0: When D X ¼ 0, the multiplex (as a whole) has two components (assuming that the networks in each layer are strongly connected), one per layer, and, therefore, Λ 2 ¼ 0 is an eigenvalue of the supra-Laplacian with multiplicity 2. Thus, from the decomposition of the determinant in Eq. (7), we conclude that each factor has a root at Λ ¼ 0 with multiplicity 1 (note that L þ has a single null eigenvector 1 and sinceL must have two, the second necessarily comes from the second factor). For any value of D X > 0, Λ ¼ 0 has multiplicity 1 because now the two layers are connected to each other. From Eq. (7), the smallest eigenvalue of jðL þ þ 2D X I − ΛIÞj now grows linearly being equal to 2D X . The zero eigenvalue coming from the second factor jðL þ − ΛIÞ − L − ðL þ þ 2D X I − ΛIÞ −1 L − j remains zero independent of the value of D X , since there must be always a zero eigenvalue forL, and this cannot be supplied by the first factor. Therefore, by continuity, we can conclude that for small D X , Λ 2 ¼ 2D X is the smallest nonzero eigenvalue ofL. (ii) Case D X → ∞: Examining for this case the two factors of the decomposition in Eq. (7), we note that for D X → ∞, in the second factor, the term L − ðL þ þ 2D X I − ΛIÞ −1 L − is negligible vs ðL þ − ΛIÞ, so effectively the second factor reduces to ðL þ − ΛIÞ. Thus, the spectrum in this case splits into two groups: (1) the roots of jðL þ − ΛIÞj, which are the eigenvalues of L þ (Laplacian of the superposition of the graphs in each layer scaled by a factor 0.5), and (2) the roots of jL þ þ 2D X I − ΛIj, which are the eigenvalues of L þ shifted by 2D X . Thus, the eigenvalues in the second group are just a shifted version of the first group of eigenvalues tending to infinity linearly when D X → ∞. Regarding Λ 2 , the above arguments show that Λ 2 ¼ 2D X for small values of D X , while it tends to the smallest nonzero eigenvalue of L þ asymptotically as D X → ∞. These results are equivalent to those reported and proven by Ref. [17]. Here, moreover, we proceed to prove that for undirected multiplex (i.e., all its edges are undirected), the second smallest eigenvalue increases monotonically with D X . The smallest nonzero eigenvalue Λ 2 of the supra-Laplacian can be expressed as follows: where the v that satisfies Eq. (8) is the eigenvector that corresponds to the Λ 2 eigenvalue. To show that Λ 2 is a monotonic increasing function of D X , let us take two values D X and D 0 X , such that D X < D 0 X , then in the sense that the difference of the right-hand side minus the left-hand side matrices is a non-negative definite matrix. Thus, for any v, and, therefore, the minima of the respective choices of v are ordered in the same way, proving that Λ 2 increases when D X increases. Note, however, that different concepts of optimality have been used in undirected multiplex within different contexts, e.g., optimal level of interconnectivity [41] or optimal couplings due to penalizing cost functions [42,43].

IV. DIRECTED MULTIPLEX: EMERGENCE OF OPTIMAL COUPLING
Next, we show that the heterogeneity introduced by the edge directionality of any of the layers of the multiplex networks modifies substantially the behavior of diffusion processes operating on them. In this case, at an intermediate degree of interlayer coupling, the multiplex may exhibit a faster diffusion (faster convergence to the steady state) than when the two layers are fully coupled (fully coupled would have been expected by extrapolation of the undirected case [17] to have the fastest rates). We refer to this region of coupling as the prime regime and show that within it there is an optimal coupling for which the multiplex achieves the fastest diffusion. To this end, we study the characteristic timescale of convergence to the steady state of a CTRW acting on different directed multiplex networks consisting of (without loss of generality) two layers, wherein the intralayer connectivity is represented by unweighted graphs (i.e., we use the combinatorial Laplacian). The supra-Laplacian is defined similarly to the undirected scenario as Notice that D 1 , D 2 , and D X can be interpreted for a CTRW as scalars that control the relative speed of a walker in each layer and between layers, respectively.
In what follows, we demonstrate the existence of the prime regime via (i) numerical and (ii) analytical exploration of simple directed multiplex providing evidence of the nonmonotonicity of Λ 2 (in contrast to the monotonicity in the undirected case) and the existence of at least one maximum of Λ 2 (optimal coupling) for intermediate values of interlayer coupling D X . For more general networks, (iii) we show that the emergence of the prime regime and the magnitude of the maximum attainable rate of convergence to the steady state for the diffusion process necessitates that one of the layers has a sufficient networkscale directionality as quantified by the NDI.
A. Numerical evidence of the existence of the prime regime We examine three prototype multiplex networks ( Fig. 1, top panels), which are chosen to exemplify different types of heterogeneities that can appear when directed networks are considered: (a) Multiplex 1 consists of a directed network (layer 1) and its undirected counterpart (layer 2); this example serves to illustrate the effect of directionality in only one of the layers. (b) Multiplex 2 is composed of the same directed network in the two layers; this example serves to show the effect of different diffusion coefficients (D 1 vs D 2 ) at the different layers. (c) Multiplex 3 where layer 1 contains the reverse directed network (opposite directionality of all the links) of the network present in layer 2; this example serves to illustrate the effect of directionality when the topology and rates of transition are the same. Note that for all three multiplexes, the interlayer links are undirected. Figure 1 (bottom panels) presents the smallest nonzero eigenvalue (in terms of its real part) ReðΛ 2 Þ of the L out as a proxy of the timescale of convergence to the steady state τ (recall that τ ∼ 1 ReðΛ 2 Þ ) when the coupling between layers (via D X ) is varied for each of the multiplexes shown. The specific values of D 1 and D 2 used to produce the results are annotated in the figure. The results from further exploration of these parameters are shown in Appendix A, and a comparison between the dynamics on the directed multiplex with respect to their undirected counterparts is shown in Appendix B. Before discussing the specifics of each panel in Fig. 1, the existence of four distinct regimes stands out. (1) Linear, where the layers are effectively decoupled and D X is the limiting factor. In this case, ReðΛ 2 Þ increases linearly as 2D X . (2) Sublinear, wherein D X is larger than the smallest nonzero eigenvalue (in terms of the real part) of the slowest layer [Reðλ slow 2 Þ, where the superscript slow refers to the slowest layer]. In this scenario, the slowest layer becomes the limiting factor (with respect to the rate of convergence of the multiplex to achieve the steady state), and, therefore, an increase in D X translates into a sublinear rate of increase for ReðΛ 2 Þ. (3) Prime, corresponding to the range of D X for which ReðΛ 2 Þ exceeds the value of ReðΛ 2 Þ for D X → ∞. Within the prime regime, we can define the optimal coupling as the value of D X for which ReðΛ 2 Þ achieves the absolute maximum. The optimal coupling occurs for values of D X in the vicinity of the smallest nonzero eigenvalue (in terms of the real part) of the fastest layer [Reðλ fast 2 Þ, where the superscript fast refers to the fastest layer]. In this case, the speed of the transport within the fastest layer competes with the transport across layers, achieving a configuration where both layers contribute significantly to the total transport but conserving a relative degree of independence in their internal dynamics (i.e., not fully synchronized or decoupled). Notice that the mathematical nature of the maximum reached at the optimal coupling can vary from smooth to nonanalytical depending (4) Asymptotic, wherein D X ≫ Reðλ slow 2 Þ, Reðλ fast 2 Þ. The two layers are completely coupled, and the counterpart nodes in the different layers are fully synchronized, behaving as a single node.
Not all four regimes are present for every multiplex configuration. In fact, the linear and the asymptotic regimes (it could be even argued for the sublinear regime as well) were observed and defined for undirected multiplex in Ref. [17]. However, the prime regime is a characteristic that solely applies to multiplex networks that have a directed network in at least one layer (see Appendix E for further evidence and proof in Sec. III, which shows that for an undirected multiplex, Λ 2 is always a monotonically increasing function of D X ). The difference in the dynamics of the directed and undirected multiplexes emerges from the fact that the directionality of the links in directed layers allows faster exploration of the nodes in the same layer (e.g., when a random walker leaves through a link, there might not be a way back through the same route). Note that this asymmetry in the path directionality is, in general, not advantageous with respect to achieving stationary states for diffusionlike processes. However, directed networks, when integrated in a multiplex, can be catalysts of diffusionlike processes when a right balance between their coupling with other layers (high enough to have access, through different layers, to shortcuts that overcome the asymmetry of the paths) and the degree of independent dynamics they preserve (take advantage of their faster exploratory capability within the layer) exists. Notice that a complete coupling of the layers in a multiplex can be interpreted as a monoplex resulting in the superposition of the connectivity of the different layers. This effect would result in the removal of the constraint on path asymmetries and, therefore, potential suboptimal times of convergence to steady state.
The diffusion in directed multiplex networks shows a richer phenomenology as compared to the undirected scenario, especially regarding the prime regime and the optimal coupling. For Multiplex 1, a prime regime is observed as a result of accessing different topological paths when the directed network is coupled to its undirected counterpart. In the case of the Multiplex 2, the networks are identical in the different layers, and, therefore, no new topological paths are created when the layers are combined in the multiplex. However, different values of the coefficients D 1 and D 2 can create distinct dynamic paths (i.e., same topology but different rates of transport) giving rise to different gradients among counterpart nodes in the separate layers, resulting in the emergence of a prime regime. For this multiplex, a transition from the linear to sublinear regime is also apparent when values of D X exceed Reðλ slow 2 Þ ¼ 0.5. The case of Multiplex 3 is particularly interesting. In this case, an optimal coupling is observed for D X ∼ 0.8 resulting from access to different topological paths when the directed network is coupled to its reverse. This configuration serves as a word of caution, since in many instances, undirected networks are assumed if transitions in both directions of a link (i → j and j → i) are possible. However, if the process is such that it operates with a certain degree of independence in each direction, the properties of the dynamics of such a system can substantially differ from those of a fully undirected network. In fact, we show how the multiplex formed by the two reverse networks can function faster than each of the layers and also faster than the equivalent undirected network (λ 2 ¼ 1).
Finally, we show that there exist two regimes of coupling wherein a superdiffusionlike behavior (i.e., the rate of convergence of the multiplex to steady state is faster than in the fastest of its individual layers) can be observed.
(1) Asymptotic is equivalent to the one described for the undirected multiplex [17], where the multiplex (at the asymptotic high-coupling regime) can be approximated by the superposition network, where the availability of paths between nodes is increased, enhancing the efficiency of diffusive processes. (2) Prime, where for an intermediate degree of coupling, the speed of the processes is higher than when the two layers are fully coupled (e.g., the asymptotic scenario). The emergence of the Prime regime can be intuitively interpreted by acknowledging that intermediate values of coupling permit maintaining a degree of independence of the intralayer dynamics and, therefore, enhancing the fastest node exploration due to the directionality of the links while also permitting interlayer shortcuts that overcome the intralayer asymmetry of the paths. Note that directed multiplex networks can exhibit superdiffusion within the prime regime, even in scenarios where the asymptotic regime is not superdiffusive [e.g., Figs. 1(a) and 1(b) bottom panels and Fig. 3 (bottom panel)].
B. Analytical evidence of the existence of the prime regime As we show numerically in the previous section, for the directed multiplex (i.e., at least one of the layers consists of a directed network), the behavior of Λ 2 is not necessarily a monotonic function of D X . The fact that directed multiplex networks include directionality implies that the matrices required for their representation are not symmetric (complex spectrum), which makes it challenging to find analytically the conditions for the existence of a prime regime as a result of the nonmonotonic behavior of Λ 2 . Furthermore, from our analysis (see Appendix A), we show that even keeping the topologies of the different layers constant, different combinations of the diffusion coefficients ½D 1 ; D 2 can result in transitions from monotonic to nonmonotonic behavior of Λ 2 . Here, we show the analytical solution for the evolution of Λ 2 as a function of D X for the case of a two-layer multiplex with layers consisting of exactly the same topology but different diffusion coefficients corresponding to each layer (e.g., multiplex 2).
Let us assume without loss of generality that D 1 ¼ 1 and D 2 ¼ α (and, therefore, L 1 ¼ L and L 2 ¼ αL), then the supra-Laplacian of the multiplex can be written as follows: Let us denote A ¼ 1 2 ð1 þ αÞL and B ¼ 1 2 ðα − 1ÞL. Since both A and B are scaled versions of L, they commute, and they are simultaneously diagonalizable.
To compute an analytical expression for the second smallest eigenvalue, we use the Schur complement and the fact that A and B commute to conclude that Since A, B, and L are all simultaneously diagonalizable, Eq. (13) reduces to simple scalar polynomial expressions. Specifically, the smallest nonzero eigenvalue Λ 2 satisfies where λ 2 is the second smallest eigenvalue of L.
Recall that the rate of convergence of the multiplex to the stationary state is controlled by the real part of Λ 2 . Thus, to prove that ReðΛ 2 Þ is not a monotonic function of D X , we need to show that Reð dΛ 2 dD X Þ becomes zero for at least one value of D X ∈ R þ . Given that u ¼ Reðλ 2 Þ and v ¼ Imðλ 2 Þ, we can write the derivative of Λ 2 with respect to D X as follows: Re Note that for λ 2 ∈ R, as is the case in undirected multiplex, ð dΛ 2 dD X Þ > 0 ∀ D X > 0. Therefore, Λ 2 is an increasing monotonic function of D X as expected (and proven in general in the previous section). For values of λ 2 ∈ C, i.e., v > 0, we analyze Eq. (16) for two different cases: u > v and v > u.
(i) Case u > v: In this case, all the terms of the radicand are strictly positive (given that α ≠ 1), and, therefore, the denominator is always greater than 2D X . Thus, the Reð dΛ 2 dD X Þ > 0 ∀ D X , which implies that ReðΛ 2 Þ is a monotonically increasing function of D X . (ii) Case v > u: We examine the asymptotic behaviors of Reð dΛ 2 dD X Þ as D X → 0 and as D X → ∞, Equation (18) shows that for D X large enough, the expression Reð dΛ 2 dD X Þ is negative. Thus, for cases where v > u, a change in the sign of the derivative occurs, which implies by continuity that Reð dΛ 2 dD X Þ ¼ 0 at least for one value of D X ∈ R þ and, therefore, that ReðΛ 2 Þ is not a monotonic function of D X .

C. Conditions for the existence of the prime regime
The three directed multiplexes illustrated in Fig. 1 are instrumental to unveil the prime regime and how the emergence of this regime relies on the faster exploration of nodes prompted by the directed layers of those multiplexes. However, the rather simple topologies (e.g., six nodes per layer) of those three multiplexes and the acknowledged fact that not all directed multiplexes can exhibit the prime regime call for the identification of the relevant topological attributes of directed multiplexes underpinning the potential emergence of the prime regime.
Here we argue that the existence of a directed layer is not a sufficient condition for the potential faster exploration of the constitutive nodes of the multiplex, since the directionality of the links does not necessarily imply the emergence of directionality at the network scale. To quantify the directionality at the network scale of a monoplex, we propose the network directionality index (NDI) as a metric, where hd ij i ¼ 1 NðN−1Þ P i;j d ij is the average path length between all pairs of nodes (note that given the directionality of the graph, for each pair of nodes i and j, two distances d ij and d ji are defined), and hjΔd ij ji is the average path asymmetry, and it is computed as follows: hjΔd ij ji ¼ Thus, for an undirected network, NDI ¼ 0 since hjΔd ij ji ¼ 0, and for a directed circle (closed chain) network wherein all the links are oriented in the same direction (i.e., clockwise or anticlockwise), as the number of nodes N increases NDI tends to 1 (see Appendix E for proof).
We hypothesize that only a directed multiplex with a directed layer characterized by a sufficiently high value of NDI can exhibit the prime regime. To explore the validity of this hypothesis, we examine a two-layer multiplex wherein one layer consists of a fully connected network (undirected), and the other layer consists of a directed layer, whose connectivity is modified to progressively decrease the NDI [see Fig. 2(a)]. To systematically explore the existence of the prime regime as a function of the NDI, we introduce the prime regime index For P > 0, the prime regime exists, and the magnitude of P quantifies the degree to which the maximum rate of transport exceeds that obtained at full coupling. In Fig. 2(b), we provide evidence supporting our hypothesis that a sufficient network-scale directionality is needed in one of the layers for the multiplex to exhibit the prime regime. Note that P as a function of the NDI does not exhibit only a threshold behavior, i.e., P ¼ 0 for NDI < 0.3, but also P exhibits an increasing trend for values of NDI > 0.3 making more prominent the effects of optimal coupling as the directionality at the network scale increases. Finally, we remark that the actual threshold value of the NDI, beyond which the prime regime emerges, depends on the topology (e.g., directionality, complementariness, etc.) of the layers in the multiplex. In this sense, it would be interesting to explore these dependencies in the future, for instance, also investigating whether it is possible to have a prime regime when the NDI of each layer is below its respective threshold but that characterizing the whole multiplex is enough to induce this regime.

D. Real-world network examples
In this section, we illustrate with two types of real-world network topologies the potential application of the results reported in this paper: (i) Delta channel networks: Delta channel networks are directional graphs which have been recently studied using graph theory to understand how topology and transport dynamics might relate to the geomorphic properties of the system, such as sediment cohesiveness [25,26], how upstream changes might affect downstream propagation of fluxes [24], and how deltas self-organize to maximize the diversity in the flux delivery to the shoreline [28]. Deltas are extremely low-relief landforms with overall slopes of the order of 10 −5 , where transport processes can have an important diffusive component. Examples of such transport processes are (i) nutrient transport (including nitrogen where the timescale of transport is of interest since it can limit the denitrification rates), (ii) sediment transport (critical for land formation and delta self-maintenance), and (iii) contaminant dispersion (with apparent environmental impacts). Thus, here we utilize the multiplex framework to analyze how advective (directional) versus diffusive (nondirectional) transport on the channel network affects the overall transport timescales. We utilize the channel network of the third largest delta on Earth, the Mekong Delta [see Fig. 3 (top) and Appendix G for network description and more information on the Mekong Delta]. We represent the advective transport by the directed (downstream) channel network (layer 1) and the diffusive transport by its undirected counterpart (layer 2). The interlayer coupling represents the rate of change from one transport regime to the other. Depending on the diffusion coefficients of the two layers (here chosen as D 1 ¼ 10.0 and D 2 ¼ 1.0) and the interlayer degree of coupling, the coupled advective-diffusive transport on the multiplex can exhibit transport timescales slower or faster than either of the individual layers. As shown in Fig. 3 (bottom panel), we observe the prime regime for intermediate values of coupling and an optimal coupling for D X ∼ 1, i.e., fastest transport timescale. According to our previous results, we explain the existence of the prime regime based on the presence of a sufficient network-scale directionality in the directed layer as measured by the network directionality index, which in this case is computed to be NDI ¼ 0.47. Although an illustrative example, it serves to demonstrate that in highly directional systems, overlying processes that lack directionality can explore pathways of transport that result in enhancing the overall system transport timescales with significant practical consequences. (ii) Author-citation networks: We utilize two authorcitation networks extracted, respectively, from publications in Physical Review Letters (PRL) and Physical Review C (PRC) [see Fig. 4(a) and Appendix H for network description]. We construct a multiplex network wherein one layer consists of a citation network extracted from PRL, and the second layer is the equivalent network obtained by extracting the citation among the same authors using PRC publications. We carry out a systematic exploration of the parameter space ½D PRL ; D PRC ; D X and find no combination of those parameters that give rise to the prime regime. These results are within the expectations obtained from the relative low values of the NDI in both layers (NDI PRL ¼ 0.24 and NDI PRC ¼ 0.18). However, if a network (layer) is engineered to have a higher value of the NDI, our results suggest that it would be possible to induce a prime regime for such a multiplex. To demonstrate this, we create a synthetic citation network to substitute the PRC layer, keeping certain statistics such as the out-degree distribution equal to those observed in PRC. As shown in Fig. 4(b), when we analyze the new multiplex, the prime regime emerges for intermediate values of coupling. The implications of the results suggest that by engineering layers with high directionality at the network scale, we can speed up the process of diffusion of information for intermediate values of coupling-a system state that is more realistic than full coupling.

V. CONCLUSIONS
In conclusion, we reveal the existence of a prime regime in which directed multiplex networks may exhibit a faster system-wide diffusion for intermediate values of coupling than in the asymptotic limit when the different layers are fully coupled. Within that regime, there is an optimal value of coupling at which the diffusion is the fastest; i.e., the timescale of convergence to steady state is minimal. Furthermore, the rates of diffusion in multiplex configurations where the prime regime is exhibited can be such that a new superdiffusion regime emerges, where the multiplex diffusivity as a whole is faster than that of the fastest of its layers for intermediate values of coupling. We argue that it is precisely the directionality of paths at the network scale that sets an anisotropic layout for the process, and combined with a balance of (1) significant connectivity across layers (making accessible paths in other layers) and (2) a degree of independence in the intralayer dynamics can catalyze the overall system transport dynamics. Using synthetic and real-world multiplex networks, we show how a sufficient degree of network-scale directionality in the directed layer appears to be a necessary condition for the emergence of the prime regime. We prove analytically that undirected multiplexes do not exhibit a prime regime for intermediate couplings, and their rate of convergence to steady state is a monotonically increasing function of the interlayer coupling. Thus, the extrapolation of the expectation from previous results on dynamics on undirected multiplex networks, wherein diffusion processes achieve the fastest rate of convergence to steady state when the different layers are fully coupled, generally leads to a wrong assessment of system dynamics. Our results open up new paths of research addressing questions such as whether natural complex systems self-organize to configurations where the optimal coupling is accessible to their dynamics. We also believe that extending this study to other real settings might have important implications. For instance, many social or technosocial networks (like online social networks) have many directed relationships, i.e., followingfollower relations are not always reciprocal or symmetric, which make them good candidates to observe the phenomenology described here. Alternatively, our results could also inspire the design of directed multilayer networks in which directionality is exploited so as to make these systems operate in a regime in which the optimal coupling is reachable for the dynamics.

APPENDIX A: EFFECT OF THE INTRALAYER DIFFUSION COEFFICIENT
We explore the characteristic rate of convergence to the stationary state solution of a continuous-time Markov chain on multiplexes 1, 2, and 3 (see Figs. 5, 6, and 7, respectively) for different pairs ðD 1 ; D 2 Þ as a function of the interlayer coupling D X . We show that depending on the relative values of the intralayer diffusion coefficients ðD 1 ; D 2 Þ for a given topology, we can observe the optimal coupling regime [e.g., Figs. 5(a), 6(c), and 7(d)] or not [e.g., Figs. 5(d) and 6(a)], and even more interesting for some triplets ðD 1 ; D 2 ; D X Þ, wherein D X lies within the optimal coupling regime, the multiplex as a whole achieves superdiffusion [e.g., Figs. 5(b) and 7(a)]; i.e., the rate of convergence to the steady-state solution is faster than the rate corresponding to the fastest of the layers.
(i) Multiplex 1.-From Fig. 5, we highlight the following important observations: (i) We first observe that for pairs of D 1 and D 2 where the rate of convergence is larger for the directed layer [Reðλ 1 2 Þ, superscript 1 refers to layer 1, blue line] than for the undirected network [Reðλ 2 2 Þ, superscript 2 refers to layer 2, red line], an optimal coupling emerges where the maximum Reðλ 2 Þ for the multiplex network (black line) is achieved exhibiting higher values than those observed for the asymptotic limit (green line) [see Figs. 5(a) and 5(b)]. (ii) In Fig. 5(b) ðD 1 ¼ 2.5; D 2 ¼ 1.0Þ, the rates of convergence in the individual layers Reðλ 1 2 Þ and Reðλ 2 2 Þ are very similar, and, therefore, each individual layer has similar transport properties. This multiplex configuration when working at optimal coupling is characterized by much faster rates of convergence to the stationary solution than the rates of the individual layers, exhibiting a superdiffusive behavior (even though no superdiffusive behavior is observed in the asymptotic limit). (iii) Figure 5(c) ðD 1 ¼ D 2 ¼ 1Þ illustrates a scenario in which although the directed network is slightly slower than the undirected network Reðλ 1 2 Þ ⪅ Reðλ 2 2 Þ for intermediate values of coupling, ReðΛ 2 Þ can exceed the asymptotic values for D X → ∞. (iv) Finally, Fig. 5(d) shows a scenario where the undirected layer is substantially faster than the directed layer [Reðλ 2 2 Þ > Reðλ 1 2 Þ]. This multiplex network, although directed, when subject to a diffusion process, presents dynamics dominated by its undirected (faster) layer. Consequently, the trends observed in Fig. 5(d), especially the lack of an optimal coupling regime, are comparable to those reported for undirected multiplex networks (e.g., Ref. [17]). (ii) Multiplex 2. -From Fig. 6, we highlight the following important observations: (i) We first report in Fig. 6(a)  similar emerging behavior: When the diffusivity of the two layers differs, an optimal coupling exists where the rate of convergence is the fastest. There are two main phenomena that when combined give rise to the emergence of this optimum transport regime: (1) The different diffusivity coefficients in the different layers form distinct dynamical paths (meaning, the same topological paths but with different celerity) that enhance the mixing. (2) For intermediate values of coupling, the two layers contribute substantially to the total transport (the random walker spends time in both layers) but keeping a certain degree of independence in their dynamics, in contrast with high-coupling scenarios (where the transport across counterpart nodes in different layers is fastest acting as a constraint in the dynamics of all the layers, since similar gradients of concentration will appear among different nodes within each of the layers). Note that in this case, we do not present results for values of D 1 exceeding D 2 for obvious symmetry in the dynamics of the multiplex since the two layers are identical. There are two main phenomena that when combined give rise to the emergence of this optimum transport regime: (1) The different directionality of the paths in the different layers enhances the mixing. This phenomenon requires a high enough value of coupling to make the alternative paths in different layers accessible (i.e., the probability of transition between layers has to be significant).
(2) For intermediate values of coupling, the two layers contribute substantially to the total transport (the random walker spends time in both layers) but keeping a certain degree of independence in their dynamics, in contrast with high-coupling scenarios (where the transport across counterpart nodes in different layers is fastest acting as a constraint in the dynamics of all the layers since similar gradients of concentration will appear among different nodes within each of the layers). (ii) Precisely when the rate of convergence of the individual layers is similar [Reðλ 1 2 Þ ≈ Reðλ 2 2 Þ], the rate of convergence in the multiplex is more enhanced with respect to the rates of the individual layers. In these scenarios, the overall system presents superdiffusionlike behavior for the optimal coupling regime [e.g., Figs. 7(a) and 7(b)]; i.e., the rates of convergence that characterize the multiplex are even faster than both of the individual layers. Note that in this case, we do not present results for values of D 1 exceeding D 2 for obvious symmetry in the dynamics of the multiplex since the two layers are topologically identical (one layer is the reverse of the other).

APPENDIX B: DIRECTED VS UNDIRECTED MULTIPLEXES
We compare the characteristic time of convergence of a diffusionlike process acting on the different multiplexes with their undirected counterparts (Figs. 8, 9, and 10 for multiplexes 1, 2, and 3, respectively) as a function of the interlayer coupling D X . The optimal coupling regime, i.e., a range of intermediate values of D X for which the rate of convergence to the stationary state is faster than when the two layers are fully coupled, is observed only for the directed version of the multiplex.
(i) Multiplex 1.-Multiplex 1 (Fig. 8 left panels) consists of a directed network (layer 1, top) and its undirected counterpart (layer 2, bottom). This example serves to illustrate the effect of directionality in only one of the layers. We compare the characteristic time of convergence of a diffusionlike process acting on this multiplex with its undirected counterpart (undirected multiplex 1). Figure 8 (right panels) shows the value of the smallest (in terms of its real part) nonzero eigenvalue (Λ 2 ) of the supra-Laplacian as a function of interlayer coupling D X for both the directed (top) and undirected (bottom) multiplex and for the same pair of diffusion coefficients ðD 1 ; D 2 Þ ¼ ð2.5; 1.0Þ. We highlight the existence of an optimal coupling that is exhibited only by the directed multiplex (D X ≈ 1.2), wherein the rate of convergence to the steady state is the fastest (even faster than for high values of coupling [asymptotic limit, green line]). (ii) Multiplex 2.-Multiplex 2 (Fig. 9, left panels) is composed of the same directed network (top) in the two layers. This example is to illustrate the effect of different rates (D 1 vs D 2 ) at the different layers. We compare the characteristic time of convergence of a diffusionlike process acting on this multiplex with its undirected counterpart (undirected multiplex 2, bottom). Figure 9 (right panels) shows the value of the smallest (in terms of its real part) nonzero eigenvalue ðΛ 2 Þ of the supra-Laplacian as a function of interlayer coupling D X for both the directed (top) and undirected (bottom) multiplex and for the same pair of diffusion coefficients ðD 1 ; D 2 Þ ¼ ð1.0; 5.0Þ. We highlight the existence of an optimal coupling that is exhibited only by the directed multiplex ðD X ≈ 2.4Þ wherein the rate of convergence to the steady state is the fastest (even faster than for high values of coupling [asymptotic limit, green line]). and undirected (bottom) multiplex and for the same pair of diffusion coefficients ðD 1 ; D 2 Þ ¼ ð1.0; 1.0Þ. We highlight the existence of an optimal coupling that is exhibited only by the directed multiplex ðD X ≈ 0.8Þ wherein the rate of convergence to the steady state is the fastest (even faster than for high values of coupling [asymptotic limit, green line]). Furthermore, in this example, the rate of convergence observed in the direct multiplex at the optimal coupling regime is even faster than the fastest rate computed for its undirected counterpart.

APPENDIX C: MATHEMATICAL NATURE OF THE OPTIMAL COUPLING MAXIMUM
We present evidence of the diverse mathematical nature of the maximum observed for the optimal coupling. This can range from a smooth analytical maxima in the case of multiplex 2 to nonanalytical for the cases of multiplexes 1 and 3 (see Fig. 11).

APPENDIX D: EDGE ASYMMETRY
We further illustrate the mechanism that gives rise to the prime regime by systematically exploring the existence of this regime as a function of the degree of directionality of the multiplex. The directionality of both layers is modulated by introducing a reverse edge by each existing edge with strength controlled by the edge directionality symmetry index (EDS) (see Fig. 12, left panels). Note that for EDS ¼ 0, both layers are completely directed (maximum asymmetry; i.e., if node i is connected to node j, node j is not connected to node i), and for EDS ¼ 1, both layers are undirected (all the edges are undirected, and, therefore, the supra-adjacency matrix of the multiplex is symmetric). For different values of the EDS, we compute the prime regime index, P ¼ max½ReðΛ 2 Þ−Re½Λ 2 ðD X →∞Þ any value of D X , and P > 0 otherwise. Figure 12 (right panels) shows that the existence of the prime regime requires different degrees of directionality (i.e., values of the EDS) depending on the topology of the networks in each layer, and more specifically, depending on the complementariness of these topologies. Thus, more complementary topologies (e.g., multiplex 3) can achieve the fastest rates of convergence for higher values of the EDS. These results also serve to show that the prime regime is not a characteristic of all directed multiplexes, but its existence requires a minimum degree of directionality, whose value depends on the complementariness of the topologies between layers.

APPENDIX E: NETWORK DIRECTIONALITY INDEX
Consider a directed strongly connected graph with N nodes, and let denote its distance matrix. That is, d ij equals the number of links of the shortest path from node i to node j. It can be easily computed via where A is the adjacency matrix of the graph since, as is well known [6], the ði; jÞ entry of A k , which is a nonnegative integer, is precisely the number of paths of length k that connect node i to node j.
We define the NDI as the ratio of the average absolute difference over the average distance between nodes, i.e., Equivalently, This index quantifies the asymmetry in the path connectivity between nodes; it clearly vanishes precisely when d ij − d ji ¼ 0 for all pairs ði; jÞ. We now prove that for a directed circle (closed chain) network consisting of N nodes and N directed links, all oriented in the same direction (clockwise or anticlockwise), the NDI tends to 1 when the number of nodes N tends to infinity.
Note that for any strongly connected directed network, d ij ≠ 0 for any pair of nodes k ¼ ði; jÞ. Define M k ¼ maxðd ij ; d ji Þ and m k ¼ minðd ij ; d ji Þ and note that M k þ m k ¼ N in a circle network. Let us also denote the total number of pairs of nodes by P ¼ NðN−1Þ 2 . From Eq. (E2), we deduce that We first show that for a directed circle network the average distance between nodes [denominator in Eq. (E3)] is equal to N=2. Indeed, Next, we examine the average absolute difference hjΔd ij ji [numerator in Eq. (E3)] separately in two cases depending on whether the number of nodes N is an odd or even number: Case N odd: For a circle network, and therefore, Case N even: In this case, for a circle network, and therefore, In view of the above arguments, we conclude that when N → ∞, NDI → 1. It is interesting to note that a number of basic questions about graphs are surprisingly challenging. For instance, conditions on the entries of D ¼ ½d ij , in addition to the obvious triangular inequality d ij þ d jk ≥ d ik that needs to hold for all triplets ði; j; kÞ, to ensure that D is a distance matrix of some directed graph is, to the best of the authors' knowledge, not known. It appears that the conjecture that NDI < 1 is equally challenging. The fact that a circular directed network, which displays extreme asymmetry, has a NDI tending to 1 as N → ∞, gives credence to this conjecture.

APPENDIX F: RANDOM NETWORKS WITH VARIABLE DIRECTIONALITY AT THE NETWORK SCALE
Here we detail a model to generate strongly connected directed networks with variable directionality at the network scale used to develop the examples presented in Fig. 2. Three model parameters N, k, and r need to be specified: N being the number of nodes in the graph, k the out-degree assumed constant for all nodes (but it can be generalized to follow a distribution), and r the reachability parameter. Networks are generated according to the following rules: (i) All nodes are first arranged in a directed circle pattern corresponding to entries of the adjacency matrix being aði; i þ 1Þ ¼ 1 for i ¼ 1; 2; …; N − 1, aðN; 1Þ ¼ 1, and aðj; kÞ ¼ 0 for all other pairs of indices. (ii) Next, each node i is linked to another k − 1 nodes chosen at random by sampling without replacement from the indexing set of nodes fi þ 1; …; i þ rg. If i þ r > N and the chosen index i þ k to establish connection exceeds N as well, i.e., i þ k > N, it is replaced by N − ði þ kÞ to satisfy the circular boundary condition. The value of r determines the maximum distance between nodes (in terms of links in the original circle network) at which links can be established according to rule (ii); large values of r have the effect of equalizing and then reversing directionality at the network scale. The model is used to generate the examples in Fig. 2(a) in the main text with parameters N ¼ 100, k ¼ 6, and increasing values of r: r ¼ 1 (NDI ¼ 0.99, circle network), r ¼ 10 (NDI ¼ 0.88), r ¼ 15 (NDI ¼ 0.77), r ¼ 20 (NDI ¼ 0.67), r ¼ 25 (NDI ¼ 0.54), r ¼ 30 (NDI ¼ 0.44), r ¼ 50 (NDI ¼ 0.29), and r ¼ 100 (NDI ¼ 0.28).
In the description of the model and in the examples provided in Fig. 2(a), k is considered uniformly distributed for simplicity. The synthetic citation network in Fig. 4(b) (main text) is also constructed using the same model, but k is chosen to follow the empirical out-degree distribution of the citation network of Physical Review C described in Appendix H. Note that in order to generalize the out-degree distribution, a secondary rule needs to be added in the case where certain node degrees k i exceed the value of r. In the synthetic citation network, to faithfully follow the chosen out-degree distribution, nodes with k i > r are allowed to establish links that reach nodes beyond the distance r. Specifically, such nodes link to all the subsequent nodes in the range ½i þ 1; i þ k − 1.

APPENDIX G: MEKONG DELTA
The Mekong Delta is located in Southern Vietnam (10.1°, 150.6°) and receives input from the Mekong River with an average water discharge of 14 770 m 3 s −1 and sediment discharge of 1.60 × 10 7 tons yr −1 . The main forcing of the delta is river and wave [45]. The Mekong Delta considered as the third largest delta with an area of 93 781 km 2 and population of 17 million [46] is one of the world's main food baskets [47].
Delta channel network: We represent the delta channel network as a graph, where the edges represent channels, and nodes correspond to the locations where one channel splits into new channels (bifurcation) or two or more channels merge into a single channel (confluence) (see Fig. 3, top panel). All the channels that connect the apex (inlet) of the delta to the shoreline, with widths equal or larger than 50 m, are considered in this network. The network consists of 253 nodes and 374 links.
Mekong Delta multiplex: We define a two-layer multiplex based on the channel network of the Mekong Delta. Layer 1 consists of the directed channel network, where the direction of the edges corresponds to the main direction of water flow, and layer 2 consists of the undirected counterpart of the network in layer 1. Note that the outlet nodes (nodes located at the shoreline that drain water directly to the ocean) of the delta have been reconnected to the apex (inlet) by adding additional edges to guarantee that the directed layer contains a strongly connected network. Similar rewiring has been used in Ref. [24], arguing conservation of mass to compute the steady-state distribution of fluxes in delta channel networks.

APPENDIX H: AUTHOR-CITATION NETWORKS
We produce two directed networks of citations among authors using publications between the years 1980 and 1989 in PRL and PRC, respectively. More specifically, each of the two networks is constructed as follows: (i) The nodes correspond to the authors who published articles during the specified period in the considered journal. (ii) A link connects author A to author B if a publication of author A in the considered journal was cited by author B (flow of information from author A to author B) in an article of the same journal and during the specified period.
Author-citation multiplex: We aim to construct a directed multiplex network, where one layer corresponds to the author-citation network in PRL, and the other layer corresponds to the author-citation network obtained from PRC publications. Given the constraints required to build the directed multiplex (i.e., a multiplex requires the same set of nodes in both layers, and in order to have nontrivial steady-state distributions of the diffusion process in both layers, we subsample the initial networks to obtain layers consisting of networks, which have a single strongly connected component) the number of nodes (authors) is N ¼ 1368 in each layer, with 33 449 links in the network extracted from PRL and 60 726 links in the layer corresponding to the PRC network.