Fractional dynamics on circulant multiplex networks: optimal coupling and long-range navigation for continuous-time random walks

TL;DR

This paper investigates fractional continuous-time random walks on multiplex networks, revealing an optimal inter-layer coupling that minimizes relaxation time and demonstrating linear mean square displacement growth despite enhanced diffusion conditions.

Contribution

It provides analytical expressions and numerical analysis of fractional random walks on circulant multiplex networks, highlighting the impact of inter- and intra-layer coefficients on dynamics.

Findings

Relaxation time has a minimum at an optimal inter-layer to intra-layer coefficient ratio.

Mean square displacement increases linearly with time even under enhanced diffusion.

Analytical solutions are derived for circulant multiplex networks with finite nodes.

Abstract

This work analyzes fractional continuous-time random walks on two-layer multiplexes. A node-centric dynamics is used, in which it is assumed a Poisson distribution of a walker to become active, while a jump to one of its neighbors depends on the connection weight. Synthetic multiplexes with well known topology are used to illustrate dynamical features obtained by numerical simulations, while exact analytical expressions are presented for multiplexes assembled by circulant layers with finite number of nodes. Special attention is given to the effect of inter- $D_x$ and intra-layer $D_i$ coefficients on the system's behavior. In opposition to usual discrete time dynamics, the relaxation time has a well defined minimum at an optimal value of $D_x/D_i$. It is found that, even for the enhanced diffusion condition, the walkers mean square displacement increases linearly with time.

Equations

$$\frac{\mathrm{d}\vec{x}(t)}{\mathrm{d}t}=-{\mathcal{L}}^{\mathcal{M}}\vec{x}(t)$$

$${\mathcal{L}}^{\mathcal{M}}={\mathcal{L}}^{\ell }+{\mathcal{L}}^x$$

$${\mathcal{L}}^{\ell }=\left(\begin{array}{cccc}{D}_1{\mathbf{L}}^1& & & \\ & D_2{\mathbf{L}}^2& & \\ & & \ddots & \\ & & & D_M{\mathbf{L}}^M\end{array}\right),$$

$${\mathcal{L}}^x=\left(\begin{array}{cccc}{\sum }_{\alpha }{D}_{1\alpha }{\mathbf{I}}_N& -D_{12}{\mathbf{I}}_N& \cdots & -D_{1M}{\mathbf{I}}_N\\ -D_{21}{\mathbf{I}}_N& {\sum }_{\alpha }{D}_{2\alpha }{\mathbf{I}}_N& \cdots & -D_{2M}{\mathbf{I}}_N\\ ⋮& ⋮& \ddots & ⋮\\ -D_{M1}{\mathbf{I}}_N& -D_{M2}{\mathbf{I}}_N& \cdots & {\sum }_{\alpha }{D}_{M\alpha }{\mathbf{I}}_N\end{array}\right).$$

$${\sigma }_i^{\alpha }={\left({\mathcal{L}}^{\mathcal{M}}\right)}_{ff},$$

$$\frac{\mathrm{d}\vec{p}(t{)}^T}{\mathrm{d}t}=-\vec{p}(t{)}^T{\mathcal{S}}^{-1}{\mathcal{L}}^{\mathcal{M}}=-\vec{p}(t{)}^T\text{mathcalboondox}L,$$

$$\frac{\mathrm{d}\vec{p}(t{)}^T}{\mathrm{d}t}=-\vec{p}(t{)}^T{\mathcal{L}}^{\mathcal{M}}.$$

$$\text{mathcalboondox}{L}_{fg}=\frac{{\mathcal{L}}_{fg}^{\mathcal{M}}}{{\mathcal{L}}_{ff}^{\mathcal{M}}}={\delta }_{fg}-\text{mathcalboondox}{T}_{fg},$$

$$\mathrm{MSD}\equiv ⟨r^2(t)⟩=⟨{\left(x(t)-x_0\right)}^2⟩\sim t^{\epsilon },$$

$$p(t{)}_{i^{\alpha }\to j^{\beta }}=\vec{p}(t{)}_g={\vec{\mathcal{C}}}_g{}^T\exp \left(-t\text{mathcalboondox}L\right){\vec{\mathcal{C}}}_f,$$

$$⟨r^2(t)⟩=\frac{1}{N^2{M}^2}\sum _{\alpha =1}^M\sum _{\beta =1}^M\sum _{i=1}^N\sum _{j=1}^N{\left(d_{i^{\alpha }\to j^{\beta }}\right)}^{2}p(t{)}_{i^{\alpha }\to j^{\beta }},$$

$$\frac{\mathrm{d}\vec{x}(t)}{\mathrm{d}t}=-{\left({\mathcal{L}}^{\mathcal{M}}\right)}^{\gamma }\vec{x}(t),$$

$${\left({\mathcal{L}}^{\mathcal{M}}\right)}^{\gamma }={\mathbf{Q\Delta }}^{\gamma }{\mathbf{Q}}^{†}=\sum _{m=1}^{NM}{\mu }_m^{\gamma }{\vec{\psi }}_j{\vec{\psi }}_j^{†},$$

$${\left({\sigma }_i^{\alpha }\right)}^{(\gamma )}={\left({\mathcal{L}}^{\mathcal{M}}\right)}_{ff}^{\gamma }=\sum _{m=1}^{NM}{\mu }_m^{\gamma }{\mathbf{Q}}_{fm}{\mathbf{Q}}_{fm}^{\ast },$$

$${\left({\mathcal{L}}^{\mathcal{M}}\right)}_{fg}^{\gamma }=\sum _{m=1}^{2N}{\mu }_m^{\gamma }{\mathbf{Q}}_{fm}{\mathbf{Q}}_{gm}^{\ast }.$$

$$\text{mathcalboondox}{L}_{fg}^{\left(\gamma \right)}=\frac{{\left({\mathcal{L}}^{\mathcal{M}}\right)}_{fg}^{\gamma }}{{\left({\mathcal{L}}^{\mathcal{M}}\right)}_{ff}^{\gamma }}={\delta }_{fg}-\text{mathcalboondox}{T}_{fg}^{\left(\gamma \right)},$$

$$p(t{)}_{i^{\alpha }\to j^{\beta }}^{\left(\gamma \right)}=\vec{p}(t{)}_g^{\left(\gamma \right)}={\vec{\mathcal{C}}}_g{}^T\exp \left(-t\text{mathcalboondox}{L}^{(\gamma )}\right){\vec{\mathcal{C}}}_f,$$

$${⟨r^2(t)⟩}^{\left(\gamma \right)}=\frac{1}{N^2{M}^2}\sum _{\alpha =1}^M\sum _{\beta =1}^M\sum _{i=1}^N\sum _{j=1}^N{\left(d_{i^{\alpha }\to j^{\beta }}\right)}^{2}p(t{)}_{i^{\alpha }\to j^{\beta }}^{\left(\gamma \right)}.$$

$${\mathcal{L}}^{\mathcal{M}}=\left(\begin{array}{cc}{\mathbf{L}}^1+D_x{\mathbf{I}}_N& -D_x{\mathbf{I}}_N\\ -D_x{\mathbf{I}}_N& {\mathbf{L}}^2+D_x{\mathbf{I}}_N\end{array}\right)=\left(\begin{array}{cc}{\mathbf{C}}^1& -D_x{\mathbf{I}}_N\\ -D_x{\mathbf{I}}_N& {\mathbf{C}}^2\end{array}\right),$$

$${\mathbf{F}}^{-1}\left(\begin{array}{cc}{\mathbf{C}}^1& -D_x{\mathbf{I}}_N\\ -D_x{\mathbf{I}}_N& {\mathbf{C}}^2\end{array}\right)\mathbf{F}=\left(\begin{array}{cc}{\mathrm{\Xi }}^1& -D_x{\mathbf{I}}_N\\ -D_x{\mathbf{I}}_N& {\mathrm{\Xi }}^2\end{array}\right)$$

$$\mathbf{F}=\left(\begin{array}{cc}\mathbf{U}& 0\\ 0& \mathbf{U}\end{array}\right)$$

$${\xi }_m^{\alpha }$$

$${\mu }_{2m-1}=\frac{{\xi }_m^1+{\xi }_m^2+\sqrt{{\left({\xi }_m^1-{\xi }_m^2\right)}^2+4D_x^2}}{2},$$

$${\mu }_{2m}=\frac{{\xi }_m^1+{\xi }_m^2-\sqrt{{\left({\xi }_m^1-{\xi }_m^2\right)}^2+4D_x^2}}{2},$$

$$\begin{array}{c}{\mathbf{Q}}_{fg}=\frac{1}{\sqrt{N\left(1+M_g^2\right)}}\exp \left(-\mathfrak{i}\frac{2\pi }{N}(f-1)\lfloor \frac{g-1}{2}\rfloor \right)\\ {\mathbf{Q}}_{(f+N)g}=\frac{M_g}{\sqrt{N\left(1+M_g^2\right)}}\exp \left(-\mathfrak{i}\frac{2\pi }{N}(f-1)\left(\lfloor \frac{g-1}{2}\rfloor +N\right)\right)\end{array}\}\mathrm{for}f,g\in \left\{1,\cdots ,N\right\},$$

$$\begin{array}{c}{\mathbf{Q}}_{(f-N)g}=\frac{1}{\sqrt{N\left(1+M_g^2\right)}}\exp \left(-\mathfrak{i}\frac{2\pi }{N}(f-1-N)\left(\lfloor \frac{g-1}{2}\rfloor +N\right)\right)\\ {\mathbf{Q}}_{fg}=\frac{M_g}{\sqrt{N\left(1+M_g^2\right)}}\exp \left(-\mathfrak{i}\frac{2\pi }{N}(f-1)\left(\lfloor \frac{g-1}{2}\rfloor +N\right)\right)\end{array}\}\mathrm{for}f,g\in \left\{N+1,\cdots ,2N\right\}$$

$$M_g=\frac{{\xi }_{\left(1+\lfloor (g-1)/2\rfloor \right)}^1-{\mu }_g}{D_x},$$

$${\sigma }_1^{(\gamma )}$$

$${\sigma }_2^{(\gamma )}$$

$${\lambda }_{2m-1}=\frac{\kappa +\sqrt{{\kappa }^2-\nu }}{2},$$

Full text

Keywords: Diffusion, Network dynamics, Diffusion in random media

I Introduction

Over the last two decades, modeling complex systems as networks has proven to be a sucessful approach to characterize the structure and dynamics of many real-world systems albert02 ; newman03 ; bocaletti06 ; barrat08 . Different dynamics have been investigated on top of networks, such as spreading processes arruda18 , percolation Dorogovtsev08 , or synchronization tang14 ; rodrigues16 . Among these dynamical processes, diffusion and random walks (RWs) have been analyzed thoroughly. Indeed, RWs models have widespread use both in the analysis of diffusion and navigability in networks as in exploring their fine-grained organization noh04 ; klafter2011first ; masuda17 . Most of the research on RWs relies on the nearest-neighbor (NN) paradigm noh04 in which the walker can only hop to one of the NNs of its current node (or position). However, other RWs definitions in discrete and continuous times, allow for both NN and long distance hops. Well known examples are, e.g., Lévy RWs riascos12 ; guo16 ; estrada17b ; nigris17 , RWs based on the $d$-path Laplacian operator estrada17b , and those defined by fractional transport in networks riascos14 ; riascos15 ; michelitsch19 . Such non-NN strategies often correspond to the better options for randomly reaching a target in an unknown environment as, for instance, the foraging of species in a given environment lomholt2008levy ; humphries2010environmental ; song2010modelling ; rhee2011levy ; ghandi2011foraging .

More recently, multiplex networks and other multilayer structures were identified as more comprehensive frameworks to describe those complex systems in which agents may interact through several channels. Here, links representing channels with different meaning and relevance are embedded in distinct layers boccaletti14 ; kivela14 ; battiston17 ; bianconi18 . As before, multiplexes have widespread use to describe, among others, social szell10 ; cozzo13 ; li15 ; arruda17 , biochemical cozzo12 ; battiston17 , or transportation systems dedomenico14 ; aleta17 . Here, layers may represent underground, bus and railway networks in large cities, each one associated to a different spatial and temporal scale. A diffusion process in a transport system like this occurs both within and across layers arruda18 ; gomez13 ; cencetti19 ; dedomenico16 ; tejedor18 , accounting for the actual displacement to reach different places as within the change of embarking platforms. Multiplexes have also been used to study related dynamical systems, such as reaction-diffusion asllani14 ; kouvaris15 ; busiello18 and synchronization processes gambuzza15 ; sevilla15 ; genio16 ; allen17 . The eigenvalue spectrum of the supra-Laplacian matrix associated to the multiplex plays a key role in the description such diffusive process. Several results indicate the presence of super-diffusive behavior in undirected and directed multiplexes gomez13 ; sole13 ; radicchi13 ; cozzo12 ; cozzo16 ; sanchez14 ; cozzob16 ; arruda18b ; serrano17 ; tejedor18 ; cencetti19 , meaning that their relaxation time to the steady state is smaller than any other observed for the isolated layer.

On the other hand, it seems that fractional diffusion on multiplex networks has not yet received due attention. In this work, we address this matter, following the framework introduced in Refs. riascos14 ; riascos15 ; michelitsch19 . Our study considers the fractional diffusion of node-centric continuous-time RWs on undirected multiplex with two layers. We present numerical results of main dynamical features for several multiplexes with well known topology, and derive exact analytical expressions for circulant layers. An important aspect is the role played by the ratio between the inter- and intra-layer coefficients. Our results indicate a nonmonotonic behavior in the rate of convergence to the steady state as the inter-layer coefficient increases. The walkers mean square displacement (MSD$\equiv\left\langle r^{2}(t)\right\rangle$) illustrates the existence of an optimal diffusive regime depending on both the inter-layer coupling and on the fractional parameter. Following the nomenclature in dedomenico14 , we show that fractional dynamics turns classic random walkers into a new type of physical random walkers, which are allowed to (i) switch layer and (ii) perform long hops to another distant vertex in the same jump. In spite of such enhanced diffusion, our results also show that the MSD still increases linearly with time when the number of multiplex nodes is finite.

The paper is organized as follows. In Sec. II, we describe the diffusion dynamics on undirected multiplex networks and MSD. Section III defines fractional dynamics on such systems. Analytical results for fractional random walks with continuous time on regular multiplex with circulant layers are presented in Sec. IV. Finally, our conclusions are summarized in Sec. V.

II Diffusion dynamics on multiplex networks

Let us consider a multiplex $\mathcal{M}$ with $N$ nodes and $M$ layers. Let $\mathbf{A}^{\alpha}=\left(\mathbf{A}_{ij}^{\alpha}\right)$ denote the adjacency matrix for the $\alpha$th layer with $1\leq\alpha\leq M$. In this work we focus on multiplexes $\mathcal{M}$ whose layers are undirected and unsigned (i.e., the edge weights are nonnegative), and contain no self-loops, i.e., $\mathbf{A}_{ij}^{\alpha}=\mathbf{A}_{ji}^{\alpha}=w_{ij}>0$ if there is a link between the nodes $i$ and $j$ in the layer $\alpha$ (and $i\neq j$), and 0 otherwise. If the layers of $\mathcal{M}$ are also unweighted, then $w_{ij}=1$. On the other hand, the strength $s_{i}^{\alpha}$ of a vertex $i$ with respect to its connections with other vertices $j$ (with $j=1,\cdots,N$) in the same layer $\alpha$ is given by $s_{i}^{\alpha}=\sum_{j=1}^{N}\mathbf{A}_{ij}^{\alpha}$.

On a discrete space, diffusive phenomena are described in terms of Laplacian matrices, which can be formally obtained as a discretized version of the Laplacian operator $\left(-\bigtriangledown^{2}\right)$ on regular lattices, and have been generalized for more complex topologies riascos14 ; riascos15 . In the case of multiplexes, we let $\vec{x}$ be a $NM\times 1$ state (column) vector whose entry $i+(\alpha-1)N$ (with $i=1,\cdots,N$) describes the concentration of a generic flowing quantity at time $t$ on node $i$ at the $\alpha$th layer, $x_{i}^{\alpha}$. Therefore, the usual diffusion equation in matrix form reads:

[TABLE]

where

[TABLE]

denotes the $NM\times NM$ (combinatorial) supra-Laplacian matrix defined in masuda17 ; gomez13 ; cencetti19 ; sole13 ; bianconi18 , and $\mathcal{L}^{\ell}$ and $\mathcal{L}^{x}$ represent the intra-layer and the inter-layer supra-Laplacian matrices, respectively, given by:

[TABLE]

and

[TABLE]

In the above equations, $\mathcal{L}^{\ell}$ is a block-diagonal matrix, $D_{\alpha}$ denotes the intra-layer diffusion constant in the $\alpha$th layer, $D_{\alpha\beta}$ (with $\alpha,\beta\in\left\{1,\cdots,M\right\}$ and $\beta\neq\alpha$) refers to the inter-layer diffusion constant between the $\alpha$th and $\beta$th layers, $\mathbf{I}_{N}$ represents the $N\times N$ identity matrix, and $\mathbf{L}^{\alpha}$ is the usual $N\times N$ (combinatorial) Laplacian matrix of the layer $\alpha$, with elements $\mathbf{L}^{\alpha}_{ij}=s_{i}^{\alpha}\delta_{ij}-\mathbf{A}_{ij}^{\alpha}$, and $\delta$ is the Kronecker delta function. Thus, the matrix $\mathbf{\mathcal{L}}^{\mathcal{M}}$ represents the generalization of the graph Laplacian to the case of linear diffusion on multiplex networks. For simplicity we will consider only diffusion processes where $D_{\alpha\beta}=D_{\beta\alpha}$, so that $\mathcal{L}^{x}$ and $\mathbf{\mathcal{L}}^{\mathcal{M}}$ are symmetric matrices.

Finally, according to Eq. (2), the elements of the main diagonal of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ represent the total strength of a given node at a given layer, i.e., the sum of (i) the strength of such vertex with respect to its connections with other vertices in the same layer and (ii) the strength of the same vertex with respect to connections to its counterparts in different layers. To denote the total strength of node $i$ in layer $\alpha$, we introduce the following short-hand notation:

[TABLE]

where $f=i+(\alpha-1)N$.

II.1 CTRWs and MSD on multiplex networks

The usual discrete-time random walk (DTRW) is a random sequence of vertices generated as follows: given a starting vertex $i$, denoted as “origin of the walk”, at each discrete time step $t$, the walker jumps to one NN of its current node masuda17 ; aldous02 ; lovasz93 ; noh04 . In the case of multiplex networks, because of its peculiar interconnected structure, DTRW can also move from one layer to another one, provided that such layers ($\alpha$ and $\beta$) are connected with each other (i.e., $D_{\alpha\beta}\neq 0$).

In the case of continuous-time random walks (CTRWs), it is assumed that the duration of the walkers waiting times between two moves obeys a given probability density function masuda17 . For that reason, the actual timing of the moves must be taken into account. For the sake of simplicity, in this work we consider that the waiting times are distributed according to a Poisson distribution with constant rate (i.e.,the exponential distribution).

Here, it becomes necessary to distinguish between two different cases of Poissonian CTRWs: Node-centric and edge-centric RWs. The Poissonian node-centric CTRWs follow the same assumption of DTRWs: when a walker becomes active, it moves from its current node to one of the neighbors with a probability proportional to the weight of the connection between such nodes. On the other hand, in the Poissonian edge-centric CTRWs, each edge (rather than a node) is activated independently according to a renewal process. Thus, if a trajectory includes many nodes with large strengths, the number $n$ of moves in the time interval $[0,t]$ tends to be larger than for trajectories that traverse many nodes with small strengths. For a wider description of the specific features of each random walk the reader is referred to Ref. masuda17 .

To generalize the fractional diffusion framework introduced in Refs. riascos14 ; riascos15 ; michelitsch19 to multiplex networks, we restrict our analyzes to CTRWs. Let $\vec{p}$ be a $NM\times 1$ vector whose entry $i+(\alpha-1)N$ (with $i=1,\cdots,N$) is the probability of finding the random walker at time $t$ on node $i$ at the $\alpha$th layer. The transition rules governing the diffusion dynamics of the node-centric random walks are determined a master equation which, in terms of suitably defined matrices, can be written as

[TABLE]

On the other hand, the dynamics of the edge-centric ones are described by

[TABLE]

In Eqs. (6) and (7), $\mathbf{X}^{T}$ stands for the transpose of matrix $\mathbf{X}$, and $\mathbf{\mathcal{S}}$ is the $NM\times NM$ diagonal matrix with elements $\mathbf{\mathcal{S}}_{ii}=\mathbf{\mathcal{L}}^{\mathcal{M}}_{ii}$. The $NM\times NM$ matrix $\mathcalboondox{L}$ denotes the “random walk normalized supra-Laplacian” masuda17 (or just “normalized supra-Laplacian” dedomenico14 ). According to the definition of $\mathcalboondox{L}$, its elements can be expressed as,

[TABLE]

where $\mathcalboondox{T}_{fg}$ are the elements of the $NM\times NM$ transition matrix $\mathcalboondox{T}$ of a discrete-time random walk, describing transition probability from one node to its NNs in the corresponding layer or to the node’s counterparts in different layers, with equal probability noh04 ; klafter2011first ; riascos14 ; riascos15 ; michelitsch19 ; masuda17 . Indeed, note that $\mathcalboondox{T}$ is a stochastic matrix, that satisfies $\mathcalboondox{T}_{ff}=0$ and $\sum_{g=1}^{NM}\mathcalboondox{T}_{fg}=1$. Following dedomenico14 ; riascos14 ; riascos15 ; michelitsch19 , heareafter we will consider only the case of node-centric CTRW (or classical random walker as in dedomenico14 ).

The MSD, defined by $\left\langle r^{2}(t)\right\rangle$, is a measure of the ensemble average distance between the position of a walker at a time $t$, $x(t)$, and a reference position, $x_{0}$. Assuming that $\left\langle r^{2}(t)\right\rangle$ has a power law dependence with respect to time, we have

[TABLE]

where the value of the parameter $\varepsilon$ classifies the type of diffusion into normal diffusion ($\varepsilon=1$), sub-diffusion ($\varepsilon<1$), or super-diffusion ($\varepsilon>1$). Although MSD is one of the used measures to analyze general stochastic data almaas03 ; gallos04 , in order to better characterize diffusion, additional measures are also required, e.g., first passage observables masuda17 . For the type of results we discuss here, $\left\langle r^{2}(t)\right\rangle$ is essential to provide a clear cut way to characterize the time dependence.

According to Eq. (6), the probability of finding the random walker at node $j$ in the $\beta$th layer (at time $t$), when the random walker was initially located at node $i$ in the $\alpha$th layer, is given by:

[TABLE]

where $g=j+(\beta-1)N$ and $f=i+(\alpha-1)N$ (with $i,j\in\left\{1,\cdots,N\right\}$ and $\alpha,\beta\in\left\{1,\cdots,M\right\}$), and $\vec{\mathcal{C}}_{\ell}$ represents the ($NM\times 1$) $\ell-$th vector of the canonical base of $\mathbb{R}^{NM}$ with components $\left(\vec{\mathcal{C}}_{g}\right)_{m}=\delta_{m\ell}$). Therefore, in the case of node-centric CTRWs, we can quantify $\left\langle r^{2}(t)\right\rangle$ at time $t$ as follows:

[TABLE]

where $d_{i^{\alpha}\rightarrow j^{\beta}}$ is the length of the shortest path distance between $i$ in the $\alpha$th layer and $j$ in the $\beta$th layer, that is, the smallest number of edges connecting those nodes.

III Fractional diffusion on multiplex networks

III.1 General Case

In this sub-section we present the general expressions for the combinatorial and normalized supra-Laplacian matrices required to study fractional diffusion in any multiplex network. Thus, following Refs. riascos14 ; riascos15 ; michelitsch19 , we generalize Eq. (1) as

[TABLE]

where $\gamma$ is a real number ($0<\gamma<1$) and $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$, the combinatorial supra-Laplacian matrix raised to a power $\gamma$, denotes here the fractional (combinatorial) supra-Laplacian matrix.

Let us briefly discuss some mathematical properties of the model defined by Eq. (12), as well as qualitative aspects of the expected behavior, limiting cases, and relations to other scenarios characterized by anomalous diffusion. An immediate consequence is that we recover Eq. (1) in the limit $\gamma\rightarrow 1$. This way of defining the fractional supra-Laplacian matrix preserves the essential features of Laplacian matrices, namely: $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$ is (i) positive semidefinite, (ii) stochastic, and (iii) all its non-diagonal elements are non-positive. On the other hand, by setting $D_{\alpha}=1$ and $D_{\alpha\beta}=1$ (with $\alpha,\beta\in\left\{1,\cdots,M\right\}$ and $\beta\neq\alpha$), $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$ is equivalent to the fractional Laplacian matrices of monolayer networks described in riascos14 ; riascos15 ; michelitsch19 . For such cases, it has been shown analytically that the continuum limits of the fractional Laplacian matrix (with $0<\gamma<1$) are connected with the operators of fractional calculus. Indeed, in the case of cycle graphs and its continuum limits, the distributional representations for fractional Laplacian matrices take the forms of Riesz fractional derivatives (see Chapter 6 of michelitsch19 for further details). Besides that, when the above definition of fractional Laplacian matrix is considered, the asymptotic behavior of node-centric CTRWs on homogeneous networks and their continuum limits (with homogeneous and isotropic node distributions) shows explicitly the convergence to a Lévy propagator associated with the emergence of Lévy flights with self-similar inverse power-law distributed long-range steps and anomalous diffusion (see Chapter 8 of michelitsch19 for further details). Alternatively, by using (non-fractional) Laplacian matrices (i.e., $\gamma=1$), Brownian motion (Rayleigh flights) and Gaussian diffusion appear. Both types of asymptotic behaviors are in good agreement with the findings presented in Ref. Metzler2000 for the CTRW model with Poisson distribution of waiting times in homogeneus, isotropic systems, when a Lévy distribution of jump lengths and a Gaussian one are considered, respectively.

We can obtain a set of eigenvalues $\mu_{j}$ and eigenvectors $\vec{\psi}_{j}$ of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ that satisfy the eigenvalue equation $\mathbf{\mathcal{L}}^{\mathcal{M}}\vec{\psi}_{j}=\mu_{j}\vec{\psi}_{j}$ for $j\in\left\{1,\cdots,NM\right\}$ and the orthonormalization condition $\vec{\psi}_{i}^{T}\vec{\psi}_{j}=\delta_{ij}$. Since $\mathbf{\mathcal{L}}^{\mathcal{M}}$ is a symmetric matrix, the eigenvalues $\mu_{j}$ are real and nonnegative. In the case of connected multiplex networks, the smallest eigenvalue is 0 and all others are positive.

Following Refs. riascos14 ; riascos15 ; michelitsch19 , we define the orthonormal matrix $\mathbf{Q}$ with elements $\mathbf{Q}_{ij}=\vec{i}\,^{T}\vec{\psi}_{j}$ and the diagonal matrix $\mathbf{\Delta}=\mathrm{diag}\left(\mu_{1},\mu_{2},\cdots,\mu_{NM}\right)$. These matrices satisfy $\mathbf{\mathcal{L}}^{\mathcal{M}}\mathbf{Q}=\mathbf{Q\Delta}$, therefore $\mathbf{\mathcal{L}}^{\mathcal{M}}=\mathbf{Q\Delta}\mathbf{Q}^{\dagger}$, where the matrix $\mathbf{Q}^{\dagger}$ is the conjugate transpose (or Hermitian transpose) of $\mathbf{Q}$. Therefore, we have:

[TABLE]

where $\mathbf{\Delta}^{\gamma}=\mathrm{diag}\left(\mu_{1}^{\gamma},\mu_{2}^{\gamma},\cdots,\mu_{NM}^{\gamma}\right)$. According to Eq. (13), $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}\vec{\psi}_{j}=\mu_{j}^{\gamma}\vec{\psi}_{j}$ for $j\in\left\{1,\cdots,NM\right\}$. Consequently, the eigenvalues of $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$ are equal to those of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ to the power of gamma, $\mu_{j}^{\gamma}$, and the eigenvectors $\vec{\psi}_{j}$ remain the same for both the supra-Laplacian and the fractional supra-Laplacian matrices.

On the other hand, the diagonal elements of the fractional supra-Laplacian matrix defined in Eq. (13) introduce a generalization of the strength $\sigma_{i}^{\alpha}=\mathbf{\mathcal{L}}^{\mathcal{M}}_{ff}$ with $f=i+(\alpha-1)N$ to the fractional case. In this way, the fractional strength of node $i$ at layer $\alpha$ is given by:

[TABLE]

where $X^{*}$ denotes the conjugate of $X$. In general, the elements of the fractional (combinatorial) supra-Laplacian matrix can be calculated as follows:

[TABLE]

Now, by analogy with the random walk normalized supra-Laplacian matrix $\mathcalboondox{L}=\mathbf{\mathcal{S}}^{-1}\mathbf{\mathcal{L}}^{\mathcal{M}}$, we introduce the normalized fractional supra-Laplacian matrix $\mathcalboondox{L}^{\left(\gamma\right)}$ with elements

[TABLE]

where $\mathcalboondox{T}_{fg}^{\left(\gamma\right)}$ denotes the elements of the $NM\times NM$ fractional transition matrix $\mathcalboondox{T}^{\left(\gamma\right)}$. Note that $\mathcalboondox{T}^{\left(\gamma\right)}$ is a stochastic matrix, that satisfies $\mathcalboondox{T}_{ff}^{\left(\gamma\right)}=0$ and $\sum_{g=1}^{NM}\mathcalboondox{T}_{fg}^{\left(\gamma\right)}=1$.

Finally, when fractional diffusion takes place for a given $\gamma$, the probability $p(t)_{i^{\alpha}\rightarrow j^{\beta}}^{\left(\gamma\right)}$ of finding a node-centric CTRW at node $j$ in the $\beta$th layer (at time $t$), when the random walker was initially located at node $i$ in the $\alpha$th layer, is expressed by:

[TABLE]

where $g=j+(\beta-1)N$ and $f=i+(\alpha-1)N$ with $i,j\in\left\{1,\cdots,N\right\}$ and $\alpha,\beta\in\left\{1,\cdots,M\right\}$. Thus, the MSD for fractional dynamics, denoted as $\left\langle r^{2}(t)\right\rangle^{\left(\gamma\right)}$, is given by:

[TABLE]

According to Eq. (18), the time evolution of $\left\langle r^{2}(t)\right\rangle^{\left(\gamma\right)}$ and the corresponding diffusive behavior for the multiplexes $\mathcal{M}$ considered in this work depend on $\gamma$, the number of nodes $N$, the total amount of layers $M$, the topology of each layer, given by $\mathbf{A}^{\alpha}$, and the intra- and inter-layer diffusion constants $D_{\alpha}$ and $D_{\alpha\beta}$ (with $\alpha,\beta\in\left\{1,\cdots,M\right\}$ and $\beta\neq\alpha$), respectively. As expected, when $\gamma\rightarrow 1$, Eqs. (14), (16), (17), and (18) reproduce the corresponding equations in the previous section.

III.2 Circulant multiplexes

In this subsection we analyze fractional diffusion on a $N$-node multiplex network in which all its $M$ layers consist of interacting cycle graphs i.e., each layer $\alpha$ (with $\alpha\in\left\{1,\cdots,M\right\}$) contains a $N-$ring topology in which each node is connected to its $J^{\alpha}$ left and $J^{\alpha}$ right nearest nodes. Thus, $J^{\alpha}$ represent the interaction parameter of the layer $\alpha$. It is easy to see that, if $N$ is odd, then $J^{\alpha}\in\left\{1,\cdots,(N-1)/2\right\}$. Note that, when $J^{\alpha}=1$, the $\alpha$-th layer contains a cycle graph whereas, if $J^{\alpha}=(N-1)/2$, it corresponds to a complete graph. For the purpose of deriving exact expressions for the eigenvalues, hearafter we only consider multiplex networks with $M=2$ layers and odd number of nodes. Besides, to emphasize the inter-layer diffusion process and simplify the notation, we choose the diffusion coefficients $D_{1}=D_{2}=1$ and $D_{12}/D_{\alpha}=D_{x}$ for $\alpha\in\left\{1,2\right\}$ gomez13 ; cencetti19 ; tejedor18 .

According to Eqs. (2)-(4), the combinatorial supra-Laplacian matrix of the multiplex is written as

[TABLE]

where both $\mathbf{C}^{1}$ and $\mathbf{C}^{2}$ are $N\times N$ circulant matrices. Since exact analytical expressions for the eigenvalues and eigenvectors of circulant matrices are well known Mieghem11 , it is also possible to obtain similar expressions for $\mu_{j}$ and $\vec{\psi}_{j}$ (for $j\in\left\{1,\cdots,2N\right\}$). So let us write

[TABLE]

where

[TABLE]

is a $2N\times 2N$ block-diagonal matrix, $\mathbf{U}$ is the $N\times N$ hermitian matrix with elements $\mathbf{U}_{ij}=\omega^{(i-1)(j-1)}/\sqrt{N}$, $\omega\equiv\exp(-\mathfrak{i}2\pi/N)$, $\mathfrak{i}=\sqrt{-1}$, $\mathrm{\Xi}^{\alpha}=\mathrm{diag}\left(\xi_{1}^{\alpha},\cdots,\xi_{N}^{\alpha}\right)$, and $\xi_{m}^{\alpha}$ are the eigenvalues of $\mathbf{C}^{\alpha}$, given by

[TABLE]

for $1<m\leq N$, and $\xi_{m}^{\alpha}=D_{x}$ for $m=1$. Since the matrices $-D_{x}\mathbf{I}_{N}$ and $\left(\mathrm{\Xi}^{2}-\mu_{m}\mathbf{I}_{N}\right)$ commute, the eigenvalues of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ can be obtained as:

[TABLE]

and

[TABLE]

for $m\in\left\{1,\cdots,N\right\}$. Note that the eigenvalues $\mu_{m}$ are not ordered from smallest to largest and vice versa (for instance, when $m=1$, $\mu_{2}=0$). Given such set of eigenvalues, the corresponding hermitian matrix of eigenvectors $\mathbf{Q}=\left(\begin{matrix}\vec{\psi}_{1}&\cdots&\vec{\psi}_{2N}\end{matrix}\right)$ has the following elements:

[TABLE]

and

[TABLE]

where

[TABLE]

and $\left\lfloor.\right\rfloor$ denotes the floor function [see Appendix A for further details on the derivation of Eqs. (23)-(27)].

Using the eigenvalue spectrum of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ [Eqs. (23) and (24)] and its eigenvectors [Eqs. (25) and (26)], the fractional strength of any node at layer 1 is

[TABLE]

whereas the fractional strength of the nodes at layer 2 is given by

[TABLE]

for $i\in\left\{1,\cdots,N\right\}$. Note that Eqs. (28) and (29) do not depend on $i$, as expected for circulant layers of interacting cycles.

The set of Eqs. (15)-(18) and (22)-(29) allows to derive the the fractional (combinatorial) supra-Laplacian matrix ($\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$), the normalized fractional supra-Laplacian matrix ($\mathcalboondox{L}^{\left(\gamma\right)}$), the fractional transition matrix ($\mathcalboondox{T}^{\left(\gamma\right)}$), the probability of finding a node-centric CTRW at a given position ($p(t)_{i^{\alpha}\rightarrow j^{\beta}}^{\left(\gamma\right)}$), and $\left\langle r^{2}(t)\right\rangle^{\left(\gamma\right)}$.

Finally, note that, according to Eq. (13), the fractional supra-Laplacian matrix $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$ is a block-matrix whose blocks are circulant and, consequently, so is the normalized fractional supra-Laplacian matrix $\mathcalboondox{L}^{\left(\gamma\right)}$. Therefore, using a strategy similar to that previously described for $\mathbf{\mathcal{L}}^{\mathcal{M}}$, it is possible to derive the eigenvalue spectrum of $\mathcalboondox{L}^{\left(\gamma\right)}$. After conducting the necessary manipulation, the resulting eigenvalues are given by:

[TABLE]

and

[TABLE]

where

[TABLE]

[TABLE]

[TABLE]

and $C^{2}=1+S^{2}$ for $m\in\left\{1,\cdots,N\right\}$ [see Appendix B for further details on the derivation of Eqs. (30)-(34)]. Note that the eigenvalues $\mu_{m}$ and $\mu_{m}^{\gamma}$ of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ and $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$ are not ordered and, consequently, neither are the $\lambda_{m}$. It is also noteworthy that $\lambda_{m}$ (with $m\in\left\{1,\cdots,2N\right\}$) depend on the value of the inter-layer diffusion constant, $D_{x}$.

By using Eqs. (30) and (31), it is possible to obtain the algebraic connectivity of $\mathcalboondox{L}^{\left(\gamma\right)}$, i.e., its second-smallest eigenvalue, denoted here as $\Lambda_{2}$. When $D_{x}\rightarrow 0$, $\Lambda_{2}=\lambda_{1}$, that is:

[TABLE]

In the case of $\Lambda_{2}\neq\lambda_{1}$, it is easy to see that $\Lambda_{2}=\lambda_{2m}<\lambda_{2m-1}$ for $m\in\left\{1,\cdots,N\right\}$ [see Eqs. (30) and (31)]. Thus, for $J^{1}\neq J^{2}$ and $D_{x}\rightarrow\infty$, the algebraic connectivity can be approximated as:

[TABLE]

According to Eqs. (35) and (36), the algebraic connectivity of the normalized supra-Laplacian matrix $\mathcalboondox{L}^{\left(\gamma\right)}$ is nonmonotonic. When $D_{x}\rightarrow 0$, the inter-layer diffusion disappears and the dynamics reduces to those for diffusion on single (isolated) layers. On the other hand, when $D_{x}\rightarrow\infty$, the strength of the vertices are approximately equal to interlayer connection, i.e., $\sigma_{i}^{(\gamma)}\approx D_{x}^{\gamma}$. Thus, the centric-node random walkers spend most of the time in switching layer, instead of jumping to other vertices. Consequently, when the interlayer diffusion coefficient is very small or very large, the diffusion is hindered. For that reason, there is an optimal range of values for $D_{x}$. Note that the previous nonmonotonic trends of node-centric CTRWs emerge from the multiplex structure itself and, consequently, they persist even in the case of $\gamma=1$.

IV Results

In this section, we present our results for the diffusion processes of Poissonian node-centric CTRWs on multiplex networks. Our discussion is mainly focused on their rate of convergence to the steady state in such systems, given by the diffusion time scale $\tau\sim 1/\Lambda_{2}$ masuda17 . To do so, we analyze the nonmonotonic dependence of $\Lambda_{2}$ on the inter-layer coupling (i.e., $D_{x}$), as well as the influence of fractional dynamics (i.e., $\gamma$). In subsection IV.1, we present our analytical results for $\Lambda_{2}$ in regular multiplex networks. By analyzing the case of regular multiplexes with $J=1$, in subsection IV.2 we show that the enhanced diffusion induced by fractional dynamics is due to the emergence of a new type long-range navigation between layers. In subsection IV.3, we provide an example of optimal convergence to the steady state of node-centric CTRWs, discussing the dependence of $\left\langle r^{2}(t)\right\rangle^{(\gamma)}$ on $D_{x}$ and $\gamma$. Finally, in subsection IV.4 we present our results for $\Lambda_{2}$ in non-regular multiplexes.

IV.1 Regular multiplexes

Regular multiplex networks meet the condition $J\equiv J^{1}=J^{2}$. Therefore, in these systems $\xi_{m}\equiv\xi_{m}^{1}=\xi_{m}^{2}=D_{x}+A_{m}$ with $A_{m}\equiv A_{m}^{1}=A_{m}^{2}$ and $m\in\left\{1,\cdots,N\right\}$. Under these circumstances, the eigenvalues of $\mathbf{\mathcal{L}}^{\mathcal{M}}$ reduce to $\mu_{2m-1}=A_{m}+2D_{x}$ and $\mu_{2m}=A_{m}$ (for $m\in\left\{1,\cdots,N\right\}$). Thus, according to Eq. (27), $(M_{m})^{2}=1$ and the fractional strength of the nodes of both layers [see Eqs. (28) and (29)] is given by:

[TABLE]

On the other hand, since the fractional strength is a constant for all nodes, the eigenvalues of the normalized supra-Laplacian matrix are given by $\lambda_{m}=\mu_{m}/\sigma^{(\gamma)}$ for $m\in\left\{1,\cdots,2N\right\}$, while the matrix of the corresponding eigenvectors is $\mathbf{Q}$ [see Eqs. (25) and (26)]. Thus, in regular multiplexes the algebraic connectivity can be calculated as:

[TABLE]

where $c\in\left\{2,N\right\}$ represents the natural number that minimizes $A_{m}$ [see Eq. (22)]. According to Eq. (38), $\Lambda_{2}$ reaches a global maximum at $D_{x}={A_{c}}/{2}$. Therefore, such value of the interlayer diffusion constant guarantees the fastest convergence to the steady state of node-centric CTRWs in these regular multiplexes (for an example see subsection IV.3). On the other hand, it is worth mentioning that $D_{x}={A_{c}}/{2}$ does not depend on $\gamma$. Therefore, for a given $J$ and $N$, the optimal value remains the same. However, the smaller the $\gamma$, the larger the optimal algebraic connectivity. Thus, for a given value of $N$ and $J$, the more intense the fractional diffusion (i.e., the smaller the paramenter $\gamma$), the larger $\Lambda_{2}$ is and, therefore, the faster the convergence to the steady state is ($\tau\sim 1/\Lambda_{2}$). In Fig. 1(a), we show the dependence of $\Lambda_{2}$ on $D_{x}$ for several regular multiplex networks with two layers. As can be observed, the numerical results are in excelent agreement with Eq. (38). Finally, note that, according to Eqs. (22) and (38), for $c=2$ and fixed values of $J$ and $\gamma$, the smaller the system size $N$, the larger $D_{x}=A_{c}/2$ and $\Lambda_{2}$, corroborating the expected results of a faster convergence to the steady state.

IV.2 Emergence of interlayer long-range navigation

In this section we explore the navigation strategy of node-centric CTRWs on regular multiplex networks that are formed by cycle graphs (i.e., $J=1$). We will explore the probability transition between two nodes $i$ and $j$ that are located in different layers. Let us suppose that $i$ is at layer 1 and $j$ is at layer 2. Following Refs. riascos14 ; riascos15 ; michelitsch19 , by using Eqs. (25)-(27), it is possible to approximate the element of the fractional supra-Laplacian that refers to $i$ and $j$ as:

[TABLE]

for $\left|D_{x}\right|>\left|A_{m}\right|$ and $D_{x}>0$, where $d\equiv d_{i^{1}\rightarrow j^{1}}$, $\theta_{m}=2\pi\left(m-1\right)/N$ represents the shortest path distance between $i$ and $j$ at layer 1, and $A_{m}=2+2\cos\left(2\theta_{m}\right)$ [see Eq.(22) for $J^{1}=J^{2}=1$ and Appendix C for further details of these derivations]. Besides that, note that $\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)_{i(j+N)}^{\gamma}=0$ for $D_{x}=0$.

On the other hand, by using Eq. (37), a similar expression can be derived for the fractional strength:

[TABLE]

Following Refs. riascos14 ; riascos15 ; michelitsch19 , Eqs. (39) and (40) can be expressed in terms of an integral in the thermodynamic limit (i.e., $N\rightarrow\infty$) which can be explored analytically (see Ref. [43] for a discussion on that integral and Appendix C). The resulting expressions are given by:

[TABLE]

and

[TABLE]

where

[TABLE]

and $\Gamma(x)$ is the $\Gamma$ function [see Appendix C for further details on the derivation of Eqs. (41)-(43)].

According to Eqs. (16), in the therodynamic limit, the elements of the transition matrix between two nodes $i$ and $j$ that belong to different layers, can be approximated by:

[TABLE]

for $d\gg 1$ (i.e., $K_{d}=0$). By using the asymptotic property $\Gamma\left(x+\gamma\right)\approx\Gamma\left(x\right)x^{\gamma}$, it is possible to express $\mathcalboondox{T}_{i(j+N)}^{(\gamma)}$ as follows:

[TABLE]

where $\varrho=1+2\gamma$. Consequently, in the simple case of regular multiplexes with $J=1$, a power-law relation emerge for the transitions between both layers when $d\gg 1$ and $N\rightarrow\infty$. Once $0<\gamma<1$, in this process the long-range transitions between different layers decay with exponent $1<\varrho<3$, in a similar way as in the case of fractional diffusion in monolayer regular networks riascos14 ; riascos15 ; michelitsch19 . In Fig. 2(a) we show the dependence of $\mathcalboondox{T}_{i(j+N)}^{(\gamma)}$ on $d$ for regular multiplex networks with two layers and $J=1$. As can be observed, the predicted exponent $\varrho$ is in excellent agreement with the results obtained from Eqs. (16). Besides, Fig. 2(a) shows that, for a given $D_{x}$, the larger the value of $\gamma$, the larger $\mathcalboondox{T}_{i(j+N)}^{(\gamma)}$ is when $d\sim 1$ (see inset). On the other hand, in Fig. 2(b) we illustrate the dependence of $\mathcalboondox{T}_{i(j+N)}^{(\gamma)}$ on $D_{x}$. As expected, the larger $D_{x}$, the smaller (larger) the element of the transition matrix when $d\gg 1$ ($d\sim 1$). In the case of $d\gg 1$, an increase on the interlayer diffusion coefficient also increases $\sigma^{(\gamma)}$, and the latter is inversely proportional to $\mathcalboondox{T}_{i(j+N)}^{(\gamma)}$ [see Eq. (16)]. The results in Fig. 2(b) also confirm that $\varrho$ does not depend on $D_{x}$ (when $d\gg 1$).

Finally, it is worth mentioning that fractional diffusion induces a novel mechanism of interlayer diffusion: fractional node-centric CTRWs are allowed to switch layer and jump to another vertex that may be very far away. For instance, Lévy RWs in guo16 are not allowed to switch layer and hop during the same jump. On the other hand, the physical RWs presented in dedomenico14 reduce to classic RWs on monolayer networks, which are subject to the NN paradigm. However, these fractional node-centric CTRWs exhibit long-range hops on top of monolayer networks (see riascos14 ; riascos15 ; michelitsch19 ).

IV.3 Gaussian enhanced diffusion: MSD on regular multiplex networks

In this section, we show an example of the nonmonotonic increase in the rate of convergence to the steady state of centric-node CTRWs diffusion. To do so, we study the dependence of the mean square displacement of the walkers, $\left\langle r^{2}(t)\right\rangle^{(\gamma)}$, on $D_{x}$ and $\gamma$ [see Eqs. (17) and (18)]. In Fig. 3 we show the results for a regular multiplex network without fractional diffusion ($\gamma=1.0$, see left panel) and with it ($\gamma=0.5$, see right panel), for several values of $D_{x}$: the optimal interlayer coefficient $D_{x}=A_{c}/2$, and two aditional values, which are very large and very small in comparison to it. As expected, the results show that, for a given value of $\gamma$, the fastest convergence to the steady state corresponds to the optimal $D_{x}$. We can observe that the differences between the results for $D_{x}=A_{c}/2$ and $D_{x}\ll A_{c}/2$ are very small. In both cases the layers are barely coupled due to the very small value of $A_{c}$ [see Eq. (22) when $c=2$]. Nonetheless, when $D_{x}=A_{c}/2$, the inter-layer connection is stronger and $\Lambda_{2}$ reaches a maximum, i.e., the diffusion is enhanced. In the case of very large values of $D_{x}$, the diffusion of the node-centric CTRWs is hindered once, as $\sigma^{(\gamma)}\approx D_{x}$, the walkers spend most of their time switching layers instead of hoping to other nodes inside the layers. On the other hand, for a given value of $D_{x}$, it can be seen that, the smaller the $\gamma$, the larger $\Lambda_{2}$, and the faster the diffusion (the diffusion time scale $\tau\sim 1/\Lambda_{2}$). Thus, the previous findings are in good agreement with Eq. (38) and the data in Fig. 1(a).

Finally, it is worth mentioning that the increased algebraic connectivity induced by fractional dynamics is reflected in the long-range navigation of node-centric walkers. As can be observed in Fig. 3(b), in the case of $\gamma=0.5$, $\left\langle r^{2}(t)\right\rangle^{(\gamma)}\approx t$ (i.e. $\varepsilon=1$). Thus, a Gaussian behavior emerges from the fractional dynamics in finite circulant multiplex networks with two layers. Other examples of circulant multiplex networks with different values of $N$, $J^{1}$, $J^{2}$ and $\gamma$ are presented in the Supplemental Material accompanying this paper, and all of them show perfect agreement with the developed analysis.

IV.4 Optimal diffusion on non-regular multiplexes

In Fig. 4 we present examples of multiplexes with circulant and noncirculant layers, when $\gamma=1$ (left panel) and $\gamma=0.5$ (right panel). As can be seen, in all the cases the nonmonotonic trend of $\Lambda_{2}$ is present. Indeed, when $D_{x}\rightarrow 0$ ($D_{x}\rightarrow\infty$), its dependence on $D_{x}$ is similar to $D_{x}^{\gamma}$ ($D_{x}^{-\gamma}$). For that reason, there is a nonmonotonic increase in the rate of convergence to the steady state of node-centric CTRWs. Besides, in all the topologies tested, the smaller the value of $\gamma$, the larger the global maximum of $\Lambda_{2}$. Therefore, there exist optimal combinations of $D_{x}$ and $\gamma$ that enhance diffusion processes of node-centric RWs and make them faster than those obtained when layers are fully coupled and $\gamma=1$. In the case of circulant multiplexes, the findings presented here are in excellent agreement with Eqs. (35) and (36). On the other hand, these results suggest that the apparent plateau observed in circulant multiplexes is not present in other topological configurations. It seems that the more random the layers are, the larger $\Lambda_{2}$ is. For that reason, finding an optimal value of $D_{x}$ is more crucial in noncirculant multiplex networks than in circulant ones.

V Conclusions

In this work, we have extended the continuous time fractional diffusion framework (for simple networks) introduced in Refs. riascos14 ; riascos15 ; michelitsch19 to multiplex networks with undirected and unsigned layers. Hence fractional diffusion is defined here in terms of the fractional supra-Laplacian matrix of the system, i.e., the combinatorial supra-Laplacian matrix of the multiplex $\mathbf{\mathcal{L}}^{\mathcal{M}}$ to a power $\gamma$, where $0<\gamma<1$. For the purpose of deriving exact analytical expressions, we have considered only diffusion processes in which $\mathbf{\mathcal{L}}^{\mathcal{M}}$ is a symmetric matrix.

We have focused our discussion on the characterization of Poissonian node-centric continuous-time random-walks on circulant multiplexes with two layers, and explored the combined effect of inter and intra layer diffusion with fractional dynamics. We have directed our attention to (i) the effect of the fractional dynamics on the nonmonotonic increase in the rate of convergence to the steady state of such process, (ii) the existence of an optimal regime that depends on both the inter-layer coupling $D_{x}$ and on the fractional parameter $\gamma$, and (iii) the emergence of a new type of long-range navigation on multiplex networks. For circulant multiplexes, analytical expressions were obtained for the main quantities involved in these dynamics, namely: the eigenvalues and eigenvectors of the combinatorial supra-Laplacian matrix and of the normalized supra-Laplacian matrix, the fractional strength of the nodes, the fractional transition matrix, the probability of finding the walkers at time $t$ on any node of a given layer, and the mean-square displacement for fractional dynamics. For other multiplex topologies some of these quantities were obtained by numerical evaluations.

We have shown that, for a given circulant multiplex network, the more intense the fractional diffusion (i.e., the smaller the paramenter $\gamma$), the larger the algebraic connectivity of the normalized supra-Laplacian matrix, denoted as $\Lambda_{2}$. Since the diffusion time scale of the Poissonian node-centric CTRWs on the multiplex $\tau$ is inversely proportional to $\Lambda_{2}$, the smaller the value of $\gamma$, the faster the convergence to the steady state is (i.e., the smaller $\tau$ is). Additionally, in multiplexes with two layers, both $\Lambda_{2}$ and $\tau$ exhibit a nonmonotonic dependence on $D_{x}$, respectively, whether or not there are fractional diffusion. Consequently, the rate of convergence to the steady state must be optimized when using fractional dynamics.

On the other hand, in the simple case of circulant (regular) multiplexes with $J=1$, we have illustrated that, once the fractional diffusion is present (i.e., $0<\gamma<1$), long-range transitions between different layers appear. Indeed, a new continuous-time random walk process appears, since here walkers are allowed to (i) switch layer and (ii) perform long hops to another distant vertex in the same jump. Additionally, in the thermodinamic limit, we have shown that the probability of long range transitions decay according to a power-law with exponent $1<\varrho<3$, in a similar way as in the case of fractional diffusion in monolayer regular networks riascos14 ; riascos15 ; michelitsch19 . We have also shown that the larger $D_{x}$, the smaller (larger) the transition probability between nodes that are very far away (very close).

Finally, the evaluation of $\equiv\left\langle r^{2}(t)\right\rangle$ indicates the existence of the optimal regime that depends on both $D_{x}$ and the fractional parameter $\gamma$. On the other hand, we have shown that the enhaced diffusion induced by fractional dynamics on finite circulant multiplexes exhibits a Gaussian behavior ($\left\langle r^{2}(t)\right\rangle\sim t$) before saturation appears.

The introduction of fractional dynamics on multiplex networks opens new possibilities for analyzing and optimizing (anomalous) diffusion on such arrangements. For instance, given the attention devoted recently to the optimal diffusion dynamics on directed multiplex networks, the generalization of this framework to such systems can be of great interest. On the other hand, the emerging long-range transitions can enhance the efficiency of the navigation on non-circulant multiplex topologies. While studing the dependence of $\Lambda_{2}$ on $D_{x}$ in circulant multiplexes, we have found an apparent plateau in $\Lambda_{2}$ that is not present in other topological configurations. Indeed, it seems that the more random the layers are, the more pronounced is the optimal $\Lambda_{2}$. Thus, finding optimal combinations of $D_{x}$ and $\gamma$ seems more crucial in noncirculant multiplex networks than in circulant ones. For that reason, an exhaustive research on the dependence of $\Lambda_{2}$ on noncirculant topologies with fractional diffusion is being conducted and it will be published elsewhere.

Acknowledgements.

This work was supported by the Brazilian agencies CAPES and CNPq through the Grants 151466/2018-1 (AA-P) and 305060/2015-5 (RFSA). RFSA also acknowledges the support of the National Institute of Science and Technology for Complex Systems (INCT-SC Brazil).

APPENDIX A: Eigenvalues and Eigenvectors of the combinatorial supra-Laplacian matrix for circulant multiplexes with two layers.

Let $\mathcal{M}$ be an undirected multiplex network with $N$ nodes and $M=2$ layers, all of which consist of interacting cycle graphs. According to Eqs. (19) and (20), the eigenvalues of $\mathbf{\mathcal{L}}^{\mathcal{M}}$, denoted as $\mu_{g}$ for $g\in\left\{1,\cdots,2N\right\}$, meet the following condition:

[TABLE]

where $\mathrm{det}\left(\mathbf{X}\right)$ refers to the determinant of a matrix $\mathbf{X}$. Taking into account that (i) the blocks of $\Upsilon$ are square matrices of the same order, and (ii) matrices $-D_{x}\mathbf{I}_{N}$ and $\mathrm{\Xi}^{2}-\mu_{g}\mathbf{I}_{N}$ commute, it is possible to show Silvester00 that Eq. (46) reduces to

[TABLE]

By definition, $\mathrm{\Xi}^{\alpha}=\mathrm{diag}\left(\xi_{1}^{\alpha},\cdots,\xi_{N}^{\alpha}\right)$, where $\xi_{m}^{\alpha}$ are the eigenvalues of $\mathbf{C}^{\alpha}$ [see Eqs. (19) and (22)], and $m\in\left\{1,\cdots,N\right\}$. Consequently, Eq. (47) is equivalent to the following $N$ equations:

[TABLE]

For a given value of $m$, we donote the two roots of Eq. (48) as $\mu_{2m-1}$ and $\mu_{2m}$, respectively. Thus, we obtain Eqs. (23) and (24).

On the other hand, the corresponding eigenvector of $\mu_{g}$, denoted by $\vec{\iota}_{g}$, can be calculated from

[TABLE]

where $\vec{v}_{g}^{1}$ and $\vec{v}_{g}^{2}$ are $N\times 1$ vectors. According to Eq. (49), the elements of $\vec{v}_{g}^{1}$ and of $\vec{v}_{g}^{2}$ should meet simultaneously the following conditions:

[TABLE]

It is possible to see that the previous restrictions are equivalent to Eq. (48). Therefore, in the case of $g\in\left\{2m-1,\,2m\right\}$, Eq. (50) requires that only the elements $\left(\vec{v}_{g}^{1}\right)_{m}$ and $\left(\vec{v}_{g}^{1}\right)_{m}$ are non-zero. Consequently, to normalize $\vec{\iota}_{g}$, we set $\left(\vec{v}_{g}^{1}\right)_{m}=1/T_{g}$ and $\left(\vec{v}_{g}^{2}\right)_{m}=(\xi_{m}^{1}-\mu_{g})/(D_{x}T_{g})$, where

[TABLE]

for $g\in\left\{2m-1,\,2m\right\}$ [see Eq. (27)]. According to the previous results, given the matrix definedby the right hand side of Eq. (20) and its corresponding eigenvectors $\vec{\iota}_{g}$, the matrix $\mathbf{\Omega}=\left(\begin{matrix}\vec{\iota}_{1}&\cdots&\vec{\iota}_{2N}\end{matrix}\right)$ has elements

[TABLE]

for $g\in\left\{1,\cdots,2N\right\}$ and $f=1+\left\lfloor(g-1)/2\right\rfloor$ (i.e., $f\in\left\{1,\cdots,N\right\}$), and zero otherwise.

Finally, the eigenvectors of the combinatorial supra-Laplacian matrix $\mathbf{\mathcal{L}}^{\mathcal{M}}$, i.e. $\vec{\psi}_{g}$, and the matrix $\mathbf{Q}=\left(\begin{matrix}\vec{\psi}_{1}&\cdots&\vec{\psi}_{2N}\end{matrix}\right)$ [in Eqs. (25) and (26)] can be obtained from

[TABLE]

as can be seen from Eq. (21).

APPENDIX B: Eigenvalues and Eigenvectors of the normalizedl supra-Laplacian matrix for circulant multiplexes with two layers.

Let $\mathcal{M}$ be an undirected $N-$node multiplex network in which all its $M=2$ layers consist of interacting cycle graphs. According to Eqs. (14)-(16), the normalized fractional supra-Laplacian matrix of the multiplex $\mathcalboondox{L}^{\left(\gamma\right)}$ is given by $\mathcalboondox{L}^{\left(\gamma\right)}=\mathcal{K}^{-1}\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)^{\gamma}$, where

[TABLE]

Taking into account Eqs. (13), (21), (52) and (53), it is possible to write

[TABLE]

where $\mathbf{\Omega}^{-1}=\mathbf{\Omega}^{T}$, since $\mathbf{Q}=\mathbf{F}\mathbf{\Omega}$, $\mathbf{Q}\mathbf{Q}^{\dagger}=\mathbf{I}_{2N}$, and $\mathbf{I}_{2N}$ represents the $2N\times 2N$ identity matrix. Thus, the eigenspectrum of $\mathcalboondox{L}^{\left(\gamma\right)}$ is equal to that of $\mathcal{K}^{-1}\mathbf{\Omega}\mathbf{\Delta}^{\gamma}\mathbf{\Omega}^{T}$.

Considering that $\mathbf{\Delta}^{\gamma}=\mathrm{diag}\left(\mu_{1}^{\gamma},\mu_{2}^{\gamma},\cdots,\mu_{2N}^{\gamma}\right)$, as well as the definition of $\Omega$ [Eq. (52)], the right hand side of Eq. (55) can be rewritten as

[TABLE]

where $\mathbf{D}^{1}$, $\mathbf{D}^{2}$ and $\mathbf{D}^{3}$ are $N\times N$ diagonal matrices, whose respective $m$ elements are given by

[TABLE]

[TABLE]

and

[TABLE]

for $m\in\left\{1,\cdots,N\right\}$, where

[TABLE]

and $C^{2}=1+S^{2}$. Notice that, by making use of Eqs. (27) and (51), it is possible to express some terms in Eqs. (57)-(59) in terms of $M_{g}$ and $T_{g}$. Being more specific, we can write $\left(M_{2m-1}\right)^{2}=(S+C)^{2}$ and $\left(M_{2m}\right)^{2}=(S-C)^{2}$, as well as $\left(T_{2m-1}\right)^{2}=1+(S+C)^{2}$ and $\left(T_{2m}\right)^{2}=1+(S-C)^{2}$ [see also the definition of $\Omega$ given by Eq. (52)].

According to Eq. (56), to calculate the eigenvalues of $\mathcalboondox{L}^{\left(\gamma\right)}$, denoted as $\lambda_{g}$ for $g\in\left\{1,\cdots,2N\right\}$, we solve the following equation:

[TABLE]

Since (i) the blocks of $\Phi$ are square matrices of the same order, and (ii) matrices $\frac{1}{\sigma_{2}^{(\gamma)}}\mathbf{D}^{3}$ and $\frac{1}{\sigma_{2}^{(\gamma)}}\mathbf{D}^{2}-\lambda_{g}\mathbf{I}_{N}$ commute, it is possible to show Silvester00 that Eq. (61) reduces to

[TABLE]

Finally, the eigenspectrum of $\mathcalboondox{L}^{\left(\gamma\right)}$ is obtained by calculating the roots of the following $N$ equations:

[TABLE]

which are equivalent to Eq. (62). For a given value of $m\in\left\{1,\cdots,N\right\}$, we denote the two roots of Eq. (63) as $\lambda_{2m-1}$ and $\lambda_{2m}$, respectively. Thus, we obtain obtain Eqs.(30)-(33).

APPENDIX C: Transitions between nodes that are located in different layers.

Let $\mathcal{M}$ be an undirected multiplex network with $M=2$ layers, both of which are cycle graphs (i.e., $J=1$). Let us consider two nodes in $\mathcal{M}$, which are located in different layers: $i$ is at layer 1 and $j$ is at layer 2. Following Refs. riascos14 ; riascos15 ; michelitsch19 , by using Eqs. (25)-(27) and conducting the necessary manipulation in Eq. (15), the element of the fractional supra-Laplacian that refers to $i$ and $j$ can be approximated as follows:

[TABLE]

for $\left|D_{x}\right|>\left|A_{m}\right|$ and $D_{x}>0$, where $d\equiv d_{i^{1}\rightarrow j^{1}}$, $\theta_{m}=2\pi\left(m-1\right)/N$, and $A_{m}=2+2\cos\left(2\theta_{m}\right)$ [see Eq.(22) for $J^{\alpha}=1$ and $\alpha\in\left\{1,2\right\}$]. If $D_{x}=0$, then $A_{m}^{\gamma}-\left(A_{m}+2D_{x}\right)^{\gamma}=0$, and this element of the fractional supra-Laplacian is equal to zero.

In the limit $N\rightarrow\infty$, the sums in Eq. (64) can be replaced by integrals, so that

[TABLE]

Taking into account that (i) the last right-hand term in Eq. (65) is equal to $-\left(2Dx\right)^{\gamma}/2$ for $d=0$ and 0 otherwise, and (ii) the analytical results

[TABLE]

and

[TABLE]

we obtain Eq.(41) [see riascos14 ; riascos15 ; michelitsch19 ; Zoia07 for further details about the derivation of Eqs. (66) and (67)].

Finally, note that Eq.(42) can be derived from Eq. (37) by using the same process presented in this Appendix and considering $d=0$ ($\left(\mathbf{\mathcal{L}}^{\mathcal{M}}\right)_{ii}^{\gamma}=\sigma^{(\gamma)}$).