# Discrete-time Quantum Walk on Multilayer Networks

M. N. Jayakody  Faculty of Engineering and the Institute of Nanotechnology and
Advanced Materials, Bar-Ilan University, Ramat Gan 5290002, Israel  Priodyuti
Pradhan  Department of Computer Science and Engineering, Indian Institute of
Information Technology Raichur, Karnataka - 584135, India  Dana Ben Porath
Faculty of Engineering and the Institute of Nanotechnology and Advanced
Materials, Bar-Ilan University, Ramat Gan 5290002, Israel  E. Cohen  Faculty
of Engineering and the Institute of Nanotechnology and Advanced Materials,
Bar-Ilan University, Ramat Gan 5290002, Israel

###### Abstract

Multilayer network is a potent platform which paves a way to study the
interactions among entities in various networks with multiple types of
relationships. In this study, the dynamics of discrete-time quantum walk on a
multilayer network are explored in detail. We derive recurrence formulae for
the coefficients of the wave function of a quantum walker on an undirected
graph with finite number of nodes. By extending these formulae 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 in relation to
Fourier and Grover walks on multilayer networks. Furthermore, we analyze the
impact of decoherence on the 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 [4, 3].

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 $\mathcal{G}=\\{V,E\\}$ where
$V=\\{v_{1},\ldots,v_{n}\\}$ is the set of vertices (nodes) and
$E=\\{(v_{i},v_{j})\ |\ v_{i},v_{j}\in 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
$\mathcal{G}$ as ${\bf A}\in\mathbb{R}^{n\times n}$ which is defined in the
following way

$a_{ij}=\begin{cases}1&\quad\text{If}\ (v_{i},v_{j})\in E\\\
0&\quad\text{Otherwise}\\\ \end{cases}$ (1)

The positive integers $n=|V|$ and $m=|E|$ represent the number of nodes and
number of edges in $\mathcal{G}$, respectively. The number of edges linked to
a particular node $v_{i}$ is referred to as its degree and denoted by
$d_{i}=\sum_{j=1}^{n}a_{ij}$. Let us now model the propagation of a QW on the
graph $\mathcal{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\rangle_{p}$ where $x=1,\ldots,n$. Then, the set of vertices
$V=\\{|1\rangle_{p},|2\rangle_{p},\ldots,|n\rangle_{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\rangle_{p}$ in $V$
as $H_{p}^{(x)}$. That is, $H_{p}^{(x)}=\text{span}\\{|x\rangle_{p}\\}$ where
$|x\rangle_{p}\in V$. Now, for each vertex $|x\rangle_{p}$, let us assign a
coin Hilbert space $H_{c}^{(x)}$ spanned by the coin basis $\\{|r\rangle_{c}\
|\ r=1,\ldots,d_{x}\\}$ where $d_{x}$ is the degree of the vertex
$|x\rangle_{p}$. Note that, the dimension of the coin Hilbert space
$H_{c}^{(x)}$ 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 $\mathcal{B}_{x}$ that comprises the labels of the
vertices adjacent to the vertex $|x\rangle_{p}$ as $\mathcal{B}_{x}=\\{y\ :\
\text{vertex}\ |x\rangle_{p}\ \text{and}\ |y\rangle_{p}\ \text{have a
connection}\\}$. Note that, $|\mathcal{B}_{x}|=d_{x}$. Now, define a function
as, $f_{x}(r)=r^{th}\text{ element of}\ \mathcal{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 $\mathcal{B}_{x}$ in the
ascending order. Now, we can denote the basis of the coin Hilbert space
$H_{c}^{(x)}$ as $\\{|f_{x}(r)\rangle_{c}\ |\ r=1,\ldots,d_{x}\\}$. 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\rangle_{p}$ at time $t$ can be written as

$|\psi(x,t)\rangle=\sum_{r=1}^{d_{x}}\alpha_{x,f_{x}(r)}(t)|x\rangle_{p}|f_{x}(r)\rangle_{c}$
(2)

where $\alpha_{x,f_{x}(r)}(t)\in\mathbb{C}$ are called the probability
amplitudes, $|x\rangle_{p}\in H_{p}^{(x)}$, $|f_{x}(r)\rangle_{c}\in
H_{c}^{(x)}$ and $|\psi(x,t)\rangle\in H_{p}^{(x)}\otimes H_{c}^{(x)}$. Note
that, $dim(H_{p}^{(x)})=1$, $dim(H_{c}^{(x)})=d_{x}$ and
$dim(H_{p}^{(x)}\otimes H_{c}^{(x)})=d_{x}$. For each position $x$, the
coefficients of $\alpha_{x,f_{x}(r)}(t)$ are defined as follows

$\alpha_{x,r}(t)=\begin{cases}\alpha_{x,r}(t)&\quad\text{If}\
r\in\mathcal{B}_{x}\\\ 0&\quad\text{Otherwise}\\\ \end{cases}$ (3)

Our next task is to write an expression for the total wave function of the
quantum walker on $\mathcal{G}$ at time $t$. For that, we need to sum the
state vectors $|\psi(x,t)\rangle$ in (2) over all the vertices. However, we
are unable to perform such a summation because the state vectors
$|\psi(x,t)\rangle$ corresponding to each vertex resides in different
composite Hilbert spaces. Note that, the size of the coin Hilbert space
$H_{c}^{(x)}$ 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_{p}^{(x)}\otimes H_{c}^{(x)}\\}_{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
$H=\bigoplus\limits_{x=1}^{n}\bigg{(}H_{p}^{(x)}\otimes H_{c}^{(x)}\bigg{)}$
where $\oplus$ denotes the external direct sum of Hilbert spaces and
$\textit{dim}(H)=\sum_{x=1}^{n}d_{x}$. According to the definition of the
global Hilbert space $H$, it is obvious that for each vertex $x$, the state
vector $|\psi(x,t)\rangle\in H$. Hence, the total wave function of the quantum
walker at time $t$ can be calculated by summing $|\psi(x,t)\rangle$ in (2)
over all the vertices. Then, the total wave function at time $t$ can be
written as

$|\psi_{t}\rangle=\sum_{x=1}^{n}\sum_{r=1}^{d_{x}}\alpha_{x,f_{x}(r)}(t)|x\rangle_{p}|f_{x}(r)\rangle_{c}$
(4)

where $|\psi_{t}\rangle\in H$ and
$\sum_{x=1}^{n}\sum_{r=1}^{d_{x}}|\alpha_{x,f_{x}(r)}(t)|^{2}=1$. Coin
operator $C^{(x)}$ which acts on the coin states associated to the vertex
$|x\rangle_{p}$ holds the transition probabilities from $|x\rangle_{p}$ to its
neighbouring vertices. Hence, $C^{(x)}$ can be defined as

$C^{(x)}=\sum_{i=1}^{d_{x}}\sum_{j=1}^{d_{x}}C^{(x)}_{ij}|f_{x}(i)\rangle\langle
f_{x}(j)|$ (5)

where $\\{|f_{x}(r)\rangle_{c}\\}_{r=1}^{d_{x}}$ are the basis elements of
$H_{c}^{(x)}$. The coin coefficients $C^{(x)}_{ij}\in\mathbb{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)})^{\dagger}C^{(x)}=C^{(x)}(C^{(x)})^{\dagger}=\mathbb{I}$. By
combining each local coin operator $C^{(x)}$, one can write the global coin
operator $C$ acting on $|\psi_{t}\rangle\in H$ as follows

$C=\sum_{x=1}^{n}\sum_{i=1}^{d_{x}}\sum_{j=1}^{d_{x}}C^{(x)}_{ij}|x\rangle\langle
x|\otimes|f_{x}(i)\rangle\langle f_{x}(j)|$ (6)

The shift operator of the QW on a graph is defined as follows

$S|x\rangle_{p}|y\rangle_{c}=|y\rangle_{p}|x\rangle_{c}$ (7)

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

$|\psi_{t+1}\rangle=U|\psi_{t}\rangle$ (8)

where $U=SC$ is the evolution operator and $S$ and $C$ are the shift and coin
operators respectively.

###  II.1 Matrix representation and simulation

The local coin operator $C^{(x)}$, associated to the vertex $|x\rangle_{p}$
holds the transition probabilities from $|x\rangle_{p}$ to its neighbouring
nodes. Suppose the vertex $|x\rangle_{p}$ is linked to $d_{x}$ number of
vertices denoted by $|y_{1}\rangle_{p},\ldots,|y_{d_{x}}\rangle_{p}$.
Moreover, suppose that $y_{1}<\ldots<y_{d_{x}}$. Hence, for each
$r\in\\{1,\ldots,d_{x}\\}$, we can write $y_{r}=f_{x}(r)$. Then, the block
matrix of $C^{(x)}$ can be written as

$C^{(x)}=\bordermatrix{&\langle f_{x}(1)|&\ldots&\langle
f_{x}(d_{x})|\cr|f_{x}(1)\rangle&C^{(x)}_{11}&\ldots&C^{(x)}_{1d_{x}}\cr\vdots&\vdots&\vdots&\vdots\cr|f_{x}(d_{x})\rangle&C^{(x)}_{d_{x}1}&\ldots&C^{(x)}_{d_{x}d_{x}}\cr}$
(9)

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

$|\psi_{t}\rangle={\bf
N}_{t}=\bordermatrix{&|1\rangle_{c}&\ldots&|n\rangle_{c}\cr|1\rangle_{p}&\alpha_{1,1}(t)&\ldots&\alpha_{1,n}(t)\cr\vdots&\vdots&\vdots&\vdots\cr|n\rangle_{p}&\alpha_{n,1}(t)&\ldots&\alpha_{n,n}(t)\cr}$
(10)

In matrix ${\bf 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 $\alpha_{x,r}(t)$, given in (3), some entries of the
matrix ${\bf 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 $\alpha_{x,r}(t)$, in other words, all the
elements of ${\bf 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 ${\bf N}_{t}$ appropriately, one can determine the matrix ${\bf
N}_{t+1}$ and hence the total wave function of the quantum walker at time
$t+1$. The elements of ${\bf N}_{t}$ are updated in two processes. First, an
intermediate matrix ${\bf\widetilde{N}}_{t}$ is generated using the following
formula

$\tilde{\alpha}_{x,f_{x}(i)}(t)=\sum_{j=1}^{d_{x}}{\alpha_{x,f_{x}(j)}(t)C^{(x)}_{ij}}$
(11)

Afterwards, the matrix ${\bf N}_{t+1}$ is determined by taking the transpose
of ${\bf\widetilde{N}}_{t}$. The update rule is given by

$\alpha_{x,f_{x}(i)}(t+1)=\tilde{\alpha}_{f_{x}(i),x}(t)$ (12)

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.

###  II.2 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

$P_{q}(x,t)=\sum_{r=1}^{d_{x}}|\alpha_{x,f_{x}(r)}(t)|^{2}$ (13)

Note that, $P_{q}(x,t)$ can be determined by summing the elements in the
$x^{th}$ row of the matrix ${\bf N}_{t}\odot{\bf N}_{t}^{*}$ where $\odot$ is
the Hadamard product of matrices and ${\bf N}_{t}^{*}$ is the complex
conjugate of the block matrix given in (10).

##  III QWs on Multilayer networks

A multilayer network is a pair $\mathcal{M}=(\mathcal{G},\mathcal{C})$ where
$\mathcal{G}=\\{\mathcal{L}_{\alpha};\;\alpha\in\\{1,2,\ldots,l\\}\\}$ is a
family of undirected graphs of
$\mathcal{L}_{\alpha}=\\{V_{\alpha},E_{\alpha}\\}$ (called layers of
$\mathcal{M}$) with
$V_{\alpha}=\\{v_{1}^{\alpha},v_{2}^{\alpha},\ldots,v_{n_{\alpha}}^{\alpha}\\}$
is the set of vertices and
$E_{\alpha}=\\{e_{1}^{\alpha},e_{2}^{\alpha},\ldots,e_{r_{\alpha}}^{\alpha}:e_{r_{\alpha}}=(v_{i}^{\alpha},v_{j}^{\alpha})\\}$
is the set of edges in the $\mathcal{L}_{\alpha}$ layer of the multilayer
network [6]. The positive integers $l$ , $n_{\alpha}$ and $r_{\alpha}$ are
termed as the number of layers in $\mathcal{M}$, number of vertices and edges
of the layer $\mathcal{L}_{\alpha}$ respectively. Moreover,
$\mathcal{C}=\\{E_{\alpha\beta}\subseteq V_{\alpha}\times
V_{\beta}:\;\alpha,\beta\in\\{1,2,\ldots,l\\},\alpha\neq\beta\\}$ is the set
of edges between $\mathcal{L}_{\alpha}$ and $\mathcal{L}_{\beta}$ layers. The
elements of $\mathcal{C}$ are called crossed layers. Further, the elements of
each $E_{\alpha}$ are called the set of intralayer edges and the elements of
each $E_{\alpha\beta}$ ($\alpha\neq\beta$) are called the interlayer edges of
$\mathcal{M}$ [23]. Let us denote the adjacency matrices corresponding to each
layer $\mathcal{L}_{\alpha}$ as
$A^{(\alpha)}=(a^{\alpha}_{ij})\in\mathbb{R}^{n_{\alpha}\times n_{\alpha}}$
which is defined by

$a^{\alpha}_{ij}=\begin{cases}1&\quad\text{If}\
(v_{i}^{\alpha},v_{j}^{\alpha})\in E_{\alpha}\\\
0&\quad\text{Otherwise}\end{cases}$ (14)

for $1\leq i,j\leq n_{\alpha}$ and $1\leq\alpha\leq l$ where $n_{\alpha}$ is
the number of nodes in layer $\mathcal{L}_{\alpha}$. The inter layer adjacency
matrix corresponding to $E_{\alpha\beta}$ is the matrix
$A^{[\alpha,\beta]}\in\mathbb{R}^{n_{\alpha}\times n_{\beta}}$ given by

$a^{\alpha\beta}_{ij}=\begin{cases}1&\quad\text{If}\
(v_{i}^{\alpha},v_{j}^{\beta})\in E_{\alpha\beta}\\\
0&\quad\text{Otherwise}\end{cases}$ (15)

By combining the adjacency matrices corresponding to each layer
$\mathcal{L}_{\alpha}$ and the inter-layer adjacency matrices appropriately,
one can derive a supra-adjacency matrix [23] which characterize the multilayer
network $\mathcal{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
$\\{\mathcal{B}_{x}\\}_{x=1}^{n}$ along with the set of functions
$\\{f_{x}(r)\\}_{x=1}^{n}$ 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

$F^{(x)}=\frac{1}{\sqrt{d_{x}}}\sum_{r=1}^{d_{x}}\sum_{s=1}^{d_{x}}e^{2i\pi(r-1)(s-1)/d_{x}}|r\rangle\langle
s|$ (16)

where $d_{x}$ is the degree of vertex $x$ [24]. Note that, the relationship of
$F^{(x)}(F^{(x)})^{\dagger}=(F^{(x)})^{\dagger}F^{(x)}=\mathbb{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

$G^{(x)}_{ij}=\begin{cases}\frac{2-d_{x}}{d_{x}}&\quad\text{If}\ i=j\\\
\frac{2}{d_{x}}&\quad\text{Otherwise}\end{cases}$ (17)

where $1\leq i,j\leq d_{x}$ and $d_{x}$ is the degree of vertex $x$. Note
that, Grover coin holds the relationship of
$G^{(x)}(G^{(x)})^{\dagger}=(G^{(x)})^{\dagger}G^{(x)}=\mathbb{I}$. 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
[26, 19]. By adopting the same concept, one can define a CRW on a multilayer
network as well [27]. Let $\Omega_{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

$P_{c}(x,t)=\sum_{r=1}^{d_{x}}\Omega_{f_{x}(r),x}P_{c}(f_{x}(r),t-1)$ (18)

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 $\Omega^{(ub)}_{x,y}$ for an unbiased CRW on any graph structure
is defined as follows

$\Omega^{(ub)}_{x,y}=\begin{cases}\frac{1}{d_{x}}&\quad\text{If}\ (x,y)\
\text{is connected}\\\ 0&\quad\text{Otherwise}\\\ \end{cases}$ (19)

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 $\mathcal{L}_{1}=\\{V_{1},E_{1}\\}$ and
$\mathcal{L}_{2}=\\{V_{2},E_{2}\\}$ respectively where $V_{1}=\\{1,2,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

$A_{sup}=\begin{pmatrix}\mbox{\Large$A^{(1)}$}&\vline&\mbox{\Large$\mathbb{I}$}\\\
\hline\cr\mbox{\Large$\mathbb{I}$}&\vline&\mbox{\Large$A^{(2)}$}\end{pmatrix}=\begin{pmatrix}\begin{matrix}0&1&1&1\\\
1&0&1&1\\\ 1&1&0&1\\\ 1&1&1&0\\\ \end{matrix}&\vline&\begin{matrix}1&0&0&0\\\
0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\ \end{matrix}\\\
\hline\cr\begin{matrix}1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\ 0&0&0&1\\\
\end{matrix}&\vline&\begin{matrix}0&1&0&1\\\ 1&0&1&0\\\ 0&1&0&1\\\ 1&0&1&0\\\
\end{matrix}\end{pmatrix}$ (20)

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 $\mathcal{L}_{1}$ and $\mathcal{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 $\\{\mathcal{B}_{x}\\}_{x=1}^{8}$ along with the set of
functions $\\{f_{x}(r)\\}_{x=1}^{8}$ and then simulate the evolution of the QW
on the multilayer network by updating the elements of the block matrix in
(10).

Figure 1: Schematic diagram of a multilayer network with four vertices and two
layers named $L_{1}$ and $L_{2}$. The set of labels $\\{x,x+4\\}$ where
$x\in\\{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 [29, 30, 31, 28].

###  V.1 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\rangle_{p}\otimes|3\rangle_{c}$ and $|1\rangle_{p}\otimes|5\rangle_{c}$.
That is, at the beginning we, place the quantum walker at the position state
$|1\rangle_{p}$ (vertex 1) attaching the coin state of $|3\rangle_{c}$ and
$|5\rangle_{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.

Figure 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\rangle_{p}\otimes|3\rangle_{c}$ and (c)
$|1\rangle_{p}\otimes|5\rangle_{c}$

.

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

$\overline{P}_{T}(x)=\frac{1}{T}\sum_{t=0}^{T-1}P_{q}(x,t)$ (21)

where $T\in\mathbb{N}$ [3]. Note that, when $T$ becomes lager,
$\overline{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
$|\phi_{1}\rangle\equiv|x\rangle_{p}\otimes\frac{1}{\sqrt{d_{x}}}\sum_{r=1}^{d_{x}}|f_{x}(r)\rangle_{c}$
where $x\in\\{1,...,8\\}$. 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
$|\phi_{2}\rangle\equiv|x\rangle_{p}\otimes\bigg{(}\frac{i}{\sqrt{d_{x}}}|f_{x}(1)\rangle_{c}+\frac{1}{\sqrt{d_{x}}}\sum_{r=2}^{d_{x}-1}|f_{x}(r)\rangle_{c}-\frac{i}{\sqrt{d_{x}}}|f_{x}(d_{x})\rangle_{c}\bigg{)}$
where $x\in\\{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).

Figure 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
$|\phi_{1}\rangle\equiv|x\rangle_{p}\otimes\frac{1}{\sqrt{d_{x}}}\sum_{r=1}^{d_{x}}|f_{x}(r)\rangle_{c}$.
The Fourier and Grover walks in (d) and (e) are initiated from the localized
position state of the form
$|\phi_{2}\rangle\equiv|x\rangle_{p}\otimes\bigg{(}\frac{i}{\sqrt{d_{x}}}|f_{x}(1)\rangle_{c}+\frac{1}{\sqrt{d_{x}}}\sum_{r=2}^{d_{x}-1}|f_{x}(r)\rangle_{c}-\frac{i}{\sqrt{d_{x}}}|f_{x}(d_{x})\rangle_{c}\bigg{)}$.
Note that, both localized initial states $|\phi_{1}\rangle$ and
$|\phi_{2}\rangle$ have a uniform superposition of coin states, yet
$|\phi_{2}\rangle$ contains some complex coefficients.

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 $|\phi_{1}\rangle$.
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
$|\phi_{2}\rangle$, we can confine the walk to a particular layer with a
higher probability, which is visible in Figure 3(d). When initiated from
$|\phi_{1}\rangle$, Grover walk exhibits no significant behaviour, but for the
initial condition of $|\phi_{2}\rangle$ 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.

###  V.2 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

$P=1-\frac{1}{\sum_{t=1}^{\infty}P(x_{0},t)}$ (22)

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

$P=1-\Pi_{t=1}^{\infty}[1-P(x_{0},t)]$ (23)

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
[35, 34]. 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\rightarrow\infty$. 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}\in\\{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,
$|1\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{c}$).
Figure 4 shows the convergence of the partial Pólya number for CRW, Grover and
Fourier walks.

Figure 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 $T_{p}\in\\{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.

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.

###  V.3 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.

Figure 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\rangle_{p}\otimes|2\rangle_{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.

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),
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
$|1\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{c}$.
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) & 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.

Figure 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\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{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.

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\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{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}\in\\{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,
$|1\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{c}$).
Figure 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.

Figure 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\rangle_{p}\otimes\frac{1}{\sqrt{d_{1}}}\sum_{r=1}^{d_{1}}|f_{1}(r)\rangle_{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.

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.

Figure 8: Convergence of the partial Pólya number for Grover (green dotted
line), Fourier (blue dash line) and Classical (red solid line) walkers 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), 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. Partial Pólya number is
calculated by choosing a set of finite time steps
$T_{p}\in\\{1,5,10,...,100\\}$ 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.

##  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.

##  VIII Acknowledgement

M. N. Jayakody acknowledges the President Scholarship Program at Bar-Ilan
University. P. Pradhan acknowledges Science and Engineering Research Board
(SERB) grant TAR/2022/000657, Govt. of India. This research was funded in part
by the Israeli Innovation Authority under Project No. 73795, by the Pazy
Foundation, by the Israeli Ministry of Science and Technology, and by the
Quantum Science and Technology Program of the Israeli Council of Higher
Education.

## References

  * Childs  _et al._ [2013] A. M. Childs, D. Gosset, and Z. Webb, Universal computation by multiparticle quantum walk, Science 339, 791 (2013).
  * Kempe [2003] J. Kempe, Quantum random walks: an introductory overview, Contemporary Physics 44, 307 (2003).
  * Wang and Manouchehri [2013] J. Wang and K. Manouchehri,  _Physical implementation of quantum walks_ (Springer, 2013).
  * Venegas-Andraca [2012] S. E. Venegas-Andraca, Quantum walks: a comprehensive review, Quantum Information Processing 11, 1015 (2012).
  * Aleta and Moreno [2019] A. Aleta and Y. Moreno, Multilayer networks in a nutshell, Annual Review of Condensed Matter Physics 10, 45 (2019).
  * Jalan and Pradhan [2018] S. Jalan and P. Pradhan, Localization of multilayer networks by optimized single-layer rewiring, Physical Review E 97, 042314 (2018).
  * Pilosof  _et al._ [2017] S. Pilosof, M. A. Porter, M. Pascual, and S. Kéfi, The multilayer nature of ecological networks, Nature Ecology & Evolution 1, 1 (2017).
  * Gosak  _et al._ [2018] M. Gosak, R. Markovič, J. Dolenšek, M. S. Rupnik, M. Marhl, A. Stožer, and M. Perc, Network science of biological systems at different scales: A review, Physics of life reviews 24, 118 (2018).
  * Gallotti and Barthelemy [2015] R. Gallotti and M. Barthelemy, The multilayer temporal network of public transport in great britain, Scientific data 2, 1 (2015).
  * Gallotti and Barthelemy [2014] R. Gallotti and M. Barthelemy, Anatomy and efficiency of urban multimodal mobility, Scientific reports 4, 1 (2014).
  * Musmeci  _et al._ [2017] N. Musmeci, V. Nicosia, T. Aste, T. Di Matteo, and V. Latora, The multiplex dependency structure of financial markets, Complexity 2017 (2017).
  * Galiceanu and Strunz [2016] M. Galiceanu and W. T. Strunz, Continuous-time quantum walks on multilayer dendrimer networks, Physical Review E 94, 022307 (2016).
  * Maquiné Batalha  _et al._ [2022] G. Maquiné Batalha, A. Volta, W. T. Strunz, and M. Galiceanu, Quantum transport on honeycomb networks, Scientific Reports 12, 1 (2022).
  * Maciel  _et al._ [2020] C. M. Maciel, C. F. Mendes, W. T. Strunz, and M. Galiceanu, Quantum transport on generalized scale-free networks, Physical Review A 102, 032219 (2020).
  * Ambainis [2003] A. Ambainis, Quantum walks and their algorithmic applications, International Journal of Quantum Information 1, 507 (2003).
  * Aharonov  _et al._ [2001] D. Aharonov, A. Ambainis, J. Kempe, and U. Vazirani, Quantum walks on graphs, in  _Proceedings of the thirty-third annual ACM symposium on Theory of computing_ (2001) pp. 50–59.
  * Godsil and Zhan [2019] C. Godsil and H. Zhan, Discrete-time quantum walks and graph structures, Journal of Combinatorial Theory, Series A 167, 181 (2019).
  * Watrous [2001] J. Watrous, Quantum simulations of classical random walks and undirected graph connectivity, Journal of computer and system sciences 62, 376 (2001).
  * Kendon [2006] V. Kendon, Quantum walks on general graphs, International Journal of Quantum Information 4, 791 (2006).
  * Feldman and Hillery [2004] E. Feldman and M. Hillery, Scattering theory and discrete-time quantum walks, Physics Letters A 324, 277 (2004).
  * Jayakody  _et al._ [2021] M. N. Jayakody, C. Meena, and P. Pradhan, One-dimensional discrete-time quantum walks with general coin, arXiv preprint arXiv:2102.07207 (2021).
  * Manouchehri and Wang [2009] K. Manouchehri and J. Wang, Quantum random walks without walking, Physical Review A 80, 060304 (2009).
  * Boccaletti  _et al._ [2014] S. Boccaletti, G. Bianconi, R. Criado, C. I. Del Genio, J. Gómez-Gardenes, M. Romance, I. Sendina-Nadal, Z. Wang, and M. Zanin, The structure and dynamics of multilayer networks, Physics reports 544, 1 (2014).
  * Mukai and Hatano [2020] K. Mukai and N. Hatano, Discrete-time quantum walk on complex networks for community detection, Physical Review Research 2, 023378 (2020).
  * Lovász [1993] L. Lovász, Random walks on graphs, Combinatorics, Paul erdos is eighty 2, 4 (1993).
  * Venegas-Andraca [2008] S. E. Venegas-Andraca, Quantum walks for computer scientists, Synthesis Lectures on Quantum Computing 1, 1 (2008).
  * Baptista  _et al._ [2022] A. Baptista, A. Gonzalez, and A. Baudot, Universal multilayer network exploration by random walk with restart, Communications Physics 5, 170 (2022).
  * Gomez  _et al._ [2013] S. Gomez, A. Diaz-Guilera, J. Gomez-Gardenes, C. J. Perez-Vicente, Y. Moreno, and A. Arenas, Diffusion dynamics on multiplex networks, Physical review letters 110, 028701 (2013).
  * Li and Patra [2010] Y. Li and J. C. Patra, Genome-wide inferring gene–phenotype relationship by walking on the heterogeneous network, Bioinformatics 26, 1219 (2010).
  * Granell  _et al._ [2013] C. Granell, S. Gómez, and A. Arenas, Dynamical interplay between awareness and epidemic spreading in multiplex networks, Physical review letters 111, 128701 (2013).
  * Granell  _et al._ [2014] C. Granell, S. Gómez, and A. Arenas, Competing spreading processes on multiplex networks: awareness and epidemics, Physical review E 90, 012808 (2014).
  * Faccin  _et al._ [2013] M. Faccin, T. Johnson, J. Biamonte, S. Kais, and P. Migdał, Degree distribution in quantum walks on complex networks, Physical Review X 3, 041007 (2013).
  * Pólya [1921] G. Pólya, Über eine Aufgabe der Wahrscheinlichkeitsrechnung betreffend die Irrfahrt im Straßennetz, Math. Ann. 84, 149 (1921).
  * Štefaňák  _et al._ [2008a] M. Štefaňák, T. Kiss, and I. Jex, Recurrence properties of unbiased coined quantum walks on infinite $d$-dimensional lattices, Phys. Rev. A 78, 032306 (2008a).
  * Štefaňák  _et al._ [2008b] M. Štefaňák, I. Jex, and T. Kiss, Recurrence and pólya number of quantum walks, Physical review letters 100, 020501 (2008b).
  * Leung  _et al._ [2010] G. Leung, P. Knott, J. Bailey, and V. Kendon, Coined quantum walks on percolation graphs, New Journal of Physics 12, 123018 (2010).
  * Kollár  _et al._ [2012] B. Kollár, T. Kiss, J. Novotnỳ, and I. Jex, Asymptotic dynamics of coined quantum walks on percolation graphs, Physical review letters 108, 230505 (2012).
  * Kiumi  _et al._ [2022] C. Kiumi, N. Konno, and S. Tamura, Return probability of quantum and correlated random walks, Entropy 24, 584 (2022).

##  Appendix A Coin flip and Shift operation

Let us represent the state vector $|\psi_{t}\rangle$ given by (4) in the block
matrix form. That is $|\psi_{t}\rangle\equiv{\bf N}_{t}$. Now let us apply the
global coin operator $C$ given in (6) on the total wave function
$|\psi_{t}\rangle$ in (4). Then we get the following expression

$C|\psi_{t}\rangle=\sum_{x=1}^{n}\sum_{i=1}^{d_{x}}\biggl{(}\sum_{j=1}^{d_{x}}\alpha_{x,f_{x}(j)}(t)C^{(x)}_{ij}\biggr{)}|x\rangle_{p}|f_{x}(i)\rangle_{c}$
(24)

Let us write
$\tilde{\alpha}_{x,f_{x}(i)}(t)=\biggl{(}\sum_{j=1}^{d_{x}}\alpha_{x,f_{x}(j)}(t)C^{(x)}_{ij}\biggr{)}$.
Then we can rewrite (24) as follows

$C|\psi_{t}\rangle=\sum_{x=1}^{n}\sum_{i=1}^{d_{x}}\tilde{\alpha}_{x,f_{x}(i)}(t)|x\rangle_{p}|f_{x}(i)\rangle_{c}$
(25)

Note that, $C|\psi_{t}\rangle$ can also be represented in the block matrix
form. Let us denote $C|\psi_{t}\rangle\equiv{\bf\widetilde{N}}_{t}$. Now let
us apply the shift operator $S$ given in (7) on $C|\psi_{t}\rangle$. The
result can be written as

$SC|\psi_{t}\rangle=\sum_{x=1}^{n}\sum_{i=1}^{d_{x}}\tilde{\alpha}_{x,f_{x}(i)}(t)|f_{x}(i)\rangle_{p}|x\rangle_{c}$
(26)

Let the block matrix form of $SC|\psi_{t}\rangle$ be ${\bf{M}}_{t}$. By
observing (25) and (26) one can argue that the elements in the $x^{th}$ row of
${\bf\widetilde{N}}_{t}$ is equal to the elements in the $x^{th}$ column of
${\bf M}_{t}$. That is, ${\bf{M}}_{t}$ is equal to the transpose of
${\bf\widetilde{N}}_{t}$. In addition, we have the following expression

$SC|\psi_{t}\rangle=|\psi_{t+1}\rangle=\sum_{x=1}^{n}\sum_{i=1}^{d_{x}}\alpha_{x,f_{x}(i)}(t+1)|x\rangle_{p}|f_{x}(i)\rangle_{c}$
(27)

Thus, for each $x\in\\{1,\ \ldots\ n\\}$ and $r\in\\{1,\ \ldots\ d_{x}\\}$ we
can derive a recurrence relationship between the elements of ${\bf N}_{t}$ and
${\bf N}_{t+1}$ as follows

$\alpha_{x,f_{x}(r)}(t+1)=\tilde{\alpha}_{f_{x}(r),x}(t)$ (28)

This completes the proof.

##  Appendix B Example of graph with four nodes

Let us consider the graph given in Figure 9. The corresponding set of vertices
and edges can be written as
$V=\\{|1\rangle_{p},|2\rangle_{p},|3\rangle_{p},|4\rangle_{p}\\}$ and
$E=\\{\\{1,2\\},\\{1,3\\},\\{2,1\\},\\{2,3\\},\\{3,1\\},\\{3,2\\},\\{3,4\\},\\{4,3\\}\\}$
respectively.

Figure 9: Schematic diagram of a graph of four vertexes and the coin operators
associated with each vertex. For each vertex $x$, the elements of the coin
operator associated to $x$ are chosen in such a way that
$C^{(x)}(C^{(x)})^{\dagger}=(C^{(x)})^{\dagger}C^{(x)}=\mathbb{I}$ where
$x\in\\{1,2,3,4\\}$

Now let us define the sets of $\mathcal{B}_{x}$ as follows;
$\mathcal{B}_{1}=\\{2,\ 3\\}$, $\mathcal{B}_{2}=\\{1,\ 3\\}$,
$\mathcal{B}_{3}=\\{1,\ 2,\ 4\\}$ and $\mathcal{B}_{4}=\\{3\\}$. The function
$f_{x}(r)$ ($r^{th}$ element of the $\mathcal{B}_{x}$ set) can be defined as
follows; $f_{1}(1)=2,\ f_{1}(2)=3$; $f_{2}(1)=1,\ f_{2}(2)=3$; $f_{3}(1)=1,\
f_{3}(2)=2,\ f_{3}(3)=4$ and $f_{4}(1)=3$. Then, the state vector of the
quantum walker at time $t$ can be written as

$|\psi_{t}\rangle=\sum_{x=1}^{4}\sum_{r=1}^{d_{x}}\alpha_{x,f_{x}(r)}(t)|x\rangle_{p}|f_{x}(r)\rangle_{c}$
(29)

The state vector given in (29) can be represented in terms of the block matrix
${\bf N}_{t}$ as

${\bf
N}_{t}=\bordermatrix{&|1\rangle_{c}&|2\rangle_{c}&|3\rangle_{c}&|4\rangle_{c}\cr|1\rangle_{p}&0&\alpha_{1,2}(t)&\alpha_{1,3}(t)&0\cr|2\rangle_{p}&\alpha_{2,1}(t)&0&\alpha_{2,3}(t)&0\cr|3\rangle_{p}&\alpha_{3,1}(t)&\alpha_{3,2}(t)&0&\alpha_{3,4}(t)\cr|4\rangle_{p}&0&0&\alpha_{4,3}(t)&0\cr}$
(30)

Note that, by updating the elements of ${\bf 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.

Figure 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.