# Time-inhomogeneous Quantum Markov Chains with Decoherence on Finite State
Spaces

###### Abstract

We introduce and study time-inhomogeneous quantum Markov chains with parameter
$\zeta\geq 0$ and decoherence parameter $0\leq p\leq 1$ on finite spaces and
their large scale equilibrium properties. Here $\zeta$ resembles the inverse
temperature in the annealing random process and $p$ is the decoherence
strength of the quantum system. Numerical evaluations show that if $\zeta$ is
small, then quantum Markov chain is ergodic for all $0<p\leq 1$ and if $\zeta$
is large, then it has multiple limiting distributions for all $0<p\leq 1$. In
this paper, we prove the ergodic property in the high temperature region
$0\leq\zeta\leq 1$. We expect that the phase transition occurs at the critical
point $\zeta_{c}=1$. For coherence case $p=0$, a critical behavior of
periodicity also appears at critical point $\zeta_{o}=2$.

Chia-Han Chou and Wei-Shih Yang Department of Mathematics, Temple University,
Philadelphia, PA 19122

Email<EMAIL_ADDRESS><EMAIL_ADDRESS>

KEY WORDS: Quantum Markov Chain; Quantum Decoherence; Time-Inhomogeneous

AMS classification Primary: 81Q10, Secondary: 82B10

PACS numbers: 03.65.Yz, 05.30.-d, 03.67.Lx, 02.50.Ga

##  1 Introduction

In order to develop more efficient algorithms for tackling a wide variety of
problems in classical computer science, researchers started utilizing
randomness techniques such as Ulam and von Nuemann’s Markov Chain Monte Carlo
(MCMC) method [14] in 1940s. This method was later refined and made well known
as the Metropolis-Hastings algorithm [10] with applications in different
areas. Even though Monte Carlos methods could sometimes return incorrect
solutions with given probability, the key idea behind the methods was that the
true solution can be approximated with high probability by repeating Monte
Carlo simulations.

More recently, the notion of quantum computation has gained popularity,
”qubit” takes a complex unit instead of ”bit” the usual binary values of zero
and one. To preserve a cohesive quantum system, the family of qubits
comprising the memory of the computer go through unitary evolution, rather
than the traditional system of gates in classical computation theory. The
state of the quantum system can be observed, and collapsing the system to one
unique state from a superposition of various state after the completion of
each algorithm. The probability of observing any given state after observation
is proportional to the absolute value squared of the amplitude of the system
at that state. So, a false solution may in fact be observed which is similar
to Monte-Carlo methods. However, if the algorithm is cleverly constructed, the
correct solution is observed with significant likelihood.

Due to the quantum mechanical nature of quantum computation, new types of
quantum algorithms have appeared. Moreover, these algorithms are more
efficient than existing classical algorithms because the run times are better.
For instance, both integer factorization and discrete logarithms undergo an
exponential speedup using Shor’s algorithm [15]. Not only an exponential
improvement, Grover’s search algorithm provides a quadratic speedup over any
known classical search algorithm [9], and on a discrete space, Grover’s
algorithm is defined by discrete-time quantum walk, which is the natural
extension of a Markov chain driven classical walk to the quantum setting.

On the other hand, if a quantum system were perfectly isolated, it would
maintain coherence indefinitely, but it would be impossible to manipulate or
investigate it. If it is not perfectly isolated, for example during a
measurement, coherence is shared with the environment and appears to be lost
with time which is called quantum decoherence. This concept was first
introduced by H. Zeh [18] in 1970, and then formulated mathematically for
quantum walks by T. Brun [6]. For both coin and position space decoherent
Hadamard walk, K. Zhang proved in [19] that with symmetric initial conditions,
it has Gaussian limiting distribution. More recently, the fact that the
limiting distribution of the rescaled position discrete-time quantum random
walks with general unitary operators subject to only coin space decoherence is
a convex combination of normal distributions under certain conditions is
proved by S. Fan, Z. Feng, S. Xiong and W. Yang [8]. The decoherent quantum
analogues of Markov chains and random walks on a finite space will be defined
and elaborated in this paper.

In fact, classical Markov chain limit theorems for the discrete time walks are
well known and have had important applications in related areas [7] and [13].
However, the primary goal of this paper is to examine the limiting behavior of
the new model, discrete time-inhomogeneous quantum walk with decoherence on
finite spaces and generalize the results from the classical theorems to the
quantum analogues. In this paper, we introduce and study time-inhomogeneous
quantum Markov chains with parameter $\zeta\geq 0$ decoherence parameter
$0\leq p\leq 1$ on finite spaces and their large scale equilibrium properties.
Here $\zeta$ resembles the inverse temperature in the annealing random process
and $p$ is the decoherence strength of the quantum system. Numerical
evaluations show that if $\zeta$ is small, then quantum Markov chain is
ergodic for all $0<p\leq 1$ and if $\zeta$ is large, then it has multiple
limiting distributions for all $0<p\leq 1$. In this paper, we prove the
ergodic property in the high temperature region $0\leq\zeta\leq 1$. We
conjecture that the phase transition is at the critical point $\zeta_{c}=1$.
For coherence case $p=0$, a critical behavior of periodicity also appears at
critical point $\zeta_{o}=2$.

For the time homogeneous case $\zeta=0$ and $o<p\leq 1$, the limiting
distribution has been obtained by Lagro, Yang and Xiong [12] under very
general conditions similar to the Perron-Frobenius type of conditions. Our
paper extends [12] to time-inhomogeneous case, but with stronger assumptions
on the transition matrices.

This paper is organized as follows. In Section 2, we set the notations,
definitions and introduce the model. In Section 3, we develop a path integral
formula for time-inhomogeneous decoherent quantum Markov chains, Theorem 3.2.
In Section 4, we obtain a limiting theorem for classical time-inhomogeneous
Markov chains, Theorem 4.1. In Section 5, we use these theorems to prove our
main result, Theorem 5.2. In Section 6, we provide numerical evaluation to
support our conjectures. In Section 7, we make our conclusions and discuss
some problems for further study.

##  2 Notations and Definitions

In classical probability, a random walk on $\mathbb{Z}$ is a Markov process
described by a stochastic transition matrix $T$. On the other hand, for a
quantum walk, instead of the transition matrix, the evolution of the system is
described by a unitary operator U acting on a Hilbert space H. Several
different models for quantum walks have been popularized. The two most
commonly used are the coined walk of Aharonov et al [1] and the quantum markov
chain of Szegedy [16]. Recently, S. Fan, Z. Feng, S. Xiong and W. Yang et al.
[8] demonstrated that under certain conditions, the limiting distribution of
the rescaled position discrete-time decoherent quantum coined walks is a
convex combination of normal distributions. All quantum walks elaborated here
in this section are based on homogeneous coined Markov chains.

We consider the quantum states $|1\rangle$, …, $|m\rangle$. If $m=2$, then it
represents the head and tail respectively when we flip a coin. Let $H$ be the
Hilbert space spanned by the orthonormal basis, $|1\rangle$, …, $|m\rangle$,

$H=span\\{|1\rangle,...,|m\rangle\\}.$

We denote the inner product of the space by
$\langle|i\rangle,|j\rangle\rangle:=\langle i|j\rangle$, and
$|i\rangle^{*}=\langle i|$ for $x=1,2,...,m$. For n=1, 2, 3, …, let
$U_{n}:H\rightarrow H$ be a unitary operator acting on the Hilbert space $H$
itself. Now, in order to define a decoherent quantum Markov chain over the
Hilbert space $H$, we consider the decoherence parameter $p\in[0,1]$, and
define for $i=1,...,m$

$A_{i}=\sqrt{p}\cdot|i\rangle\langle i|,$

and

$A_{0}=\sqrt{1-p}\cdot I$

Note that $\\{A_{0},A_{1},...,A_{m}\\}$ has the property that
$\sum_{i}{A_{i}}^{*}A_{i}=I$, and is called a measurement over the space $H$.
This notion will be generally defined later.

######  Definition 2.1.

Let $L(H)$ be the set of linear operators on $H$. For $\rho\in L(H)$, we
define the $n$-th step time-inhomogeneous decoherent operator $\Phi_{n}$ such
as

$\Phi_{n}(\rho)=\sum_{i=0}^{m}A_{i}U_{n}\rho U_{n}^{*}A_{i}^{*}$

Now, suppose that $x\in\\{1,...,m\\}$ is the initial position, and
$\rho_{0}=|x\rangle\langle x|$, and the $n$ step time-inhomogeneous quantum
Markov chain is defined as

$\displaystyle\rho_{n}=\underset{n\text{
times}}{\Phi_{n}\cdots\Phi_{1}}(\rho_{0}),$ (2.1)

and, the probability of getting $y\in\\{1,...,m\\}$ from the initial position
$x\in\\{1,...,m\\}$ after $n$ steps is defined as the trace of
$|y\rangle\langle y|\rho_{n}$, we denote it by

$\displaystyle P_{n}(x,y):=Tr\Big{(}|y\rangle\langle y|\rho_{n}\Big{)}.$ (2.2)

######  Example 2.1.

Let $U_{n}=U$ such as,

$\displaystyle U=\left[\begin{array}[]{ c c
}\frac{1}{\sqrt{2}}&\frac{1}{\sqrt{2}}\\\
\frac{1}{\sqrt{2}}&-\frac{1}{\sqrt{2}}\end{array}\right]$ (2.5)

Then the quantum Markov chain is called the time-homogeneous Hadamard (fair
coined) quantum Markov chain.

######  Example 2.2.

Let $\displaystyle U_{n}:H\rightarrow H$, be unitary defined by

$U_{n}=\begin{bmatrix}\sqrt{1-\frac{\lambda}{n^{\zeta}}}&\sqrt{\frac{\lambda}{n^{\zeta}}}\\\
\sqrt{\frac{\lambda}{n^{\zeta}}}&-\sqrt{1-\frac{\lambda}{n^{\zeta}}}\end{bmatrix},$

time inhomogeneous unitary operators, where $\lambda$ and $\zeta$ are non
negative real numbers. Note that we can easily obtain the homogeneous case
from this setting by letting $\zeta=0$. For example, if $\zeta=0$ and
$\lambda=\frac{1}{2}$, we have the fair coined quantum Markov chain from
Example 2.1.

######  Example 2.3.

Let $H=span\\{|1\rangle,...,|m\rangle\\}$ with $m\in\mathbb{N}$, and
$U_{n}=e^{i\frac{G}{\sqrt{n^{\zeta}}}}$ where $G$ is an $m\times m$ self-
adjoint matrix, where $\zeta$ is a non-negative real number.

##  3 Compound Markov chain representation and basic properties

Let’s start by proving some basic properties before we go to the equilibrium
convergence theorem. First, if we suppose that the quantum Markov chain we
defined is completely decoherent, $p=1$, we have

$\Phi_{n}(\rho)=\sum_{i=1}^{m}A_{i}U_{n}\rho U_{n}^{*}A_{i}^{*}$

and

$P_{n}(i,j)=Tr\Big{(}|j\rangle\langle j|\Phi_{n}|i\rangle\langle i|\Big{)}$

where $\rho$ is any $m$ by $m$ density matrix, and $i,j=1,...,m$. Using the
fact that ${|k\rangle,k=1,2,...}$ are orthonormal basis, we obtain

$P_{n}(i,j)=\sum_{k=1}^{m}Tr\Big{(}|j\rangle\langle j|k\rangle\langle
k|U_{n}|i\rangle\langle i|U_{n}^{*}|k\rangle\langle k|\Big{)}$
$=\sum_{k=1}^{m}Tr\Big{(}\langle j|k\rangle\langle k|U_{n}|i\rangle\langle
i|U_{n}^{*}|k\rangle\langle k|j\rangle)\Big{)}$ $=Tr\Big{(}\langle
j|U_{n}|i\rangle\langle
i|U_{n}^{*}|j\rangle\Big{)}=U_{n}(j,i)U_{n}^{*}(i,j)=|U_{n}(j,i)|^{2}$

If $\zeta=0$ and $p=1$, then it reduces to time-homogeneous Markov chain and
it is well known that $\rho_{n}\rightarrow\rho_{\infty}$ when
$n\rightarrow\infty$ by the Ergodic Theorem for finite-state Markov chains.
Our focus now will be for the case of $\zeta>0$ or $0<p<1$.

By Kolomogorov $0$-$1$ law, for $0<p\leq 1$, there exists an infinite sequence
of measurements when $t\rightarrow\infty$. Let $X_{1},X_{2},...$ be the
outcomes of the measurements. In the following Proposition 3.1, we will show
that given the measurement times, $\\{X_{n}\\}_{n=1}^{\infty}$ is a time-
inhomogeneous Markov chain.

######  Definition 3.1.

Let $\\{T_{n}\\}_{n=1}^{\infty}$ be i.i.d. geometric random variables with
mean $\frac{1}{p}$. Let $\sigma_{0}=0$ and $\sigma_{n}=T_{1}+\cdots+T_{n}$.
For $n\geq 1$, we define a Markov transition matrix

$\displaystyle Q_{\sigma_{n-1}}(i,j):=\big{|}\langle
j|U_{\sigma_{n-1}+T_{n}}\cdots U_{\sigma_{n-1}+1}|i\rangle\big{|}^{2}$ (3.6)

######  Proposition 3.1.

Let $\\{T_{n}\\}_{n=1}^{\infty}$ be i.i.d. geometric random variables with
mean $\frac{1}{p}$, and let $\sigma_{0}=0$ and
$\sigma_{n}=T_{1}+\cdots+T_{n}$. Let $X_{1},X_{2},...$ be the outcomes of the
measurements at measurement times $\sigma_{0},\sigma_{1},...$. If
$\rho_{0}=|i\rangle\langle i|$, then

  1. (a)

$\displaystyle
P_{\sigma}\big{(}X_{1}=i_{1},X_{2}=i_{2},...,X_{n}=i_{n}\big{)}=\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots$

$\cdots\big{|}\langle i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$

  2. (b)

$\displaystyle
P_{\sigma}\big{(}X_{n}=i_{n}|X_{1}=i_{1},X_{2}=i_{2},...,X_{n-1}=i_{n-1}\big{)}=P_{\sigma}\big{(}X_{n}=i_{n}|X_{n-1}=i_{n-1}\big{)}$

$=Q_{\sigma_{n-1}}\big{(}i_{n-1},i_{n}\big{)}$

Here $i_{0}=i$ and $P_{\sigma}$ denotes the conditional probability given
$\\{\sigma_{1},\sigma_{2},...\\}$.

###### Proof of Proposition 3.1.

Let $\displaystyle\Phi_{k}(\rho)=\sum_{i=0}^{m}A_{i}U_{k}\rho U_{k}^{*}A_{i}$,
$\rho\in\mathcal{E}$, where $\mathcal{E}=$the set of all density operators on
$H$. The quantum operation of the partially decoherent quantum process at step
$k$ is given by $\Phi_{k}$.

Suppose that the initial state is $\rho_{0}=|i\rangle\langle i|$. Then the
state at time $t$ is

$\rho_{t}=\Phi_{t}\cdots\Phi_{1}(\rho_{0})$

The probability at time $t$, the system is found at state $|j\rangle$ is given
by

$Tr(|j\rangle\langle j|\rho_{t})=\langle
j|\sum_{j_{t}=0}^{m}\cdots\sum_{j_{1}=0}^{m}A_{j_{t}}U_{t}\cdots
A_{j_{2}}U_{2}A_{j_{1}}U_{1}|i\rangle\langle i|U_{1}^{*}A_{j_{1}}^{*}\cdots
U_{t}^{*}A_{j_{t}}^{*}|j\rangle$
$=\sum_{j_{t}=0}^{m}\cdots\sum_{j_{1}=0}^{m}\langle j|A_{j_{t}}U_{t}\cdots
A_{j_{2}}U_{2}A_{j_{1}}U_{1}|i\rangle\langle i|U_{1}^{*}A_{j_{1}}^{*}\cdots
U_{t}^{*}A_{j_{t}}^{*}|j\rangle$ $=\sum_{j_{t},...,j_{1}=0}^{m}\big{|}\langle
j|A_{j_{t}}U_{t}\cdots A_{j_{2}}U_{2}A_{j_{1}}U_{1}|i\rangle\big{|}^{2}$

In each term of above sum, let $0<\sigma_{1}<\sigma_{2}<...<\sigma_{n}\leq t$
be exactly the times that $j_{\sigma_{k}}=1$, $2$, …, or $m$, and $j_{s}=0$
for all $s\neq\sigma_{1},...,\sigma_{n}$.
$0<\sigma_{1}<\sigma_{2}<...<\sigma_{n}$ are called the decoherence times, and
we put $i_{k}=j_{\sigma_{k}}$. So, the sum can be written as,

$Tr(|j\rangle\langle j|\rho_{t})=\sum_{n=0}^{\infty}\text{
}\sum_{0=\sigma_{0}<...<\sigma_{n}\leq
t}p^{n}q^{t-n}\sum_{i_{n}=1}^{m}\cdots\sum_{i_{1}=1}^{m}\big{|}\langle
j|U_{t}\cdots U_{\sigma_{n}+1}|i_{n}\rangle\big{|}^{2}\cdot$
$\displaystyle\cdot\big{|}\langle i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$ (3.7)

Now, we prove Part (a) of the proposition.

$P\big{(}X_{1}=i_{1},...,X_{n}=i_{n}\big{)}=\sum_{t=1}^{\infty}P\big{(}X_{1}=i_{1},...,X_{n}=i_{n},\sigma_{n}=t\big{)}$

By Equation (3.7), we have

$=\sum_{t=1}^{\infty}\text{
}\sum_{0<\sigma_{1}<\cdots<\sigma_{n}=t}p^{n}q^{t-n}\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$
$=\sum_{0<\sigma_{1}<\cdots<\sigma_{n}<\infty}pq^{\sigma_{1}-\sigma_{0}-1}pq^{\sigma_{2}-\sigma_{1}-1}...pq^{\sigma_{n}-\sigma_{n-1}-1}\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}.$

On the other hand,

$P\big{(}X_{1}=i_{1},...,X_{n}=i_{n}\big{)}=E\Big{[}P_{\sigma}(X_{1}=i_{1},...,X_{n}=i_{n})\Big{]},$

and

$P(T_{1}=\sigma_{1},T_{1}+T_{2}=\sigma_{2},...,T_{1}+...+T_{n}=\sigma_{n})=pq^{\sigma_{1}-\sigma_{0}-1}pq^{\sigma_{2}-\sigma_{1}-1}...pq^{\sigma_{n}-\sigma_{n-1}-1}.$

Therefore, we have

$P_{\sigma}(X_{1}=i_{1},...,X_{n}=i_{n})=\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$

To prove (b), by (a) and definition, we have

$P_{\sigma}\big{(}X_{1}=i_{1},X_{2}=i_{2},...,X_{n}=i_{n}\big{)}=Q_{\sigma_{0}}(i,i_{1})\cdots
Q_{\sigma_{n-1}}(i_{n-1},i_{n})$

Then,

$P_{\sigma}\big{(}X_{n}=i_{n}|X_{1}=i_{1},...,X_{n-1}=i_{n-1}\big{)}=\frac{P_{\sigma}\big{(}X_{1}=i_{1},...,X_{n-1}=i_{n-1},X_{n}=i_{n}\big{)}}{P_{\sigma}\big{(}X_{1}=i_{1},...,X_{n-1}=i_{n-1}\big{)}}$
$=Q_{\sigma_{n-1}}(i_{n-1},i_{n})$

We also note that

$P_{\sigma}\big{(}X_{n}=i_{n}|X_{n-1}=i_{n-1}\big{)}=\frac{P_{\sigma}\big{(}X_{n}=i_{n},X_{n-1}=i_{n-1}\big{)}}{P_{\sigma}\big{(}X_{n-1}=i_{n-1}\big{)}}$
$=\frac{\sum_{i_{1},...,i_{n-2}=1}^{m}Q_{\sigma_{0}}(i,i_{1})\cdots
Q_{\sigma_{n-2}}(i_{n-2},i_{n-1})Q_{\sigma_{n-1}}(i_{n-1},i_{n})}{\sum_{i_{1},...,i_{n-2}=1}^{m}Q_{\sigma_{0}}(i,i_{1})\cdots
Q_{\sigma_{n-2}}(i_{n-2},i_{n-1})}=Q_{\sigma_{n-1}}(i_{n-1},i_{n})$

∎

######  Remark 3.1.1.

  1. (a)

It follows from Proposition 3.1 (b) that given measurement times $\sigma$,
$X_{1},X_{2},...$ is a time in-homogeneous Markov chain with transition
probability

$P_{\sigma}\big{(}X_{n}=j|X_{n-1}=i\big{)}=Q_{\sigma_{n-1}}(i,j),\ 1\leq
i,j\leq m,$

and this is illustrated in the following diagram

$0$$1$…$\sigma_{1}$$\uparrow$$X_{1}$…$\sigma_{2}$$\uparrow$$X_{2}$…$\sigma_{n-1}$$\uparrow$$X_{n-1}$$\sigma_{n}$$\uparrow$$X_{n}$…$t$$\sigma_{n+1}$$\uparrow$$X_{n+1}$time$T_{1}$$T_{2}$$T_{n}$$T_{n+1}$

  2. (b)

This proposition gives a probabilistic interpretation of decoherence time
$0<\sigma_{0}<\sigma_{1}<...<\sigma_{n}<...$ as a sequence of arrival times of
independence Bernoulli trials with success probability $p$, and
$\\{X_{n}\\}_{n=1}^{\infty}$ can be viewed as discrete version of a compound
Poisson process.

For $\sigma_{n}\leq t$, and from Equation (3.7), we obtain an expression for
the probability that the system is found at state $j$ with exact decoherent
time $\sigma_{1},...\sigma_{n}$ and outcomes $i_{1},i_{2},...,i_{n}$. We will
formualte this as our main representation theorem as follows.

######  Definition 3.2.

Let $t\geq 0$ and $n_{t}=\max\\{n\geq 0;\sigma_{n}\leq t\\}$ be the number of
occurences before or at time $t$. Let

$W_{\sigma_{n_{t}}}(i,j)=\big{|}\langle j|U_{t}\cdots
U_{\sigma_{n_{t}-1}+1}|i\rangle\big{|}^{2}$

From Equation (3.7) and Definition 3.2, we have the following

######  Theorem 3.2 (Compound Markov Chain Representation).

Suppose $\rho_{0}=|i\rangle\langle i|$. Then

$Tr(|j\rangle\langle
j|\rho_{t})=E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}(i,j)]$

Note that in the above theorem, the matrices $Q$ and $W$ are pure quantum
random walks. In between decohenernces times, the successive processes are
coherent quantum walks. The Compound Markov Chain Representation Theorem 3.2
not only gives a very intuitive expression for the $0<p\leq 1$ decoherent
quantum walks with the probability of a specific path from the initial state
to the target state, but also provides a very useful tool for generating
numerical simulations, as we shall see in Section 6.

###### Proof of Theorem 3.2.

The right hand side of the theorem,

$E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}(i,j)]$
$=\sum_{n=0}^{\infty}\text{ }\sum_{0=\sigma_{0}<...<\sigma_{n}\leq
t<\sigma_{n+1}}P(T_{1}=\sigma_{1},T_{1}+T_{2}=\sigma_{2},...,T_{1}+...+T_{n}=\sigma_{n},T_{1}+...+T_{n+1}=\sigma_{n+1})$
$\sum_{i_{n}=1}^{m}\cdots\sum_{i_{1}=1}^{m}\big{|}\langle j|U_{t}\cdots
U_{\sigma_{n}+1}|i_{n}\rangle\big{|}^{2}\cdot\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$
$=\sum_{n=0}^{\infty}\text{ }\sum_{0=\sigma_{0}<...<\sigma_{n}\leq
t<\sigma_{n+1}}p^{n}q^{t-n}q^{\sigma_{n+1}-t-1}p\sum_{i_{n}=1}^{m}\cdots\sum_{i_{1}=1}^{m}\big{|}\langle
j|U_{t}\cdots U_{\sigma_{n}+1}|i_{n}\rangle\big{|}^{2}\cdot$ $\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$
$=\sum_{n=0}^{\infty}\text{ }\sum_{0=\sigma_{0}<...<\sigma_{n}\leq
t}p^{n}q^{t-n}\sum_{i_{n}=1}^{m}\cdots\sum_{i_{1}=1}^{m}\big{|}\langle
j|U_{t}\cdots U_{\sigma_{n}+1}|i_{n}\rangle\big{|}^{2}\cdot$ $\big{|}\langle
i_{n}|U_{\sigma_{n}}\cdots
U_{\sigma_{n-1}+1}|i_{n-1}\rangle\big{|}^{2}\cdots\big{|}\langle
i_{1}|U_{\sigma_{1}}\cdots U_{1}|i\rangle\big{|}^{2}$ $=Tr(|j\rangle\langle
j|\rho_{t}),$

by Equation (3.7).

∎

##  4 Convergence to equilibrium

We first show a limiting theorem for time-inhomogeneous Markov chains in
classical probability as follows.

######  Theorem 4.1.

Let $\Pi_{ij}=\Pi_{1j}\geq 0$ for all $i,j$ and $\sum_{j}\Pi_{ij}=1$. Let
$P_{k}$ be the Markov transition matrices on a finite states space
$\Sigma=\\{1,2,...,m\\}$. Suppose there is a $k_{0}$ such that for all $k\geq
k_{0}$, $\Pi P_{k}=\Pi$ and $P_{k}\geq\delta_{k}\Pi$, where
$0\leq\delta_{k}\leq 1$. If
$\displaystyle\prod_{k=k_{0}}^{\infty}(1-\delta_{k})=0$, then $P_{s}\cdots
P_{n}\rightarrow\Pi$ as $n\rightarrow\infty$, for all $1\leq s\leq k_{0}$.

###### Proof of Theorem 4.1.

Note that since $P_{k}$ are Markov transition matrices, if
$P_{k}\geq\delta_{k}\Pi$, for some $\delta_{k}>0$, then $\delta_{k}\leq 1$.
Moreover, if $\delta_{k}=1$, then $P_{k}=\Pi$. So if there is a $k\geq k_{0}$
such that $\delta_{k}=1$, then $P_{1}\cdots P_{n}=\Pi$, since $P_{k}\Pi=\Pi$
and by assumption $\Pi P_{k}=\Pi$. Therefore it is sufficient to consider the
case $0\leq\delta_{k}<1$.

We also note that since $P_{k}$ are Markov transition matrices, if

$P_{k_{0}}P_{k_{0}+1}\cdots P_{n}\rightarrow\Pi,$

then

$P_{s}\cdots P_{n}\rightarrow\Pi,$

as $n\rightarrow\infty$, for all $1\leq s\leq k_{0}$.

For $n\geq k_{0}$, we have

$P_{n}=\alpha_{n}\Big{(}\frac{P_{n}-\delta_{n}\Pi}{\alpha_{n}}\Big{)}+\delta_{n}\Pi=\alpha_{n}\tilde{P}_{n}+\delta_{n}\Pi$

Since $\alpha_{n}=1-\delta_{n}$ and $\tilde{P}_{n}$ is a Markov transition
matrix, and by assumption, we have the following

$\displaystyle\alpha_{n}+\delta_{n}=1$ (4.8) $\displaystyle\Pi\cdot\Pi=\Pi$
(4.9)
$\displaystyle\displaystyle\tilde{P_{n}}\cdot\Pi=\frac{P_{n}-\delta_{n}\Pi}{\alpha_{n}}\Pi=\frac{P_{n}\Pi-\delta_{n}\Pi^{2}}{\alpha_{n}}=\frac{\Pi-\delta_{n}\Pi}{\alpha_{n}}=\frac{1-\delta_{n}}{\alpha_{n}}\Pi=\Pi$
(4.10) $\displaystyle\Pi\cdot\tilde{P_{n}}=\Pi$ (4.11)

Therefore, we have

$P_{k_{0}}\cdots
P_{n}=\prod_{k=k_{0}}^{n}\big{[}\alpha_{n}\tilde{P}_{n}+\delta_{n}\Pi\big{]}=\prod_{k=k_{0}}^{n}\alpha_{k}\tilde{P}_{k}+\big{(}1-\prod_{k=k_{0}}^{n}\alpha_{k}\big{)}\Pi$

The last equality can be proved by induction, if $n=k_{0}+1$

$\big{[}\alpha_{k_{0}}\tilde{P}_{k_{0}}+(1-\alpha_{k_{0}})\Pi\big{]}\cdot\big{[}\alpha_{k_{0}+1}\tilde{P}_{k_{0}+1}+(1-\alpha_{k_{0}+1})\Pi\big{]}$
$=\alpha_{k_{0}}\alpha_{k_{0}+1}\tilde{P}_{k_{0}}\tilde{P}_{k_{0}+1}+(1-\alpha_{k_{0}})\alpha_{k_{0}+1}\Pi\tilde{P}_{k_{0}+1}+\alpha_{k_{0}}(1-\alpha_{k_{0}+1})\tilde{P}_{k_{0}}\Pi+(1-\alpha_{k_{0}})(1-\alpha_{k_{0}+1})\Pi^{2}$
$=\alpha_{k_{0}}\alpha_{k_{0}+1}\tilde{P}_{k_{0}}\tilde{P}_{k_{0}+1}+\Pi-\alpha_{k_{0}}\alpha_{k_{0}+1}\Pi+\alpha_{k_{0}}\Pi-\alpha_{k_{0}}\alpha_{k_{0}+1}\Pi-\Pi+\alpha_{k_{0}}\alpha_{k_{0}+1}\Pi+\Pi-\alpha_{k_{0}}\Pi=$
$\alpha_{k_{0}}\alpha_{k_{0}+1}\tilde{P}_{k_{0}}\tilde{P}_{k_{0}+1}+(1-\alpha_{k_{0}}\alpha_{k_{0}+1})\Pi$

And now, let’s assume that the formula is true for $n$, and prove it for
$n+1$,

$\prod_{k=k_{0}}^{n+1}\big{[}\alpha_{n}\tilde{P}_{n}+\delta_{n}\Pi\big{]}=\Big{[}\prod_{k=k_{0}}^{n}\alpha_{k}\tilde{P}_{k}+\big{(}1-\prod_{k=k_{0}}^{n}\alpha_{k}\big{)}\Pi\Big{]}\Big{[}\alpha_{n+1}\tilde{P}_{n+1}+(1-\alpha_{n+1})\Pi\Big{]}$
$=\prod_{k=k_{0}}^{n+1}\alpha_{k}\tilde{P}_{k}+\prod_{k=k_{0}}^{n}\alpha_{k}(1-\alpha_{n+1})\Pi+\big{(}1-\prod_{k=k_{0}}^{n}\alpha_{k}\big{)}\alpha_{n+1}\Pi+\big{(}1-\prod_{k=k_{0}}^{n}\alpha_{k}\big{)}(1-\alpha_{n+1})\Pi$
$=\prod_{k=k_{0}}^{n+1}\alpha_{k}\tilde{P}_{k}+\prod_{k=k_{0}}^{n}\alpha_{k}\Pi-\prod_{k=1}^{n+1}\alpha_{k}\Pi+\alpha_{n+1}\Pi-\prod_{k=k_{0}}^{n+1}\alpha_{k}\Pi+\Pi-\alpha_{n+1}\Pi-\prod_{k=k_{0}}^{n}\alpha_{k}\Pi+\prod_{k=1}^{n+1}\alpha_{k}\Pi$
$=\prod_{k=k_{0}}^{n+1}\alpha_{k}\tilde{P}_{k}+\big{(}1-\prod_{k=k_{0}}^{n+1}\alpha_{k}\big{)}\Pi$

For an $m\times m$ matrix $B$, we denote by $||B||_{\infty}$ the operator norm
from $(\mathbb{C}^{m},||\cdot||_{\infty})$ with sup norm to itself. Note that
if $P$ is a Markov transition matrix, then $||P||_{\infty}=1$.

Then we have

$||\prod_{k=k_{0}}^{n}\alpha_{n}\tilde{P}_{n}||_{\infty}\leq\prod_{k=k_{0}}^{n}\alpha_{n}$

We observe that if $\displaystyle\prod_{k=k_{0}}^{\infty}\alpha_{k}=0$, then
$\displaystyle P_{k_{0}}P_{k_{0}+1}\cdots P_{n}\rightarrow\Pi$, and we have
obtained the conclusion. ∎

We will need the following estimates for a unitary semigroup.

######  Lemma 4.2.

Let $G$ be an $m\times m$ self-adjoint matrix, and $U=e^{i\theta G}$ where
$\theta$ is a non-negative real number. Suppose there exists a constant
$\epsilon_{0}>0$ such that $|G_{jk}|\geq\epsilon_{0}$, for all $j,k$. If
$\theta\leq\frac{\epsilon_{0}}{4||G||_{\infty}^{2}}\wedge\frac{1}{4||G||_{\infty}}$,
then

(a) for all $j\neq k$, we have

$\frac{1}{2}\theta\epsilon_{0}\leq|U_{jk}|\leq 2\theta||G||_{\infty},$

(b) for all $j$, we have

$\frac{1}{2}\leq|U_{jj}|.$

###### Proof of Lemma 4.2.

(a) Let $j\neq k$. For the first inequality, we have

$|U_{jk}|=|(e^{i\theta G})_{jk}|=|[e^{i\theta G}-I]_{jk}|=|i\theta
G_{jk}+\sum_{n=2}^{\infty}\frac{(i\theta G)^{n}_{jk}}{n!}|$ $\geq|i\theta
G_{jk}|-|\sum_{n=2}^{\infty}\frac{(i\theta G)^{n}_{jk}}{n!}|$
$\displaystyle\geq\theta\epsilon_{0}-\sum_{n=2}^{\infty}\frac{\theta^{n}||G||_{\infty}^{n}}{n!}$
(4.12)

The second term of the above is bounded by

$\displaystyle\sum_{n=2}^{\infty}\frac{\theta^{n}||G||_{\infty}^{n}}{n!}\leq\sum_{n=2}^{\infty}\theta^{n}||G||_{\infty}^{n}=\frac{\theta^{2}||G||_{\infty}^{2}}{1-\theta||G||_{\infty}}\leq
2\theta^{2}||G||_{\infty}^{2},$ (4.13)

since $\theta\leq\frac{1}{4||G||_{\infty}}$. Therefore, we have

$\displaystyle|U_{jk}|\geq\theta\epsilon_{0}-2\theta^{2}||G||_{\infty}^{2}\geq\frac{1}{2}\theta\epsilon_{0},$
(4.14)

since $\theta\leq\frac{\epsilon_{0}}{4||G||_{\infty}^{2}}$.

For the second inequality in (a), we have

$|U_{jk}|=|(e^{i\theta G})_{jk}|=|[e^{i\theta
G}-I]_{jk}|=|\sum_{n=1}^{\infty}\frac{(i\theta G)^{n}_{jk}}{n!}|$
$\displaystyle\leq\sum_{n=1}^{\infty}\frac{\theta^{n}||G||_{\infty}^{n}}{n!}$
$\displaystyle\leq\frac{\theta||G||_{\infty}}{1-\theta||G||_{\infty}}\leq
2\theta||G||_{\infty},$ (4.15)

since $\theta\leq\frac{1}{4||G||_{\infty}}$.

(b) For any $j$, we have

$|U_{jj}|=|(e^{i\theta G})_{jj}|=|1+\sum_{n=1}^{\infty}\frac{(i\theta
G)^{n}_{jj}}{n!}|$ $\geq 1-|\sum_{n=1}^{\infty}\frac{(i\theta
G)^{n}_{jk}}{n!}|$ $\displaystyle\geq
1-\sum_{n=1}^{\infty}\frac{\theta^{n}||G||_{\infty}^{n}}{n!}$ (4.16)

The second term of the above is bounded by

$\displaystyle\sum_{n=1}^{\infty}\frac{\theta^{n}||G||_{\infty}^{n}}{n!}\leq\sum_{n=1}^{\infty}\theta^{n}||G||_{\infty}^{n}=\frac{\theta||G||_{\infty}}{1-\theta||G||_{\infty}}\leq
2\theta||G||_{\infty},$ (4.17)

since $\theta\leq\frac{1}{4||G||_{\infty}}$. Therefore, we have

$\displaystyle|U_{jj}|\geq 1-2\theta||G||_{\infty}\geq\frac{1}{2},$ (4.18)

since $\theta\leq\frac{1}{4||G||_{\infty}}$.

∎

Before we prove the convergence, let’s first observe that $Q_{\sigma_{n}}$
defined in Definition 3.1 is doubly stochastic which will be useful later in
the proofs.

######  Proposition 4.3.

$Q_{\sigma_{n-1}}$ is doubly stochastic for all $n\geq 1$.

###### Proof of Proposition 4.3.

Note that a matrix $A=(a_{ij})$, $a_{ij}\geq 0$, is doubly stochastic if

$\sum_{i}a_{ij}=1\ \text{ for all }j$

and

$\sum_{j}a_{ij}=1\ \text{ for all }i$

So, we have for all $j=1,2,...,m$,

$\sum_{i=1}^{m}Q_{\sigma_{n-1}}(i,j)=\sum_{i=1}^{m}\langle
j|U_{\sigma_{n-1}+T_{n}}\cdots U_{\sigma_{n-1}+1}|i\rangle\langle
i|U_{\sigma_{n-1}+1}^{*}\cdots U_{\sigma_{n-1}+T_{n}}^{*}|j\rangle$
$=\sum_{i=1}^{m}Tr\big{[}|j\rangle\langle j|U_{\sigma_{n-1}+T_{n}}\cdots
U_{\sigma_{n-1}+1}|i\rangle\langle i|U_{\sigma_{n-1}+1}^{*}\cdots
U_{\sigma_{n-1}+T_{n}}^{*}\big{]}$ $=Tr\big{(}|j\rangle\langle j|\big{)}=1$

And also, $\displaystyle\sum_{j=1}^{m}Q_{\sigma_{n-1}}(i,j)=1$ for all
$i=1,2,...,m$ by similar argument. ∎

Finally, we have the equilibrium property

######  Proposition 4.4.

Let $H=span\\{|1\rangle,...,|m\rangle\\}$ with $m\in\mathbb{N}$, and
$U_{n}=e^{i\frac{G}{\sqrt{n^{\zeta}}}}$ where $G$ is an $m\times m$ self-
adjoint matrix, where $\zeta$ is a non-negative real number. Suppose there
exists a constant $\epsilon_{0}>0$ such that $|G_{ij}|\geq\epsilon_{0}$, for
all $i,j$. Let $\Pi$ be an $m\times m$ matrix with $\Pi_{ij}=\frac{1}{m}$, for
all $i,j$. If $0<\zeta\leq 1$ and $0<p\leq 1$, then for all $0\leq s<\infty$,

$Q_{\sigma_{s}}Q_{\sigma_{s+1}}Q_{\sigma_{s+2}}\cdots
Q_{\sigma_{n-1}}\rightarrow\Pi,$

almost surely, as $n\rightarrow\infty$.

###### Proof of Proposition 4.4 .

Step 1. We first show that with probability 1,
$Q_{\sigma_{n-1}}\geq\delta_{n-1}\Pi$, where
$\delta_{n-1}=\frac{m\epsilon_{0}^{2}T_{n}^{2}}{4(\sigma_{n})^{\zeta}}$, for
sufficeintly large $n$.

Let

$\theta_{0}=\frac{\epsilon_{0}}{4||G||_{\infty}^{2}}\wedge\frac{1}{4||G||_{\infty}}$

Let

$\theta=\frac{1}{(\sigma_{n-1}+1)^{\zeta/2}}+...+\frac{1}{(\sigma_{n-1}+T_{n})^{\zeta/2}}.$

Then we have

$\displaystyle\frac{T_{n}}{\sigma_{n}^{\zeta/2}}\leq\theta\leq\frac{T_{n}}{\sigma_{n-1}^{\zeta/2}}$
(4.19)

This implies that with probability one, $\theta\to 0$ as $n\to\infty$.
Therefore, with probability 1, there is an $n_{0}>s$ such that
$\theta<\theta_{0}$, for all $n\geq n_{0}>s$.

By definition,

$\displaystyle Q_{\sigma_{n-1}}(i,j)=\langle j|U_{\sigma_{n-1}+T_{n}}\cdots
U_{\sigma_{n-1}+1}|i\rangle\langle i|U_{\sigma_{n-1}+1}^{*}\cdots
U_{\sigma_{n-1}+T_{n}}^{*}|j\rangle$ (4.20) $\displaystyle=|\langle
j|U_{\sigma_{n-1}+T_{n}}\cdots U_{\sigma_{n-1}+1}|i\rangle|^{2}=|e^{i\theta
G}(i,j)|^{2}$ (4.21) $\displaystyle\geq\frac{\epsilon_{0}^{2}\theta^{2}}{4},$
(4.22)

for $i\neq j$, by Lemma 4.2, if $\theta\leq\theta_{0}$. For $i=j$, by Lemma
4.2, we have a lower bound $1/4$. Combining (4.19), we have
$Q_{\sigma_{n-1}}\geq\delta_{n-1}\Pi$, where
$\delta_{n-1}=\frac{m\epsilon_{0}^{2}T_{n}^{2}}{4(\sigma_{n})^{\zeta}}$, for
some $n_{0}>s$ with $n\geq n_{0}$.

Now

$\prod_{k=n_{0}}^{n}\alpha_{k}=\prod_{k=1}^{n}(1-\delta_{k})=\prod_{k=n_{0}}^{n}\Big{[}1-\frac{m\epsilon_{0}^{2}T_{k}^{2}}{4(\sigma_{k})^{\zeta}}\Big{]}$
$=e^{\sum_{k=n_{0}}^{n}\ln\Big{[}1-\frac{m\epsilon_{0}^{2}T_{k}^{2}}{4(\sigma_{k})^{\zeta}}\Big{]}}\leq
e^{-\sum_{k=n_{0}}^{n}\frac{m\epsilon_{0}^{2}T_{n}^{2}}{4(\sigma_{n})^{\zeta}}}$

since $\ln(1-x)\leq-x$ for $0<x<1$.

By the Strong Law of Large Numbers, with probability one,

$\frac{k^{\zeta}}{\sigma_{k}^{\zeta}}=\frac{1}{\big{(}\frac{T_{1}+T_{2}+\cdots+T_{k}+1}{k}\big{)}^{\zeta}}\rightarrow\frac{1}{(E(T_{1}))^{\zeta}}={p^{\zeta}}$

as $k\rightarrow\infty$. Therefore,

$\prod_{k=n_{0}}^{n}\alpha_{k}\leq
e^{-\sum_{k=n_{0}}^{n}\frac{m\epsilon_{0}^{2}}{4\sigma_{n}^{\zeta}}}\leq\exp\Big{(}-\frac{m\epsilon_{0}^{2}}{4}\sum_{k=n_{0}}^{n}\frac{p^{\zeta}}{k^{\zeta}}\frac{\frac{1}{\sigma_{k}^{\zeta}}}{\frac{p^{\zeta}}{k^{\zeta}}}\Big{)}$

Since
$\displaystyle\frac{\frac{1}{\sigma_{k}^{\zeta}}}{\frac{p^{\zeta}}{k^{\zeta}}}\rightarrow
1$, almost surely, as $k\rightarrow\infty$, we have that
$\prod_{k=n_{0}}^{n}\alpha_{k}\rightarrow 0$ a.s., if $\zeta\leq 1$. And, we
conclude that, $Q_{\sigma_{s}}\cdots Q_{\sigma_{n-1}}\rightarrow\Pi$, a.s.,
for all $0<\zeta\leq 1$ and $0<p\leq 1$.

∎

Note that for the special homogeneous case $\zeta=0$, the result is already
known, and it was proved by Lagro et al. [12] that the quantum Markov chain is
convergent. [12] did not use the compound Markov chain representation. It used
the spectral theory of the density operators. So it does not contain this type
of results for $\zeta=0$

##  5 Convergence of density operators

In this section we will prove convergence of density operators for $\zeta\leq
1$ and $0<p\leq 1$. The homogeneous case $\zeta=0$ has been done in Lagro,
Yang and Xiong [12] using spectral analysis of density operators. So we will
just focus on the case when $0<\zeta\leq 1$. Let us recall a few definitions
that will be needed for this section.

We consider $H=span\\{|1\rangle,...,|m\rangle\\}$ with $m\in\mathbb{N}$, and
$U_{n}=e^{i\frac{G}{\sqrt{n^{\zeta}}}}$ where $G$ is an $m\times m$ self-
adjoint matrix such that there exists $\epsilon_{0}>0$,
$|G_{ij}|>\epsilon_{0}$ for all $i,j$. Let

$\Phi_{n}(\rho)=\sum_{i=0}^{m}A_{i}U_{n}\rho U_{n}^{*}A_{i}^{*}$

and

$\displaystyle\rho_{n}=\Phi_{n}\cdots\Phi_{1}(\rho)$ (5.23)

where $A_{0}=\sqrt{1-p}I$ and $A_{i}=\sqrt{p}|i\rangle\langle i|$, $1\leq
i\leq m$, and let $\rho_{0}=|i\rangle\langle i|$ fixed, and

$P_{n}(i,j)=Tr\big{(}|j\rangle\langle j|\rho_{n}\big{)}$

Let $T_{1},T_{2},...$ be i.i.d. with distribution Geo($p$), and
$\sigma_{n}=T_{1}+\cdots+T_{n}$. For a fixed $t$, we let
$n_{t}=\max\\{n:\sigma_{n}\leq t\\}$. By Theorem 3.2, we have

$\displaystyle P_{t}(i,j)=Tr\big{(}|j\rangle\langle
j|\rho_{t}\big{)}=E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}(i,j)]$
(5.24)

where $Q_{\sigma_{n-1}}(i,j)=|\langle j|U_{\sigma_{n-1}+T_{n}}\cdots
U_{\sigma_{n-1}+1}|i\rangle|^{2}$ and $W_{\sigma_{n_{t}}}(i,j)=|\langle
j|U_{t}\cdots U_{\sigma_{n}+1}|i\rangle|^{2}$.

We first prove convergence of probabilities as follows.

######  Proposition 5.1.

Suppose there exists $\epsilon_{0}>0$ such that $|G_{ij}|>\epsilon_{0}$ for
all $i,j$. For all probability distribution $V$ and $j=1,2,...,m$, if
$0<\zeta\leq 1$ and $0<p\leq 1$, then

$\sum_{i}V_{i}P_{t}(i,j)\rightarrow\pi_{j}=\frac{1}{m},$

as $t\rightarrow\infty$.

Using the same argument as in the proof, under the same conditions of
Proposition 5.1, we have

######  Corollary 5.1.1.

Suppose there exists $\epsilon_{0}>0$ such that $|G_{ij}|>\epsilon_{0}$ for
all $i,j$. If $0<\zeta\leq 1$ and $0<p\leq 1$, then for all $0\leq s<\infty$,
we have

$Q_{\sigma_{s}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}\to\Pi,$

almost surely, as $t\to\infty$.

###### Proof of Proposition 5.1.

Let $\Pi(i,j)=\pi_{j}=1/m$. Observe that we can write $P_{t}(i,j)-\Pi(i,j)=$

$E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}-\Pi]=E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}W_{\sigma_{n_{t}}}-Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}+Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}-\Pi]$
$=E[Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}(W_{\sigma_{n_{t}}}-I\big{)}]+E\big{[}Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}-\Pi\big{]}$

Since $0<\zeta\leq 1$, the second term tends to $0$ as $t\rightarrow\infty$ by
Proposition 4.4. Since $Q_{\sigma_{0}}Q_{\sigma_{1}}...Q_{\sigma_{n_{t}-1}}$
and $W_{\sigma_{n_{t}}}-I$ are bounded, by the Bounded Convergence Theorem, it
is sufficient to show that

$\displaystyle\lim_{t\to\infty}W_{\sigma_{n_{t}}}=I,a.s.$ (5.25)

$W_{\sigma_{n_{t}}}(i,j)=|\langle j|U_{t}\cdots
U_{\sigma_{n_{t}}+1}|i\rangle|^{2}=\big{|}\langle
j|e^{\frac{iG}{\sqrt{t^{\zeta}}}}\cdots
e^{\frac{iG}{\sqrt{(\sigma_{n_{t}}+1)^{\zeta}}}}|i\rangle\big{|}^{2}$
$=\big{|}\langle
j|e^{iG\sum_{k=\sigma_{n_{t}}+1}^{t}\frac{1}{\sqrt{k^{\zeta}}}}|i\rangle\big{|}^{2}$

Note that if we show that $|\langle
j|e^{iG\sum_{k=\sigma_{n_{t}}+1}^{t}\frac{1}{\sqrt{k^{\zeta}}}}|i\rangle\big{|}^{2}\rightarrow
0$ a.s. for $i\neq j$, then $W_{\sigma_{n_{t}}}(i,j)\rightarrow\delta_{ij}$,
a.s. and we will obtain the result.

Let’s write $U_{n}=U_{n}^{F}+U_{n}^{D}$ for every $n$ where $U_{n}^{F}$ is the
off-diagonal part and $U_{n}^{D}$ is the diagonal parts of $U_{n}$. So, for
$i\neq j$,

$|\langle
j|e^{iG\sum_{k=\sigma_{n_{t}}+1}^{t}\frac{1}{\sqrt{k^{\zeta}}}}|i\rangle\big{|}^{2}=\Big{|}\sum_{l=\sigma_{N_{t}}+1}^{t}\langle
j|U_{t}\cdots U_{l+1}U_{l}^{F}U_{l-1}^{D}\cdots
U_{\sigma_{N_{t}}+1}^{D}|i\rangle\Big{|}^{2}$

The above expansion is made by applying $U\cdots U(U^{F}+U^{D})=U\cdots
UU^{F}+U\cdots U(U^{F}+U^{D})U^{D}=\cdots$, successively. By the Cauchy-
Schwarz Inequality,

$\leq(t-\sigma_{N_{t}})\sum_{l=\sigma_{N_{t}}+1}^{t}\big{|}\langle
j|U_{t}\cdots U_{l+1}U_{l}^{F}U_{l-1}^{D}\cdots
U_{\sigma_{N_{t}}+1}^{D}|i\rangle\Big{|}^{2}$

By definition, we have $\displaystyle U_{n}=e^{\frac{iG}{\sqrt{n^{\zeta}}}}$,
where $G$ is self-adjoint $m$ by $m$ matrix and $|U_{n}(i,j)|\leq 1$ for all
$i,j$, the above is

$\leq(t-\sigma_{N_{t}})\sum_{l=\sigma_{N_{t}}+1}^{t}\big{|}\langle
j|U_{t}\cdots U_{l+1}U_{l}^{F}|i\rangle U_{l-1}(i,i)\cdots
U_{\sigma_{N_{t}}+1}(i,i)\big{|}^{2}$
$\leq(t-\sigma_{N_{t}})\sum_{l=\sigma_{N_{t}}+1}^{t}\big{|}\sum_{k\neq
i}\langle j|U_{t}\cdots U_{l+1}|k\rangle U_{l}(k,i)U_{l-1}(i,i)\cdots
U_{\sigma_{N_{t}}+1}(i,i)\big{|}^{2}$

By the Cauchy-Schwarz Inequality again,

$\leq(t-\sigma_{N_{t}})\sum_{l=\sigma_{N_{t}}+1}^{t}(m-1)\sum_{k\neq
i}\big{|}\langle j|U_{t}\cdots
U_{l+1}|k\rangle\big{|}^{2}\big{|}U_{l}(k,i)\big{|}^{2}\big{|}U_{l-1}(i,i)\cdots
U_{\sigma_{N_{t}}+1}(i,i)\big{|}^{2}$

Using the fact that $\big{|}\langle j|U_{t}\cdots
U_{l+1}|k\rangle\big{|}^{2}\leq 1$, and by Lemma 4.2,
$\big{|}U_{l}(k,i)\big{|}^{2}\leq 4||G||_{\infty}^{2}(\frac{1}{l^{\zeta}})$,
for sufficiently large $l$, and

$\displaystyle\big{|}U_{l-1}(i,i)\cdots
U_{\sigma_{N_{t}}+1}(i,i)\big{|}^{2}\leq 1,$

we conclude

$(t-\sigma_{N_{t}})\sum_{l=\sigma_{N_{t}}+1}^{t}(m-1)^{2}4||G||_{\infty}^{2}\big{(}\frac{1}{l^{\zeta}}\big{)}\leq(m-1)^{2}(t-\sigma_{N_{t}})^{2}4||G||_{\infty}^{2}\Big{(}\frac{1}{{(\sigma_{N_{t}}+1)}^{\zeta}}\Big{)}$

Now, it’s enough to prove that
$\displaystyle\frac{(t-\sigma_{N_{t}})^{2}}{(\sigma_{N_{t}})^{\zeta}}\leq\frac{T_{N_{t}+1}^{2}}{(T_{1}+\cdots+T_{N_{t}})^{\zeta}}\rightarrow
0$ a.s. The first inequality is trivial, and note that

$\sum_{n=1}^{\infty}P\Big{(}\Big{|}\frac{T_{n}^{2}}{n^{\zeta}}\Big{|}>\epsilon\Big{)}=\sum_{n=1}^{\infty}P\Big{(}T_{n}^{2}>n^{\zeta}\epsilon\Big{)}=\sum_{n=1}^{\infty}P\Big{(}\frac{T_{1}^{\frac{2}{\zeta}}}{\epsilon^{\frac{1}{\zeta}}}>n\Big{)}\leq
E\big{(}\frac{T_{1}^{\frac{2}{\zeta}}}{\epsilon^{\frac{1}{\zeta}}}\big{)}+1<\infty.$

By the Borel-Cantelli Lemma,

$P\Big{(}\Big{|}\frac{T_{n}^{2}}{n^{\zeta}}\Big{|}>\epsilon\text{ i.o.
}\Big{)}=0$

Therefore

$\frac{T_{n}^{2}}{n^{\zeta}}\rightarrow 0\text{ a.s.}$

So, using the Strong Law of Large Numbers, we can conclude that

$\frac{T_{n+1}^{2}}{(T_{1}+\cdots
T_{n})^{\zeta}}=\frac{\frac{T_{n+1}^{2}}{n^{\zeta}}}{(\frac{T_{1}+\cdots
T_{n}}{n})^{\zeta}}\rightarrow\frac{0}{E(T_{1})^{\zeta}}\text{ a.s. }$

∎

Now we show our main result for the convergence of density operators. For
$\zeta=0$, the following result has been obtained in [12] with more general
conditions. Here we prove for the case $0<\zeta\leq 1$.

######  Theorem 5.2.

Suppose that $|G_{ij}|>\epsilon_{0}>0$ for all $i,j$. Suppose $0<\zeta\leq 1$
and $0<p\leq 1$. Then for any initial density matrix $\rho_{0}$, we have

$\displaystyle\rho_{t}\rightarrow\sum_{i=1}^{m}\pi_{i}\cdot|i\rangle\langle
i|$ (5.26)

where $\pi_{i}=\frac{1}{m}$, for $i=1,...,m$.

###### Proof of Theorem 5.2.

Let $0<\zeta\leq 1$ and $0<p\leq 1$.

Step 1. We first consider the case $\rho_{0}=|i\rangle\langle i|$. By
Proposition 5.1, we have

$\displaystyle\lim_{t\to\infty}Tr\big{(}|j\rangle\langle
j|\rho_{t}\big{)}=\lim_{t\to\infty}P_{t}(i,j)=\pi_{j}.$ (5.27)

To prove (5.26), it is then sufficient to show

$\displaystyle\lim_{t\to\infty}Tr\big{(}|j\rangle\langle k|\rho_{t}\big{)}=0,$
(5.28)

for all $k\neq j$.

$Tr\big{(}|k\rangle\langle j|\rho_{t}\big{)}=\langle
j|\rho_{t}|k\rangle=\langle j|(\Phi_{t}\cdots\Phi_{2}\Phi_{1}|i\rangle\langle
i|)|k\rangle$ $=\sum_{i_{1},i_{2},...,i_{t}=0}^{m}\langle
j|A_{i_{t}}U_{t}\cdots A_{i_{1}}U_{1}|i\rangle\langle
i|U_{1}^{*}A_{i_{1}}^{*}\cdots U_{t}^{*}A_{i_{t}}^{*}|k\rangle$
$=E\Big{[}\sum_{j_{1},...,j_{N_{t}}=1}^{m}\langle j|U_{t}\cdots
U_{\sigma_{N_{t}}+1}|j_{N_{t}}\rangle\langle
j_{N_{t}}|U_{\sigma_{N_{t}}}\cdots
U_{\sigma_{N_{t}-1}+1}|j_{N_{t}-1}\rangle\cdots$ $\cdots\langle
j_{1}|U_{\sigma_{1}}\cdots U_{2}U_{1}|i\rangle\langle
i|U_{1}^{*}U_{2}^{*}\cdots U_{\sigma_{1}}^{*}|j_{1}\rangle\cdots\langle
j_{N_{t}-1}|U_{\sigma_{N_{t}-1}+1}^{*}\cdots
U_{\sigma_{N_{t}}}^{*}|j_{N_{t}}\rangle\times$ $\times\langle
j_{N_{t}}|U_{\sigma_{N_{t}}+1}^{*}\cdots U_{t}^{*}|k\rangle\Big{]}$
$=E\Big{[}\sum_{j_{N_{t}}=1}^{m}\langle j|U_{t}\cdots
U_{\sigma_{N_{t}}+1}|j_{N_{t}}\rangle(Q_{\sigma_{0}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}})(i,j_{N_{t}})\times$ $\times\langle
j_{N_{t}}|U_{\sigma_{N_{t}}+1}^{*}\cdots U_{t}^{*}|k\rangle\Big{]}$

Let

$W(i,j,k):=\langle j|U_{t}\cdots U_{\sigma_{N_{t}}+1}|i\rangle\langle
i|U_{\sigma_{N_{t}}+1}^{*}\cdots U_{t}^{*}|k\rangle$ $=\langle j|U_{t}\cdots
U_{\sigma_{N_{t}}+1}|i\rangle\overline{\langle k|U_{t}\cdots
U_{\sigma_{N_{t}}+1}|i\rangle}$

So, we have that

$Tr\big{(}|k\rangle\langle
j|\rho_{t}\big{)}=E\big{[}Q_{\sigma_{0}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W\big{]}(i,j,k),$

here for any matrix $Q$, we define

$QW(i,j,l)=\sum_{l}Q(i,l)W(l,j,k)$

Then

$|Q_{\sigma_{0}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W(i,j,k)|=\big{|}\sum_{l}Q_{\sigma_{0}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}(i,l)W(l,j,k)\big{|}$
$\leq\sup_{l}||w(l,j,k)|\sum_{l}Q_{\sigma_{0}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}(i,l)=\sup_{l}|w(l,j,k)|.$

Since either $j$ or $k\neq l$, and $\zeta>0$, we may without loss of
generality, assume that $j\neq l$. Then

$\big{|}\langle j|U_{t}\cdots
U_{\sigma_{N_{t}+1}}|l\rangle\big{|}^{2}\rightarrow 0,a.s.$

by (5.25). Since

$\big{|}\langle k|U_{t}\cdots U_{\sigma_{N_{t}}}|l\rangle\big{|}^{2}\leq 1,$

we have

$\lim_{t\to\infty}\sup_{l}|w(l,j,k)|\rightarrow 0,a.s.$

Therefore, (5.28) holds by the Dominated Convergence Theorem.

Step 2. Now we consider for general $\rho_{0}$. For any $j,k$, by definition,

$Tr\big{(}|k\rangle\langle j|\rho_{t}\big{)}=\langle
j|\rho_{t}|k\rangle=\langle
j|(\Phi_{t}\cdots\Phi_{2}\Phi_{1}\rho_{0})|k\rangle$
$=\sum_{i_{1},i_{2},...,i_{t}=0}^{m}\langle j|A_{i_{t}}U_{t}\cdots
A_{i_{1}}U_{1}\rho_{0}U_{1}^{*}A_{i_{1}}^{*}\cdots
U_{t}^{*}A_{i_{t}}^{*}|k\rangle$
$=E\Big{[}\sum_{j_{1},...,j_{N_{t}}=1}^{m}\langle j|U_{t}\cdots
U_{\sigma_{N_{t}}+1}|j_{N_{t}}\rangle\langle
j_{N_{t}}|U_{\sigma_{N_{t}}}\cdots
U_{\sigma_{N_{t}-1}+1}|j_{N_{t}-1}\rangle\cdots$ $\cdots\langle
j_{1}|U_{\sigma_{1}}\cdots U_{2}U_{1}\rho_{0}U_{1}^{*}U_{2}^{*}\cdots
U_{\sigma_{1}}^{*}|j_{1}\rangle\cdots\langle
j_{N_{t}-1}|U_{\sigma_{N_{t}-1}+1}^{*}\cdots
U_{\sigma_{N_{t}}}^{*}|j_{N_{t}}\rangle\times$ $\times\langle
j_{N_{t}}|U_{\sigma_{N_{t}}+1}^{*}\cdots U_{t}^{*}|k\rangle\Big{]}$
$=E\Big{[}\sum_{j_{1}=1}^{m}\sum_{j_{N_{t}}=1}^{m}v_{j_{1}}\langle
j|U_{t}\cdots U_{\sigma_{N_{t}}+1}|j_{N_{t}}\rangle(Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}})(j_{1},j_{N_{t}})\times$ $\times\langle
j_{N_{t}}|U_{\sigma_{N_{t}}+1}^{*}\cdots U_{t}^{*}|k\rangle\Big{]}$

where $v_{j_{1}}=\langle j_{1}|U_{\sigma_{1}}\cdots
U_{2}U_{1}\rho_{0}U_{1}^{*}U_{2}^{*}\cdots U_{\sigma_{1}}^{*}|j_{1}\rangle$.

So, we have that

$Tr\big{(}|k\rangle\langle
j|\rho_{t}\big{)}=E\big{[}\sum_{j_{1}=1}^{m}v_{j_{1}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W\big{]}(j_{1},j,k).$

If $j\neq k$, then

$|Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W(j_{1},j,k)|=\big{|}\sum_{l}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}(i,l)W(l,j,k)\big{|}$
$\leq\sup_{l}||w(l,j,k)|\sum_{l}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}(i,l)=\sup_{l}|w(l,j,k)|.$

Since either $j$ or $k\neq l$, and $\zeta>0$, suppose that $j\neq l$

$\big{|}\langle j|U_{t}\cdots
U_{\sigma_{N_{t}+1}}|l\rangle\big{|}^{2}\rightarrow 0$

by (5.25) and

$\big{|}\langle k|U_{t}\cdots U_{\sigma_{N_{t}}}|l\rangle\big{|}^{2}\leq 1,$

we have

$\lim_{t\to\infty}\sup_{l}|w(l,j,k)|\rightarrow 0,a.s.$

If $j=k$, then

$\lim_{t\to\infty}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W(j_{1},j,k)=\pi_{j},a.s.,$

by Corollary 5.1.1. Since $v_{j}$ is a probability distribution, by the
Dominated Convergence Theorem, we have

$Tr\big{(}|k\rangle\langle
j|\rho_{t}\big{)}=E\big{[}\sum_{j_{1}=1}^{m}v_{j_{1}}Q_{\sigma_{1}}\cdots
Q_{\sigma_{N_{t}-1}}W\big{]}(j_{1},j,k)\to\pi_{j}.$

We have thus proved the theorem

∎

######  Remark 5.2.1.

Using a similar argument, we can show that the conclusion of Theorem 5.2 also
holds for Example 2.2, i.e., if $0<\zeta\leq 1$ and $0<p\leq 1$, then

$\rho_{t}\rightarrow\sum_{i=1}^{2}\pi_{i}\cdot|i\rangle\langle i|$

where $\pi_{i}=\frac{1}{2}$, for $i=1,2$.

##  6 Critical behaviors and simulations

In this section, we consider critical behaviors and discuss their numerical
results. Our main observations are the following,

  1. 1.

For $p=0$, our quantum Markov chain exhibits a time-inhomogeneous periodic
behavior for $\zeta\leq 2$ and non-periodic for $\zeta>2$. At time $n$, let
$T_{n}$ denote the first period after time $n$ for fixed $\lambda$ and
$\zeta$, we have for large $n$

     * •

If $\zeta<2$,

$\displaystyle
T_{n}\approx\Big{[}\big{(}1-\frac{\zeta}{2}\big{)}\frac{2\pi}{\lambda_{1}-\lambda_{2}}+n^{1-\frac{\zeta}{2}}\Big{]}^{\frac{1}{1-\frac{\zeta}{2}}}-n,$
(6.29)

     * •

If $\zeta=2$,

$\displaystyle T_{n}\approx
n\cdot(e^{\frac{2\pi}{\lambda_{1}-\lambda_{2}}}-1),$ (6.30)

where $\lambda_{1}>\lambda_{2}$ are two eigenvalues of the matrix
$G=\left[\begin{array}[]{ c c }\lambda&\lambda\\\
\lambda&\lambda\end{array}\right]$. Moreover, let $T_{n}^{(\lambda)}$ denote
the first period after time $n$ for fixed $\zeta$ depending on $\lambda$, the
relation

$\displaystyle T_{n}^{(\lambda)}\sim\frac{T_{n}^{(1)}}{\lambda}$ (6.31)

can be directly deduced by Formulas (6.29) and (6.30) using Taylor series
expansion.

  2. 2.

For small $p>0$, the time-inhomogeneous periodicity behaves the same manner as
the case $p=0$.

  3. 3.

For small $p>0$, the amplitudes of each period decays exponentially as $n$
gets large. The exponential decay rate is approximately $p/2$, independent of
the choices of $\zeta\leq 2$ and $\lambda>0$.

  4. 4.

For $p\approx 1$, the decay to limiting distribution is exponentially fast if
$\zeta$ is small, and it is power law decay if $\zeta$ is large. The critical
point $\zeta_{c}$ appears in the interval $[0.6,0.7]$.

  5. 5.

For $p>0$, the critical phase transition occurs when $\zeta_{c}=1$ for the
decoherent Markov chain.

  6. 6.

Instead of fully connected graph we considered, simulation results show that
the decoherent system converges equilibrium limit with $N$ dimensional cyclic
graph $G$. For $p=0$, time-inhomogeneous periodicity is not predictable
although some kind of periodicity was observed.

In what follows, we will discuss the numerical results for the above
observations assuming that all the entries of the matrix $G$ are $\lambda>0$
for simplicity.

###  6.1 Pure quantum case, $p=0$, critical behavior at $\zeta_{o}=2$

Let us consider the pure quantum case when $p=0$. For $2$ dimensional case
$i,j=1,2$, recall Definition 2.1 and Equation 2.2, the probability starting
from the state $|i\rangle$ and getting the state $|j\rangle$ after $n$ steps
can be deduced to

$\displaystyle P(i,j)=\big{|}\langle j|U_{n}\cdots U_{1}|i\rangle\big{|}^{2},$
(6.32)

using the fact $U_{n}=e^{i\frac{G}{\sqrt{n^{\zeta}}}}$, where $G=$ is a self-
adjoint matrix with non-zero entries,

$P(i,j)=\big{|}\langle
j|e^{iG\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}|i\rangle\big{|}^{2}=\big{|}\langle
j|B\left[\begin{array}[]{ c c
}e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}&0\\\
0&e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\end{array}\right]B^{*}|i\rangle\big{|}^{2},$

where $\lambda_{1}$ and $\lambda_{2}$ are eigenvalues of
$G=\left[\begin{array}[]{ c c }\lambda&\lambda\\\
\lambda&\lambda\end{array}\right]$ with $\lambda_{1}>\lambda_{2}$, and B the
associated orthogonal matrix. Let

$\Lambda_{n}=\left[\begin{array}[]{ c c
}e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}&0\\\
0&e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\end{array}\right],$

and we have

$P_{n}(i,j)=\Big{|}\sum_{l=1}^{2}\langle j|B|l\rangle\langle
l|\Lambda_{n}|l\rangle\langle
l|B^{*}|i\rangle\Big{|}^{2}=\Big{|}\sum_{l=1}^{2}\langle
j|B|l\rangle\cdot\Lambda_{n}(l,l)\cdot\langle l|B^{*}|i\rangle\Big{|}^{2}$
$=\Big{|}\langle j|B|1\rangle\cdot\Lambda_{n}(1,1)\cdot\langle
1|B^{*}|i\rangle+\langle j|B|2\rangle\cdot\Lambda_{n}(2,2)\cdot\langle
2|B^{*}|i\rangle\Big{|}^{2}$ $=\Big{|}\langle j|B|1\rangle\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\langle
1|B^{*}|i\rangle+\langle j|B|2\rangle\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\langle
2|B^{*}|i\rangle\Big{|}^{2}.$

Suppose that $i=j=1$, with $B=\left[\begin{array}[]{ c c }b_{11}&b_{12}\\\
b_{21}&b_{22}\end{array}\right]$, we obtain

$P_{n}(1,1)=\Big{|}\langle 1|B|1\rangle\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\langle
1|B^{*}|1\rangle+\langle 1|B|2\rangle\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\langle
2|B^{*}|1\rangle\Big{|}^{2}$ $=\Big{|}B(1,1)\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot
B^{*}(1,1)+B(1,2)\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot
B^{*}(2,1)\Big{|}^{2}$ $=\Big{|}b_{11}\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\bar{b}_{11}+b_{12}\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\cdot\bar{b}_{12}\Big{|}^{2}=\Big{|}|b_{11}|^{2}\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}+|b_{12}|^{2}\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\Big{|}^{2}$
$=\Big{(}|b_{11}|^{2}\cdot
e^{i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}+|b_{12}|^{2}\cdot
e^{i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\Big{)}\Big{(}|b_{11}|^{2}\cdot
e^{-i\lambda_{1}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}+|b_{12}|^{2}\cdot
e^{-i\lambda_{2}\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}\Big{)}$
$=|b_{11}|^{4}+|b_{11}|^{2}|b_{12}|^{2}e^{i(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}+|b_{12}|^{4}+|b_{12}|^{2}|b_{11}|^{2}e^{i(\lambda_{2}-\lambda_{1})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}}$
$=|b_{11}|^{4}+|b_{12}|^{4}+2|b_{12}|^{2}|b_{11}|^{2}\cos\big{(}(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}\big{)}$
$=1-2|b_{12}|^{2}|b_{11}|^{2}\Big{[}1-\cos\big{(}(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}\big{)}\Big{]}$

Finally, we conclude

$P_{n}(1,1)=1-2|b_{12}|^{2}|b_{11}|^{2}\Big{[}1-\cos\big{(}(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}\big{)}\Big{]},$

and

$P_{n}(2,1)=1-\Big{(}1-2|b_{12}|^{2}|b_{11}|^{2}\Big{[}1-\cos\big{(}(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}\big{)}\Big{]}\Big{)}$
$=2|b_{12}|^{2}|b_{11}|^{2}\Big{[}1-\cos\big{(}(\lambda_{1}-\lambda_{2})\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}\big{)}\Big{]}.$

We observe from above deduction that the probability at $n$ depends on the
eigenvalues and eigenvectors of $G$ and $\zeta$. In our case, $B$ has non-zero
entries with $\lambda_{1}>\lambda_{2}$, and the probabilities in $P_{n}$ has
periodicity property if $\zeta\leq 2$, the series
$\displaystyle\sum_{k=1}^{n}\frac{1}{\sqrt{k^{\zeta}}}$ is divergent as $n$
tends to infinity. On the other hand, if $\zeta>2$, $P_{n}$ converges, and the
limit depends on the matrix $G$ and the value of $\zeta$.

(a) $\zeta=1$

(b) $\zeta=1.5$

(c) $\zeta=2$

Figure 1: Simulation results for $P_{n}$ fixing $p=0$ and $\lambda=1$ with
$t=50000$

(a) $\lambda=0.5$

(b) $\lambda=1.5$

(c) $\lambda=2$

Figure 2: Simulation results for $P_{n}$ fixing $p=0$ and $\zeta=1$ with
$t=50000$

Figures 1 and 2 show the periodicity behavior of the $2\times 2$ matrix
$P_{n}$ by fixing $\zeta$ and $\lambda$ with respectively. We observe from
Figure 1 that with fixed $\lambda$, the periodicity of $P_{n}$ decreases if
$\zeta$ increases. On the other hand, Figure 2 illustrates the phenomenon that
with fixed $\zeta$, the periodicity of $P_{n}$ increases if $\lambda$
increases.

To analyze the periodicity of $P_{n}$ when $\zeta\leq 2$, we use the fact that
$\cos(\cdot)$ has period $2\pi$, and let $T_{n}$ be the period for fixed $n$.

$(\lambda_{1}-\lambda_{2})\sum_{k=n+1}^{n+T_{n}}\frac{1}{\sqrt{k^{\zeta}}}=2\pi.$

Asymptotically, we have for large $n$

$\int_{n}^{n+T_{n}}\frac{1}{\sqrt{x^{\zeta}}}dx\approx\frac{2\pi}{\lambda_{1}-\lambda_{2}}.$

If $\zeta<2$, it can be deduced to

$\Big{[}\frac{x^{1-\frac{\zeta}{2}}}{1-\frac{\zeta}{2}}\Big{]}_{x=n}^{x=n+T_{n}}\approx\frac{2\pi}{\lambda_{1}-\lambda_{2}}\implies\frac{(n+T_{n})^{1-\frac{\zeta}{2}}}{1-\frac{\zeta}{2}}-\frac{n^{1-\frac{\zeta}{2}}}{1-\frac{\zeta}{2}}\approx\frac{2\pi}{\lambda_{1}-\lambda_{2}},$

then, we have

$(n+T_{n})^{1-\frac{\zeta}{2}}-n^{1-\frac{\zeta}{2}}\approx(1-\frac{\zeta}{2})\frac{2\pi}{\lambda_{1}-\lambda_{2}},$

and

$\displaystyle
T_{n}\approx\Big{[}\big{(}1-\frac{\zeta}{2}\big{)}\frac{2\pi}{\lambda_{1}-\lambda_{2}}+n^{1-\frac{\zeta}{2}}\Big{]}^{\frac{1}{1-\frac{\zeta}{2}}}-n.$
(6.33)

If $\zeta=2$, we have

$\Big{[}\ln|x|\Big{]}_{x=n}^{x=n+T_{n}}\approx\frac{2\pi}{\lambda_{1}-\lambda_{2}}\implies\ln|n+T_{n}|-\ln|n|\approx\frac{2\pi}{\lambda_{1}-\lambda_{2}},$

then,

$\displaystyle T_{n}\approx
n\cdot(e^{\frac{2\pi}{\lambda_{1}-\lambda_{2}}}-1).$ (6.34)

Figure 3 compares the periods generated by simulation and Formula (6.33) with
different values of $\lambda$ and $\zeta$. As results, we observe that the
ratios between the simulated periods and theoretical periods are $1$ as $n$ is
large, which numerically proved the asymptotic formulas discussed.

(a) $\lambda=1$ and $\zeta=0.8$

(b) $\lambda=1$ and $\zeta=1.2$

(c) $\lambda=0.5$ and $\zeta=1.2$

(d) $\lambda=1.5$ and $\zeta=1.2$

Figure 3: Simulated vs theoretical periods for $P_{n}$ with $t=50000$

###  6.2 Periodicity and decay rates for decoherent quantum Markov chain,
$p>0$

Previously, we have numerically shown the periodic behavior for the pure
quantum Markov chains using Formulas (6.33), (6.34). For $p>0$, despite the
equilibrium convergence property proved in Theorem 5.2, we numerically analyze
the periodicity and convergence rates to equilibrium for decoherent Markov
chains.

(a) $\lambda=0.5$ and $\zeta=0.4$

(b) $\lambda=0.5$ and $\zeta=0.6$

(c) $\lambda=0.5$ and $\zeta=0.8$

Figure 4: Probabilities for $P_{n}(1,1)$ with $t=1000$

Regarding the periodicity of the decoherent system, if $p$ is close to $0$,
the probability to make measurements is low, the system maintains more quantum
behavior by intuition. Conversely, the system tends to have less quantum
coherence while $p$ increases. Simulation results not only accord with above
observation, but also demonstrate that only small $p$’s conserve notable
quantum periodic behavior (see Figure 4). Furthermore, the results show that
for larger $\zeta$’s, fixed $p$’s and $\lambda$’s, the periodicity decreases
similarly as the coherent case.

(a) $\lambda=0.5$ and $\zeta=0.4$

(b) $\lambda=0.5$ and $\zeta=0.6$

(c) $\lambda=0.5$ and $\zeta=0.8$

Figure 5: Simulated vs pure quantum theoretical periods for $P_{n}$ with
$p=0.01$ with $t=2000$

(a) $\lambda=0.5$ and $\zeta=0.4$

(b) $\lambda=0.5$ and $\zeta=0.6$

(c) $\lambda=0.5$ and $\zeta=0.8$

Figure 6: Simulated vs pure quantum theoretical periods for $P_{n}$ with
$p=0.03$ with $t=2000$

We also calculate the ratios between the simulated periodicity and the
theoretical periodicity from Formula (6.33) fixing different values of $p$ and
$\zeta$. As results, Figures 5 and 6 show that the ratios are asymptotically
$1$. In other words, even for $p>0$, the decoherent system maintains the same
periodicity as in pure quantum case for small $p$.

In order to analyze the convergence rates for small $p$ and different values
of $\zeta$’s and $\lambda$’s, we fit the local maximums of the $P_{n}(1,1)$ in
non-linear regression with exponential decay model: $ce^{-rt}-\frac{1}{2}$,
and rational decay model: $ct^{-r}-\frac{1}{2}$ for $p$ near $0$. Since for
every choice of $\lambda$ and $\zeta$, the number of local maximums are
different, we consider the adjusted $R^{2}$ coefficients to perform the
comparison. Consequently, for all parameters mentioned above, not only the
adjusted $R^{2}$ of exponential models with decay rate $r$ are greater than
rational models’ adjusted $R^{2}$, but also they are all greater than $0.9$.
Even more, we notice that the linear relationship between $r$ and $p$ when $p$
tends to $0$ for all $0\leq\zeta\leq 1$, which is, the convergence rates are
proportional to $p$ and they are independent of $\zeta$ and $\lambda$ as $p$
tends to $0$.

More specifically, Tables 1, 2 and 3 show the convergence rates obtained from
exponential decay non-linear regression model for $\lambda=0.2$, $0.35$ and
$0.5$, and we numerically conclude that the decay rates $r\sim\frac{p}{2}$ as
$p$ approaches $0$. In other words, when environmental interaction is
extremely low, the convergence to equilibrium does not depend on the
parameters $\zeta$ and $\lambda$.

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.005$ | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025
$0.01$ | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.0051
$0.015$ | 0.0075 | 0.0075 | 0.0075 | 0.0077 | 0.0077 | 0.0077 | 0.0077 | 0.0077 | 0.0076 | 0.0077
$0.02$ | 0.0101 | 0.0101 | 0.0101 | 0.0101 | 0.0101 | 0.0102 | 0.0102 | 0.0101 | 0.0102 | 0.0103
$0.025$ | 0.0126 | 0.0126 | 0.0127 | 0.0127 | 0.0127 | 0.0127 | 0.0128 | 0.0127 | 0.0128 | 0.0129
Table 1: Estimated convergence rates $r$ for $\lambda=0.2$ by non-linear regression model: $ce^{-rt}$ $p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.005$ | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025
$0.01$ | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005
$0.015$ | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0076 | 0.0076 | 0.0076 | 0.0076
$0.02$ | 0.0101 | 0.01 | 0.01 | 0.01 | 0.0101 | 0.0101 | 0.0101 | 0.0101 | 0.0102 | 0.0102
$0.025$ | 0.0126 | 0.0124 | 0.0125 | 0.0125 | 0.0126 | 0.0126 | 0.0127 | 0.0127 | 0.0128 | 0.0126
Table 2: Estimated convergence rates $r$ for $\lambda=0.35$ by non-linear regression model: $ce^{-rt}$ $p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.005$ | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025
$0.01$ | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005 | 0.005
$0.015$ | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0075 | 0.0076 | 0.0075
$0.02$ | 0.0101 | 0.0101 | 0.01 | 0.0101 | 0.0101 | 0.01 | 0.01 | 0.01 | 0.0101 | 0.0101
$0.025$ | 0.0126 | 0.0127 | 0.0125 | 0.0127 | 0.0127 | 0.0125 | 0.0125 | 0.0126 | 0.0127 | 0.0125
Table 3: Estimated convergence rates $r$ for $\lambda=0.5$ by non-linear
regression model: $ce^{-rt}$

On the other hand, we similarly fit the probabilities $P_{n}(1,1)$ with
$n=200$ in non-linear regression with exponential decay model:
$ce^{-rt}-\frac{1}{2}$, and rational decay model: $ct^{-r}-\frac{1}{2}$ for
different values of $\lambda$, $\zeta$ and $p$ close to $1$ where the
probability of measurement is high. In this case, we compare the $R^{2}$
coefficients for above parameters since there have no local maximums for $p$
close to $1$. Tables 4, 5 and 6 show that for $\lambda=0.2$, $0.35$, $0.5$ and
$p$ close to $1$, exponential decay model fits better for $0<\zeta\leq 0.6$
with $R^{2}$ coefficients greater than $0.9$ while rational decay model fits
better for $0.7<\zeta\leq 1$ with $R^{2}>0.9$. Therefore, not only we notice
that the decay behavior is very different than when $p$ is close to 0, but
also there is a phase of transition $\zeta_{o}$ such that $0.6\leq\zeta\leq
0.7$ for $p\sim 1$.

Indeed, we also observe from these tables that the convergence rates decrease
in function of $p$ when $p$ tends to $0$, and the convergence rates increase
when $p$ tends to $1$.

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.7$ |  |  | 1 |  |  |  |  |  | 1 | 2
$0.8$ | 1 |  | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.9$ | 1 | 1 |  | 1 | 1 | 2 | 2 | 2 | 2 | 2
$1$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2

Model Fitness: 1 for Exponential decay model and 2 for Rational decay model

$\zeta\ \backslash\ p$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---
$0.1$ |  | 1.33 | 0.84 | 0.54
$0.2$ |  |  | 0.74 | 0.47
$0.3$ | 1.65 | 1.06 |  | 0.40
$0.4$ |  | 0.94 | 0.56 | 0.34
$0.5$ |  | 0.83 | 0.47 | 0.27
$0.6$ |  | 0.72 |  | 0.21
$0.7$ |  |  |  |
$0.8$ |  |  |  |
$0.9$ | 0.84 |  |  |
$1$ |  |  |  |

Exponential Decay Rates

$\zeta\ \backslash\ p$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---
$0.1$ |  |  |  |
$0.2$ |  |  |  |
$0.3$ |  |  |  |
$0.4$ |  |  |  |
$0.5$ |  |  |  |
$0.6$ |  |  | 1.11 |
$0.7$ |  | 1.32 | 1.01 | 0.83
$0.8$ |  | 1.18 | 0.89 | 0.72
$0.9$ |  | 1.02 | 0.75 | 0.6
$1$ | 1.4 | 0.84 | 0.61 | 0.47

Rational Decay Rates

Table 4: Estimated convergence rates $r$ for $\lambda=0.5$ by non-linear
regression model

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.6$ |  |  | 1 |  | 1 | 1 | 1 | 2 | 2 | 2
$0.7$ |  | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.8$ |  |  | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.9$ |  | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$1$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2

Model Fitness: 1 for Exponential decay model and 2 for Rational decay model

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ |  |  |  |  | 0.23
$0.2$ |  |  |  | 0.25 | 0.19
$0.3$ | 0.64 | 0.45 | 0.3 | 0.21 | 0.15
$0.4$ |  | 0.38 | 0.24 | 0.16 | 0.11
$0.5$ | 0.49 | 0.31 | 0.19 | 0.12 | 0.08
$0.6$ | 0.42 | 0.24 | 0.14 | 0.08 | 0.05
$0.7$ | 0.34 |  |  |  |
$0.8$ |  |  |  |  |
$0.9$ |  |  |  |  |
$1$ |  |  |  |  |

Exponential Decay Rates

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ |  |  |  |  |
$0.2$ |  |  |  |  |
$0.3$ |  |  |  |  |
$0.4$ |  |  |  |  |
$0.5$ |  |  |  |  |
$0.6$ |  |  |  |  |
$0.7$ |  | 0.85 | 0.71 | 0.61 | 0.52
$0.8$ | 0.93 | 0.73 | 0.6 | 0.49 | 0.41
$0.9$ | 0.79 | 0.6 | 0.48 | 0.39 | 0.32
$1$ | 0.64 | 0.47 | 0.37 | 0.29 | 0.23

Rational Decay Rates

Table 5: Estimated convergence rates $r$ for $\lambda=0.35$ by non-linear
regression model

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.6$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.7$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.8$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.9$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$1$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2

Model Fitness: 1 for Exponential decay model and 2 for Rational decay model

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ | 0.19 | 0.13 | 0.1 | 0.08 | 0.06
$0.2$ | 0.15 | 0.11 | 0.08 | 0.06 | 0.05
$0.3$ | 0.11 | 0.08 | 0.06 | 0.04 | 0.03
$0.4$ | 0.08 | 0.05 | 0.04 | 0.03 | 0.02
$0.5$ | 0.05 | 0.03 | 0.02 | 0.02 | 0.01
$0.6$ | 0.03 | 0.02 | 0.01 | 0.01 | 0.01
$0.7$ |  |  |  |  |
$0.8$ |  |  |  |  |
$0.9$ |  |  |  |  |
$1$ |  |  |  |  |

Rxponential Decay Rates

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ |  |  |  |  |
$0.2$ |  |  |  |  |
$0.3$ |  |  |  |  |
$0.4$ |  |  |  |  |
$0.5$ |  |  |  |  |
$0.6$ |  |  |  |  |
$0.7$ | 0.44 | 0.36 | 0.3 | 0.25 | 0.21
$0.8$ | 0.34 | 0.27 | 0.22 | 0.18 | 0.15
$0.9$ | 0.26 | 0.2 | 0.16 | 0.13 | 0.11
$1$ | 0.19 | 0.15 | 0.12 | 0.1 | 0.08

Rational Decay Rates

Table 6: Estimated convergence rates $r$ for $\lambda=0.2$ by non-linear
regression model

###  6.3 High dimensional quantum Markov chain for dimension $N>2$

It would be interesting to see if our observations above also hold for higher
dimensional quantum Markov chain.

For the pure quantum case, we first observe from the simulation results that
periodicity behavior remains and the probability will concentrates on the
diagonal of the matrix $P_{n}$ varying from $\frac{1}{N}$ to $1$ where $N$ is
the dimension of the space and $P_{n}$ is a $N\times N$ matrix, and the off-
diagonal entries vary from $0$ to $\frac{1}{N}$ due to the symmetry property
(See Figure 7).

On the other hand, the periodicity for high dimensional Markov chains also
follows the theoretical formula (6.33) for 2 dimensional case. Figure 8 shows
that the ratio between the periods simulated and the theoretical period from
(6.33) is asymptotically $1$ for $\lambda=1$, $\zeta=1$ and space dimension
$N=3,4,5$.

(a) $N=3$

(b) $N=5$

Figure 7: Simulated probabilities for $P_{n}(1,:)$ for $p=0$, $\lambda=1$,
$\zeta=1$, dimension $N$ and $t=500$

(a) $N=3$

(b) $N=4$

(c) $N=5$

Figure 8: Simulated vs pure quantum theoretical periods for $P_{n}$,
$\lambda=1$ and $\zeta=1$, dimension $N$ with $t=2000$

Consider now the decoherent case when $p$ near $0$. Unlike the symmetrical
decay to the equilibrium of the $2$ dimensional case, the local maximums of
the diagonal decrease to the equilibrium limit from $1$ to $\frac{1}{N}$ while
the local minimums maintain near $\frac{1}{N}$. On the other hand, for the
off-diagonal entries, the local minimums increase to the equilibrium limit
from $0$ to $\frac{1}{N}$. Figure 9 illustrates the $4$ dimensional decoherent
Markov chain convergence to the equilibrium of the first row of the matrix
$P_{n}$ when $p=0.005$, $\lambda=1$, $\zeta=0.5$ and $t=1000$.

Figure 9: Simulated probabilities for $P_{n}(1,:)$ for $p=0.005$, $\lambda=1$,
$\zeta=0.5$, dimension $N=4$ and $t=1000$

Similarly, we fit the local maximums of the probabilities at the diagonal
entries of the matrix $P_{n}$ in non-linear regression with exponential decay
model: $ce^{-rt}-\frac{1}{N}$ and rational decay model: $ct^{-rt}-\frac{1}{N}$
to estimate the decay rates, and compare their adjusted $R^{2}$ coefficients.
As results, the exponential decay model has better fit than rational model
with adjusted $R^{2}$ coefficients greater than $0.9$ for different values of
$\zeta$ and $\lambda$ for $p\sim 0$. Even further, the convergence rates near
$p\sim 0$ has the similar behavior as $2$ dimensional case, they are
proportional to $p$ by the relation $r=\frac{p}{2}$, and independent of the
other parameters (See Table 7).

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.005$ | 0.0024 | 0.0024 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025 | 0.0025
$0.01$ | 0.0047 | 0.0048 | 0.0049 | 0.005 | 0.005 | 0.0051 | 0.005 | 0.0051 | 0.0051 | 0.005
$0.015$ | 0.0071 | 0.0072 | 0.0074 | 0.0075 | 0.0076 | 0.0075 | 0.0076 | 0.0075 | 0.0076 | 0.0075
$0.02$ | 0.0095 | 0.0095 | 0.0099 | 0.0099 | 0.01 | 0.0102 | 0.0101 | 0.0101 | 0.0101 | 0.01
$0.025$ | 0.0119 | 0.0119 | 0.0124 | 0.0125 | 0.0126 | 0.0128 | 0.0126 | 0.0127 | 0.0126 | 0.0125
Table 7: Estimated $5$ dimensional convergence rates $r$ for $\lambda=0.5$ by
non-linear regression model: $ce^{-rt}$

For $p$ close to $1$, we also fit the diagonal entries of the matrix $P_{n}$
in the same non-linear regression model, the same phase of transition point we
obtained in the $2$ dimensional case also occurs in high dimensional case.
Table 8 shows the best fit non-linear regression results: exponential decay
model or rational decay model with $R^{2}$ coefficients greater than $0.9$ and
their estimated rates $r$ with respectively as before, and we numerically
obtained that the phase of transition is between $0.6<\zeta<0.7$ (See Table
8).

$p\ \backslash\ \zeta$ | $0.1$ | $0.2$ | $0.3$ | $0.4$ | $0.5$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---|---|---|---|---|---
$0.6$ | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2
$0.7$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.8$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$0.9$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2
$1$ | 1 | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 2 | 2

Model Fitness: 1 for Exponential decay model and 2 for Rational decay model

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ | 0.29 | 0.27 | 0.24 | 0.2 | 0.17
$0.2$ | 0.26 | 0.24 | 0.2 | 0.17 | 0.14
$0.3$ | 0.23 | 0.2 | 0.17 | 0.14 | 0.11
$0.4$ | 0.2 | 0.17 | 0.14 | 0.1 | 0.08
$0.5$ | 0.17 | 0.14 | 0.1 | 0.07 | 0.05
$0.6$ | 0.13 | 0.1 | 0.07 | 0.05 | 0.03
$0.7$ |  |  |  |  | 0.01
$0.8$ |  |  |  |  |
$0.9$ |  |  |  |  |
$1$ |  |  |  |  |

Exponential Decay Rates

$\zeta\ \backslash\ p$ | $0.6$ | $0.7$ | $0.8$ | $0.9$ | $1$
---|---|---|---|---|---
$0.1$ |  |  |  |  |
$0.2$ |  |  |  |  |
$0.3$ |  |  |  |  |
$0.4$ |  |  |  |  |
$0.5$ |  |  |  |  |
$0.6$ |  |  |  |  |
$0.7$ | 0.4 | 0.47 | 0.55 | 0.63 |
$0.8$ | 0.62 | 0.53 | 0.45 | 0.38 | 0.32
$0.9$ | 0.52 | 0.43 | 0.36 | 0.3 | 0.25
$1$ | 0.43 | 0.35 | 0.28 | 0.23 | 0.19

Rational Decay Rates

Table 8: Estimated $5$ dimensional convergence rates $r$ for $\lambda=0.2$ by
non-linear regression model

###  6.4 Decoherent case, $p>0$, critical behavior $\zeta_{c}=1$

Recall our mainly result (Theorem 5.2) states that for $0<\zeta\leq 1$, the
$N$ dimensional time-inhomogeneous decoherent Markov chains converge to theirs
equilibrium limit, $\frac{1}{N}$. However, we expect that the critical
transition occurs exactly when $\zeta_{c}=1$ which means that if $\zeta>1$,
the large time scale limit is not $\frac{1}{N}$ and depends on the the values
of $\lambda$ and the decoherence parameter $p$.

Figure 10 support our expectation and shows the simulation results for $2$
dimensional $P_{n}(1,1)$ with $\zeta=1.1$ and different values of $\lambda$
and the decoherent parameter $p$. We observe from these figures that for fixed
$\zeta$, the probabilities converge, and they don’t converge to $\frac{1}{2}$
which is uniform distribution. In other words, the large time scale limits
depend on the values of $\zeta$, $\lambda$ and the decoherence parameter $p$.
More specifically, the large time scale limits depend on the initial
conditions.

(a) $\lambda=0.3$

(b) $\lambda=0.5$

(c) $\lambda=0.7$

Figure 10: Convergence of $2$ dimensional $P_{n}(1,1)$ fixing $\zeta=1.1$ for
different $\lambda$’s and $p$’s

###  6.5 Time-inhomogeneous quantum Markov chain with cyclic graph

Instead of the matrix $G$ with strictly positive entries we considered
previously, let us consider the quantum Markov chain with cyclic graphs.
Consider the same decoherent quantum Markov chain model with the matrix

$G=\left[\begin{array}[]{ c c c c c}0&\lambda&0&0&\lambda\\\
\lambda&0&\lambda&0&0\\\ 0&\lambda&0&\lambda&0\\\ 0&0&\lambda&0&\lambda\\\
\lambda&0&0&\lambda&0\end{array}\right]$

where $\lambda>0$.

For the pure quantum $p=0$ case, simulation results show that the time-
inhomogeneous periodic behavior exists. Moreover, probabilities still
concentrate more at the diagonal entries of $P_{n}$ varying near $\frac{1}{N}$
to 1 while the probabilities on the off-diagonal entries vary from $0$ to
values near $\frac{1}{N}$ due to the symmetry.

On the other hand, the probabilities at each entry also numerically converge
to $\frac{1}{N}$ for relatively small $\zeta$’s and $\lambda$’s.

Figure 11(a) illustrates the periodicity with local maximums and minimums of
the first column of the matrix $P_{n}$ for $5$ dimensional pure cyclic quantum
Markov chain with $\lambda=\frac{1}{2}$ and $\zeta=\frac{1}{2}$, Figures 11(b)
and 11(c) show the convergence to equilibrium property for $5$ dimensional
decoherent cyclic quantum Markov chain to $\frac{1}{5}$ with the same
parameters when $p=0.01$ and the completely decoherent case $p=1$. Although G
does not have all positive entries, we still believe that the transition
points are $\zeta_{c}=1$ and $\zeta_{o}=2$ for decoherent case and pure
quantum case respectively by similar argument in section 6.1.

Furthermore, we believe that the phase of transition occurs at $\zeta_{c}=1$
and $\zeta_{o}=2$ for any irreducible $G$.

(a) $p=0$, $\lambda=0.5$, $\zeta=0.5$

(b) $p=0.01$, $\lambda=0.5$, $\zeta=0.5$

(c) $p=1$, $\lambda=0.5$, $\zeta=0.5$

Figure 11: Simulated probabilities for $P_{n}(1,:)$ using cyclic graph $G$
with dimension $5$ and $t=500$

##  7 Conclusion and future work

In this paper, we mainly considered the time-inhomogeneous unitary operators,
and defined the time-inhomogeneous quantum analogue of the classical Markov
chain with decoherence parameter on finite dimensional state space and we
interpreted the decoherent parameter as the probability to perform a
measurement, that means that at each step, we perform a measurement with a
certain probability. We proved the Markov properties at the geometric
measurement times using the path integral representation.

We also made the conclusion that the time-inhomogeneous quantum Markov chain
on finite dimensional state spaces with non zero probability of measurement
converges to an equilibrium limit when $0\leq\zeta\leq 1$ as time approaches
to infinity.

Additionally, we analyzed the periodicity, the decay rate and the convergence
of our quantum Markov chains. We numerically concluded that our quantum Markov
chain exhibits a time-inhomogeneous periodic and non-periodic behaviors with
phase of transition $\zeta_{o}=2$ for the probability of measurement close to
zero, and a large time scale equilibrium behavior with phase of transition
$\zeta_{c}=1$ for the positive probability of measurement. Moreover, the
decoherent Markov chain decays exponentially for small $\zeta$ and decays with
power rule for large $\zeta$ with critical point $\zeta_{d}$ between $0.6$ and
$0.7$.

Instead of fully connected graph, numerical simulations were also done
considering the cyclic graph, and illustrated similar periodic behavior for
coherent quantum Markov chain and convergence to equilibrium for decoherent
quantum Markov chain for small $\lambda$’s $\zeta$’s. As result, we expect
that the phase of transition occurs at $\zeta_{c}=1$ and $\zeta_{o}=2$ for
cyclic graphs.

Not to mention, our work is not only based on finite quantum state space, also
called ”Qudit” in high dimensional quantum computation which researchers have
developed many applications recently including quantum circuit building,
quantum algorithm design and quantum experimental methods (see [17]), but
might also be related to one of the fundamental aspects of quantum annealing
algorithm which was first proposed by B. Apolloni, N. Cesa Bianchi and D. De
Falco in [2] and [3] and later formulated by T. Kadowaki and H. Nishimori in
[11]. With this motivation, our future work will focus on applications of
quantum search algorithms with time-inhomogeneous quantum Markov chain in high
dimensional spaces or general graph.

For $0<p\leq 1$ and $\zeta=0$, our model is related to open quantum random
walks in [4] and [5] which has a lot of potential applications in quantum
computing. It would be very interesting to see how open quantum walk can be
generalized to time-inhomogeneous case.

Furthermore, our numerical results will also lead us challenging analytical
problems, e.g., generalizing our results to time-inhomogeneous quantum Markov
chain associated with any connected graphs.

## References

  * [1] D. Aharonov, A. Ambainis, J. Kempe and U. Vazirani, Quantum walks on graphs, Proceedings of the 33rd annual ACM symposium on Theory of computing, 50-59, (2001).
  * [2] B. Apolloni, N. Cesa-Bianchi and D. De Falco, A numerical implementation of quantum annealing, Stochastic Processes, Physics and Geometry, Proceedings of the Ascona-Locarno Conference., 97-111, (1988).
  * [3] S. Attal, N. Guillotin-Plantard and C. Sabot, Central limit theorems for open quantum random walks and quantum measurement records, Preprint: arXiv:1206.1472, (2013).
  * [4] S. Attal, F. Petruccione, C. Sabot and I. Sinayskiy, Open quantum random walks, Journal of Statistical Physics. 147: 832–852, (2012).
  * [5] B. Apolloni, M. Carvlho and D. De Falco, Quantum stochastic optimization, Stoc. Proc. Appl. 33 (2): 233–244, (1989).
  * [6] T. Brun, H. Carteret and A. Ambainis, Quantum random walks with decoherent coins, Phys. Rev., A 67, 032304, (2003).
  * [7] R. Durrett, Probability: theory and examples, 5th edition, Cambridge, (2019).
  * [8] S. Fan, Z. Feng, S. Xiong and W. Yang, Convergence of quantum random walks with decoherence, Phys. Rev., A 84, 042317, (2011).
  * [9] L. Grover, A fast quantum mechanical algorithm for database search, Proceeding of 28th annual ACM symposium Theory of computing, 212-219, (1996).
  * [10] W. Hastings and Monte Carlo, Sampling methods using Markov chains and their applications, Biometrika 57 1, 97-109, (1970).
  * [11] T. Kadowaki and H. Nishimori, Quantum annealing in the transverse Ising model. Phys. Rev. E. 58 (5): 5355, (1998).
  * [12] M. Lagro, W. Yang and S. Xiong, A Perron–Frobenius Type of Theorem for Quantum Operations, Journal of Statistical Physics, 169, 38-62, (2017).
  * [13] D. Levin, Y. Peres and E. Wilmer, Markov Chains and Mixing Times, AMS, (2008).
  * [14] N. Metropolis and S. Ulam, The Monte Carlo method, Journal of the American statistical association, 44 247, 124-134, (1949).
  * [15] P. Shor, Algorithms for quantum computation: Discrete logarithms and factoring, Foundations of Computer Science, 1994 Proceedings 35th Annual Symposium, 335-341, (1994).
  * [16] M. Szegedy, Quantum speed-up of Markov chain based algorithms, Foundations of Computer Science Proceedings 45th Annual IEEE Symposium, 32-41, (2004).
  * [17] Y. Wang, Z. Hu, B Sanders and S. Kais, Qudits and high-dimensional quantum computing, Front. Phys., 8, 479. 10.3389, (2020).
  * [18] H. Zeh, On the Interpretation of Measurement in Quantum Theory, Foundations of Physics, 1, 69-79, (1970).
  * [19] K. Zhang, Limiting distribution of decoherent quantum random walks, Phys. Rev., A 77, 062302, (2008).