Discrete-Time Quantum Walk on Multilayer Networks

A Multilayer network is a potent platform that paves the way for the study of the interactions among entities in various networks with multiple types of relationships. This study explores the dynamics of discrete-time quantum walks on a multilayer network. We derive a recurrence formula for the coefficients of the wave function of a quantum walker on an undirected graph with a finite number of nodes. By extending this formula to include extra layers, we develop a simulation model to describe the time evolution of the quantum walker on a multilayer network. The time-averaged probability and the return probability of the quantum walker are studied with Fourier, and Grover walks on multilayer networks. Furthermore, we analyze the impact of decoherence on quantum transport, shedding light on how environmental interactions may impact the behavior of quantum walkers on multilayer network structures.


I. INTRODUCTION
Quantum walks (QWs) are the quantum analogues of classical random walks (CRWs).Importantly, QWs contribute to theoretical and applied studies of quantum computing [1] in general and quantum algorithms [2] in particular.There are two broad classes of QWs known as discrete-time QWs (DTQW) and continuous-time QWs (CTQW), each of which has significant distinctions in their mathematical formalism [3].Considerable body of work can be found in literature which explore the dynamics of linear and cyclic QWs in two and higher dimensional spaces, as well as on specific graphs [3,4].
In the study of complex systems, multilayer networks play a crucial role as a modeling tool [5,6].A multilayer network consists of nodes and edges, yet the nodes exist in separate layers representing different forms of interactions.Multilayer networks are used to describe the behavior patterns and evolution of ecological systems [7], complex interactions across multiple layers of biological systems [8], public transportation systems [9,10] and the structure of financial markets [11].In literature, one can find several studies which utilize the framework of CTQWs to study the transport properties of multilayer dendrimer networks [12], honeycomb networks [13] and scale-free networks [14] in detail.Nonetheless, it is clear that there has been insufficient exploration of the dynamics of DTQWs on multilayer networks.Therefore, to address this gap, we present a comprehensive study of DTQWs on a multilayer network.The most general form of networks and multilayer networks are mathematically represented using the notion of a graph [5].Hence, its worth exploring the research work on DTQWs that incorporate graphs.Several formulations of the DTQWs on specific graph structures can be found in literature.However, defining a DTQWs on an arbitrary graph is more difficult than that of a CTQWs [15].Ref. [16] has proposed a framework for defining QWs on regular graphs.Sometimes this framework is termed as Shunt-Decomposition model [17].In [2], it is suggested to add one or more self-loops to each vertex with the purpose of obtaining a regular graphs when modeling the QWs on irregular graphs.Extending the idea given in [18], Kendon [19] has presented a method to simulate DTQWs on general undirected graphs.Sometimes this method is identified as Arc-Reversal model [17].Ref. [20] has proposed another framework for QWs on a graph in which the motion of the quantum walker takes place on the edges of the graph rather than the vertices.
In our work, we give a block matrix representation for the state of the quantum walker instead of the usual column matrix representation.Then, we derive recurrence formulae that imitate the coin toss and the shifting operation of the QW.We adopt the shift operation given in Arc-Reversal model to develop our framework.Using the recurrence formulae and the block matrix representation, we develop a simulation to mimic the progression of QWs on an undirected graph.Later we extend our framework to mimic QWs on multilayer networks.The paper is organized in the following way.In section II, we present our mathematical model for QWs on a graph.Section III is dedicated to extend our mathematical model to include QWs on multilayer networks.Moreover, for the sake of comparison purposes, we model a CRW on a multilayer network as well in section IV.Numerical implementation of QWs on a toy model and some synthetic multilayer networks are given in section V and VI, respectively, along with a detailed analysis on time-average probability, return probability and decoherence.

II. QWS ON A GRAPH
Consider a finite undirected graph G = {V, E} where V = {v 1 , . . ., v n } is the set of vertices (nodes) and E = {(v i , v j ) | v i , v j ∈ V } is the set of edges (connections).Note that, the graphs which are studied here are finite as opposed to the unrestricted line or the integer line we used to define the QWs on a line.In our study, we consider simple graphs (i.e.graphs without self-loops or parallel edges).We denote the adjacency matrix corresponding to G as A ∈ R n×n which is defined in the following way The is d x .For the purpose of gaining an advantage in simulation, we adopt the following strategy to modify the labeling of the coin states.We define a set B x that comprises the labels of the vertices adjacent to the vertex |x⟩ p as B x = {y : vertex |x⟩ p and |y⟩ p have a connection}.Note that, |B x | = d x .Now, define a function as, f x (r) = r th element of B x .Since parallel edges are excluded in our study, the function f x (r) is a bijection.Without the loss of generality, we always arrange the elements of B x in the ascending order.Now, we can denote the basis of the coin Hilbert space Note that, by using the function f x (r), we have labeled the coin states of each vertex using the edges connected to it.Such an approach can be found in the study [18].Now, the state vector of the quantum walker at position |x⟩ p at time t can be written as where α x,fx(r) (t) ∈ C are called the probability amplitudes, For each position x, the coefficients of α x,fx(r) (t) are defined as follows Our next task is to write an expression for the total wave function of the quantum walker on G at time t.For that, we need to sum the state vectors |ψ(x, t)⟩ in (2) over all the vertices.However, we are unable to perform such a summation because the state vectors |ψ(x, t)⟩ corresponding to each vertex resides in different composite Hilbert spaces.Note that, the size of the coin Hilbert space changes with the degree of the node x.Therefore, to perform such a summation, one needs to combine the set of composite Hilbert spaces {H c } x in a reasonable manner, with the purpose of forming a bigger Hilbert space that includes all the state vectors.The operation of direct sum of vector spaces paves a way to combine the composite Hilbert spaces to cater our demand.Let us define where ⊕ denotes the external direct sum of Hilbert spaces and dim(H) = n x=1 d x .According to the definition of the global Hilbert space H, it is obvious that for each vertex x, the state vector |ψ(x, t)⟩ ∈ H. Hence, the total wave function of the quantum walker at time t can be calculated by summing |ψ(x, t)⟩ in (2) over all the vertices.Then, the total wave function at time t can be written as where |ψ t ⟩ ∈ H and x) which acts on the coin states associated to the vertex |x⟩ p holds the transition probabilities from |x⟩ p to its neighbouring vertices.Hence, C (x) can be defined as where {|f x (r)⟩ c } dx r=1 are the basis elements of H ij ∈ C are chosen in such a way that the condition of unitarity of the QW is preserved, i.e. the total probability is unity at all time steps.Hence, (C (x) ) † C (x) = C (x) (C (x) ) † = I.By combining each local coin operator C (x) , one can write the global coin operator C acting on |ψ t ⟩ ∈ H as follows The shift operator of the QW on a graph is defined as follows Note that, both C and S are unitary operators associated to the Hilbert space H. Hence, a single step progression of the quantum walker on the graph is given by where U = SC is the evolution operator and S and C are the shift and coin operators respectively.

A. Matrix representation and simulation
The local coin operator C (x) , associated to the vertex |x⟩ p holds the transition probabilities from |x⟩ p to its neighbouring nodes.Suppose the vertex |x⟩ p is linked to d x number of vertices denoted by |y 1 ⟩ p , . . ., |y dx ⟩ p .Moreover, suppose that y 1 < . . .< y dx .Hence, for each r ∈ {1, . . ., d x }, we can write y r = f x (r).Then, the block matrix of C (x) can be written as As an example, the local coin operators associated to a four vertex graph is shown in Appendix B. In standard mathematical formalism of QWs [21], the state of the quantum walker at time t is represented by a column vector.
However, alternatively, one can give a convenient block matrix representation for the total wave function of the quantum walker on a graph given in (4) as follows In matrix N t , the rows represent the position states and columns represent the coin states.Single row holds the coefficients corresponding to the coin states associated to a single position.According to the definition of the coefficients of α x,r (t), given in (3), some entries of the matrix N t become zero (see Appendix B).Such a block matrix representation of the total wave function can be found in the study [22].Note that, one can view the block matrix representation given in (10) as an adjacency matrix of a weighted graph having time dependent complex weights.By knowing all the coefficients of α x,r (t), in other words, all the elements of N t , we can uniquely determine the total wave function of the quantum walker on the graph at time t.In addition, by updating the elements of N t appropriately, one can determine the matrix N t+1 and hence the total wave function of the quantum walker at time t + 1.The elements of N t are updated in two processes.First, an intermediate matrix N t is generated using the following formula Afterwards, the matrix N t+1 is determined by taking the transpose of N t .The update rule is given by Note that, the formulas given in ( 11) and ( 12) correspond to the coin flip and shift operation of the QW on the graph.Proof is given in the Appendix A.

B. Probability Calculation
The probability P q (x, t) of finding the quantum walker at vertex x at time t can be calculated using (2) as follows Note that, P q (x, t) can be determined by summing the elements in the x th row of the matrix N t ⊙ N * t where ⊙ is the Hadamard product of matrices and N * t is the complex conjugate of the block matrix given in (10).

III. QWS ON MULTILAYER NETWORKS
A multilayer network is a pair } is the set of edges in the L α layer of the multilayer network [6].The positive integers l , n α and r α are termed as the number of layers in M, number of vertices and edges of the layer L α respectively.Moreover, C = {E αβ ⊆ V α × V β : α, β ∈ {1, 2, . . ., l}, α ̸ = β} is the set of edges between L α and L β layers.The elements of C are called crossed layers.Further, the elements of each E α are called the set of intralayer edges and the elements of each E αβ (α ̸ = β) are called the interlayer edges of M [23].Let us denote the adjacency matrices corresponding to each layer L α as A (α) = (a α ij ) ∈ R nα×nα which is defined by for 1 ≤ i, j ≤ n α and 1 ≤ α ≤ l where n α is the number of nodes in layer L α .The inter layer adjacency matrix corresponding to E αβ is the matrix By combining the adjacency matrices corresponding to each layer L α and the inter-layer adjacency matrices appropriately, one can derive a supra-adjacency matrix [23] which characterize the multilayer network M. Now let us define a QW on the multilayer network structure.Recall that, in section II, we developed a mathematical model to mimic the propagation of a QW on any given undirected graph.In our model, we acquire the information of the graphical structure by using the adjacency matrix of the graph.According to our model, when the adjacency matrix is given, we define the sets of {B x } n x=1 along with the set of functions {f x (r)} n x=1 and then simulate the evolution of the QW on the graph by updating the elements of the block matrix in (10).Likewise, one can use the same mathematical model to mimic the propagation of a quantum walker on a multilayer network just by following the same procedure given in section II with the supra-adjacency matrix of the multilayer network.In a QW, the transition probabilities from vertex x to its neighbouring vertices are given by the coin operator C (x) attached to the vertex x.One can choose suitable coin operators according to the simulation to control the probability flow from one vertex to the neighbouring vertices in the multilayer network.When the transition probabilities from a vertex x to its neighbouring vertices are same, we say that the QW on the multilayer network is unbiased.To model such an unbiased QW, we can attach a Fourier coin F (x) to each vertex x given by where d x is the degree of vertex x [24].Note that, the relationship of F (x) (F (x) ) † = (F (x) ) † F (x) = I is preserved by the Fourier coin.In the studies of QWs on graphs, it has been shown that a QW with Grover coin tends to be localized around the initial vertex [24].Hence, it is worth exploring Grover walk on a multilayer network as well.The ij th element of the Grover coin G (x) attached to the vertex x on a multilayer network can be written as where 1 ≤ i, j ≤ d x and d x is the degree of vertex x.Note that, Grover coin holds the relationship of In section V and VI, we analyze the dynamics of some QWs on specific multilayer networks with different choices of coin operators.

IV. CLASSICAL RANDOM WALK ON MULTILAYER NETWORKS
The propagation of the classical random walk (CRW) on different graph structures is a topic that has been studied extensively [25].In general, the propagation of a random walker on any given network structure is modeled using the transition probabilities from a vertex x to its neighbouring vertices [19,26].By adopting the same concept, one can define a CRW on a multilayer network as well [27].Let Ω x,y be the transition probability from vertex x to y.Then, the probability P c (x, t) of finding the random walker at position x at time t is given by the following recurrence relations where d x is the degree of vertex x and for each r, the function f x (r) gives the labels of the neighbouring vertices of x.When the transition probabilities from a vertex x to its neighbouring vertices are same, we say that the CRW is unbiased.In usual practice [26], transition probability Ω (ub) x,y for an unbiased CRW on any graph structure is defined as follows Ω (ub)  x,y = where d x is the degree of vertex x.One can adopt the same definition given in (19) to model an unbiased CRW on a multilayer network.

V. NUMERICAL IMPLEMENTATION ON A TOY MODEL
In this section, we perform a CRW and a QW on an illustrated toy multilayer network structure (Figure 1) and examine the flow of probability through various layers.Our intention is to explore the fundamental differences between classical and quantum dynamics on a mulitlayer network.From the definition of the mulitlayer network given in section III, one can consider diverse configurations of structures with multilayers.Nonetheless, a two-layer network, which consists of two distinct graphs can be understood as the simplest multilayer network.Hence, here we consider one such simplest multilayer structure to perform a CRW and a QW.Since the inter-layer edges of the toy model in Figure 1 link only the vertices representing the same entity in different layers, this network can be classified as a multiplex network which is a special class of multilayer networks [28].Let the top and bottom layers be 3, 4} and V 2 = {5, 6, 7, 8} are the set of vertices in each layer and E 1 , E 2 are the set of edges corresponding to each layer.The supra-adjacency matrix [23] of the multilayer network in Figure 1 can be written as follows where the first (A (1) ) and second (A (2) ) block matrices along the diagonal (i.e.top left corner and bottom right corner) represent the adjacency matrices of the layer L 1 and L 2 respectively.Top right and bottom left block matrices represent the connection between the layers.Using the supra-adjacency matrix in (20), one can define the sets of {B x } 8 x=1 along with the set of functions {f x (r)} 8 x=1 and then simulate the evolution of the QW on the multilayer network by updating the elements of the block matrix in (10).

A. Probability Distribution
In this section, we perform an unbiased QW and CRW on the toy multilayer network structure and investigate the probability of finding the walker on each layer after 100 steps.Recall that, to perform an unbiased QW, one needs to use the Fourier coin given in (16).For the case of unbiased CRW, one needs to use the transition probabilities given in (19).We initialize the CRW from the vertex 1.For the QW, we consider two localized initial states of the form |1⟩ p ⊗ |3⟩ c and |1⟩ p ⊗ |5⟩ c .That is, at the beginning we, place the quantum walker at the position state |1⟩ p (vertex 1) attaching the coin state of |3⟩ c and |5⟩ c to the walker.We select these localized initial states to replicate conditions similar to those in CRW, which enables a meaningful comparison between CRWs and QWs.After 100 steps, we calculate the probability of finding the walker at each layer by summing the probability of finding the walker at each node corresponding to a given layer.
An interesting distinction between the unbiased classical and quantum walkers is illustrated in Figure 2.For the case of CRW, the probabilities of finding the walker on the top and bottom layers eventually stabilize to a steady state as time progresses.This can be seen in Figure 2(a).Conversely, in QW, the quantum walker displays dynamic changes in probability over time (see Figure 2(b) and 2(c)).It is important to note that, in the case of CRWs, the probability of locating the walker on the top layer is consistently higher than that of the bottom layer, which is expected since we initialize the walk from the top layer.However, in QWs, despite the initial placement on the top layer, there exists certain time steps where the probability of locating the walker on the bottom layer surpasses that of the top layer.For instead, in Figure 2(c) one can see certain time steps where the probability of finding the walker on the bottom layer is higher than that of the top layer.In addition, from Figure 2(b) and 2(c), one can identify that even though the initial position state is the same, different initial coin states of the QWs can control the temporarily transition of the quantum walker from top layer to bottom layer.Such a behaviour has no analogy to CRWs.Hence, unlike a classical walker who tend to stay within the initial layer, the ability of a quantum walker to temporarily transition to the bottom layer with higher probability implies that the QW could explore a broader portion of the multilayer structure.This enhanced exploration could be useful when searching through large, complex databases represented as multilayer networks.In section VI, we further examine this behaviour by employing different types of lager mulilayer networks.
In Figure 2(b) and 2(c) we have seen that the quantum probability is fluctuating.It is a widely acknowledged fact that, the unitary characteristic of QWs prevents the quantum walker from achieving a steady state [32].Hence, to obtain an idea of the static picture, one can calculate the time-averaged probability of finding the quantum walker at vertex x defined by where T ∈ N [3].Note that, when T becomes lager, P T (x) becomes a better measure that depicts the static picture.
In Figure 3 we illustrate the time-averaged probabilities for a time period of T = 100 for Fourier and Grover walks.The vertical axis of the heatmaps shows the initial node from which the walker starts the walk and the horizontal axis shows the target node where the walker ends the walk.For both Fourier and Grover walk, we consider two cases.In the first case, we initiate both Fourier and Grover walker from the localized initial state of the form That is, we start both QWs from each node by attaching a uniform superposition of coin states.Then, for each node, the time-averaged probability for a time period of T = 100 is calculated using (21).The results are given in Figure 3(b) and 3(c).In the second case, we repeat the same procedure by using the localized initial state of the form ⟩ c where x ∈ {1, ..., 8}.The results are given in Figure 3(d) and 3(e).For the unbiased CRW, we initiate the walker from each node and calculate the probability of finding the walker at each node after 100 time steps.The result is given in Figure 3(a).
According to the Figure 3, the classical walker tends to stay on the top layer of the toy model after 100 time steps irrespective of its initial node.Since the top layer comprises a complete graph, we can expect this result.Time-average probability profile of the Fourier walk also exhibits a similar behavior like the classical walker when initiated from |ϕ 1 ⟩.That is, the time-average probability of finding the Fourier walker on the top layer is relatively higher than that of the bottom layer irrespective of the initial node.However, when the Fourier walker is initiated from |ϕ 2 ⟩, we can confine the walk to a particular layer with a higher probability, which is visible in Figure 3(d).When initiated from |ϕ 1 ⟩, Grover walk exhibits no significant behaviour, but for the initial condition of |ϕ 2 ⟩ Grover walker tends to stay, with higher probability, at the initial position and at the corresponding position on the other layer.As a result, one can see a sharp line along the anti-diagonal of the grid in Figure 3(e) and two lines parallel to this sharp anti-diagonal line.Hence, one can control the quantum dynamics on the multilayer network by changing the initial state, which has no direct analogy in CRW.

B. Return Probability
Another interesting question one could ask related to CRWs or QWs on a multilayer network is, how long it would take for the walker to return to it's initial position.This could be understood as the recurrence of the walk on the multilayer network.Recurrence in a CRW is characterized by the Pólya number [33], which can be written as where P (x 0 , t) is the probability of finding the walker at the initial node x 0 at time step t.For the condition of P = 1, the walk is identified as recurrent.Otherwise, the walk is called transient.Moreover, the expression can also be used as a definition for the Pólya number as it provides the same criteria for the classification of the walk [34].With the formula given in (23), notion of Pólya number can be extended to the study of recurrence in QWs [34,35].For practical purposes, one can calculate the partial Pólya number using either (22) or ( 23) for a finite number of time steps t = T p rater than extending t → ∞.We calculate the partial Pólya number for CRW, Grover and Fourier walks on the toy multilayer network by choosing a set of finite time steps T p ∈ {1, 5, 10, ..., 200} with a gap of 5 units.Our purpose is to make an estimation of the convergence of the P'olya number for each walk.We initialize the walker from node 1 and for both Grover and Fourier walks, initial coin state is chosen as the uniform superposition of coin states.(That is, 4 shows the convergence of the partial Pólya number for CRW, Grover and Fourier walks.According to Figure 4, Grover and Fourier walkers exhibit recurrence during the progression of the walk on the toy multilayer structure.However, Grover walker returns to it's initial position faster than both Fourier and Classical walker.

C. Impact of decoherence
In this section we study the impact of decoherence on the quantum dynamics of the walker that propagates on the toy mulitlayer network model.With regards to this, we study the impact of decoherence that arises from randomly broken links.In the context of QWs on graphs, broken link decoherence specifically relates to the loss of coherence in a QW due to the imperfections or disruptions in the graph structure [36,37].This can occur when edges in the graph are altered or removed, creating discontinuities in the QW.These disruptions can be caused by various factors, including physical imperfections, noise, or intentional modifications to the graph.We perform a Fourier walk for 100 time steps on the toy mulitlayer network model while breaking the connection between node 1 and 3 with a 0.5 probability at each time step.Afterwards, the mean probability distribution was calculated by averaging over 1000 trials.Figure 5(c) displays the mean probability distribution of the Fourier walk over 100 time steps when subjected to the broken link decoherence model.Additionally, for comparison purposes, the probability distributions of the unbiased CRW and the standard Fourier walk are also presented in Figure 5(a) and 5(c) respectively.
According to the Figure 5(a), the classical walker tends to stay on the top layer with higher probability after 100 time steps, a behaviour which we have consistently seen in Figure 2 and 3. On the other hand, the standard Fourier walk exhibits a very different probability profile compered to a CRW as shown in Figure 5(b).However, when subjected to broken link decoherence model, the classical signature emerges in the average probability distribution of the Fourier walk as depicted in Figure 5(c).Recall that, in this study, we have realized the broken link decoherence model by breaking the connection between the node 1 and 3 of the toy multilayer network model with a 0.5 probability at each time step.That is, we are allowing to alter only a single edge in the toy multilayer network.Yet, the affect of decoherence make a substantial impact on the dynamics of walk eventually converging it to the classical distribution.This implies that QWs on a multilayer network could be very sensitive to decoherence models like broken links.

VI. NUMERICAL IMPLEMENTATION ON SYNTHETIC MULTILAYER NETWORKS
We next apply our model to perform QWs on six different two-layered multiplex networks, each consists of 100 nodes.The top and bottom layers of each multiplex network are constructed from combinations of scale-free (SF), The mean probability distribution of the Fourier walk subjected to broken link decoherence model is given in (c).The mean probability distribution is generated by breaking the connection between the node 1 and 3 of the toy multilayer network model with a 0.5 probability at each time step of the Fourier walk and by averaging over 1000 trials.complete (CP) and star networks with 50 nodes.For the case of SF-SF multiplex network, two different scale-free networks are chosen for top and bottom layers.In addition, the hub node of the star network is taken as the first node.First we perform a Fourier walk on each of the six different two-layered multiplex networks.The walk is initiated from the localized state of For each time step up to 100 time steps, we calculate the probability of finding the walker on each layer by summing the probabilities of finding the walker at each node corresponding to that layer and the results are given in Figure 6.For the comparison purposes, we perform an unbiased CRW on the same two-layered multiplex network structures.The classical walker is initiated from node 1.The results are given in Figure 10.By comparing the results given in Figures 6 and 10, one could conclude that the probability of finding the Fourier and classical walkers on most of the multiplex networks we study here follow a similarly trend.For instead, probability profiles for the combinations of SF-SF, SF-CP and SF-STAR given in Figures 6(a The probability profiles of the Grover walk on the six different multiplex networks are given in Figure 7.The Grover walker is initiated from the localized state of |1⟩ p ⊗ 1 √ d1 d1 r=1 |f 1 (r)⟩ c and for each time step up to 100 time steps, we calculate the probability of finding the walker on each layer by summing the probabilities of finding the walker at each node corresponding to that layer.On average, for the cases of SF-SF and SF-STAR, the Grover walker tends to be on both layers with an equal probability.This can be identified from Figures 7 (a) and 7(c).Grover walker on the CP-STAR multiplex network behaves in a very similar way like a classical walker (see Figure 7 (e) and 10 (e)).A periodic behaviour of the probability of finding the Grover walker on top and bottom layers can be seen on CP-CP and STAR-STAR multiplex networks (see Figure 7 (d) and 7(f)).The periodicity reflects the Grover walker's coherent oscillations between the network layers and could be influenced by the topology and connectivity of the CP-CP and STAR-STAR multiplex network.Further research and analysis are essential to unlock the full potential of these insights for practical quantum applications.
In addition to the exploration on probability profiles, we have studied the recurrence probability of the Fourier, Grover and classical walkers on the six different multiplex networks.We calculate the partial Pólya number for CRW, Grover and Fourier walks by choosing a set of finite time steps T p ∈ {1, 5, 10, ..., 100} with a gap of 5 units.Our purpose is to make an estimation of the convergence of the Pólya number for each walk.We initialize each walker from node 1 and for both Grover and Fourier walks, initial coin state is chosen as the uniform superposition of coin states.(That is, 8 shows the convergence of the partial Pólya number for CRW, Grover and Fourier walks on the six different multiplex networks.From Figure 8, one can identify that the Grover walk exhibits recurrence on most of the network structures studied here.On the other hand, Fourier and classical walkers show no recurrence within 100 time steps.For all the plots in Figure 8, initially, the convergence of the partial Pólya number of the Fourier walk is low compared to that of a classical walker.However, as time elapses, the convergence of the Fourier walk surpasses that of CRW except for the case of STAR-STAR multiplex network.
We apply the broken link decoherence model for the six different multiplex networks as well by following the same procedure described in section 5. To realize the broken link decoherence model, we have removed some edges of the multiplex networks randomly and have calculated the average probability distribution of the QW after 100 time steps by averaging over 1000 trials.In section 5, we observed the emergence of classical signature in the probability distribution even for a single broken link on the toy multilayer network.However, since the number of nodes in the six different multiplex networks are relatively large, we couldn't observe a fast convergence to the classical behaviour when a single edge is broken.However, when the number of broken links increases, the convergence to the classical distribution becomes faster.Hence, the impact of decoherence depends upon the number of broken edges.While our investigation has shed some light on decoherence on multilayer network, it is important to acknowledge that the scope of this research is not fully comprehensive.There remains significant potential for future studies to expand upon these findings to gain a deeper understanding of how decoherence impact on the propagation of the quantum walker on multilayer networks.

VII. DISCUSSION
In this paper we studied the dynamics of discrete-time quantum walks on multilayer networks.We derived recurrence formulae for the coefficients of the wave function of a quantum walker on an undirected graph with finite nodes.Then, by extending these formulae to include extra layers, we developed a simulation to mimic the evolution of the quantum walker on a multilayer network.While multilayer networks have been studied in the context of CTQWs, to the best of our knowledge, there is a lack of literature related to DTQWs on multilayer networks.Hence, the prime objective of this study was to present a comprehensive mathematical framework to model DTQWs on multilayer networks with the aim of bridging this gap.In this regard, we employed our mathematical model to analyze the time-averaged probability and the return probability of the quantum walker on multilayer networks in relation to Fourier and Grover walks.Moreover, we studied the impact of decoherence on the progression of Fourier walk on mulilayer networks.For sake of clarity and readability, first we used a toy muliayer network to conduct our analysis.Later we extended our analysis to much larger synthetic multilayer networks.Our study reveled that the Grover walk on a multilayer network exhibits rich dynamics.For instance, the Grover walker displays a periodic behaviour of occupying the top and bottom layers of a two-layered multiplex network constructed from a complete graph or a star graph.Further research and analysis are essential to unlock the full potential of these insights for practical quantum applications, e.g. for quantum computation and quantum communication.Moreover, in relation to the recurrence probability, the Grover walker returns to the initial position faster than both Fourier and Classical walkers.In the context of QWs, the recurrence probability has a deep link to the localization property, playing a pivotal role in diverse applications, including quantum search algorithms and topological insulators [38].Hence, there seems to be significant potential for future studies related to the return probability on multilayer networks.Another finding of this study is that the QWs on multilayer networks are vulnerable to decoherence arsing from randomly broken links.Multilayer networks with a smaller number of nodes are sensitive to defects even in a single edge.However, tolerance of the QWs for the decoherence arising from a defect in single edge increases as the number of nodes increases.Nonetheless, when there exists more defects in the edges, the probability distribution of the QWs converges to the classical distribution very quickly.As a future extension of this study, one could explore how other forms of decoherence models, like Pauli channels as well as amplitude and phase damping in the coin degree of freedom, make an impact on the QWs on multilayer networks.In summary, we anticipate that the mathematical analysis we have performed here may have a profound influence on a broad spectrum of problems which can be modelled or assisted by DTQWs on multilayer networks.
The state vector given in (B1) can be represented in terms of the block matrix N t as Note that, by updating the elements of N t using the relationships given in (11) and ( 12) one can mimic the evolution of a QW on the graph in Figure 9.
Appendix C: Probability of finding the classical walker on top and bottom layers of six different multiplex networks The following set of graphs shows the probability of finding the classical walker on six different synthetic two-layered multiplex networks with 100 nodes.For the case of SF-SF, two different scale-free networks are chosen.Moreover, the hub node of the star network is taken as the first node.CRW is initiated from vertex 1 and for each time step up to 100 steps, the probability of finding the walker on a given layer is calculated by summing the probabilities of finding the walker at each node corresponding to that layer.
positive integers n = |V | and m = |E| represent the number of nodes and number of edges in G, respectively.The number of edges linked to a particular node v i is referred to as its degree and denoted by d i = n j=1 a ij .Let us now model the propagation of a QW on the graph G.For the sake of convenience, let us relabel the vertices of the graph or in QWs' terminology, the position states of the QW as |x⟩ p where x = 1, . . ., n.Then, the set of vertices V = {|1⟩ p , |2⟩ p , . . ., |n⟩ p } becomes the position basis set of the QW which spans the position Hilbert space H p .Let us denote the subspace spanned by the basis element |x⟩ p in V as H (x) p .That is, H (x) p = span{|x⟩ p } where |x⟩ p ∈ V .Now, for each vertex |x⟩ p , let us assign a coin Hilbert space H (x) c spanned by the coin basis {|r⟩ c | r = 1, . . ., d x } where d x is the degree of the vertex |x⟩ p .Note that, the dimension of the coin Hilbert space H (x) c

FIG. 1 .
FIG. 1. Schematic diagram of a multilayer network with four vertices and two layers named L1 and L2.The set of labels {x, x + 4} where x ∈ {1, 2, 3, 4} represent the same entity in different layers.The solid lines connect the vertices within each layer and the dotted lines connect inter-layer.Since the inter-layer edges link only the vertices representing the same entity in different layers, this network can be classified as a multiplex network which is a special class of multilayer networks.Top layer represents an undirected regular graph and the bottom layer represents an undirected connected graph.Parallel edges and self-loops are excluded here.This schematic diagram is inspired by [28-31].

FIG. 2 .
FIG. 2. Probability of finding the walker on top layer (red solid line) and bottom layer (blue dotted line) for each time step up to 100 steps (a) unbiased CRW is initiated from vertex 1. Fourier walk is initiated from (b) |1⟩p ⊗ |3⟩c and (c) |1⟩p ⊗ |5⟩c.

FIG. 3 .
FIG. 3. Heatmaps are depicted for (a) unbiased CRW, (b,d) Fourier walk and (c,e) Grover walk.The vertical axis shows the initial node from which the walker starts the walk and the horizontal axis shows the target node where the walker ends the walk.Each square corresponding to the Fourier and Grover walks indicates the value of time-averaged probabilities for a time period of T = 100.For the unbiased CRW, each square corresponds to the probability after 100 time steps.Both Fourier and Grover walks in (b) and (c) are initiated from the localized position state of the form |ϕ1⟩ ≡ |x⟩p ⊗ 1 √ dx dx r=1 |fx(r)⟩c.The Fourier and Grover walks in (d) and (e) are initiated from the localized position state of the form |ϕ2⟩ ≡ |x⟩p ⊗ i √ dx |fx(1)⟩c +

FIG. 4 .
FIG.4.Convergence of the partial Pólya number for Grover (green dotted line), Fourier (blue dash line) and Classical (red solid line) walkers on the toy multilayer network.Partial Pólya number is calculated by choosing a set of finite time steps Tp ∈ {1, 5, 10, ..., 200} with a gap of 5 units.Each walk is initiated from the node 1 and for both Grover and Fourier walks, initial coin state is chosen as the uniform superposition of coin states.

FIG. 5 .
FIG. 5. Probability distributions of the (a) unbiased CRW and (b) the standard Fourier walk after 100 time steps on the toy multilayer network.The CRW is initiated from node 1 and Fourier walk is initiated from the localized state of |1⟩p ⊗ |2⟩c.The mean probability distribution of the Fourier walk subjected to broken link decoherence model is given in (c).The mean probability distribution is generated by breaking the connection between the node 1 and 3 of the toy multilayer network model with a 0.5 probability at each time step of the Fourier walk and by averaging over 1000 trials.
) & 10 (a), 6(b) & 10 (b) and 6(c) & 10 (c) exhibit a similar trend.However, for the cases of CP-CP, CP-STAR and STAR-STAR, Fourier walk shows slight differences in its probability profiles in the first few time steps compared to the CRW.Nonetheless, as time elapses, overall trend of the probability profiles of the Fourier walker becomes similar to that of CRW.

FIG. 6 . 1 √ d 1 d 1 r=1
FIG. 6.This figure illustrates the probability of finding the Fourier walker on the top layer (red solid line) and the bottom layer (blue dotted line) of six different two-layered multiplex networks, each consists of 100 nodes.The top and bottom layers of each multiplex network are constructed from combinations of scale-free (SF), complete (CP) and star networks with 50 nodes.(a) SF-SF (b) SF-CP (c) SF-STAR (d) CP-CP (e) CP-STAR and (f) STAR-STAR.For the case of SF-SF, two different scale-free networks are chosen.Moreover, the hub node of the star network is taken as the first node.Fourier walk is initiated from the localized position state of the form |1⟩p ⊗ 1 √ d 1 d 1r=1 |f1(r)⟩c and for each time step up to 100 steps, the probability of finding the walker on a given layer is calculated by summing the probabilities of finding the walker at each node corresponding to that layer.

FIG. 7 . 1 √ d 1 d 1 r=1
FIG. 7.This figure illustrates the probability of finding the Grover walker on the top layer (red solid line) and the bottom layer (blue dotted line) of six different two-layered multiplex networks, each consists of 100 nodes.The top and bottom layers of each multiplex network are constructed from combinations of scale-free (SF), complete (CP) and star networks with 50 nodes.(a) SF-SF (b) SF-CP (c) SF-STAR (d) CP-CP (e) CP-STAR and (f) STAR-STAR.For the case of SF-SF, two different scale-free networks are chosen.Moreover, the hub node of the star network is taken as the first node.Grover walk is initiated from the localized position state of the form |1⟩p ⊗ 1 √ d 1 d 1 r=1 |f1(r)⟩c and for each time step up to 100 steps, the probability of finding the walker on a given layer is calculated by summing the probabilities of finding the walker at each node corresponding to that layer.In cases (d) and (f), the time step has been extended to 200 steps to improve the clarity and visibility of the plot shapes.

FIG. 10 .
FIG. 10.This figure illustrates the probability of finding the classical walker on the top layer (red solid line) and the bottom layer (blue dotted line) of six different two-layered multiplex networks, each consists of 100 nodes.The top and bottom layers of each multiplex network are constructed from combinations of scale-free (SF), complete (CP) and star networks with 50 nodes.(a) SF-SF (b) SF-CP (c) SF-STAR (d) CP-CP (e) CP-STAR and (f) STAR-STAR.For the case of SF-SF, two different scale-free networks are chosen.Moreover, the hub node of the star network is taken as the first node.CRW is initiated from vertex 1 and for each time step up to 100 steps, the probability of finding the walker on a given layer is calculated by summing the probabilities of finding the walker at each node corresponding to that layer.