# Checking Continuous Stochastic Logic against Quantum Continuous-Time Markov
Chains

Jingyi Mei  Shanghai Key Lab of Trustworthy Computing, & MoE Engineering
Research Center of Software/Hardware Co-design Technology and Application,
East China Normal University, China<EMAIL_ADDRESS>,  Ming Xu  Shanghai Key
Lab of Trustworthy Computing, & MoE Engineering Research Center of
Software/Hardware Co-design Technology and Application, East China Normal
University, China<EMAIL_ADDRESS>,  Ji Guan  State Key Lab of Computer
Science, Institute of Software, Chinese Academy of Sciences, China
<EMAIL_ADDRESS>,  Yuxin Deng  Shanghai Key Lab of Trustworthy Computing,
& MoE Engineering Research Center of Software/Hardware Co-design Technology
and Application, East China Normal University, China<EMAIL_ADDRESS>and  Nengkun Yu  Centre for Quantum Software and Information, Faculty of
Engineering and Information Technology, University of Technology, Sydney,
Australia<EMAIL_ADDRESS>
###### Abstract.

Verifying quantum systems has attracted a lot of interest in the last decades.
In this paper, we study the quantitative model-checking of quantum continuous-
time Markov chains (quantum CTMCs). The branching-time properties of quantum
CTMCs are specified by continuous stochastic logic (CSL), which is famous for
verifying real-time systems, including classical CTMCs. The core of checking
the CSL formulas lies in tackling multiphase until formulas. We develop an
algebraic method using proper projection, matrix exponentiation, and definite
integration to symbolically calculate the probability measures of path
formulas. Thus the decidability of CSL is established. To be efficient,
numerical methods are incorporated to guarantee that the time complexity is
polynomial in the encoding size of the input model and linear in the size of
the input formula. A running example of Apollonian networks is further
provided to demonstrate our method.

###### Key words and phrases:

model checking, quantitative reasoning, formal logic, quantum computing

##  1\. Introduction

Rapid development has taken place on quantum hardware in very recent years. In
2020, the quantum computer Jiuzhang implemented a type of Boson sampling on 76
photonic qubits, in which case the quantum computer spent less than 20 seconds
while a classical supercomputer would require 600 million years [ZWD+20]. One
year later, IBM unveiled its 127-qubit quantum processor Eagle, which broke
the 100-qubit processor barrier [CDG21].

In the meantime, quantum software will be crucial in implementing quantum
algorithms and harnessing the power of quantum computers, e.g. Shor’s
algorithm with an exponential speed-up for integer factorization [Sho94],
Grover’s algorithm with a quadratic speed-up for unstructured search [Gro96],
and the HHL algorithm with an exponential speed-up for solving proper (sparse
and well-conditioned) systems of linear equations [HHL09]. The first practical
quantum programming language QCL appeared in Ömer’s work [Öme98]. The quantum
guarded command language (qGCL) was presented to program a “universal” quantum
computer [SZ00]. Selinger [Sel04] proposed functional programming languages
QFC with high-level features. Detailed survey on programming languages can be
found in [GLR+13, VRS+15]. To ensure reliability, all kinds of verification
techniques [DP06, YF21] need to be developed for various quantum algorithms
and protocols.

Model-checking [CGP99, BK08] is one of the most successful techniques for
verifying classical hardware and software systems. There are also many
explorations on applying model-checking to discrete-time quantum systems
[YF21], such as quantum programs and quantum (discrete-time) Markov chains.
However, as suggested by Feynman [Fey86] that continuous-time operations are
physically more realistic than discrete-time operations, it is worth studying
continuous-time models in quantum computing. In this paper, we consider a
novel quantum continuous-time Markov chain (quantum CTMC). It is a composite
system consisting of a classical subsystem and a quantum subsystem. In a
classical CTMC, the instantaneous descriptions (IDs) are given by probability
distributions over a finite set $S$ of classical states. Whereas, in a quantum
CTMC, the IDs are given by positive semi-definite matrices, which bear quantum
information, distributed over those classical states $S$. The ID’s evolution
of a particular class of quantum CTMCs that have only one classical state is a
_closed_ system in which there are no other classical states (its environment)
to interact with. In comparison, the ID’s evolution of a quantum CTMC with
more than one classical state is an _open_ system. Let us imagine a swinging
pendulum in a mechanical clock. In an ideal situation, the pendulum does not
suffer any unwanted interaction, such as air friction, from its environment,
which is a closed system. In practice, however, air friction indeed exists,
which would eventually make the pendulum cease. As interactions happen between
the pendulum and its environment, it is an open system. It is evident that
open systems are more common in the real world, which motivates us to study
the open system — quantum CTMC. The dynamical system of quantum CTMC can be
characterized by the Lindblad master equation:

$\frac{\mathrm{d}\,\rho(t)}{\mathrm{d}\,t}=\mathcal{L}(\rho(t))$ (1)

where $\rho(t)$ is the ID of the system at time $t$, and $\mathcal{L}$ is a
linear function of $\rho(t)$ (to be described in Subsection 2.2).

Under the model of quantum CTMC, we can develop the cylinder set of paths with
a computable probability measure, which is obtained by proper projection on
the ID $\rho$ in Eq. (1) for ruling out dissatisfying paths and matrix
exponentiation of the linear function $\mathcal{L}$ for computing the
probability distribution as time goes by. Thereby we establish the probability
space to do the formal quantitative reasoning. As a kind of real-time system,
the branching-time properties are popularly specified by continuous stochastic
logic (CSL). The core of checking the CSL formulas lies in tackling the
multiphase until formula. It requires the temporal property that any
satisfying path should meet a sequence of state constraints in turn, and
requires the timing constraint that the switch time-stamps should be in some
appointed time intervals. Specifically, the switch time-stamps are in the
absolute timing that are measured from the start. It guides us to design a
nested definite integral concerning the time variables interpreted as the
switch time-stamps to collect the measure of the candidate paths phase by
phase. By the mathematical tools of projection, matrix exponentiation and
definite integration, we can decide the CSL formula with time complexity being
polynomial in the encoding size of the input quantum CTMC and linear in the
size of the input CSL formula, besides the complexity of the query — how
sufficiently the Euler constant $\mathrm{e}$ should be approached. For the
consideration of practical efficiency, we also incorporate numerical methods —
the scaling and squaring method for matrix exponentiation and the Riemann sum
for definite integration — to yield a polynomial-time procedure for the
satisfaction probability. Furthermore, we provide a running example of
Apollonian networks to demonstrate our method, and show the usefulness of our
method via checking the percolation performance of those networks.

###  1.1. Related Work

#### Checking on CTMC

The IDs for a continuous-time Markov chain (CTMC) are represented by
probability distributions. Aziz et al. initialized the model-checking on CTMC
in 1996 [ASSB96]. They introduced the continuous stochastic logic (CSL) to
specify temporal properties of the CTMC. Roughly speaking, the syntax of CSL
amounts to that of computation tree logic (CTL) plus the multiphase until
formula
$\Phi_{0}\mathrm{U}\,^{\mathcal{I}_{0}}\Phi_{1}\mathrm{U}\,^{\mathcal{I}_{1}}\Phi_{2}\cdots\mathrm{U}\,^{\mathcal{I}_{K-1}}\Phi_{K}$
and the probability quantifier formula $\Pr_{>\texttt{c}}(\,\cdot\,)$ defined
on those IDs. They proved the decidability of CSL by number-theoretic
analysis. An approximate model-checking algorithm for a reduced version of CSL
was provided by Baier et al. [BKH99], in which multiphase until formulas are
restricted to binary until formulas
$\Phi_{1}\mathrm{U}\,^{\mathcal{I}}\Phi_{2}$. Under this logic, they applied
the efficient numerical technique — _uniformization_ — for _transient
analysis_ [BHHK00]. The approximate algorithms have been extended for
multiphase until formulas using _stratification_ [ZJNH11, ZJNH12]. After that,
Xu et al. considered the multiphase until formulas over the CTMC equipped with
a cost/reward structure [XZJ+16]. An algebraic algorithm was proposed to
attack this problem, whose effectiveness is ensured by number-theoretic
results and algebraic manipulation. Most of the above algorithms have been
implemented in probabilistic model checkers, like PRISM [KNP11], Storm
[DJKV17] and EPMC [FHL+22]. However, the IDs of the quantum CTMC are
represented by positive semi-definite matrices, which generalize the diagonal
matrices encoded for probability distributions. It leads to the fact that
quantum CTMCs generalize classical CTMCs, so new formalism must be developed
for this situation, as done in the current paper.

#### Checking on quantum DTMC

A quantum discrete-time Markov chain (quantum DTMC) is a composite model on
classical state space (a finite set) and quantum state space (a continuum), on
which the evolution is discrete-time given by quantum operations [LF15]. Gay
et al. [GNP08] restricted the quantum operations to Clifford group gates
(i.e., Hadamard, CNOT and the phase gates) and the whole state space as a
finite set of describable states, named stabilizer states, that are closed
under those Clifford group gates. They applied the model checker PRISM to
verify the protocols of quantum superdense coding, quantum teleportation, and
quantum error correction. Whereas, Feng et al. proposed Markov chain with
transitions given by the general quantum operations [FYY13]. Under the model,
the authors considered the reachability probability [YFYY13], the repeated
reachability probability [FHTY17], and the model-checking of a quantum analogy
of computation tree logic [FYY13]. A key step in their work is decomposing the
whole state space (known as a Hilbert space) into a direct-sum of some bottom
strongly connected component (BSCC) subspaces plus a maximal transient
subspace concerning a given quantum operation describing the evolution of the
quantum DTMC [YFYY13, GFY18]. After the decomposition, all the problems above
were shown to be computable/decidable in polynomial time concerning the
dimension of the Hilbert space. In contrast to the quantum operations in
discrete-time evolution, when considering quantum CTMC, the general linear
operators should be admitted to describe the rates of state transitions in the
continuous-time evolution. For this purpose, we will characterize the
continuous-time evolution by the Lindblad master equation in the present work.

#### Linear-time logic vs branching-time logic

All of the above literature are concerned with branching-time logic that
classifies all paths into satisfying ones and dissatisfying ones, and compares
the measure of satisfying paths with a predefined threshold. Fewer works
concerning linear-time logics impose constraints on the measure of all paths
at some critical time periods. Guan and Yu introduced continuous linear logic
(CLL) over CTMCs [GY22], and proved the decidability by real root isolation of
exponential polynomials. Xu et al. initialized the model-checking on quantum
CTMCs, in which the decidability of signal temporal logic (STL) was
established, again, by real root isolation [XMGY21]. Both CLL and STL are
linear-time logics. However, checking the more popular branching-time logic
over quantum CTMCs, as far as we know, remains open.

#### Contribution

The contributions of the present paper can be summarised as:

  1. (1)

We consider the novel model of quantum CTMC that has a more precise
description of a composite system — consisting of a classical subsystem and a
quantum subsystem — than the existing work [XMGY21].

  2. (2)

We prove the decidability of CSL over quantum CTMCs and give a numerical
approach in polynomial time to compute the probability measure.

  3. (3)

A running example of Apollonian networks is provided to show the utility of
our method.

#### Organization

The rest of this paper is organized as follows. Section 2 recalls some notions
and notations from quantum computing and number theory. In Section 3 we
introduce the model of quantum CTMC, and establish the probability space.
Based on that, we will check the CSL formula for specifying branching-time
properties over quantum CTMCs in Section 4. We further prove the decidability
of CSL and give an efficient numerical approach. Section 5 is the conclusion.
The implementation is delivered in Appendix A.

##  2\. Preliminaries

###  2.1. Quantum Computing

Here we recall some notions and notations from quantum computing [NC00], which
will be widely used in this paper. Let $\mathcal{H}$ be a Hilbert space with
finite dimension $d$ that is a vector space over complex numbers $\mathbb{C}$
equipped with an inner product. We would employ the Dirac notations as
follows:

  * $\triangleright$

$\ket{\psi}$ is a unit column vector labelled with $\psi$;

  * $\triangleright$

$\bra{\psi}\coloneqq\ket{\psi}^{\dagger}$ is the Hermitian adjoint (transpose
and entry-wise complex conjugate) of $\ket{\psi}$;

  * $\triangleright$

$\ip{\psi_{1}}{\psi_{2}}\coloneqq\bra{\psi_{1}}\ket{\psi_{2}}$ is the inner
product of $\ket{\psi_{1}}$ and $\ket{\psi_{2}}$;

  * $\triangleright$

$\op{\psi_{1}}{\psi_{2}}\coloneqq\ket{\psi_{1}}\otimes\bra{\psi_{2}}$ is the
outer product where $\otimes$ denotes tensor product;

  * $\triangleright$

$\ket{\psi_{1},\psi_{2}}\coloneqq\ket{\psi_{1}}\ket{\psi_{2}}$ is a shorthand
of the tensor product $\ket{\psi_{1}}\otimes\ket{\psi_{2}}$.

Fixing an orthonormal basis $\\{\ket{\psi_{i}}\colon 1\leq i\leq d\\}$ of
$\mathcal{H}$, $\ket{\psi}\in\mathcal{H}$ can be linearly expressed as
$\sum_{i=1}^{d}c_{i}\ket{\psi_{i}}$ with $c_{i}\in\mathbb{C}$ and
$\sum_{i=1}^{d}|c_{i}|^{2}=1$. This $\ket{\psi}$ is called a _superposition_
if it has two or more nonzero terms $c_{i}\ket{\psi_{i}}$ under that basis.
For example, the most familiar elements
$\ket{\pm}\coloneqq(\ket{0}\pm\ket{1})/\sqrt{2}$ are superpositions over the
basis $\\{\ket{0},\ket{1}\\}$.

#### Linear operators

We will mainly consider the linear operators on $\mathcal{H}$. Here such a
parameter $\mathcal{H}$ can be omitted when it is clear from the context. A
linear operator $\gamma$ is _Hermitian_ if $\gamma=\gamma^{\dagger}$; it is
_positive_ if $\bra{\psi}\gamma\ket{\psi}\geq 0$ holds for all
$\ket{\psi}\in\mathcal{H}$. A _projector_ $\mathbf{P}$ is a positive operator
of the form $\sum_{i=1}^{k}\op{\psi_{i}}{\psi_{i}}$ ($k\leq d$) for some
orthonormal elements $\ket{\psi_{i}}\in\mathcal{H}$ ($1\leq i\leq k$).
Clearly, there is a bijective map between projectors
$\mathbf{P}=\sum_{i=1}^{k}\op{\psi_{i}}{\psi_{i}}$ and subspaces of
$\mathcal{H}$ that are spanned by $\\{\ket{\psi_{i}}\colon 1\leq i\leq k\\}$.
In sum, positive operators are Hermitian operators whose eigenvalues are all
nonnegative; projectors are positive operators whose all nonzero eigenvalues
are $1$. The _identity operator_ $\mathbf{I}$ is the linear operator
$\sum_{i=1}^{d}\op{\psi_{i}}{\psi_{i}}$ for some orthonormal basis
$\\{\ket{\psi_{i}}\colon 1\leq i\leq d\\}$ of $\mathcal{H}$. A linear operator
$\mathbf{U}$ is _unitary_ if
$\mathbf{U}\mathbf{U}^{\dagger}=\mathbf{U}^{\dagger}\mathbf{U}=\mathbf{I}$.
The _trace_ of a linear operator $\gamma$ is defined as
$\mathrm{tr}(\gamma)\coloneqq\sum_{i=1}^{d}\bra{\psi_{i}}\gamma\ket{\psi_{i}}$
for some orthonormal basis $\\{\ket{\psi_{i}}\colon 1\leq i\leq d\\}$ of
$\mathcal{H}$. It is worth noting that the trace function is actually
independent of the orthonormal basis selected. For a linear operator
$\gamma_{A,B}$ on the composite Hilbert space
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$, we can define the _partial traces_
as:

$\displaystyle\mathrm{tr}_{A}(\gamma_{A,B})$
$\displaystyle\coloneqq\sum_{i=1}^{d_{A}}(\bra{\psi_{i}^{A}}\otimes\mathbf{I}_{B})\gamma_{A,B}(\ket{\psi_{i}^{A}}\otimes\mathbf{I}_{B})$
(2) $\displaystyle\mathrm{tr}_{B}(\gamma_{A,B})$
$\displaystyle\coloneqq\sum_{i=1}^{d_{B}}(\mathbf{I}_{A}\otimes\bra{\psi_{i}^{B}})\gamma_{A,B}(\mathbf{I}_{A}\otimes\ket{\psi_{i}^{B}}),$

where $\mathbf{I}_{A}$ (resp. $\mathbf{I}_{B}$) are the identity operators on
$\mathcal{H}_{A}$ (resp. $\mathcal{H}_{B}$), $d_{A}$ (resp. $d_{B}$) are the
dimensions of $\mathcal{H}_{A}$ (resp. $\mathcal{H}_{B}$), and
$\\{\ket{\psi_{i}^{A}}\colon 1\leq i\leq d_{A}\\}$ (resp.
$\\{\ket{\psi_{i}^{B}}\colon 1\leq i\leq d_{B}\\}$) are some orthonormal basis
of $\mathcal{H}_{A}$ (resp. $\mathcal{H}_{B}$).

#### Quantum states

The quantum states are given by the form of probabilistic ensembles
$\\{(p_{i},\ket{\psi_{i}})\colon 1\leq i\leq k\\}$ with $p_{i}>0$ and
$\sum_{i=1}^{k}p_{i}=1$. An alternative representation is using positive
operators $\rho=\sum_{i=1}^{k}p_{i}\op{\psi_{i}}{\psi_{i}}$, named _density
operators_ , whose traces are unit as
$\mathrm{tr}(\rho)=\sum_{i=1}^{k}p_{i}\mathrm{tr}(\op{\psi_{i}}{\psi_{i}})=1$.
Here those $\ket{\psi_{i}}$ are not necessary to be orthonormal. To be more
explicit, we resort to the spectral decomposition [NC00, Box 2.2] that
$\rho=\sum_{i=1}^{d}\lambda_{i}\op{\lambda_{i}}{\lambda_{i}}$, where
$\ket{\lambda_{i}}$ are eigenvectors interpreted as the _eigenstates_ of
$\rho$ and $\lambda_{i}$ are eigenvalues interpreted as the _probabilities_ of
taking the eigenstates $\ket{\lambda_{i}}$. It is worth noting that the
spectral decomposition of $\rho$ is not unique, but the number of nonzero
eigenvalues is unique since it is exactly the rank of $\rho$. When there is
only one eigenstate with positive probability, $\rho$ is said to be a _pure_
state; otherwise it is a _mixed_ state. In other words, a pure state indicates
the system state which we completely know; a mixed state gives all possible
system states with a total probability of $1$. We denote by $\mathcal{D}$ the
set of density operators, and by $\mathcal{D}^{\leq 1}$ the set of partial
density operators that are positive operators with trace bounded by $1$.
Sometimes, we need to extract local information from a composite system. It
can be achieved by the partial traces mentioned in Eq. (2), which result in
the _reduced_ density operators $\rho_{A}=\mathrm{tr}_{B}(\rho_{A,B})$ on
$\mathcal{H}_{A}$ and $\rho_{B}=\mathrm{tr}_{A}(\rho_{A,B})$ on
$\mathcal{H}_{B}$, respectively. For example, considering the Bell state
$\rho_{A,B}=\tfrac{1}{2}(\ket{0,0}+\ket{1,1})(\bra{0,0}+\bra{1,1})$, we can
get the reduced density operators
$\rho_{A}=\mathrm{tr}_{B}(\rho_{A,B})=\tfrac{1}{2}(\op{0}{0}+\op{1}{1})$ and
$\rho_{B}=\mathrm{tr}_{A}(\rho_{A,B})=\tfrac{1}{2}(\op{0}{0}+\op{1}{1})$. By
the partial trace, the pure state $\rho_{A,B}$ on
$\mathcal{H}_{A}\otimes\mathcal{H}_{B}$ would be mixed states $\rho_{A}$ on
$\mathcal{H}_{A}$ and $\rho_{B}$ on $\mathcal{H}_{B}$, respectively. It also
reveals that the general quantum operations are non-unitary.

#### State transitions

A _closed_ system is an ideal system that does not suffer from any unwanted
interaction from outside environment. State transitions in a closed system are
characterized by some unitary transformation $\ket{\psi(t+\Delta
t)}=\mathbf{U}(\Delta t)\ket{\psi(t)}$, where $\ket{\psi(t)}$ is the state of
the system at time $t$, and $\mathbf{U}(\Delta t)=\exp(-\imath\mathbf{H}\Delta
t)$ comes from the _Hamiltonian_ (a Hermitian operator) $\mathbf{H}$ in the
Schrödinger equation:

$\frac{\mathrm{d}\,\ket{\psi(t)}}{\mathrm{d}\,t}=-\imath\mathbf{H}\ket{\psi(t)}.$
(3)

An _open_ system interacts with its environment. Composed with the
environment, it would be a larger closed system. In details, for a mixed state
$\rho=\sum_{i=1}^{k}p_{i}\op{\psi_{i}}{\psi_{i}}$ in an open system, there is
a pure state in the larger closed system

$\left(\sum_{i=1}^{k}\sqrt{p_{i}}\ket{\psi_{i},env_{i}}\right)\left(\sum_{i=1}^{k}\sqrt{p_{i}}\bra{\psi_{i},env_{i}}\right)$

for some orthonormal elements $\ket{env_{i}}$ of the environment. State
transitions in the open system can be obtained by tracing out the environment
in the larger closed system, which leads to the quantum operation
$\rho(t+\Delta t)=\sum_{j=0}^{m}\mathbf{E}_{j}(\Delta
t)\rho(t)\mathbf{E}_{j}^{\dagger}(\Delta t)$, where $\mathbf{E}_{j}$ are
linear operators satisfying the trace preservation condition
$\sum_{j=0}^{m}\mathbf{E}_{j}^{\dagger}\mathbf{E}_{j}=\mathbf{I}$. Here those
linear operators $\mathbf{E}_{j}$ come from the Lindblad master equation to be
described in the coming subsection.

###  2.2. Lindblad Master Equation

The quantum continuous-time Markov chain considered in this paper is a kind of
open system satisfying the _memoryless_ property, i.e., the dynamic system
depends only on the current state regardless of the history. To characterize
such state evolution, we employ the Lindblad master equation [Lin76, GKS76].
It is of the following form

$\displaystyle\rho^{\prime}$
$\displaystyle=-\imath[\mathbf{H},\rho]+\sum_{j=1}^{m}\left(\mathbf{L}_{j}\rho\mathbf{L}_{j}^{\dagger}-\tfrac{1}{2}\\{\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j},\rho\\}\right)$
(4)
$\displaystyle=-\imath\mathbf{H}\rho+\imath\rho\mathbf{H}+\sum_{j=1}^{m}\left(\mathbf{L}_{j}\rho\mathbf{L}_{j}^{\dagger}-\tfrac{1}{2}\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j}\rho-\tfrac{1}{2}\rho\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j}\right),$

where $\mathbf{H}$ is a Hermitian operator, $\mathbf{L}_{j}$ are linear
operators,
$[\mathbf{A},\mathbf{B}]\coloneqq\mathbf{A}\mathbf{B}-\mathbf{B}\mathbf{A}$ is
the _commutator_ between two linear operators $\mathbf{A}$ and $\mathbf{B}$,
and
$\\{\mathbf{A},\mathbf{B}\\}\coloneqq\mathbf{A}\mathbf{B}+\mathbf{B}\mathbf{A}$
is the _anti-commutator_. The term $-\imath[\mathbf{H},\rho]$ describes the
evolution of the internal system; the term
$\mathbf{L}_{j}\rho\mathbf{L}_{j}^{\dagger}-\tfrac{1}{2}\\{\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j},\rho\\}$
represents the interaction between the system and its environment. In other
words, to characterize the evolution of an open system, it is necessary to use
those linear operators $\mathbf{L}_{j}$ besides the Hermitian operator
$\mathbf{H}$.

Considering the density operator $\rho$ at time $t$ given by Eq. (4), we can
get $\rho(t+\Delta t)=\sum_{j=0}^{m}\mathbf{E}_{j}(\Delta
t)\rho(t)\mathbf{E}_{j}^{\dagger}(\Delta t)$ where

$\mathbf{E}_{0}(\Delta
t)=\mathbf{I}-\left(\frac{1}{2}\sum_{j=1}^{m}\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j}+\imath\mathbf{H}\right)\Delta
t\qquad\mbox{and}\qquad\mathbf{E}_{j}(\Delta t)=\mathbf{L}_{j}\sqrt{\Delta
t}\quad\mbox{for }j>0,$

which satisfies the trace preservation condition as $\lim_{\Delta t\to
0}\sum_{j=0}^{m}\mathbf{E}_{j}^{\dagger}(\Delta t)\mathbf{E}_{j}(\Delta
t)=\mathbf{I}$. Conversely, for a given quantum operation, we can construct a
Lindblad master equation induced by it. It entails that the Lindblad master
equation is the most general type of Markovian and time-homogeneous master
equation describing (generally non-unitary) state evolution that preserves the
laws of quantum mechanics, i.e., complete positivity and trace-preservation.

We now turn to derive the solution of Eq. (4). Two useful functions are
defined as:

  * $\triangleright$

$\mathrm{L2V}(\gamma)\coloneqq\sum_{i=1}^{d}\sum_{j=1}^{d}\bra{i}\gamma\ket{j}\ket{i,j}$
that rearranges entries of the linear operator $\gamma$ on the Hilbert space
$\mathcal{H}$ with dimension $d$ as a column vector; and

  * $\triangleright$

$\mathrm{V2L}(\mathbf{v})\coloneqq\sum_{i=1}^{d}\sum_{j=1}^{d}\bra{i,j}\mathbf{v}\op{i}{j}$
that rearranges entries of the column vector $\mathbf{v}$ as a linear
operator.

Here, $\mathrm{L2V}$ and $\mathrm{V2L}$ are pronounced “linear operator to
vector” and “vector to linear operator”, respectively. They are mutually
inverse functions, so that if a linear operator (resp. its vectorization) is
determined, its vectorization (resp. the original linear operator) is
determined. Hence, we can freely choose one of the two representations for
convenience. Using the following transformation that for any linear operators
$\mathbf{A},\mathbf{B},\mathbf{C}$,

$\mathbf{D}=\mathbf{A}\mathbf{B}\mathbf{C}\Longleftrightarrow\mathrm{L2V}(\mathbf{D})=(\mathbf{A}\otimes\mathbf{C}^{\mathrm{T}})\mathrm{L2V}(\mathbf{B})$

where $\mathrm{T}$ denotes transpose, we can reformulate Eq. (4) as the linear
ordinary differential equation

$\mathrm{L2V}(\rho^{\prime})=\mathds{M}\cdot\mathrm{L2V}(\rho),$ (5)

where
$\mathds{M}=-\imath\mathbf{H}\otimes\mathbf{I}+\imath\mathbf{I}\otimes\mathbf{H}^{\mathrm{T}}+\sum_{j=1}^{m}\big{(}\mathbf{L}_{j}\otimes\mathbf{L}_{j}^{*}-\tfrac{1}{2}\mathbf{L}_{j}^{\dagger}\mathbf{L}_{j}\otimes\mathbf{I}-\tfrac{1}{2}\mathbf{I}\otimes\mathbf{L}_{j}^{\mathrm{T}}\mathbf{L}_{j}^{*}\big{)}$
is the _governing matrix_ with $*$ denoting entry-wise complex conjugate. As a
result, we get the desired solution
$\mathrm{L2V}(\rho(t))=\exp(\mathds{M}\cdot t)\cdot\mathrm{L2V}(\rho(0))$ or
equivalently $\rho(t)=\mathrm{V2L}(\exp(\mathds{M}\cdot
t)\cdot\mathrm{L2V}(\rho(0)))$ in a closed form. It is obtained in polynomial
time by the standard method [Kai80].

###  2.3. Number Theory

We will show the decidability of a temporal logic, which is based on the
following essential facts in number theory. A number $\alpha$ is _algebraic_ ,
denoted by $\alpha\in\mathbb{A}$, if there is a $\mathbb{Q}$-polynomial
$f_{\alpha}(z)$ of least positive degree, satisfying $f_{\alpha}(\alpha)=0$;
otherwise $\alpha$ is _transcendental_. In the above definition, such a
polynomial $f_{\alpha}(z)$ is called the _minimal polynomial_ of $\alpha$. The
_degree_ $D$ of $\alpha$ is $\deg_{z}(f_{\alpha})$. The standard encoding of
$\alpha$ is the minimal polynomial $f_{\alpha}$ plus an isolation disk in the
complex plane that distinguishes $\alpha$ from other roots of $f_{\alpha}$. So
the encoding size $\|\alpha\|$ is the number of bits used to store the minimal
polynomial and the isolation disk.

Let $\mu_{1},\ldots,\mu_{m}$ be irrational complex numbers. Then the _field
extension_ $\mathbb{Q}(\mu_{1},\ldots,\mu_{m}):\mathbb{Q}$ is the smallest set
that contains $\mu_{1},\ldots,\mu_{m}$ and is closed under arithmetic
operations, i.e. addition, subtraction, multiplication and division. Here
those irrational complex numbers $\mu_{1},\ldots,\mu_{m}$ are called the
generators of the field extension. A field extension is _simple_ if it has
only one generator. For instance, the field extension
$\mathbb{Q}(\sqrt{2}):\mathbb{Q}$ is exactly the set $\\{a+b\sqrt{2}\colon
a,b\in\mathbb{Q}\\}$.

[[Loo83, Algorithm 2]] Let $\alpha_{1}$ and $\alpha_{2}$ be two algebraic
numbers of degrees $D_{1}$ and $D_{2}$, respectively. There is an algebraic
number $\mu$ of degree at most $D_{1}D_{2}$, such that the field extension
$\mathbb{Q}(\mu):\mathbb{Q}$ is exactly
$\mathbb{Q}(\alpha_{1},\alpha_{2}):\mathbb{Q}$. For a collection of algebraic
numbers $\alpha_{1},\ldots,\alpha_{m}$ appearing in the input instance, by
repeatedly applying Lemma 2.3, we can obtain a simple field extension
$\mathbb{Q}(\mu):\mathbb{Q}$ that can generate all
$\alpha_{1},\ldots,\alpha_{m}$.

[[Coh93, Corollary 4.1.5]] Let $\alpha$ be an algebraic number of degree $D$,
and $g(z)$ an $\mathbb{A}$-polynomial with degree $D_{g}$ and coefficients
taken from $\mathbb{Q}(\alpha):\mathbb{Q}$. There is a $\mathbb{Q}$-polynomial
$f(z)$ of degree at most $DD_{g}$, such that roots of $g(z)$ are those of
$f(z)$. The above lemma entails that roots of any $\mathbb{A}$-polynomial are
also algebraic.

[(Lindemann 1882) [Bak75, Theorem 1.4]] For any nonzero algebraic numbers
$\alpha_{1},\ldots,\alpha_{m}$ and any distinct algebraic numbers
$\lambda_{1},\ldots,\lambda_{m}$, the sum
$\sum_{i=1}^{m}\alpha_{i}\mathrm{e}^{\lambda_{i}}$ with $m\geq 1$ is nonzero.
Hence the sign of a given real number with the form
$\sum_{i=1}^{m}\alpha_{i}\mathrm{e}^{\lambda_{i}}$ is decidable, which can be
achieved by sufficiently approaching the Euler constant $\mathrm{e}$. But the
complexity of the procedure is unknown as in the existing literature [COW16,
HLXL18]. It becomes a query (also known as oracle) in our work.

##  3\. Quantum Continuous-Time Markov Chain

In this section, we first introduce the model of quantum continuous-time
Markov chain, which subsumes the classical continuous-time Markov chain. Based
on that, we will establish a probability space over paths to reason about
real-time properties. An algebraic approach using matrix exponentiation is
presented to calculate those measurable events in the probability space.

Let $AP$ be a set of atomic propositions throughout this paper. A labelled
_quantum continuous-time Markov chain_ (quantum CTMC for short) $\mathfrak{Q}$
over the Hilbert space $\mathcal{H}$ is a triple $(S,Q,L)$, in which

  * •

$S$ is a finite set of classical states;

  * •

$Q$ is the transition generator function given by

    * –

a Hermitian operator $\mathbf{H}$ with the form $\sum_{s\in
S}\op{s}{s}\otimes\mathbf{H}_{s}$ where $\mathbf{H}_{s}$ are Hermitian
operators on $\mathcal{H}$,

    * –

a finite set of linear operators $\mathbf{L}_{j}$ with the form
$\op{s_{2}}{s_{1}}\otimes\mathbf{L}_{j,s_{1},s_{2}}$ for some $s_{1},s_{2}\in
S$, $s_{1}\neq s_{2}$, where $\mathbf{L}_{j,s_{1},s_{2}}$ are linear operators
on $\mathcal{H}$; and

  * •

$L\colon S\to 2^{AP}$ is a labelling function.

Usually, a density operator of the form $\rho(0)=\sum_{s\in
S}\op{s}{s}\otimes\rho_{s}$ is appointed as the initial instantaneous
description (ID) of $\mathfrak{Q}$, where
$\rho_{s}\in\mathcal{D}_{\mathcal{H}}^{\leq 1}$ satisfy $\sum_{s\in
S}\mathrm{tr}(\rho_{s})=1$.

Every classical state $s\in S$ can be quantized by an orthonormal state
$\ket{s}$. Let $\mathcal{C}\coloneqq\mathrm{span}(\\{\ket{s}\colon s\in S\\})$
be the Hilbert space corresponding to the classical system, and
$\mathcal{H}_{\mathrm{cq}}\coloneqq\mathcal{C}\otimes\mathcal{H}$ the enlarged
Hilbert space corresponding to the whole classical–quantum system. The
dimension of $\mathcal{H}_{\mathrm{cq}}$ is $N\coloneqq nd$ where $n=|S|$ and
$d=\dim(\mathcal{H})$. This parameter $N$ is used to reflect the size of the
quantum CTMC $\mathfrak{Q}$.

The IDs $\rho$ for a quantum CTMC $\mathfrak{Q}$ are represented by density
operators on the enlarged Hilbert space $\mathcal{H}_{\mathrm{cq}}$ with the
mixed structure $\sum_{s\in S}\op{s}{s}\otimes\rho_{s}$ where
$\rho_{s}\in\mathcal{D}_{\mathcal{H}}^{\leq 1}$ are partial density operators,
satisfying $\sum_{s\in S}\mathrm{tr}(\rho_{s})=1$. In mathematical, such a
mixed structure entails that the entries
$(\bra{s_{1}}\otimes\mathbf{I})\rho(\ket{s_{2}}\otimes\mathbf{I})$ are zero
whenever $s_{1}\neq s_{2}$. It underlies the key difference between the
classical state space $\mathcal{C}$ and the quantum state space $\mathcal{H}$
that:

  * $\star$

The classical state space $\mathcal{C}$ has finitely many base states
$\ket{s}$ with $s\in S$, which are known and fixed in prior. After tracing out
the quantum state space $\mathcal{H}$, the reduced density operator on
$\mathcal{C}$ is a probability distribution over these classical states, and
can not be a _superposition_.

  * $\star$

The quantum state system $\mathcal{H}$ has uncountably many base states
$\ket{\psi}$, almost all of which are superpositions on the orthonormal basis
of $\mathcal{H}$, due to the postulate of quantum mechanics [NC00] that the
space of quantum base states is a _continuum_. It would bring extra technical
hardness to verifying quantum models.

We can also regard these IDs $\rho$ as positive-operator valued distributions
over classical states $S$, which generalize probability distributions
representing the IDs of a classical CTMC. Under the mixed structure, the
current model of quantum CTMC is more precise than the one proposed in
[XMGY21] that does _not_ clarify the difference between the classical state
space and the quantum state space.

The transition generator function $Q$ is functionally similar to the
transition rate matrix in the classical CTMC, but employs the Lindblad master
equation. It gives rise to a _universal_ way to describe the behavior of a
quantum CTMC that keeps the ID $\rho$ in the mixed structure $\sum_{s\in
S}\op{s}{s}\otimes\rho_{s}$, following from the generality of the Lindblad
master equation. In more detail, each $\mathbf{H}_{s}$ determines the internal
evolution at classical state $s$ and thereby $\mathbf{H}$ gives all such
individual internal evolution, while
$\mathbf{L}_{j}=\op{s_{2}}{s_{1}}\otimes\mathbf{L}_{j,s_{1},s_{2}}$ determines
the external evolution from classical state $s_{1}$ to $s_{2}$. It is notable
that, as required in Eq. (4), $\mathbf{H}_{s}$ can be chosen as _arbitrary_
Hermitian operators and $\mathbf{L}_{j,s_{1},s_{2}}$ can be _arbitrary_ linear
operators, which makes the model much expressive. A state $s$ is _absorbing_
if all transitions from it are removed, which is achieved by some projectors
to be defined later.

In fact, the model of quantum CTMC $\mathfrak{Q}=(S,Q,L)$ extends that of
classical CTMC $\mathfrak{C}=(S,\mathbf{Q},L)$, where $\mathbf{Q}$ is the
transition rate matrix. It can be seen from the following lemma with a trivial
labelling function $L$: [[XMGY21, Lemma 10]] Given a CTMC $\mathfrak{Q}$, it
can be modelled by a quantum CTMC $\mathfrak{C}$ over one-dimensional Hilbert
space $\mathcal{H}$.

For the sake of computability, the entries of $\rho_{s}$, $\mathbf{H}_{s}$ and
$\mathbf{L}_{j,s_{1},s_{2}}$ in the model of quantum CTMC are supposed to be
algebraic numbers in the simple field extension $\mathbb{Q}(\mu):\mathbb{Q}$
for some algebraic number $\mu$. The encoding size $\|\mu\|$ is used to
reflect the encoding size of operands in $\mathfrak{Q}$. Together with the
dimension $N$ of $\mathcal{H}_{\mathrm{cq}}$, they are two key parameters in
the encoding size $\|\mathfrak{Q}\|$, which will be used to analyze the
complexity of our method afterward.

Here we will study the open quantum walk (OQW) along Apollonian networks (ANs)
[AHAdS05]. ANs have small-world and scale-free properties. Small-world
property means that there is a relatively short path between two nodes despite
the often large size of networks, and scale-free property implies that most of
the nodes in the network are connected to only a few nodes while few nodes are
connected to a large number of nodes. With the two properties, ANs provide an
excellent facility to analyze dynamical processes taking place on networked
systems, including percolation and electrical conduction. In the meanwhile,
ANs serve as a valuable tool for studying properties of various walk models,
including (discrete- and continuous-time) random walks and quantum walks
[XLL08].

ANs are generated in such a way: At the initial generation, the network
consists of three nodes marked as 0, 1 and 2, respectively. At each subsequent
generation, a new node is added inside each triangle and linked to the three
vertices of that triangle, as shown in Figure 1.

$0$$1$$2$

(a) 0th generation

$0$$1$$2$$3$$\leftarrow A$$B\rightarrow$$\leftarrow
B$$A\rightarrow$$\leftarrow
A$$B\rightarrow$$A+C/\sqrt{3}\rightarrow$$\leftarrow C$$\leftarrow
C$$\tfrac{1}{\sqrt{3}}C\rightarrow$$\leftarrow B+C/\sqrt{3}$$C\rightarrow$

(b) 1st generation

$0$$1$$2$$3$$4$$5$$6$

(c) 2nd generation

Figure 1. The first three generations of ANs

Consider a qutrit, a quantum state in a 3-dimensional Hilbert space, is
walking along the first generation of AN [PGMS15], which is depicted in Figure
1(b). As the qutrit may interact with its environment, some noise from the
network channel will affect the walking of the qutrit which can be described
by linear operators, say, $A=\op{x}{x}$, $B=\op{y}{y}$ and $C=\op{z}{z}$,
where $\ket{x}=(1,1,1)^{\mathrm{T}}/\sqrt{3}$,
$\ket{y}=(1,w,w^{2})^{\mathrm{T}}/\sqrt{3}$ and
$\ket{z}=(1,w^{2},w)^{\mathrm{T}}/\sqrt{3}$ with $w=\exp(2\pi\imath/3)$. The
OQW can be modelled by the quantum CTMC $\mathfrak{Q}_{1}=(S,Q,L)$ over the
Hilbert space $\mathcal{H}=\mathrm{span}(\\{\ket{0},\ket{1},\ket{2}\\})$,
where $S=\\{0,1,2,3\\}$ represents the four nodes; Node $3$ has a label
“center” and other nodes have none; the transition function $Q$ is given by
the Hermitian operator $\mathbf{H}=0$ and the following 12 linear operators:

$\displaystyle\mathbf{L}_{1}=\op{1}{0}\otimes A$
$\displaystyle\mathbf{L}_{2}=\op{2}{0}\otimes B$
$\displaystyle\mathbf{L}_{3}=\op{3}{0}\otimes C$
$\displaystyle\mathbf{L}_{4}=\op{2}{1}\otimes A$
$\displaystyle\mathbf{L}_{5}=\op{0}{1}\otimes B$
$\displaystyle\mathbf{L}_{6}=\op{3}{1}\otimes C$
$\displaystyle\mathbf{L}_{7}=\op{0}{2}\otimes A$
$\displaystyle\mathbf{L}_{8}=\op{1}{2}\otimes B$
$\displaystyle\mathbf{L}_{9}=\op{3}{2}\otimes C$
$\displaystyle\mathbf{L}_{10}=\op{0}{3}\otimes C/\sqrt{3}$
$\displaystyle\mathbf{L}_{11}=\op{1}{3}\otimes(B+C/\sqrt{3})$
$\displaystyle\mathbf{L}_{12}=\op{2}{3}\otimes(A+C/\sqrt{3}).$

An initial ID $\rho(0)$ could be, say,
$\op{3}{3}\otimes\tfrac{1}{3}\mathbf{I}$. Here, the left operand $\op{3}{3}$
in the tensor product means that we are definitely in Node $3$ in the
classical system, the right operand $\tfrac{1}{3}\mathbf{I}$ means that we are
in one of three base states $\ket{0}$ through $\ket{2}$ with equal probability
$\tfrac{1}{3}$ in the quantum system, and they together determine the ID in
the whole classical–quantum system. ∎

A path $\omega$ of a quantum CTMC $\mathfrak{Q}$ is an infinite sequence
$s_{0}\xrightarrow{\tau_{0}}s_{1}\xrightarrow{\tau_{1}}s_{2}\xrightarrow{\tau_{2}}\cdots$.
It means that $\omega$ starts at $s_{0}$, and sojourns there up to time
$t_{0}=\tau_{0}$ (but excluding $\tau_{0}$); then $\omega$ jumps to $s_{1}$ at
time $t_{0}$, and sojourns there up to time $t_{1}=\tau_{0}+\tau_{1}$; then
$\omega$ jumps to $s_{2}$ at time $t_{1}$, and sojourns there up to time
$t_{2}=t_{1}+\tau_{2}$; and so on. These time $\tau_{k}$ ($k\geq 0$) are
called _sojourn_ time that are measured from the instants when the prior
transitions $s_{k-1}\rightarrow s_{k}$ take place; while these time $t_{k}$
($k\geq 0$) are _switch_ time-stamps that are measured from the start of the
path. A finite path $\hat{\omega}$ is a finite fragment of such a path
$\omega$, e.g.
$s_{0}\xrightarrow{\tau_{0}}s_{1}\xrightarrow{\tau_{1}}s_{2}\xrightarrow{\tau_{2}}\cdots\xrightarrow{\tau_{K-1}}s_{K}$.
Namely, the state $\omega(t)$ of $\omega$ at time $t$ is $s_{k}$ for
$t\in[t_{k-1},t_{k})$, where $t_{-1}$ is defined to be zero for convenience.
Let $Path$ be the set of paths of $\mathfrak{Q}$.

To formally reason about quantitative properties of quantum CTMC, we will
establish the probability space over paths as follows: A _measurable space_ is
a pair $(\Omega,\Sigma)$, where $\Omega$ is a nonempty set and $\Sigma$ is a
$\sigma$-algebra on $\Omega$ that is a collection of subsets of $\Omega$,
satisfying:

  * •

$\Omega\in\Sigma$, and

  * •

$\Sigma$ is closed under countable union and complement.

In addition, a _probability space_ is a triple $(\Omega,\Sigma,\Pr)$, where
$(\Omega,\Sigma)$ is a measurable space and $\Pr\colon\Sigma\to[0,1]$ is a
probability measure, satisfying:

  * •

$\Pr(\Omega)=1$,

  * •

$\Pr(A)\geq 0$ for any $A\in\Sigma$, and

  * •

$\Pr(\biguplus_{i}A_{i})=\sum_{i}\Pr(A_{i})$ for any pairwise disjoint
$A_{i}\in\Sigma$.

Note that any finite path
$\hat{\omega}=s_{0}\xrightarrow{\tau_{0}}s_{1}\xrightarrow{\tau_{1}}s_{2}\xrightarrow{\tau_{2}}\cdots\xrightarrow{\tau_{K-1}}s_{K}$
is not measurable, since it requires the transitions from $s_{k-1}$ to $s_{k}$
($k>0$) take place instantaneously. To fix it, we use the _cylinder set_ of
paths, which is defined as

$\displaystyle
Cyl(s_{0}\xrightarrow{\mathcal{J}_{0}}s_{1}\xrightarrow{\mathcal{J}_{1}}s_{2}\xrightarrow{\mathcal{J}_{2}}\cdots\xrightarrow{\mathcal{J}_{K-1}}s_{K})\coloneqq$
(6) $\displaystyle\qquad\left\\{\omega\in
Path\colon\omega=s_{0}\xrightarrow{\tau_{0}}s_{1}\xrightarrow{\tau_{1}}s_{2}\xrightarrow{\tau_{2}}\cdots\xrightarrow{\tau_{K-1}}s_{K}\xrightarrow{\tau_{K}}\cdots\wedge\bigwedge_{k=0}^{K-1}\tau_{k}\in\mathcal{J}_{k}\right\\}.$

Here, these time intervals $\mathcal{J}_{k}$ ($0\leq k<K$) are _relatively_
timing, as they are measured from the instants when the prior transitions
$s_{k-1}\rightarrow s_{k}$ take place. Let $\Omega=Path$, and $\Pi\subseteq
2^{\Omega}$ be countably many cylinder sets $Cyl$ of paths plus the empty set
$\emptyset$. By [BK08, Chapter 10], there is a smallest $\sigma$-algebra
$\Sigma$ of $\Pi$ that contains $\Pi$ and is closed under countable union and
complement. The pair $(\Omega,\Sigma)$ forms a measurable space.

Next, for a cylinder set
$Cyl=s_{0}\xrightarrow{\mathcal{J}_{0}}s_{1}\xrightarrow{\mathcal{J}_{1}}s_{2}\xrightarrow{\mathcal{J}_{2}}\cdots\xrightarrow{\mathcal{J}_{K-1}}s_{K}$,
we define the probability measure along $Cyl$ as
$\Pr(Cyl)=\mathrm{tr}(\rho^{(K)})$, where $\rho^{(K)}$ is the partial density
operator in the $K$th phase. It can be calculated inductively by

$\left\\{\begin{aligned} \rho^{(0)}=&\ \mathbf{P}_{s_{0}}\cdot\rho(0)\\\
\rho^{(k)}=&\
\mathbf{P}_{s_{k}}\cdot\mathrm{V2L}(\exp(\mathds{M}_{s_{k-1}}|\mathcal{J}_{k-1}|)\
\cdot\\\ &\
\mathrm{L2V}(\mathbf{P}_{s_{k-1}}\cdot\mathrm{V2L}(\exp(\mathds{M}_{s_{k-1}}\inf\mathcal{J}_{k-1})\cdot\mathrm{L2V}(\rho^{(k-1)}))))&\mbox{if
}k>0,\\\ \end{aligned}\right.$ (7)

where $|\mathcal{J}_{k-1}|=\sup\mathcal{J}_{k-1}-\inf\mathcal{J}_{k-1}$ is the
length of $\mathcal{J}_{k-1}$, and $\mathds{M}_{s}$ is the governing matrix
adapted for the transition generator function $Q_{s}$ consisting of the
Hermitian operator $\mathbf{H}\cdot\mathbf{P}_{s}$ and finitely many linear
operators $\mathbf{L}_{j}\mathbf{P}_{s}$. It is preferable to calculate those
matrix multiplication in Eq. (7) with right associativity, according to the
forward semantics. Specifically, $\mathbf{P}_{s_{k-1}}\cdot\rho$ extracts the
term $\op{s_{k-1}}{s_{k-1}}\otimes\rho_{s_{k-1}}$ from the ID $\rho=\sum_{s\in
S}\op{s}{s}\otimes\rho_{s}$; $\exp(\mathds{M}_{s_{k-1}}\tau)$ gives rise to
the state transition function, which carries the density operator of state
$s_{k-1}$ to all states $s$ ($s\in S$) for some sojourn time $\tau$ and
preserves the density operator of state $s$ ($s\neq s_{k-1}$). Totally, the
induction is computing the final partial density operator $\rho^{(k)}$ from
$\rho^{(k-1)}$, satisfying a sojourn at $s_{k-1}$ for some time
$\tau_{k-1}\in\mathcal{J}_{k-1}$. Similar to Vardi’s work for classical Markov
chains [Var85], the domain of $\Pr$ can be extended to $\Sigma$, i.e.,
$\Pr\colon\Sigma\to[0,1]$, which is well-defined under the countable union
$\bigcup_{i}A_{i}$ for any $A_{i}\in\Sigma$ and under the complement
$A^{\textup{c}}$ for any $A\in\Sigma$. Hence the triple $(\Omega,\Sigma,\Pr)$
forms a probability space.

Let us continue to consider Example 3. The cylinder set
$Cyl_{1}=3\xrightarrow{(0,1)}1\xrightarrow{(1,2)}3$ contains all paths that
take the transitions from Node $3$ to Node $1$ and then back to Node $3$ with
sojourn time in the intervals $\mathcal{J}_{1}=(0,1)$ and
$\mathcal{J}_{2}=(1,2)$ respectively. We will compute its probability with the
initial ID $\rho(0)=\op{3}{3}\otimes\op{0}{0}$. Firstly, we have the
projectors $\mathbf{P}_{i}=\op{i}{i}\otimes\mathbf{I}$ for $i=0,1,2,3$ and the
corresponding adapted governing matrices

$\mathds{M}_{i}=\sum_{j=1}^{12}[(\mathbf{L}_{j}\mathbf{P}_{i})\otimes(\mathbf{L}_{j}\mathbf{P}_{i})^{*}-\tfrac{1}{2}(\mathbf{L}_{j}\mathbf{P}_{i})^{\dagger}(\mathbf{L}_{j}\mathbf{P}_{i})\otimes\mathbf{I}-\tfrac{1}{2}\mathbf{I}\otimes(\mathbf{L}_{j}\mathbf{P}_{i})^{\mathrm{T}}(\mathbf{L}_{j}\mathbf{P}_{i})^{*}].$

Then, we inductively calculate the partial density operator in the $k$th
phase:

$\displaystyle\rho^{(0)}$
$\displaystyle=\mathbf{P}_{3}\cdot\rho(0)=\op{3}{3}\otimes\op{0}{0},$
$\displaystyle\rho^{(1)}$
$\displaystyle=\mathbf{P}_{1}\cdot\mathrm{V2L}(\exp(\mathds{M}_{3}|\mathcal{J}_{1}|)\cdot\mathrm{L2V}(\mathbf{P}_{3}\cdot\mathrm{V2L}(\exp(\mathds{M}_{3}\inf\mathcal{J}_{1})\cdot\mathrm{L2V}(\rho^{(0)}))))=\op{1}{1}\otimes\rho_{1},$
$\displaystyle\rho^{(2)}$
$\displaystyle=\mathbf{P}_{3}\cdot\mathrm{V2L}(\exp(\mathds{M}_{1}|\mathcal{J}_{2}|)\cdot\mathrm{L2V}(\mathbf{P}_{1}\cdot\mathrm{V2L}(\exp(\mathds{M}_{1}\inf\mathcal{J}_{2})\cdot\mathrm{L2V}(\rho^{(1)}))))=\op{3}{3}\otimes\rho_{3},$

in which $\rho_{1}$ and $\rho_{3}$ are computable partial density operators on
$\mathcal{H}$ with the following numerical approximations obtained by the
scaling and squaring method embedded with Padé approximations [ML03, AMH09]
for matrix exponentials. That is,
$\exp(\mathds{M})=\big{(}\exp(\mathds{M}/m)\big{)}^{m}$ and
$\exp(\mathds{M}/m)\approx\mathds{P}_{k}$ for some $k,m$, where
$\mathds{P}_{k}$ is the $k$th Padé approximation to be described in Eq. (8),

$\displaystyle\rho_{1}\approx$ $\displaystyle
0.184318\op{0}{0}+(-0.092159-0.037415\imath)[\op{0}{1}+\op{2}{0}]\ +$
$\displaystyle(-0.092159+0.037415\imath)[\op{0}{2}+\op{1}{0}]+0.053744\op{1}{1}\
+$
$\displaystyle(0.03904-0.03820\imath)\op{1}{2}+(0.03904+0.03820\imath)\op{2}{1}+0.053744\op{2}{2}],$
$\displaystyle\rho_{3}\approx$ $\displaystyle
0.016333\op{0}{0}+(-0.008166+0.014145\imath)[\op{0}{1}+\op{1}{2}]\ +$
$\displaystyle(-0.008166-0.014145\imath)[\op{0}{2}+\op{1}{0}]+0.016333\op{1}{1}\
+$
$\displaystyle(-0.008166+0.014145\imath)\op{2}{0}+(-0.008166-0.014145\imath)\op{2}{1}+0.016333\op{2}{2}.$

Finally we get the probability
$\Pr(Cyl_{1})=\mathrm{tr}(\rho^{(2)})=\mathrm{tr}(\rho_{3})\approx 0.048999$.
∎

#### Complexity

By Eq. (7), the procedure of exactly calculating the probability $\Pr(Cyl)$
requires $K$ times of induction phases, for each of which the matrix
exponentiation $\exp(\mathds{M})$ dominates the computational cost. The matrix
exponential $\exp(\mathds{M})$ can be computed by Jordan decomposition as
follows. Since $\mathds{M}$ is a matrix of dimension $N^{2}$ with entries
being algebraic numbers taken from the simple field extension
$\mathbb{Q}(\mu):\mathbb{Q}$, the characteristic polynomial $g(x)$ of
$\mathds{M}$ is an $\mathbb{A}$-polynomial with degree $N^{2}$, which can be
obtained in $\mathcal{O}((N^{2})^{4})=\mathcal{O}(N^{8})$ by [BPR06, Algorithm
8.17]. To further determine the eigenvalues of $\mathds{M}$, we first convert
$g(x)$ to a $\mathbb{Q}$-polynomial $f(x)$. It is achieved by Sylvester
resultant [BPR06, Notation 4.12] — that is a determinant of dimension
$N^{2}+\deg(\mu)$ — to eliminate the common occurrence of $\mu$ in $g(x)$ and
$f_{\mu}$ (the minimal polynomial of $\mu$). It results in the desired
$\mathbb{Q}$-polynomial $f(x)$ with degree at most $D\coloneqq N^{2}\deg(\mu)$
has all roots of $g(x)$. The variable elimination is in $\mathcal{O}(D^{3})$
by the standard determinant computation. The roots of $f(x)$, including the
eigenvalues of $\mathds{M}$, can be determined in
$\mathcal{O}(D^{5}\log(D)^{2})$ by [BPR06, Algorithm 10.4]. So we get the
Jordan decomposition
$\mathds{M}=\mathds{T}\cdot\mathds{J}\cdot\mathds{T}^{-1}$ in
$\mathcal{O}((N^{2})^{3})=\mathcal{O}(N^{6})$ by computing the eigenvectors
with respect to the known eigenvalues, where $\mathds{J}$ is the Jordan
canonical form of $\mathds{M}$ and $\mathds{T}$ is the transformation matrix.
The same complexity holds for computing
$\exp(\mathds{M}\tau)=\mathds{T}\cdot\exp(\mathds{J}\tau)\cdot\mathds{T}^{-1}$.
Hence the exactly calculating procedure is in $\mathcal{O}(K\cdot
D^{5}\log(D)^{2})$ that is bounded polynomially in the encoding size
$\|\mathfrak{Q}\|$ and linearly in the length $K=|Cyl|$.

#### Numerical speed-up

To improve efficiency, we have to address the problem that the uniformization
method [BHHK00] is not directly applicable to compute $\exp(\mathds{M})$,
since the governing matrix $\mathds{M}$ has non-structure of the generator
matrix in a classical CTMC that could be discretized as a Poisson process. For
the general matrix exponentiation, we would employ the state-of-the-art
scaling and squaring method embedded with Padé approximations [ML03, AMH09].
The technique is _reliable_ , as it gives some warning whenever it introduces
excessive errors. The rationale goes as follows. We first _scale_
$\exp(\mathds{M})$ by choosing $m$ to be a power of two satisfying
$\|\exp(\mathds{M}/m)\|\leq 1$ where $\|\,\cdot\,\|$ denotes Frobenius norm.
This $m$ can be chosen to be linear in the size $N^{4}$ of $\mathds{M}$, i.e.
$m\in\mathcal{O}(N^{4})$. Then, as a subroutine, we compute the Padé
approximations $\mathds{P}_{k}$ of $\exp(\mathds{M}/m)$, which is reliable if
$\|\exp(\mathds{M}/m)\|$ is small as in our setting. The $k$th Padé
approximation $\mathds{P}_{k}$ is the fraction

$\left(\sum\limits_{i=0}^{k}\dfrac{(2k-i)!\cdot
k!\cdot(\mathds{M}/m)^{i}}{(2k)!\cdot(k-i)!\cdot
i!}\right)\Big{/}\left(\sum\limits_{i=0}^{k}\dfrac{(2k-i)!\cdot
k!\cdot(-\mathds{M}/m)^{i}}{(2k)!\cdot(k-i)!\cdot i!}\right),$ (8)

and its computational cost lies in $k$ times of matrix multiplication, for
each of which is in $\mathcal{O}((N^{2})^{3})=\mathcal{O}(N^{6})$. The
subroutine terminates whenever the difference
$\|\mathds{P}_{k+1}-\mathds{P}_{k}\|$ is bounded by a predefined tolerance
$\epsilon$. As the sequence $\\{\mathds{P}_{k}\\}_{k\in\mathbb{N}}$ is
exponentially convergent to $\exp(\mathds{M}/m)$, i.e., the error is
exponentially convergent to $0$, we get the desired index
$k\in\mathcal{O}(\log(1/\epsilon))$ or equivalently
$k\in\mathcal{O}(\|\epsilon\|)$ that yields the approximation $\mathds{P}_{k}$
within error $\epsilon$ of $\exp(\mathds{M}/m)$. Finally we _square_
$\exp(\mathds{M}/m)$ by computing the powers $\mathds{P}_{k}^{2^{l}}$ for
$l=1,2,\ldots,\log(m)$ in turn, and get the desired approximation
$\mathds{P}_{k}^{m}$ of $\exp(\mathds{M})$, which costs $\log(m)$ times of
matrix multiplication. Hence the approximately calculating procedure is in
$\mathcal{O}(K\cdot(k+\log(m))\cdot
N^{6})=\mathcal{O}(K\cdot(\|\epsilon\|+\log(N))\cdot N^{6})$. At least a
factor $N^{4}$ in the complexity of the exactly calculating procedure is saved
here.

##  4\. Checking Continuous Stochastic Logic

We introduce the continuous stochastic logic (CSL for short) that consists of
three layers: i) state formulas (which can be true or false in a specific
classical state), ii) path formulas (which can be true or false along a
specific path), and iii) model formulas (true or false for a specific model).
Then we present an algebraic algorithm for checking the CSL formulas over
quantum CTMCs.

The syntax of the CSL formulas, consisting of state formulas $\Phi$, path
formulas $\phi$ and model formulas $\chi$, are defined as follows:

$\displaystyle\Phi$
$\displaystyle\coloneqq\textup{a}\mid\neg\Phi\mid\Phi_{1}\wedge\Phi_{2}$
$\displaystyle\phi$
$\displaystyle\coloneqq\Phi_{0}\mathrm{U}\,^{\mathcal{I}_{0}}\Phi_{1}\mathrm{U}\,^{\mathcal{I}_{1}}\Phi_{2}\cdots\mathrm{U}\,^{\mathcal{I}_{K-1}}\Phi_{K}$
$\displaystyle\chi$ $\displaystyle\coloneqq\mathcal{P}_{\sim\textup{c}}(\phi)$

where $\textup{a}\in AP$ is an atomic proposition,
$\mathcal{I}_{k}\subseteq\mathbb{R}_{\geq 0}$ ($0\leq k<K$) are time intervals
with rational endpoints,
$\mathord{\sim}\in\\{\mathord{=},\mathord{<},\mathord{>}\\}$ is a comparison
operator, and $\textup{c}\in[0,1]\cap\mathbb{Q}$ is a probability threshold.
In particular, the CSL formula is referred to as the model formula $\chi$. For
convenience, we assume that those time intervals $\mathcal{I}_{k}$ are
disjoint left-open and right-closed intervals $(a_{k},b_{k}]$, and $b_{K-1}$
is allowed to be $\infty$. The path formula $\phi$ is called the _multiphase
until_ formula to be interpreted below.

Following [ASSB96, BHHK00], we give the formal semantics of CSL. The semantics
of CSL interpreted over a quantum CTMC $\mathfrak{Q}=(S,Q,L)$ is given by the
satisfaction relation $\models$:

$s$ | $\models\textup{a}$ |  | if $\textup{a}\in L(s)$;
---|---|---|---
$s$ | $\models\neg\Phi$ |  | if $s\not\models\Phi$;
$s$ | $\models\Phi_{1}\wedge\Phi_{2}$ |  | if $s\models\Phi_{1}$ and $s\models\Phi_{2}$;
$\omega$ | $\models\Phi_{0}\mathrm{U}\,^{\mathcal{I}_{0}}\Phi_{1}\mathrm{U}\,^{\mathcal{I}_{1}}\Phi_{2}\cdots\mathrm{U}\,^{\mathcal{I}_{K-1}}\Phi_{K}$ |  | if there exist $t_{0}<t_{1}<\cdots<t_{K-1}$
|  such that for each integer $0\leq k<K$, we have $t_{k}\in\mathcal{I}_{k}$
|  and
$\forall\,t^{\prime}\in[t_{k-1},t_{k})\,:\,\omega(t^{\prime})\models\Phi_{k}$
where $t_{-1}$ is defined to be zero,
|  and additionally $\omega(t_{K-1})\models\Phi_{K}$;
$\mathfrak{Q}$ | $\models\mathcal{P}_{\sim\textup{c}}(\phi)$ |  | if $\Pr(\\{\omega\,:\,\omega\models\phi\\})\sim\textup{c}$.

Here, the real numbers $t_{1},\ldots,t_{K-1}$ in the interpretation of the
multiphase until formula $\phi$ are switch time-stamps, and the time intervals
$\mathcal{I}_{k}$ ($0\leq k<K$) in $\phi$ are _absolute_ timing that are
measured from the start of the path, which differs from the LTL interpretation
of the nested until formula with the _relative_ timing that are measured from
the instants when the prior transitions take place.

The core of the model-checking algorithm lies in measuring all the paths
satisfying the multiphase until formula $\phi$, which can be tackled by the
following lemma. Before that, we would define the following four notions.

  1. (1)

Let $\mathbf{P}_{\Phi}\coloneqq\sum_{s\models\Phi}\op{s}{s}\otimes\mathbf{I}$
be the projector for some state formula $\Phi$, and
$\mathds{P}_{\Phi}=\mathbf{P}_{\Phi}\otimes\mathbf{P}_{\Phi}$.

  2. (2)

The transition generator function $Q_{\Phi}$ is given by the Hermitian
operator $\mathbf{H}\cdot\mathbf{P}_{\Phi}$ and the finite set of linear
operators $\mathbf{L}_{j}\mathbf{P}_{\Phi}$. Let $\mathds{M}_{\Phi}$ be the
governing matrix adapted for $Q_{\Phi}$. The transition generator function
$Q_{\Phi}$ and the corresponding governing matrix $\mathds{M}_{\Phi}$ keep the
behaviors of the $\Phi$-states and absorb the other states.

  3. (3)

Recall in probability theory that for an exponentially distributed random
variable with rate parameter $\kappa$, the _cumulative distribution function_
$CDF(t)$ is $1-\exp(\kappa t)$ that determines the probability of the random
variable bounded from above by $t$. So the change in the probability for the
random variable falling into a time interval $(a,b]$ is
$CDF(b)-CDF(a)=\exp(\kappa a)-\exp(\kappa b)$. Correspondingly, the matrix
exponential $\exp(\mathds{M}t)$ is the state transition function that
determines the change in the density operator representing ID in a quantum
CTMC.

  4. (4)

Again, in probability theory, the _probability density function_ $PDF(t)$ is
$\kappa\exp(\kappa t)$ that determines the change rate in the probability of
the random variable, since $CDF(b)-CDF(a)=\int_{a}^{b}PDF(t)\ \mathrm{d}\,t$.
Correspondingly, $\mathds{M}\cdot\exp(\mathds{M}t)$ determines the change rate
in the density operator.

They are powerful tools to develop the algorithm for checking the CSL formula
over quantum CTMCs, together with definite integrals as those for classical
CTMCs in [XZJ+16].

######  Lemma 1.

For a multiphase until formula
$\phi=\Phi_{0}\mathrm{U}\,^{\mathcal{I}_{1}}\Phi_{1}\mathrm{U}\,^{\mathcal{I}_{2}}\Phi_{2}\cdots\mathrm{U}\,^{\mathcal{I}_{K}}\Phi_{K}$,
the final partial density operator $\rho^{(K)}$ measuring all the paths that
satisfy $\phi$ can be calculated inductively by

$\mathrm{L2V}(\rho^{(k)})=\begin{cases}\mathds{P}_{\Phi_{0}}\cdot\mathrm{L2V}(\rho(0))\qquad
k=0,\\\
\mathds{P}_{\Phi_{k}}\cdot[\exp(\mathds{M}_{\Phi_{k-1}}(b_{k}-b_{k-1}))+\int_{a_{k}}^{b_{k}}\exp(\mathds{M}_{\Phi_{k}}(b_{k}-t_{k}))\cdot\mathds{P}_{\Phi_{k-1}}\
\cdot\\\
\qquad\mathds{P}_{\neg\Phi_{k-1}}\cdot\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(t_{k}-a_{k}))\
\mathrm{d}\,t_{k}\ \cdot\\\
\qquad\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))]\cdot\mathrm{L2V}(\rho^{(k-1)})\qquad\mbox{if
}0<k<K,\\\
\mathds{P}_{\Phi_{k}}\cdot[\exp(\mathds{M}_{\Phi_{k-1}\wedge\neg\Phi_{k}}(b_{k}-a_{k}))\cdot\mathds{P}_{\Phi_{k-1}}+\mathds{P}_{\Phi_{k-1}}]\
\cdot\\\
\qquad\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))\cdot\mathrm{L2V}(\rho^{(k-1)})\qquad\mbox{if
}k=K,\end{cases}$ (9)

where $a_{k}=\inf\mathcal{I}_{k}$, $b_{k}=\sup\mathcal{I}_{k}$ and $b_{0}$ is
defined to be zero for convenience; and thus the satisfaction probability
$\Pr(\\{\omega\colon\omega\models\phi\\})$ is $\mathrm{tr}(\rho^{(K)})$.

######  Proof 4.1.

We prove it by an induction on the phase index $k$.

Initially, $\rho^{(0)}$ is the partial density operator
$\sum_{s\models\Phi_{0}}\op{s}{s}\otimes\rho_{s}$ extracted from
$\rho(0)=\sum_{s\in S}\op{s}{s}\otimes\rho_{s}$ by the projector
$\mathbf{P}_{\Phi_{0}}$ (see Notion 1), which measures the paths that satisfy
$\phi$ at time $0$.

Let $\rho^{(k-1)}$ ($k<K$) be the partial density operator measuring the paths
that satisfy $\phi$ up to time $b_{k-1}$. We proceed to construct the partial
density operator $\rho^{(k)}$ measuring the paths that satisfy $\phi$ up to
time $b_{k}$ as follows. Based on the semantics of the $k$th phase of $\phi$,
all the paths that satisfy $\phi$ can be classified into two disjoint path
sets: one contains the paths leaving $\Phi_{k-1}$-states for $\Phi_{k}$-ones
during $\mathcal{I}_{k}$, i.e., for some switch time-stamp
$t_{k}\in\mathcal{I}_{k}=(a_{k},b_{k}]$,

$\\{\omega\colon[\forall\,t^{\prime}\in(b_{k-1},t_{k})\colon\omega(t^{\prime})\models\Phi_{k-1}]\wedge[\omega(t_{k})\models\neg\Phi_{k-1}]\wedge[\forall\,t^{\prime\prime}\in[t_{k},b_{k}]\colon\omega(t^{\prime\prime})\models\Phi_{k}]\\},$

and the other contains the paths that do not have to leave $\Phi_{k-1}$-states
during $\mathcal{I}_{k}$, i.e.,
$\\{\omega\colon[\forall\,t^{\prime}\in(b_{k-1},b_{k}]\colon\omega(t^{\prime})\models\Phi_{k-1}]\wedge[\omega(b_{k})\models\Phi_{k}]\\}$.
Here, “leaving $\Phi_{k-1}$-states” means that the path $\omega$ satisfies
$\phi$ but at some $t\in\mathcal{I}_{k}$, $\omega(t)$ dissatisfies
$\Phi_{k-1}$ and it must enter into $\Phi_{k}$-states then. Correspondingly,
the former path set is measured by the partial density operator

$\mathrm{V2L}\left(\begin{aligned}
&\mathds{P}_{\Phi_{k}}\cdot\int_{a_{k}}^{b_{k}}\exp(\mathds{M}_{\Phi_{k}}(b_{k}-t_{k}))\cdot\mathds{P}_{\neg\Phi_{k-1}}\cdot\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(t_{k}-a_{k}))\
\mathrm{d}\,t_{k}\ \cdot\\\
&\mathds{P}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))\cdot\mathrm{L2V}(\rho^{(k-1)})\end{aligned}\right),$

where

  * •

$\mathds{M}_{\Phi_{k-1}}$ and $\mathds{M}_{\Phi_{k}}$ are the governing
matrices respectively adapted to $Q_{\Phi_{k-1}}$ and $Q_{\Phi_{k}}$ (see
Notion 2),

  * •

$\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))$ and
$\exp(\mathds{M}_{\Phi_{k}}(b_{k}-t_{k}))$ are the state transition functions
(see Notion 3), and

  * •

$\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(t_{k}-a_{k}))$ gives
the rate of state transition (see Notion 4).

The latter path set is measured by the partial density operator
$\mathrm{V2L}(\mathds{P}_{\Phi_{k}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(b_{k}-b_{k-1}))\cdot\mathrm{L2V}(\rho^{(k-1)}))$.
Since the two path sets are disjoint, $\rho^{(k)}$ is the sum of the above two
partial density operators.

Let $\rho^{(K-1)}$ be the partial density operator measuring the paths that
satisfy $\phi$ up to time $b_{K-1}$. We proceed to construct the partial
density operator $\rho^{(K)}$ measuring the paths that satisfy $\phi$ up to
time $b_{K}$ as follows. Based on the semantics of the $K$th phase of $\phi$,
all the paths that satisfy $\phi$ are classified into two disjoint sets: one
is not immediately reaching $\Phi_{K}$-states during $\mathcal{I}_{K}$, i.e.,
for some switch time-stamp $t_{K}\in\mathcal{I}_{K}=(a_{K},b_{K}]$,

$\\{\omega\colon[\forall\,t^{\prime}\in[b_{K-1},t_{K})\colon\omega(t^{\prime})\models\Phi_{K-1}]\wedge[\forall\,t^{\prime\prime}\in[a_{K},t_{K})\colon\omega(t^{\prime\prime})\models\neg\Phi_{K}]\wedge[\omega(t_{K})\models\Phi_{K}]\\},$

and the other is immediately reaching $\Phi_{K}$-states at $a_{K}$,
$\\{\omega\colon[\forall\,t^{\prime}\in[b_{K-1},a_{K})\colon\omega(t^{\prime})\models\Phi_{K-1}]\wedge[\omega(a_{K})\models\Phi_{K}]\\}$.
Here, “not immediately reaching $\Phi_{K}$-states” means that the satisfying
path $\omega$ dissatisfies $\Phi_{K}$ at $a_{K}$, and it enters in
$\Phi_{K}$-states for some $t>a_{K}$. Correspondingly, the two disjoint path
sets are measured by the partial density operators:

$\displaystyle\mathrm{V2L}\big{(}\mathds{P}_{\Phi_{k}}\cdot\exp(\mathds{M}_{\Phi_{k-1}\wedge\neg\Phi_{k}}(b_{k}-a_{k}))\cdot\mathds{P}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))\cdot\mathrm{L2V}(\rho^{(k-1)})\big{)}$
$\displaystyle\mathrm{V2L}(\mathds{P}_{\Phi_{k}\wedge\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(a_{k}-b_{k-1}))\cdot\mathrm{L2V}(\rho^{(k-1)})),$

whose sum is exactly $\rho^{(K)}$. Hence $\rho^{(K)}$ is the final partial
density operator measuring the paths that satisfy the whole path formula
$\phi$.

When tackling $b_{K}=\infty$, we do not suffer from the trouble that the
bottom strongly connected component (BSCC) subspaces would make the orbit of
density operators periodical, since the probability $\mathrm{tr}(\rho^{(K)})$
of satisfying $\phi$ is monotonically increasing along $b_{K}$. So it can be
achieved by taking the limit.

The percolation performance is an essential property of networks that
describes transitions between nodes and their neighborhood. Reconsidering the
first generation of AN in Example 3, the return probability of the center node
reflects the percolation performance. We aim to compute the probability of the
following return event.

> The qutrit starts at the center node, then transits to non-center nodes in
> one unit of time and sojourns there up to (absolute) time at least $1$, and
> finally returns to the center node within time $2$.

This could be formally specified by the CSL formula
$\phi_{1}\equiv\Phi_{0}\mathrm{U}\,^{\mathcal{I}_{1}}\Phi_{1}\mathrm{U}\,^{\mathcal{I}_{2}}\Phi_{2}$,
where $\Phi_{0}\equiv\Phi_{2}\equiv\texttt{center}$,
$\Phi_{1}\equiv\neg\texttt{center}$, and time intervals
$\mathcal{I}_{1}=(0,1]$ and $\mathcal{I}_{2}=(1,2]$. Firstly, we take
$\rho(0)=\op{3}{3}\otimes\tfrac{1}{3}\mathbf{I}$ and define

$\displaystyle\mathds{P}_{\Phi_{0}}$
$\displaystyle=\mathds{P}_{\Phi_{2}}=\mathbf{P}_{\Phi_{0}}\otimes\mathbf{P}_{\Phi_{0}}=(\op{3}{3}\otimes\mathbf{I})\otimes(\op{3}{3}\otimes\mathbf{I}),$
$\displaystyle\mathds{P}_{\Phi_{1}}$
$\displaystyle=\mathds{P}_{\neg\Phi_{0}}=\mathds{P}_{\Phi_{1}\wedge\neg\Phi_{2}}=\mathbf{P}_{\Phi_{1}}\otimes\mathbf{P}_{\Phi_{1}}$
$\displaystyle=[(\op{0}{0}+\op{1}{1}+\op{2}{2})\otimes\mathbf{I}]\otimes[(\op{0}{0}+\op{1}{1}+\op{2}{2})\otimes\mathbf{I}],$
$\displaystyle\mathds{M}_{\Phi_{i}}$
$\displaystyle=\sum_{j=1}^{12}[(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})\otimes(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})^{*}-\tfrac{1}{2}(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})^{\dagger}(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})\otimes\mathbf{I}-\tfrac{1}{2}\mathbf{I}\otimes(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})^{\mathrm{T}}(\mathbf{L}_{j}\mathbf{P}_{\Phi_{i}})^{*}],$

for $i=0$ or $1$. Then, we calculate $\rho^{(2)}$ in an inductive fashion:

$\displaystyle\mathrm{L2V}(\rho^{(0)})=$ $\displaystyle\
\mathds{P}_{\Phi_{0}}\cdot\mathrm{L2V}(\rho(0)),$
$\displaystyle\mathrm{L2V}(\rho^{(1)})=$ $\displaystyle\
\mathds{P}_{\Phi_{1}}\cdot[\exp(\mathds{M}_{\Phi_{0}}(1-0))+\int_{0}^{1}\exp(\mathds{M}_{\Phi_{1}}(1-t_{1}))\cdot\mathds{P}_{\neg\Phi_{0}}\cdot\mathds{M}_{\Phi_{0}}\
\cdot$ $\displaystyle\qquad\exp(\mathds{M}_{\Phi_{0}}(t_{1}-0))\
\mathrm{d}\,t_{1}\cdot\mathds{P}_{\Phi_{0}}\cdot\exp(\mathds{M}_{\Phi_{0}}(0-0))]\cdot\mathrm{L2V}(\rho^{(0)}),$
$\displaystyle\mathrm{L2V}(\rho^{(2)})=$ $\displaystyle\
[\mathds{P}_{\Phi_{2}}\cdot\exp(\mathds{M}_{\Phi_{1}\wedge\neg\Phi_{2}}(2-1))+\mathds{P}_{\Phi_{2}}]\cdot\mathds{P}_{\Phi_{1}}\cdot\exp(\mathds{M}_{\Phi_{1}}(1-1))\cdot\mathrm{L2V}(\rho^{(1)}).$

The result is $\rho^{(2)}=\op{3}{3}\otimes\rho_{3}$ with computable partial
density operator $\rho_{3}$ on $\mathcal{H}$. We replace the definite integral
with the Riemann sum under a predefined step length, saying $\tfrac{1}{100}$,

$\frac{1}{100}\sum_{i=1}^{100}\exp(\mathds{M}_{\Phi_{1}}\tfrac{100-i}{100})\cdot\mathds{P}_{\neg\Phi_{0}}\cdot\mathds{M}_{\Phi_{0}}\cdot\exp(\mathds{M}_{\Phi_{0}}\tfrac{i}{100}),$

and get the numerical approximation

$\displaystyle\rho_{3}\approx$ $\displaystyle
0.105266[\op{0}{0}+\op{1}{1}+\op{2}{2}]+(-0.052633+0.091163\imath)[\op{0}{1}+\op{1}{2}+\op{2}{0}]\
+$ $\displaystyle(-0.052633-0.091163\imath)[\op{0}{2}+\op{1}{0}+\op{2}{1}].$

Thus we obtain the return probability
$\Pr(\phi)=\mathrm{tr}(\rho^{(2)})=\mathrm{tr}(\rho_{3})\approx 0.315798$. ∎

#### Complexity

The procedure of exactly calculating the final density operator by Eq. (9) is
based on the exact matrix exponentiation $\exp(\mathds{M}\tau)$ that costs
$\mathcal{O}(D^{5}\log(D)^{2})$. Additionally, it requires $K-1$ times of
induction phases, for each of which the definite integration

$\int_{a_{k}}^{b_{k}}\exp(\mathds{M}_{\Phi_{k}}(b_{k}-t_{k}))\cdot\mathds{P}_{\neg\Phi_{k-1}}\cdot\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(t_{k}-a_{k}))\
\mathrm{d}\,t_{k}$

dominates the computational cost. The definite integral can be computed by at
most $N^{2}-1$ times of integration-by-parts, whose integrands can be obtained
in $\mathcal{O}(N^{6})$ as products of a few matrices with dimension $N^{2}$.
Hence the exactly calculating procedure is in
$\mathcal{O}(K(D^{5}\log(D)^{2}+N^{8}))=\mathcal{O}(K\cdot D^{5}\log(D)^{2})$
that is bounded polynomially in the encoding size $\|\mathfrak{Q}\|$ and
linearly in the size $\|\phi\|$. However, if we adopt the Riemann sum to do
the approximate definite integration, after setting a step length $\delta$ or
equivalently the number of samples $M=(b_{k}-a_{k})/\delta$, we could get an
approximate value
$\frac{1}{M}\sum_{i=1}^{M}\exp(\mathds{M}_{\Phi_{k}}(M-i)\delta)\cdot\mathds{P}_{\neg\Phi_{k-1}}\cdot\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}i\delta)$,
each term is obtained in $\mathcal{O}(N^{6})$ as a product of a few matrices
with dimension $N^{2}$. Then the approximately calculating procedure is in
$\mathcal{O}(K\cdot(\|\epsilon\|+\log(N)+M)\cdot N^{6})$.

######  Theorem 2.

The CSL formula is decidable over quantum CTMCs.

######  Proof 4.2.

Using Jordan decomposition, matrix exponentiation and integration-by-parts, we
first claim that the satisfaction probability can be expressed as a number in
the explicit form

$p=\alpha_{1}\mathrm{e}^{\lambda_{1}}+\alpha_{2}\mathrm{e}^{\lambda_{2}}+\cdots+\alpha_{m}\mathrm{e}^{\lambda_{m}},$
(10)

where $\alpha_{1},\ldots,\alpha_{m}$ are algebraic numbers and
$\lambda_{1},\ldots,\lambda_{m}$ are distinct algebraic numbers. We observe
the facts:

  1. (1)

The governing matrix $\mathds{M}_{\Phi}$ appearing in Eq. (9) takes algebraic
numbers as entries, since the entries of $\mathbf{H}_{s}$ and
$\mathbf{L}_{j,s_{1},s_{2}}$ are algebraic and the endpoints of
$\mathcal{I}_{k}$ are rational.

  2. (2)

The characteristic polynomial of $\mathds{M}_{\Phi}$ is an
$\mathbb{A}$-polynomial. The eigenvalues $\lambda_{1},\ldots,\lambda_{m}$ of
$\mathds{M}_{\Phi}$ are algebraic, as they are roots of that
$\mathbb{A}$-polynomial by Lemma 2.3. Those eigenvalues make up all exponents
in Eq. (10).

  3. (3)

The entries of the matrix exponential $\exp(\mathds{M}_{\Phi}\cdot a)$ are in
the form (10).

  4. (4)

The same structure holds for entries of

$\int_{a_{k}}^{b_{k}}\exp(\mathds{M}_{\Phi_{k}}(b_{k}-t_{k}))\cdot\mathds{P}_{\neg\Phi_{k-1}}\cdot\mathds{M}_{\Phi_{k-1}}\cdot\exp(\mathds{M}_{\Phi_{k-1}}(t_{k}-a_{k}))\
\mathrm{d}\,t_{k}.$

By Theorem 2.3, we have that $p=\textup{c}$ if and only if $m=1$,
$\lambda_{1}=0$ and $\alpha_{1}=\textup{c}$. If not, we further compute the
value $p$ up to any precision by sufficiently approaching the Euler constant
$\mathrm{e}$, and thereby compare $p$ with c, which decides the truth of
$\mathfrak{Q}\models\mathcal{P}_{\sim\textup{c}}(\phi)$.

#### Complexity

The decision algorithm is in time polynomial in the encoding size of the input
quantum CTMC and linear in the size of the input CSL formula, besides the
query/oracle — how sufficiently the Euler constant $\mathrm{e}$ should be
approached. Whereas the satisfaction probability can be computed in polynomial
time using the numerical methods — the scaling and squaring method embedded
with Padé approximations and the Riemann sum.

##  5\. Concluding Remarks

This paper introduced a novel model of quantum CTMC consisting of a classical
subsystem and a quantum subsystem. The essential difference of these two
systems lies in that a quantum system allows superpositions while a classical
system not. The classical subsystem in the current model forbids
superpositions, making it more precise than the model of quantum CTMC proposed
in [XMGY21], which does not take this point into consideration. Then we
checked the famous branching logic — CSL — against the quantum CTMC. The
decidability is established by an algebraic approach to tackling multiphase
until formulas, which lie at the core of CSL. To be more efficient, numerical
methods could be incorporated to yield a polynomial-time procedure.
Furthermore, we provided the Apollonian network as a running example to
demonstrate the method mentioned above.

For future work, an interesting direction is model-checking the quantum
continuous-time Markov decision processes (CTMDP) [YY18], in which
probabilistic and nondeterministic behaviors coexist. To resolve
nondeterminism, following the verification work [Mil68, MSS20] in classical
CTMDPs, an optimal scheduler is expected to maximize the probability of a
given property over almost all time-stamps at each classical state. Here, it
follows that local optimization yields global optimization. In quantum
systems, however, there are additionally quantum states that are described by
density operators. The known way to compare density operators $\rho$ are based
on the Löwner partial order $\sqsubseteq$, i.e., $\rho_{1}\sqsubseteq\rho_{2}$
if $\rho_{2}-\rho_{1}$ is positive. Thus local optimization does not yield
global optimization generally, which would bring nontrivial technical hardness
to be resolved.

## References

  * [AHAdS05] J. S. Andrade, H. J. Herrmann, R. F. S. Andrade, and L. R. da Silva.  Apollonian networks: Simultaneously scale-free, small world, Euclidean, space filling, and with matching graphs.  Physical Review Letters, 94(1):article no. 018702, 2005.  doi:10.1103/PhysRevLett.94.018702.
  * [AMH09] A. H. Al-Mohy and N. J. Higham.  A new scaling and squaring algorithm for the matrix exponential.  SIAM Journal on Matrix Analysis and Applications, 31(3):970–989, 2009.  doi:10.1137/09074721X.
  * [ASSB96] A. Aziz, K. Sanwal, V. Singhal, and R. Brayton.  Verifying continuous time Markov chains.  In R. Alur and T. A. Henzinger, editors, Computer Aided Verification: 8th International Conference, CAV’96, volume 1102 of LNCS, pages 269–276. Springer, 1996.  doi:10.1007/3540614745\\_75.
  * [Bak75] A. Baker.  Transcendental Number Theory.  Cambridge University Press, 1975.
  * [BHHK00] C. Baier, B. R. Haverkort, H. Hermanns, and J.-P. Katoen.  Model checking continuous-time Markov chains by transient analysis.  In E. Allen Emerson and A. Prasad Sistla, editors, Computer Aided Verification, 12th International Conference, CAV 2000, volume 1855 of LNCS, pages 358–372. Springer, 2000.  doi:10.1007/10722167\\_28.
  * [BK08] C. Baier and J.-P. Katoen.  Principles of Model Checking.  MIT Press, 2008.
  * [BKH99] C. Baier, J.-P. Katoen, and H. Hermanns.  Approximate symbolic model checking of continuous-time Markov chains.  In J. C. M. Baeten and S. Mauw, editors, CONCUR’99 Concurrency Theory, volume 1664 of LNCS, pages 146–161. Springer, 1999.  doi:10.1007/3540483209\\_12.
  * [BPR06] S. Basu, R. Pollack, and M.-F. Roy.  Algorithms in Real Algebraic Geometry.  Springer, 2 edition, 2006.  doi:10.1007/3540330992.
  * [CDG21] J. Chow, O. Dial, and J. Gambetta.  IBM Quantum breaks the 100‑qubit processor barrier, 2021.  URL: https://research.ibm.com/blog/127-qubit-quantum-processor-eagle.
  * [CGP99] E. M. Clarke, O. Grumberg, and D. A. Peled.  Model Checking.  MIT Press, 1999.
  * [Coh93] H. Cohen.  A Course in Computational Algebraic Number Theory.  Springer, 1993.  doi:10.1007/9783662029459.
  * [COW16] V. Chonev, J. Ouaknine, and J. Worrell.  On the Skolem problem for continuous linear dynamical systems.  In I. Chatzigiannakis, M. Mitzenmacher, Y. Rabani, and D. Sangiorgi, editors, 43rd International Colloquium on Automata, Languages, and Programming, ICALP 2016, volume 55 of LIPIcs, pages 100:1–100:13. Schloss Dagstuhl, 2016.  doi:10.4230/LIPIcs.ICALP.2016.100.
  * [DJKV17] C. Dehnert, S. Junges, J.-P. Katoen, and M. Volk.  A Storm is coming: A modern probabilistic model checker.  In R. Majumdar and V. Kuncak, editors, Computer Aided Verification - 29th International Conference, CAV 2017, Part II, volume 10427 of LNCS, pages 592–600. Springer, 2017.  doi:10.1007/9783319633909\\_31.
  * [DP06] E. D’Hondt and P. Panangaden.  Quantum weakest preconditions.  Mathematical Structures in Computer Science, 16(3):429–451, 2006\.  doi:10.1017/S0960129506005251.
  * [Fey86] R. P. Feynman.  Quantum mechanical computers.  Foundations of Physics, 16(6):507–531, 1986.  doi:10.1007/BF01886518.
  * [FHL+22] C. Fu, E. M. Hahn, Y. Li, S. Schewe, M. Sun, A. Turrini, and L. Zhang.  EPMC gets knowledge in multi-agent systems.  In B. Finkbeiner and T. Wies, editors, Verification, Model Checking, and Abstract Interpretation - 22nd International Conference, VMCAI 2022, volume 13182 of LNCS, pages 93–107. Springer, 2022.  doi:10.1007/9783030945831\\_5.
  * [FHTY17] Y. Feng, E. M. Hahn, A. Turrini, and S. Ying.  Model checking $\omega$-regular properties for quantum Markov chains.  In R. Meyer and U. Nestmann, editors, 28th International Conference on Concurrency Theory, CONCUR 2017, volume 85 of LIPIcs, pages 35:1–35:16. Schloss Dagstuhl, 2017.  doi:10.4230/LIPIcs.CONCUR.2017.35.
  * [FYY13] Y. Feng, N. Yu, and M. Ying.  Model checking quantum Markov chains.  Journal of Computer and System Sciences, 79(7):1181–1198, 2013\.  doi:10.1016/j.jcss.2013.04.002.
  * [GFY18] J. Guan, Y. Feng, and M. Ying.  Decomposition of quantum Markov chains and its applications.  Journal of Computer and System Sciences, 95:55–68, 2018.  doi:10.1016/j.jcss.2018.01.005.
  * [GKS76] V. Gorini, A. Kossakowski, and E. C. G. Sudarshan.  Completely positive dynamical semigroups of $n$-level systems.  Journal of Mathematical Physics, 17(5):821–825, 1976.  doi:10.1063/1.522979.
  * [GLR+13] A. S. Green, P. L. Lumsdaine, N. J. Ross, P. Selinger, and B. Valiron.  Quipper: a scalable quantum programming language.  In H.-J. Boehm and C. Flanagan, editors, ACM SIGPLAN Conference on Programming Language Design and Implementation, PLDI 2013, pages 333–342. ACM, 2013.  doi:10.1145/2491956.2462177.
  * [GNP08] S. J. Gay, R. Nagarajan, and N. Papanikolaou.  QMC: A model checker for quantum systems.  In A. Gupta and S. Malik, editors, Computer Aided Verification, 20th International Conference, CAV 2008, volume 5123 of LNCS, pages 543–547. Springer, 2008.  doi:10.1007/9783540705451\\_51.
  * [Gro96] L. K. Grover.  A fast quantum mechanical algorithm for database search.  In Proc. 28th Annual ACM Symposium on the Theory of Computing, pages 212–219. ACM, 1996.  doi:10.1145/237814.237866.
  * [GY22] J. Guan and N. Yu.  A probabilistic logic for verifying continuous-time Markov chains.  In D. Fisman and G. Rosu, editors, Tools and Algorithms for the Construction and Analysis of Systems: 28th International Conference, TACAS, LNCS. Springer, 2022.  to appear.  URL: https://arxiv.org/abs/2004.08059.
  * [HHL09] A. W. Harrow, A. Hassidim, and S. Lloyd.  Quantum algorithm for linear systems of equations.  Physical Review Letters, 103(15):article no. 150502, 2009.  doi:10.1103/PhysRevLett.103.150502.
  * [HLXL18] C.-C. Huang, J.-C. Li, M. Xu, and Z.-B. Li.  Positive root isolation for poly-powers by exclusion and differentiation.  Journal of Symbol Computation, 85:148–169, 2018.  doi:10.1016/j.jsc.2017.07.007.
  * [Kai80] T. Kailath.  Linear Systems.  Prentice Hall, 1980.
  * [KNP11] M. Kwiatkowska, G. Norman, and D. Parker.  PRISM 4.0: Verification of probabilistic real-time systems.  In G. Gopalakrishnan and S. Qadeer, editors, Computer Aided Verification: 23rd International Conference, CAV 2011, volume 6806 of LNCS, pages 585–591. Springer, 2011.  doi:10.1007/9783642221101\\_47.
  * [LF15] L. Li and Y. Feng.  Quantum Markov chains: Description of hybrid systems, decidability of equivalence, and model checking linear-time properties.  Information and Computation, 244:229–244, 2015.  doi:10.1016/j.ic.2015.07.001.
  * [Lin76] G. Lindblad.  On the generators of quantum dynamical semigroups.  Communications in Mathematical Physics, 48(2):119–130, 1976.  doi:10.1007/BF01608499.
  * [Loo83] R. Loos.  Computing in algebraic extensions.  In Computer Algebra: Symbolic and Algebraic Computation, pages 173–187. Springer, 1983.  doi:10.1007/9783709175514\\_12.
  * [Mil68] B. L. Millor.  Finite state continuous time Markov decision processes with a finite planning horizon.  SIAM Journal on Control, 6(2):266–280, 1968.  doi:10.1137/0306020.
  * [ML03] C. B. Moler and C. Van Loan.  Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later.  SIAM Review, 45(1):3–49, 2003.  doi:10.1137/S00361445024180.
  * [MSS20] R. Majumdar, M. Salamati, and S. Soudjani.  On decidability of time-bounded reachability in CTMDPs.  In A. Czumaj, A. Dawar, and E. Merelli, editors, 47th International Colloquium on Automata, Languages, and Programming, ICALP 2020, volume 168 of LIPIcs, pages 133:1–133:19. Schloss Dagstuhl, 2020\.  doi:10.4230/LIPIcs.ICALP.2020.133.
  * [NC00] M. A. Nielsen and I. L. Chuang.  Quantum Computation and Quantum Information.  Cambridge University Press, 2000.
  * [Öme98] B. Ömer.  A procedural formalism for quantum computing.  Technical report, Technical University of Vienna, 1998.  doi:10.1.1.121.5017.
  * [PGMS15] L. Pawela, P. Gawron, J. A. Miszczak, and P. Sadowski.  Generalized open quantum walks on Apollonian networks.  PLOS ONE, 10(7):1–23, 2015.  doi:10.1371/journal.pone.0130967.
  * [Sel04] P. Selinger.  Towards a quantum programming language.  Mathematical Structures in Computer Science, 14(4):527–586, 2004\.  doi:10.1017/S0960129504004256.
  * [Sho94] P. W. Shor.  Algorithms for quantum computation: Discrete logarithms and factoring.  In Proc. 35th Annual Symposium on Foundations of Computer Science, pages 124–134. IEEE Computer Society, 1994.  doi:10.1109/SFCS.1994.365700.
  * [SZ00] J. W. Sanders and P. Zuliani.  Quantum programming.  In R. Backhouse and J. N. Oliveira, editors, Mathematics of Program Construction, pages 80–99. Springer, 2000.  doi:10.1007/10722010\\_6.
  * [Var85] M. Y. Vardi.  Automatic verification of probabilistic concurrent finite-state programs.  In Proc. 26th Annual Symposium on Foundations of Computer Science, pages 327–338. IEEE Computer Society, 1985.  doi:10.1109/SFCS.1985.12.
  * [VRS+15] B. Valiron, N. J. Ross, P. Selinger, D. Scott Alexander, and J. M. Smith.  Programming the quantum future.  Communications of the ACM, 58(8):52–61, 2015.  doi:10.1145/2699415.
  * [XLL08] X.-P. Xu, W. Li, and F. Liu.  Coherent transport on Apollonian networks and continuous-time quantum walks.  Physical Review E, 78:article no. 052103, 2008.  doi:10.1103/PhysRevE.78.052103.
  * [XMGY21] M. Xu, J. Mei, J. Guan, and N. Yu.  Model checking quantum continuous-time Markov chains.  In S. Haddad and D. Varacca, editors, 32nd International Conference on Concurrency Theory, CONCUR 2021, volume 203 of LIPIcs, pages 13:1–13:17. Schloss Dagstuhl, 2021.  doi:10.4230/LIPIcs.CONCUR.2021.13.
  * [XZJ+16] M. Xu, L. Zhang, D. N. Jansen, H. Zhu, and Z. Yang.  Multiphase until formulas over Markov reward models: An algebraic approach.  Theoretical Computer Science, 611:116–135, 2016.  doi:10.1016/j.tcs.2015.07.047.
  * [YF21] M. Ying and Y. Feng.  Model Checking Quantum Systems: Principles and Algorithms.  Cambridge University Press, 2021.  doi:10.1017/9781108613323.
  * [YFYY13] S. Ying, Y. Feng, N. Yu, and M. Ying.  Reachability probabilities of quantum Markov chains.  In P. R. D’Argenio and H. C. Melgratti, editors, CONCUR 2013: Concurrency Theory - 24th International Conference, volume 8052 of LNCS, pages 334–348. Springer, 2013.  doi:10.1007/9783642401848\\_24.
  * [YY18] S. Ying and M. Ying.  Reachability analysis of quantum Markov decision processes.  Information and Computation, 263:31–51, 2018.  doi:10.1016/j.ic.2018.09.001.
  * [ZJNH11] L. Zhang, D. N. Jansen, F. Nielson, and H. Hermanns.  Automata-based CSL model checking.  In L. Aceto, M. Henzinger, and J. Sgall, editors, Automata, Languages and Programming: 38th International Colloquium, ICALP 2011, Part II, volume 6756 of LNCS, pages 271–282. Springer, 2011.  doi:10.1007/9783642220128\\_21.
  * [ZJNH12] L. Zhang, D. N. Jansen, F. Nielson, and H. Hermanns.  Efficient CSL model checking using stratification.  Logical Methods in Computer Science, 8(2:17):1–18, 2012.  doi:10.2168/LMCS-8(2:17)2012.
  * [ZWD+20] H. Zhong, H. Wang, Y. Deng, M. Chen, L. Peng, Y. Luo, J. Qin, D. Wu, X. Ding, Y. Hu, P. Hu, X. Yang, W. Zhang, H. Li, Y. Li, X. Jiang, L. Gan, G. Yang, L. You, Z. Wang, L. Li, N. Liu, C. Lu, and J. Pan.  Quantum computational advantage using photons.  Science, 370(6523):1460–1463, 2020.  doi:10.1126/science.abe8770.

##  Appendix A Implementation

The presented method has been implemented in the wolfram language on
Mathematica 12.1 with Intel Core i7-10700 CPU at 2.90GHz. We deliver it as
user-friendly functions in the file Functions.nb available at
https://github.com/meijingyi/CheckQCTMC. The main functions are listed below.

  * •

ExpSS computes the matrix exponential, as a subroutine from MATLAB, using the
scaling and squaring method [AMH09].

  * •

GovernMat computes the governing matrix $\mathds{M}_{\Phi}$ which keeps the
behavior of $\Phi$-states and absorbs others.

  * •

Cyl gets the final partial density operator in Eq. (7) while computing the
probability measure along a given cylinder set.

  * •

Multiphase gets the final partial density operator in Eq. (9) while computing
the probability measure of a multiphase until formula.

For the function ExpSS, there is a useful tool MATLink available to transfer
data between Mathematica and MATLAB, thus we can use the embedded function
expm in MATLAB.

After specifying a quantum CTMC model $\mathfrak{Q}$ attached with the Hilbert
space and a multiphase until formula $\phi$, one can perform the model-
checking by calling the above functions respectively. The running examples are
validated in the file Appollonian Network.nb. Notably, using the scaling and
squaring method for matrix exponentiation, the time and the space consumption
of computing the probabilities of the cylinder set in Example 3 and the
multiphase until formula in Example 4 are summarized in Table 1.

Table 1. Performance using the scaling and squaring method method | Cyl | Multiphase | Overall
---|---|---|---
time | memory | time | memory | time | memory
ExpSS | 0.44s | 4.60MB | 4.80s | 36.54MB | 5.375s | 176.531MB