# Epsilon regularity for the Navier-Stokes equations via weak-strong
uniqueness

Dallas Albritton  Department of Mathematics, Princeton University, Princeton,
NJ 08544, USA<EMAIL_ADDRESS>,  Tobias Barker  Department of
Mathematical Sciences, University of Bath, Bath BA2 7AY, UK
<EMAIL_ADDRESS>and  Christophe Prange  Cergy Paris Université,
Laboratoire de Mathématiques AGM, UMR CNRS 8088, France
<EMAIL_ADDRESS>In honor of Olga Ladyženskaja

###### Abstract.

We give a new concise proof of a certain one-scale epsilon regularity
criterion using weak-strong uniqueness for solutions of the Navier-Stokes
equations with non-zero boundary conditions. It is inspired by an analogous
approach for the stationary system due to Struwe.

##  1\. Introduction

In the regularity theory of the three-dimensional non-stationary Navier-Stokes
equations, so-called epsilon regularity theory remains the state-of-the-art.
As well as being of independent interest, epsilon regularity at one scale
(‘one-scale epsilon regularity’) has been a crucial tool for proving some of
the best results in the field. See [3], [22], [27] and [5], for example. Thus
far, there have been three main methods for proving one-scale epsilon
regularity results:

  * •

direct iteration arguments (see [3]),

  * •

compactness arguments (see [21], [19]),

  * •

De Giorgi techniques (see [28]).

This short article is devoted to a new concise proof of a certain one-scale
epsilon regularity criterion for the three-dimensional non-stationary Navier-
Stokes equations away from boundaries. Theorem B is based on slicing
techniques and a comparison to a solution for which we have improved
integrability. To the best of our knowledge, related arguments in the Navier-
Stokes context were introduced by Struwe in [26] for the five-dimensional
stationary Navier-Stokes equations. Our paper is, in spirit, an analogue for
the non-stationary Navier-Stokes equations.

The following localized weak-strong uniqueness result is the keystone of our
paper.

######  Theorem A (first version of weak-strong uniqueness with boundary
conditions).

Let $\Omega\subset\mathbb{R}^{3}$ be a smooth bounded domain111We rely on
linear results of [10], which require $\partial\Omega$ to be of class
$C^{2,1}$. and $T>0$. There exists
$\bar{\kappa}=\bar{\kappa}_{\Omega,T}\in(0,\infty)$ and
$C(\Omega,T)\in(0,\infty)$ such that the following holds.222See (3.9) for an
estimate of $\bar{\kappa}_{\Omega,T}$. One can take
$C(\Omega,T):=K_{\Omega,T}(1+C_{0}C_{2})$. Here $C_{0}$ is defined in (3.4),
$C_{2}$ is defined in (3.7) and $K_{\Omega,T}\in(0,\infty)$ is defined in
(3.1). There exists a unique very weak solution

$U\in L^{4}(\Omega\times(0,T))$ (1.1)

to the Navier-Stokes equations on $\Omega\times(0,T)$ in the sense of
Definition 2.1 with boundary data $a$ and initial data $b$ satisfying the
integrability conditions

$a\in L^{4}(\partial\Omega\times(0,T)),\qquad b\in L^{4}(\Omega),$ (1.2)

the compatibility conditions

$\int\limits_{\partial\Omega}a(\cdot,t)\cdot
n=0,\qquad\int\limits_{\Omega}b\cdot\nabla q=0\quad\mbox{for all}\quad q\in
C^{\infty}(\mathbb{R}^{3})$ (1.3)

and the smallness condition

$\kappa:=\|a\|_{L^{4}(\partial\Omega\times(0,T))}+\|b\|_{L^{4}(\Omega)}<\bar{\kappa}.$
(1.4)

Moreover, $U\in L^{4}(0,T;L^{6}(\Omega))$, and

$\|U\|_{L^{4}(0,T;L^{6}(\Omega))}\leq C(\Omega,T)\kappa.$ (1.5)

The critical integrability in (1.5) is the main practical outcome of the
theorem.

Below, we give a second version of this localized uniqueness result, see
Theorem A’, that can be directly used to prove the epsilon regularity result
stated in Theorem B below.

######  Remark 1.1 (structure result).

We actually prove a strengthened version of the uniqueness result. Namely, we
prove that any very weak solution $U$ in the sense of Definition 2.1 with data
satisfying (1.2), (1.3) and the smallness condition (1.4) can be written as
$U=\bar{U}+V$, where $\bar{U}$ is the unique very weak solution to the Stokes
equation with boundary data $a$ and initial data $b$ and $V$ is the unique
mild solution to the perturbed Navier-Stokes equations around $\bar{U}$ in
$\Omega\times(0,T)$ with homogeneous boundary and initial data. We refer to
Section 3 below for more details.

######  Remark 1.2 (smallness assumption (1.4)).

The smallness condition (1.4) is required in Section 3 (Step 1-b) below in
order to construct a mild solution to the perturbed Navier-Stokes equations
around the solution $\bar{U}$ to the linear Stokes equations lifting the
boundary data $a$ and the initial data $b$.

######  Remark 1.3 (qualitative integrability assumption $L^{4}$).

The integrability condition (1.1) on the solution is essential to prove that
the perturbation in Section 3 (Step 2) has finite energy. The integrability
condition (1.2) on the boundary data $a$ is assumed in order to apply the
linear results of Farwig, Kozono and Sohr [10] (see Section 2 below). Notice
that according to [10], the boundary data can be taken in the larger space
$L^{4}(0,T;W^{-\frac{1}{6},6}(\partial\Omega))$. However, we do not state such
an improved result, because in the application to epsilon regularity that we
have in mind, the data $a$ satisfies (1.2).

In the second version of the localized uniqueness stated below, the main
change with respect to the previous theorem is in the conditions on the
initial data $b$. Indeed the requirement (1.3) is strong and imposes not only
that $b$ is divergence-free, but also that $b\cdot n=0$ on $\partial\Omega$.
This condition cannot be satisfied in general for data $b$ arising from
slicing in the proof of the epsilon regularity result, Theorem B below; for
details see Step 1 in Section 4. The assumptions in Theorem A’, contrary to
those of Theorem A, are immediately satisfied for the data stemming from the
slicing in the proof of Theorem B. In order to fit into the framework of
Theorem B, we take $\Omega=B_{r_{0}}$ for some
$r_{0}\in(\frac{5}{8},\frac{7}{8})$. This specific choice is only for
convenience. It is easy to extend the result to more general domains, the key
point being that one has to assume that $b$ is defined and divergence-free on
a larger domain than $\Omega$, so that it is possible to cut-off.

######  Theorem A’ (second version of weak-strong uniqueness with boundary
conditions).

Let $\Omega=B_{r_{0}}$ for some $r_{0}\in(\frac{5}{8},\frac{7}{8})$ and
$t_{0}\in(-1,-\frac{3}{4})$. There exists universal constants
$\bar{\kappa}\in(0,\infty)$ and $C\in(0,\infty)$ such that the following
holds.333See (3.14) for an estimate of $\bar{\kappa}$. One can take
$C=K(1+C_{0}C_{2}$). Here $C_{0}$ is defined in (3.4), $C_{2}$ is defined in
(3.7) and $K\in(0,\infty)$ is a universal constant (see (3.13)). There exists
a unique very weak solution

$U\in L^{4}(B_{r_{0}}\times(t_{0},0))$ (1.6)

to the Navier-Stokes equations on $B_{r_{0}}\times(t_{0},0)$ in the sense of
Definition 2.1 with boundary data $a$ and initial data $b$ satisfying the
integrability conditions

$a\in L^{4}(\partial B_{r_{0}}\times(t_{0},0)),\qquad b\in L^{4}(B_{1}),$
(1.7)

the compatibility condition

$\int\limits_{\partial B_{r_{0}}}a(\cdot,t)\cdot n=0,$ (1.8)

the incompressibility condition $\nabla\cdot b=0$ in $B_{1}$ in the sense of
distributions and the smallness condition

$\kappa:=\|a\|_{L^{4}(\partial
B_{r_{0}}\times(t_{0},0))}+\|b\|_{L^{4}(B_{1})}<\bar{\kappa}.$ (1.9)

Moreover, $U\in L^{4}(t_{0},0;L^{6}(B_{r_{0}}))$, and

$\|U\|_{L^{4}(t_{0},0;L^{6}(B_{r_{0}}))}\leq C\kappa.$ (1.10)

Again, the critical integrability is the main practical outcome of the
theorem.

We now state the epsilon regularity result. Let $(U,P)$ be a finite-energy
weak solution to the three-dimensional Navier-Stokes equations in
$Q_{1}=B_{1}(0)\times(-1,0)$ i.e.

$\displaystyle\partial_{t}U-\Delta U+U\cdot\nabla U+\nabla
P=0,\quad\nabla\cdot U=0\qquad\mbox{in}\quad Q_{1},$ in the sense of
distributions, and (1.11)
$\displaystyle\Bigg{(}\sup_{t\in(-1,0)}\int\limits_{B_{1}}|U(\cdot,t)|^{2}+\int\limits_{-1}^{0}\int\limits_{B_{1}}|\nabla
U|^{2}\Bigg{)}^{\frac{1}{2}}\leq M<\infty.$ (1.12)

######  Theorem B (epsilon regularity).

There exists $\bar{\varepsilon}\in(0,1)$ such that for any $M\in(0,\infty)$,
for any finite-energy weak solution $U$ to the Navier-Stokes equations in the
sense of (1.11)-(1.12) belonging to $C^{\infty}(B_{1}\times(-1,T))$ for all
$T\in(-1,0)$,555This assumption is there to keep technicalities to a minimum.
It makes our result applicable to rule out first-time singularities. the
following result holds.

Qualitative statement:

Assume that $U$ satisfies the smallness condition in Theorem A’, i.e.

$\|U\|_{L^{4}(Q_{1})}<\bar{\varepsilon}.$ (1.13)

Then $U\in L^{\infty}(Q_{\frac{1}{4}})$.

Quantitative statement:

Let $0<\varepsilon<\bar{\varepsilon}$. Assume that

$\|U\|_{L^{4}(Q_{1})}\leq\varepsilon.$ (1.14)

Then $U\in L^{\infty}(Q_{\frac{1}{4}})$ and in addition we have the
quantitative estimate

$\|U\|_{L^{\infty}(Q_{\frac{1}{4}})}\leq P(\varepsilon,M),$ (1.15)

where $P(\varepsilon,M)$ is a positive polynomial in $0<\varepsilon,\,M$.

This result is a mere corollary of Theorem A’.

######  Remark 1.4 (smallness condition (1.13)).

Ultimately, the reason for the space $L^{4}(Q_{1})$ in Theorem B is that the
solution of the linear Stokes equation with boundary data in
$L^{4}(\partial\Omega\times(0,T))$ and zero initial condition belongs to the
critical space $L^{4}(0,T;L^{6}(\Omega))$ (see Theorem 2.2). Notably, the
proof we present below apparently does not work in $L^{4-\varepsilon}(Q_{1})$.

######  Remark 1.5 (nonlocality and pressure).

Theorem B is in the spirit of some other epsilon regularity results that do
not involve the pressure, such as [29] and [18].

### Outline of the paper

In Section 2, we review the concept of very weak solution with $L^{p}$ initial
and boundary data. In Section 3, we prove weak-strong uniqueness, Theorems A
and A’. In Section 4, we prove the epsilon regularity criterion, Theorem B.

### Notations

For $r>0$ and $p\in(1,\infty)$, we denote by $\mathbf{A}$ the Stokes operator
realized on $L^{p}_{\sigma}(B_{r})$. Notice that
$\mathbf{A}=-\mathbb{P}\Delta_{D}$, where $\mathbb{P}$ is the Helmholtz-Leray
projection and $\Delta_{D}$ is the realization of the Laplace operator under
the Dirichlet boundary condition on $\partial B_{r}$. We also use the notation
$(e^{-t\mathbf{A}})_{t\in(0,\infty)}$ for the Stokes semigroup on $B_{r}$.
Notice that these notations do not involve explicitly the parameter $r$ in
order to lighten the notation and because the parameter $r$ will be fixed. For
further details concerning the Stokes operator and the Stokes semigroup, we
refer to [25, Chapter III.2 and IV.1] and [16]. For definitions of Sobolev
spaces on $\partial\Omega$, we refer to [25, Chapter I.3]. We define
$C_{0,\sigma}^{2}(\bar{\Omega}):=\\{f\in
C^{2}(\bar{\Omega};\mathbb{R}^{3}):\textrm{div}\,f=0,\,f|_{\partial\Omega}=0\\}.$

As is usual, the notation $C$ denotes a numerical constant possibly depending
on parameters that we do not track. This constant may change from line to
line. When a constant depends on parameters that we track, we denote this by
$C_{\alpha,\beta,\ldots}$.

##  2\. Preliminaries

Let $\Omega\subset\mathbb{R}^{3}$ be a smooth bounded domain.666Below, we
apply the results for the fixed smooth domain $\Omega=B_{r_{0}}$. Hence, we
pay no attention here to how the constants in the estimates will depend
quantitatively on the regularity of the domain. This dependence is not tracked
in [10]. See also Footnote 1. Let $-\infty<T_{1}<T_{2}<\infty$. Let $a\in
L^{1}(\partial\Omega\times(T_{1},T_{2});\mathbb{R}^{3})$, $b\in
L^{1}(\Omega;\mathbb{R}^{3})$ and $F\in
L^{1}(\Omega\times(T_{1},T_{2});\mathbb{R}^{3\times 3})$ satisfying the
compatibility conditions (1.3).

We consider the following Stokes problem:

$\displaystyle\begin{split}&\partial_{t}U-\Delta U+\nabla P=\nabla\cdot
F,\quad\nabla\cdot U=0\qquad\mbox{in}\quad\Omega\times(T_{1},T_{2}),\\\
&\bar{U}=a\qquad\mbox{on}\quad\partial\Omega\times(T_{1},T_{2}),\\\
&\bar{U}(\cdot,T_{1})=b.\end{split}$ (2.1)

This problem includes in particular the standard Navier-Stokes problem taking
$F=-U\otimes U$.

######  Definition 2.1 (very weak solution).

For $a,\,b$ and $F$ given as above, we say that $U\in
L^{1}(\Omega\times(T_{1},T_{2}))$ is a very weak solution of (2.1) if for all
$\Phi\in C^{1}_{0}([T_{1},T_{2});C_{0,\sigma}^{2}(\bar{\Omega})$ and for all
$q\in C^{\infty}(\overline{\Omega}\times[T_{1},T_{2}];\mathbb{R})$, we have

$\displaystyle\begin{split}&-\int\limits_{T_{1}}^{T_{2}}\int\limits_{\Omega}U\cdot(\partial_{t}\Phi+\Delta\Phi+\nabla
q)=-\int\limits_{T_{1}}^{T_{2}}\int\limits_{\Omega}F:\nabla\Phi\\\
&\qquad\qquad+\int\limits_{\Omega}b\cdot\Phi(\cdot,T_{1})-\int\limits_{T_{1}}^{T_{2}}\int\limits_{\partial\Omega}a\cdot(\nabla\Phi\cdot
n)-\int\limits_{T_{1}}^{T_{2}}\int\limits_{\partial\Omega}(a\cdot
n)q.\end{split}$ (2.2)

We state here an existence and uniqueness result of very weak solutions to the
Stokes system. Results in this direction were established in [9], though the
version we use below is from [10].

######  Theorem 2.2 (very weak solutions for the Stokes problem, [10, Lemma
1.2]).

Let $4\leq s,q<\infty$ such that $\frac{2}{s}+\frac{3}{q}=1$, and let
$r:=\frac{2}{3}q$. Let $a\in L^{s}(T_{1},T_{2};L^{r}(\Omega))$ and $F=b=0$.
Assume that $a$ satisfies the compatibility condition (1.3).

Then, there exists a unique very weak solution $U\in
L^{s}(T_{1},T_{2};L^{q}(\Omega))$ to the Stokes problem (2.1) in the sense of
Definition 2.1. Morever,

$\|U\|_{L^{s}(T_{1},T_{2};L^{q}(\Omega))}\leq
C_{\Omega,T_{1},T_{2},q}\|a\|_{L^{s}(T_{1},T_{2};L^{r}(\Omega))}.$ (2.3)

Notice that our boundary data $a$ is slightly less general (but general enough
for our purposes) than the data considered in [10, Lemma 1.2]. It follows that
any term in the very weak formulation of the equation (2.2) makes sense in a
classical integral sense and duality brackets are not needed. The statement of
Theorem 2.2 follows directly from [10, Lemma 1.2]. Indeed,
$W^{-\frac{1}{q},q}(\partial\Omega)$ is the dual space of
$W^{1-\frac{1}{q^{\prime}},q^{\prime}}(\partial\Omega)$ with
$\frac{1}{q}+\frac{1}{q^{\prime}}=1$, and by Sobolev embedding [4, Theorem
4.1.3], $W^{1-\frac{1}{q^{\prime}},q^{\prime}}(\partial\Omega)$ embeds into
$L^{r^{\prime}}(\partial\Omega)$ for $r=\frac{2}{3}q$ and
$\frac{1}{r}+\frac{1}{r^{\prime}}=1$. Therefore, $a\in
L^{s}(T_{1},T_{2};L^{r}(\partial\Omega))$ embeds into
$L^{s}(T_{1},T_{2};W^{-\frac{1}{q},q}(\partial\Omega))$.

######  Remark 2.3 (non-zero initial data).

We handle non-zero initial data $b$ as a separate issue. Indeed this point is
more classical than the case of rough boundary data which is treated in the
result above. For non-zero initial data, we rely on the Stokes semigroup
estimates in Lebesgue spaces, see for instance [16].

######  Remark 2.4 (on alternative linear results with rough boundary data).

It is also possible to rely on the results of Fabes, Lewis and Rivière for
boundary value problems with $L^{p}$ data obtained in [20] for the half-space
and in [8, 7] for bounded smooth domains. However, due to an additional
integrability condition on $a\cdot n$ on $\partial\Omega$ in [7, Theorem
(IV.3.3), page 643], these results require to work in $L^{4+\varepsilon}$,
$\varepsilon>0$, rather than $L^{4}$. Moreover, notice that there is a typo in
the statement of [7, Theorem (IV.3.1)]. The space for $a\cdot n$ is
$L^{\bar{q}}_{t}L^{\frac{2}{3}\bar{p}}_{x}$, not
$L^{\bar{q}}_{t}L^{\frac{3}{2}\bar{p}}_{x}$ as written in [7, Theorem IV.3.1].
Indeed, to get the estimate of the term involving the normal data, namely
$\nabla H$, one relies on Theorem IV.2.3.

######  Lemma 2.5 (Uniqueness of square integrable very weak solutions).

Let $U\in L^{2}(T_{1},T_{2};L^{2}(\Omega))$ be a very weak solution to the
Stokes problem (2.1), in the sense of Definition 2.1, with $F=a=b=0$. Then
$U\equiv 0$ on $\Omega\times(T_{1},T_{2})$.

###### Proof.

The proof relies on classical duality arguments, which we include for
completeness. Without loss of generality, let $T_{1}=-1$ and $T_{2}=0.$ Let
$f\in C^{\infty}_{0}(\Omega\times(-1,0);\mathbb{R}^{3})$ be arbitrary. Define
$\widetilde{f}(x,t)=f(x,-t)$ and let $\widetilde{\Phi}$ solve

$\partial_{t}\widetilde{\Phi}-\Delta\widetilde{\Phi}+\nabla\widetilde{q}=\widetilde{f}\quad\textrm{and}\quad\textrm{div}\,\widetilde{\Phi}=0\quad\textrm{in}\quad\Omega\times(0,\infty),\quad\widetilde{\Phi}|_{\partial{\Omega}}=0,\quad\widetilde{\Phi}(\cdot,0)=0.$

From [25, Theorem 2.7.3, IV], there is a classical solution to the above
linear problem satisfying $\widetilde{\Phi},\widetilde{q}\in
C^{\infty}(\overline{\Omega\times(\delta,2)})$ for all $\delta\in(0,2)$.
Furthermore, since $\widetilde{f}$ is supported in time away from zero we can
apply [25, Lemma 2.4.2, Chapter IV] to infer that $\widetilde{\Phi}$ is
supported in time away from zero.

We define

$\Phi(x,t)=-\widetilde{\Phi}(x,-t)\quad\textrm{and}\quad
q(x,t)=\widetilde{q}(x,-t).$

Inserting $\Phi$ and $q$ into (2.2) (where $F$, $a$ and $b$ are all zero)
gives

$\int\limits_{-1}^{0}\int\limits_{\Omega}U\cdot fdxds=0\qquad\forall f\in
C^{\infty}_{0}(\Omega\times(-1,0);\mathbb{R}^{3}).$

This implies the desired conclusion. ∎

######  Remark 2.6 (Very weak solutions and the energy equality).

Let $U\in L^{2}(T_{1},T_{2};L^{2}(\Omega))$ be a very weak solution to the
Stokes problem (2.1) in the sense of Definition 2.1, with $a=b=0$ and $F\in
L^{2}(\Omega\times(T_{1},T_{2});\mathbb{R}^{3\times 3})$. Applying [25,
Theorem 2.3.1 and Theorem 2.4.1, IV], together with Lemma 2.5, gives that
$U\in C([T_{1},T_{2}];L^{2}_{\sigma}(\Omega))\cap
L^{2}(T_{1},T_{2};W^{1,2}_{0,\sigma}(\Omega))$ with

$\|U(\cdot,t)\|_{L^{2}(\Omega)}^{2}+2\int\limits_{T_{1}}^{t}\int\limits_{\Omega}|\nabla
U|^{2}dxds=-2\int\limits_{T_{1}}^{t}\int\limits_{\Omega}F:\nabla
Udxds\quad\forall t\in[T_{1},T_{2}].$

This can also be used to show that very weak solutions to the Navier-Stokes
equations (with $a=0$) satisfy the energy equality. For results in this
direction, see [13].

##  3\. Proof of the localized weak-strong uniqueness results

###  3.1. Proof of Theorem A

Step 1: existence of a ‘strong’ solution to the Navier-Stokes system777After
completion of this paper, we became aware that in [9] the existence of a
‘strong solution’ with small non-zero boundary data is also shown. To make
this paper self-contained, we include such arguments in Step 1 below. in
$L^{4}(0,T;L^{6}(\Omega))$.

Step 1-a: linear problem.
Our goal is to lift the boundary data $a$ and the initial data $b$ satisfying
(1.2) and (1.3) by constructing a very weak solution $\bar{U}$ to the Stokes
problem (2.1) in $\Omega\times(0,T)$ with source $F=0$. By Theorem 2.2 of
Farwig, Kozono and Sohr, with $s=r=4$ and $q=6$, there exists a unique very
weak solution $U_{a}\in L^{4}(0,T;L^{6}(\Omega))$ in the sense of Definition
2.1 to the Stokes system with $F=b=0$. Hence we construct $\bar{U}$ the unique
very weak solution in the sense of Definition 2.1 to the Stokes system with
$F=0$ as follows:

$\bar{U}=U_{a}+e^{-t\mathbf{A}}b.$

Thanks to the semigroup estimates of [16], we have for $t\in(0,\infty)$,

$\|e^{-t\mathbf{A}}b\|_{L^{6}(\Omega)}\leq
C_{\Omega}t^{-\frac{1}{8}}\|b\|_{L^{4}(\Omega)}.$

Hence,

$\|\bar{U}\|_{L^{4}(0,T;L^{6}(\Omega))}\leq
K_{\Omega,T}\big{(}\|a\|_{L^{4}(\partial\Omega\times(t_{0},0))}+\|b\|_{L^{4}(\Omega)}\big{)},$
(3.1)

for $K_{\Omega,T}\in(0,\infty)$.

Step 1-b: existence of a unique strong solution $V$ to the perturbed system.
In this step, we construct a mild solution to the perturbed Navier-Stokes
system

$\displaystyle\begin{split}\partial_{t}V-\Delta V+\nabla
Q&=-\nabla\cdot(V\otimes
V)-\nabla\cdot(V\otimes\bar{U})-\nabla\cdot(\bar{U}\otimes
V)-\nabla\cdot(\bar{U}\otimes\bar{U}),\\\ \nabla\cdot
V&=0\qquad\qquad\qquad\qquad\qquad\qquad\mbox{in}\quad\Omega\times(0,T),\\\
V&=0\qquad\mbox{on}\quad\partial\Omega\times(0,T),\\\
V(\cdot,0)&=0\qquad\mbox{on}\quad\Omega.\end{split}$ (3.2)

for a small drift $\bar{U}$. We show that the Duhamel equation associated to
(3.2)

$V(\cdot,t)=-\int\limits_{0}^{t}e^{-(t-s)\mathbf{A}}\left(\nabla\cdot(V\otimes
V)+\nabla\cdot(V\otimes\bar{U})+\nabla\cdot(\bar{U}\otimes
V)+\nabla\cdot(\bar{U}\otimes\bar{U})\right)\,ds$ (3.3)

has a unique fixed point in the critical space $L^{4}(0,T;L^{6}(\Omega))$.

For $t\in(0,T)$, define

$B(D,E)(\cdot,t):=-\int\limits_{0}^{t}e^{-(t-s)\mathbf{A}}\nabla\cdot(D\otimes
E)\,ds.$

Using [17, Proposition 20], we see that

$\left\|e^{-(t-s)\mathbf{A}}\nabla\cdot(D\otimes
E)\right\|_{L^{6}(\Omega)}\leq\frac{C}{(t-s)^{\frac{3}{4}}}\|D(\cdot,t)\|_{L^{6}(\Omega)}\|E(\cdot,t)\|_{L^{6}(\Omega)},$

and hence, by Hardy-Littlewood-Sobolev’s theorem [15, Theorem 7.25],

$\left\|B(D,E)\right\|_{L^{4}(0,T;L^{6}(\Omega))}\leq
C_{0}\|D\|_{L^{4}(0,T;L^{6}(\Omega))}\|E\|_{L^{4}(0,T;L^{6}(\Omega))}.$ (3.4)

Define the linear operator

$L(D)(\cdot,t):=-\int\limits_{0}^{t}e^{-(t-s)\mathbf{A}}\nabla\cdot(D\otimes\bar{U}+\bar{U}\otimes
D)\,ds.$

Then by the same reasoning as above

$\left\|L(D)\right\|_{L^{4}(0,T;L^{6}(\Omega))}\leq
C_{1}\|D\|_{L^{4}(0,T;L^{6}(\Omega))}\|\bar{U}\|_{L^{4}(0,T;L^{6}(\Omega))},$
(3.5)

and for the source term quadratic in $\bar{U}$,

$\left\|B(\bar{U},\bar{U})\right\|_{L^{4}(0,T;L^{6}(\Omega))}\leq
C_{0}\|\bar{U}\|_{L^{4}(0,T;L^{6}(\Omega))}^{2}.$ (3.6)

Let

$C_{2}:=\min\Big{(}\frac{1}{C_{0}},\frac{1}{C_{1}}\Big{)}.$ (3.7)

Using (3.4)-(3.6), we can apply [14, Lemma 4.1]. This gives the existence of a
fixed point/strong solution $V\in L^{4}(0,T;L^{6}(\Omega))$ provided that

$\|\bar{U}\|_{L^{4}(0,T;L^{6}(\Omega))}<\frac{C_{2}}{4}.$ (3.8)

Moreover,

$\|V\|_{L^{4}(0,T;L^{6}(\Omega))}\leq
4C_{0}\|\bar{U}\|^{2}_{L^{4}(0,T;L^{6}(\Omega))}.$

In view of the linear estimate (3.1), this is achieved whenever

$\displaystyle\|\bar{U}\|_{L^{4}(0,T;L^{6}(\Omega))}\leq\ $ $\displaystyle
K_{\Omega,T}\big{(}\|a\|_{L^{4}(0,T;L^{4}(\Omega))}+\|b\|_{L^{4}(\Omega)}\big{)}$
$\displaystyle=\ $ $\displaystyle K_{\Omega,T}\kappa$ $\displaystyle<\ $
$\displaystyle K_{\Omega,T}\bar{\kappa}=\frac{C_{2}}{4},$

i.e.

$\bar{\kappa}:=\frac{C_{2}}{4K_{\Omega,T}}.$ (3.9)

Therefore, there is a strong solution such that

$\|V\|_{L^{4}(0,T;L^{6}(\Omega))}<4C_{0}K_{\Omega,T}^{2}\kappa^{2}=C_{0}C_{2}K_{\Omega,T}\kappa.$
(3.10)

The fact that $V$ is the only strong solution in $L^{4}(0,t;L^{6}(\Omega))$
follows from (3.4) and arguments in [6, Theorem (3.3)].

Step 2: weak-strong uniqueness result.

Let $W:=U-\bar{U}\in L^{4}(\Omega\times(0,T))$. Let $V\in
L^{4}(0,T;L^{6}(\Omega))$ be the strong solution constructed in Step 1-b
above. Then $V-W$ is a very weak solution to the Stokes system with $a=b=0$
and $F:=U\otimes U-(V+\bar{U})\otimes(V+\bar{U})\in L^{2}(\Omega\times(0,T))$.
Applying Remark 2.6 gives that $V-W$ has finite energy on $\Omega\times(0,T)$
with zero initial data and satisfies the energy equality

$\frac{1}{2}\int\limits_{\Omega}|V-W|^{2}(x,t)dx+\int\limits_{0}^{t}\int\limits_{\Omega}|\nabla(V-W)|^{2}dxds\\\
=\int\limits_{0}^{t}\int\limits_{\Omega}(V-W)\otimes(V+\bar{U}):\nabla(V-W)dxds\quad\forall
t\in[0,T).$

For $t\in[0,T)$, we define

$\displaystyle\mathcal{E}(t):=\sup_{s\in[0,t]}\frac{1}{2}\int\limits_{\Omega}|V-W|^{2}(x,s)dx+\int\limits_{0}^{t}\int\limits_{\Omega}|\nabla(V-W)|^{2}(x,s)dxds.$

Using the above energy equality (combined with Hölder’s inequality,
interpolation of Lebesgue spaces, the Sobolev embedding theorem and Young’s
inequality) yields that for $t\in[0,T)$ and some positive universal constant
$C$,

$\mathcal{E}(t)\leq
C\mathcal{E}(t)\int\limits_{0}^{t}\|(V+\bar{U})(\cdot,s)\|_{L^{6}(\Omega)}^{4}ds.$

Using this, together with the fact that $V+\bar{U}$ is in
$L^{4}(0,T;L^{6}(\Omega))$, allows us to use an absorbing argument to conclude
that $\mathcal{E}(t)=0$ for all $t\in(0,T)$. This concludes the proof of
Theorem A and of the structure result mentioned in Remark 1.1. Notice that the
quantitative estimate (1.5) directly follows from the linear estimate (3.1),
the estimate (3.10) for the mild solution and the definition of $\bar{\kappa}$
in (3.9).

###  3.2. Proof of Theorem A’

The proof of this result only differs from the proof of Theorem A in the
treatment of the linear evolution of the initial data $b$. Let us outline the
changes, which concern only Step 1-a of Section 3.1.

We take $\varphi\in C^{\infty}_{c}(B_{\frac{15}{16}})$ a cut-off function such
that $\varphi\equiv 1$ on $B_{\frac{7}{8}}$, $0\leq\varphi\leq 1$ and
$|\nabla\varphi|\leq 32$. Thanks to the Bogovskii operator [12] on the annulus
$B_{\frac{15}{16}}\setminus B_{\frac{7}{8}}$, there exists a divergence-free
extension $E(b)$ of $\varphi b$ defined on $\mathbb{R}^{3}$ that is compactly
supported in $B_{\frac{15}{16}}$, $E(b)=b$ on $B_{\frac{7}{8}}\subset
B_{r_{0}}$ and such that

$\|E(b)\|_{L^{4}(B_{1})}\leq C\|b\|_{L^{4}(B_{1})}.$

Let $\Gamma$ be the heat kernel on $\mathbb{R}^{3}$. We note that
$\Gamma(\cdot-t_{0})\star E(b)$ is divergence-free, which implies

$\int\limits_{\partial B_{r_{0}}}\big{(}\Gamma(\cdot-t_{0})\star
E(b)\big{)}\cdot n=0.$

Moreover, it follows from [7, Lemma IV.3.2] that

$\big{\|}\Gamma(\cdot-t_{0})\star E(b)\big{\|}_{L^{4}(\partial
B_{r_{0}}\times(t_{0},0))}\leq
C\big{(}(-t_{0})^{\frac{1}{8}}+(-t_{0})^{\frac{1}{4}}\big{)}\|E(b)\|_{L^{4}(\mathbb{R}^{3})}\leq
C\|b\|_{L^{4}(B_{1})},$ (3.11)

and

$\big{\|}\Gamma(\cdot-t_{0})\star
E(b)\big{\|}_{L^{4}(t_{0},0;L^{6}(B_{r_{0}}))}\leq
C(-t_{0})^{\frac{1}{8}}\|E(b)\|_{L^{4}(\mathbb{R}^{3})}\leq
C\|b\|_{L^{4}(B_{1})},$ (3.12)

where $C\in(0,\infty)$ denotes as usual a universal constant. We can now
construct the linear solution $\bar{U}$ as follows. Let

$\widetilde{a}=a-\Gamma(\cdot-t_{0})\star E(b)|_{\partial B_{r_{0}}}.$

By Theorem 2.2 of Farwig, Kozono and Sohr [10], with $s=r=4$ and $q=6$, there
exists a unique very weak solution $U_{\widetilde{a}}\in
L^{4}(t_{0},0;L^{6}(B_{r_{0}}))$ in the sense of Definition 2.1 to the Stokes
system with $F=b=0$ and boundary data $\widetilde{a}$. Hence we construct
$\bar{U}$ the unique very weak solution in the sense of Definition 2.1 to the
Stokes system with $F=0$ as follows:

$\bar{U}=U_{\widetilde{a}}+\Gamma(\cdot-t_{0})\star E(b).$

Hence, the estimates (2.3) and (3.12) lead to

$\|\bar{U}\|_{L^{4}(t_{0},0;L^{6}(B_{r_{0}}))}\leq
K\big{(}\|a\|_{L^{4}(\partial
B_{r_{0}}\times(t_{0},0))}+\|b\|_{L^{4}(B_{r_{0}})}\big{)},$ (3.13)

where $K\in(0,\infty)$ denotes a universal constant.

The rest of the proof, Step 1-b to Step 2 of Section 3.1 are identical,
replacing the constant $K_{\Omega,T}$ by $K$. The definition of $\bar{\kappa}$
becomes

$\bar{\kappa}:=\frac{C_{2}}{4K}.$ (3.14)

######  Remark 3.1 (on the compatibility conditions).

We emphasize that in the compatibility conditions (1.3), the condition888We
stress that this condition on $a$ comes from the linear result Theorem 2.2
taken from [10]; these conditions are also present in the work [7]. on the
boundary data $a$ is much weaker than the condition on $b$. Indeed we only ask
that $a\cdot n$ has mean zero on $\partial\Omega$, while $b\cdot n$ is
required to vanish identically on $\partial\Omega$. This fact is the essential
redeeming feature that allows the above argument to work. Notice that owing to
the fact that $\Gamma(\cdot-t_{0})\star E(b)$ is divergence-free,
$\big{(}\Gamma(\cdot-t_{0})\star E(b)\big{)}\cdot n$ has mean zero on
$\partial B_{r_{0}}$ but is not necessarily zero identically.

##  4\. Proof of the epsilon regularity result

This section is devoted to the proof of Theorem B. We directly prove the
quantitative version, i.e. the bound (1.15). We assume that $U$ is a finite-
energy weak solution to the Navier-Stokes equations, i.e. (1.11) and (1.12)
hold. In addition, we assume that $U$ belongs to
$C^{\infty}(B_{1}\times(-1,T))$ for all $T\in(-1,0)$, see Footnote 5, and
satisfies the assumption (1.14) with $0<\varepsilon\leq\bar{\varepsilon}$ and
$\bar{\varepsilon}$ defined in (4.4).

Step 1: finding good space and time scales.
We need to select a space slice and a time slice. The choices are completely
independent, so the order in which we select the slices does not matter. Let
us first select a space slice $r_{0}\in(\frac{5}{8},\frac{7}{8})$. By the
coarea formula and Fubini’s theorem, we have

$\int\limits_{-1}^{0}\int\limits_{B_{1}}|U|^{4}\,dxdt=\int\limits_{-1}^{0}\int\limits_{0}^{1}\int\limits_{\partial
B_{r}}|U|^{4}\,dS_{r}drdt=\int\limits_{0}^{1}\int\limits_{-1}^{0}\int\limits_{\partial
B_{r}}|U|^{4}\,dS_{r}dtdr.$

Therefore, the pigeonhole principle implies that there exists
$r_{0}\in(\frac{5}{8},\frac{7}{8})$ such that

$\int\limits_{-1}^{0}\int\limits_{\partial B_{r_{0}}}|U|^{4}\,dS_{r_{0}}dt\leq
4\int\limits_{-1}^{0}\int\limits_{B_{1}}|U|^{4}\,dxdt.$ (4.1)

We now select a time slice. By the pigeonhole principle, there exists
$t_{0}\in(-1,-\frac{3}{4})$ such that

$\int\limits_{B_{1}}|U|^{4}(x,t_{0})\,dx\leq
4\int\limits_{-1}^{0}\int\limits_{B_{1}}|U|^{4}\,dxdt.$ (4.2)

From now on, we call $a:=U|_{\partial B_{r_{0}}\times(-1,0)}$ and
$b:=U(\cdot,t_{0})$. By (4.1) we have $a\in L^{4}(\partial
B_{r_{0}}\times(-1,0))$ and by (4.2) we have $b\in L^{4}(B_{1})$. Moreover,

$\|a\|_{L^{4}(\partial B_{r_{0}}\times(t_{0},0))}+\|b\|_{L^{4}(B_{1})}\leq
2\sqrt{2}\|U\|_{L^{4}(Q_{1})}\leq
2\sqrt{2}\varepsilon<2\sqrt{2}\bar{\varepsilon}\leq\bar{\kappa},$

with $\bar{\kappa}$ defined by (3.14), on condition that
$\bar{\varepsilon}\leq\frac{\bar{\kappa}}{2\sqrt{2}}$.

Step 2: applying weak-strong uniqueness.
We now apply the weak-strong uniqueness result of Theorem A’ to the solution
$U$ on $B_{r_{0}}\times(t_{0},0)$ with the data $a$ and $b$ from Step 1 above.
By the quantitative estimate (1.5) it follows that

$\|U\|_{L^{4}(t_{0},0;L^{6}(B_{r_{0}}))}\leq 2\sqrt{2}\varepsilon
K(1+C_{0}C_{2}).$ (4.3)

Step 3: conclusion via Ladyženskaja-Prodi-Serrin.
We can now argue in a similar spirit to [23, 2], except we use the bound (4.3)
to set up a contraction mapping999For such a contraction mapping, see for
instance [2, Lemma 12]. related to the localized vorticity equation rather
than the velocity. In particular, there exists $\bar{\eta}$ such that for any

$0<\varepsilon<\frac{\bar{\eta}}{2\sqrt{2}K(1+C_{0}C_{2})}=:\bar{\varepsilon}$
(4.4)

we have

$\|U\|_{L^{6}(Q_{\frac{1}{2}})}\leq M(1+\varepsilon).$

With this, we can bootstrap in the same way as in [24] to obtain the estimate
(1.15). This concludes the proof of Theorem B.

######  Remark 4.1 (on the linear solution).

In this proof, notice that we are not able to assert that the solution
$\bar{U}$ to the linear Stokes problem with rough boundary data $a\in
L^{4}(B_{r_{0}}\times(t_{0},0))$ (see the structure result in Remark 1.1) is
smooth in space. However, the linear result is pivotal in order to establish
that the original solution $U$ has improved critical integrability, hence is
smooth in space. Notice that we obtain $L^{\infty}$ time integrability by
appealing to the $L^{\infty}_{t}L^{2}_{x}$ control granted by being finite-
energy. We cannot bootstrap further in time without controlling the pressure
$p$ up to time 0.

######  Remark 4.2 (on the half-space).

It remains an open problem as to whether such a proof can be done for
establishing epsilon regularity results near a boundary. Indeed, the linear
results of Farwig, Kozono and Sohr [10] and of Fabes, Lewis and Rivière [7]
ask for smoothness of the domain $\Omega$, see for instance Theorem 2.2 above.
However, we are unable to carry out the slicing procedure of Step 1 above near
a smooth boundary.

######  Remark 4.3 (comparison with Struwe [26]).

Struwe’s approach to the five-dimensional stationary Navier-Stokes equations
is technically different; it involves iterating the energy on a sequence of
balls, which we avoid. A version of our proof should work in the stationary
setting of [26] as well. There, the smallness condition in $L^{4}$ could be
replaced by smallness in $H^{1}$, since in five dimensions, the slicing
procedure yields boundary data in $H^{1}(\partial B_{r_{0}})\subset
L^{4}(\partial B_{r_{0}})$ (four-dimensional Sobolev embedding). The slicing
technique has proven useful for the six-dimensional stationary Navier-Stokes
equations [11] and advection-diffusion equation with rough drifts [1], among
others.

### Acknowledgement

DA was supported by NSF Postdoctoral Fellowship Grant No. 2002023. DA is also
grateful to ENS Paris for supporting his academic visit to Paris during which
this research was initiated. DA thanks Vladimír Šverák, who contributed an
offhanded comment in 2017 which played a role in the genesis of this paper. CP
is partially supported by the Agence Nationale de la Recherche, project BORDS,
grant ANR-16-CE40-0027-01, project SINGFLOWS, grant ANR- 18-CE40-0027-01,
project CRISIS, grant ANR-20-CE40-0020-01, by the CY Initiative of Excellence,
project CYNA (CY Nonlinear Analysis) and project CYFI (CYngular Fluids and
Interfaces).

## References

  * [1] D. Albritton and H. Dong.  Regularity properties of passive scalars with rough divergence-free drifts.  arXiv preprint arXiv:2107.12511, 2021.
  * [2] T. Barker and C. Prange.  Localized Smoothing for the Navier–Stokes Equations and Concentration of Critical Norms Near Singularities.  Arch. Ration. Mech. Anal., 236(3):1487–1541, 2020.
  * [3] L. Caffarelli, R. Kohn, and L. Nirenberg.  Partial regularity of suitable weak solutions of the Navier-Stokes equations.  Comm. Pure Appl. Math., 35(6):771–831, 1982.
  * [4] J. Droniou.  Quelques Résultats sur les Espaces de Sobolev.  working paper or preprint, Apr. 2001.
  * [5] L. Escauriaza, G. A. Seregin, and V. Šverák.  $L_{3,\infty}$-solutions of Navier-Stokes equations and backward uniqueness.  Uspekhi Mat. Nauk, 58(2(350)):3–44, 2003.
  * [6] E. B. Fabes, B. F. Jones, and N. M. Rivière.  The initial value problem for the Navier-Stokes equations with data in $L^{p}$.  Arch. Rational Mech. Anal., 45:222–240, 1972.
  * [7] E. B. Fabes, J. E. Lewis, and N. M. Rivière.  Boundary value problems for the Navier-Stokes equations.  Amer. J. Math., 99(3):626–668, 1977.
  * [8] E. B. Fabes, J. E. Lewis, and N. M. Rivière.  Singular integrals and hydrodynamic potentials.  Amer. J. Math., 99(3):601–625, 1977.
  * [9] R. Farwig, G. P. Galdi, and H. Sohr.  A new class of weak solutions of the Navier-Stokes equations with nonhomogeneous data.  J. Math. Fluid Mech., 8(3):423–444, 2006.
  * [10] R. Farwig, H. Kozono, and H. Sohr.  Global weak solutions of the Navier-Stokes equations with nonhomogeneous boundary data and divergence.  Rend. Semin. Mat. Univ. Padova, 125:51–70, 2011.
  * [11] J. Frehse and M. Ružička.  Existence of regular solutions to the steady Navier-Stokes equations in bounded six-dimensional domains.  Annali della Scuola Normale Superiore di Pisa-Classe di Scienze, 23(4):701–719, 1996.
  * [12] G. P. Galdi.  An introduction to the mathematical theory of the Navier-Stokes equations.  Springer Monographs in Mathematics. Springer, New York, second edition, 2011.  Steady-state problems.
  * [13] G. P. Galdi.  On the energy equality for distributional solutions to Navier-Stokes equations.  Proc. Amer. Math. Soc., 147(2):785–792, 2019.
  * [14] I. Gallagher, D. Iftimie, and F. Planchon.  Asymptotics and stability for global solutions to the Navier-Stokes equations.  Ann. Inst. Fourier (Grenoble), 53(5):1387–1424, 2003.
  * [15] M. Giaquinta and L. Martinazzi.  An introduction to the regularity theory for elliptic systems, harmonic maps and minimal graphs, volume 11 of Appunti. Scuola Normale Superiore di Pisa (Nuova Serie) [Lecture Notes. Scuola Normale Superiore di Pisa (New Series)].  Edizioni della Normale, Pisa, second edition, 2012.
  * [16] Y. Giga.  Solutions for semilinear parabolic equations in $L^{p}$ and regularity of weak solutions of the Navier-Stokes system.  J. Differential Equations, 62(2):186–212, 1986.
  * [17] M. Hieber and J. Saal.  The Stokes equation in the $L^{p}$-setting: well-posedness and regularity properties.  In Handbook of mathematical analysis in mechanics of viscous fluids, pages 117–206. Springer, Cham, 2018.
  * [18] H. Kwon.  The role of the pressure in the regularity theory for the Navier-Stokes equations.  arXiv preprint arXiv:2104.03160, 2021.
  * [19] O. A. Ladyzhenskaya and G. A. Seregin.  On partial regularity of suitable weak solutions to the three-dimensional navier—stokes equations.  Journal of Mathematical Fluid Mechanics, 1(4):356–387, 1999.
  * [20] J. E. Lewis.  The initial-boundary value problem for the Navier-Stokes equations with data in $L^{p}$.  Indiana Univ. Math. J., 22:739–761, 1972/73.
  * [21] F. Lin.  A new proof of the Caffarelli-Kohn-Nirenberg theorem.  Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 51(3):241–257, 1998.
  * [22] J. Nečas, M. Růžička, and V. Šverák.  On Leray’s self-similar solutions of the Navier-Stokes equations.  Acta Math., 176(2):283–294, 1996.
  * [23] G. A. Seregin.  A new version of the Ladyzhenskaya-Prodi-Serrin condition.  Algebra i Analiz, 18(1):124–143, 2006.
  * [24] J. Serrin.  On the interior regularity of weak solutions of the Navier-Stokes equations.  Arch. Rational Mech. Anal., 9:187–195, 1962.
  * [25] H. Sohr.  The Navier-Stokes equations.  Modern Birkhäuser Classics. Birkhäuser/Springer Basel AG, Basel, 2001.  An elementary functional analytic approach, [2013 reprint of the 2001 original] [MR1928881].
  * [26] M. Struwe.  On partial regularity results for the Navier-Stokes equations.  Communications on Pure and Applied Mathematics, 41(4):437–458, 1988\.
  * [27] T.-P. Tsai.  On Leray’s self-similar solutions of the Navier-Stokes equations satisfying local energy estimates.  Archive for rational mechanics and analysis, 143(1):29–51, 1998\.
  * [28] A. F. Vasseur.  A new proof of partial regularity of solutions to Navier-Stokes equations.  Nonlinear Differential Equations and Applications NoDEA, 14(5):753–785, 2007.
  * [29] J. Wolf.  On the local regularity of suitable weak solutions to the generalized Navier-Stokes equations.  Ann. Univ. Ferrara Sez. VII Sci. Mat., 61(1):149–171, 2015.


# QTM: computational package using MPI protocol for quantum trajectories
method

Marek Sawerwain  Institute of Control & Computation Engineering University of
Zielona Góra, Licealna 9, Zielona Góra 65-417, Poland;
<EMAIL_ADDRESS>Joanna Wiśniewska  Institute of Information
Systems, Faculty of Cybernetics, Military University of Technology,
Urbanowicza 2, 00-908 Warsaw, Poland<EMAIL_ADDRESS>
###### Abstract

The Quantum Trajectories Method (QTM) is one of the frequently used methods
for studying open quantum systems. The main idea of this method is the
evolution of wave functions which describe the system (as functions of time).
Then, so-called quantum jumps are applied at a randomly selected point in
time. The obtained system state is called as a trajectory. After averaging
many single trajectories, we obtain the approximation of the behavior of a
quantum system. This fact also allows us to use parallel computation methods.
In the article, we discuss the QTM package which is supported by the MPI
technology. Using MPI allowed utilizing the parallel computing for calculating
the trajectories and averaging them – as the effect of these actions, the time
taken by calculations is shorter. In spite of using the C++ programming
language, the presented solution is easy to utilize and does not need any
advanced programming techniques. At the same time, it offers a higher
performance than other packages realizing the QTM. It is especially important
in the case of harder computational tasks, and the use of MPI allows improving
the performance of particular problems which can be solved in the field of
open quantum systems.

##  1 Introduction

The Quantum Trajectories Method (QTM) is an important method actively applied
for investigation in the field of open quantum systems [30], [15], [3], [29].
It was implemented as packages in a few programming languages and tools. The
first implementation was a package written by Sze M. Tan for the Matlab
environment [25]. There is also a package for the C++ language [23], [20]. The
latest implementations are QuTIP [11], [12] for the Python Ecosystem and
QuantumOptics.jl [13] prepared for the Julia language.

The aforementioned solutions, especially QuTIP and QuantumOptics.jl, allow
utilizing parallel computing inside the environments of Python and Julia.
However, the packages are not intended for High-Performance Computing (HPC)
[17], [26] where Message Process Interface, termed as MPI [2], plays a
significant role. Using MPI allows us to re-implement the QTM for HPC systems,
regardless of their scale (the QTM is scale-free because each trajectory may
be calculated separately). The scale-free character of the QTM will allow
utilizing more computing power, and that will result in shorter time of
calculations, which is especially important for cases in which 20,000 and more
trajectories are generated.

Previously, we prepared the implementations of the QTM for CPUs and GPUs [27],
[28]. The results presented in this work relate to a brand new implementation
of the QTM for the MPI protocol (the actual source code of the QTM can be
found at [21]). The version 1.0a of the QTM package is also available at [22].

As far as the QTM is concerned, a proper selection method for solving systems
of Ordinary Differential Equations (ODEs) must be considered. This is a basic
action taken while numerical computations are realized. More precisely, a very
important issue is an Initial Value Problem (IVP). Due to that, the selection
of an appropriate method for solving IVP (in general ODEs) is crucial.
Especially when the system of equations is difficult to solve, e.g. so-called
systems of stiff ODEs. The stiff ODEs constitute a significant type of ODEs
and the correct solving of these equations is pivotal for numerical
computations in many cases – especially for the QTM where calculating many
trajectories requires solving many ODEs.

A few groups of methods [8], [9] used in the context of stiff ODEs may be
recalled, and these are the Backward Differentiation Formula (BDF) methods.
Another approach, which is commonly used in numerical computations for solving
ODEs, is Livermore Solver for ODE (LSODE) method [10]. In this work, we reuse
the LSODE variant called ZVODE (LSODE for complex numbers) for the
implementation of the QTM in the MPI environment.

The ZVODE package is a stable solution offering high accuracy, but this is not
a reentrant solution which may be directly utilized in a parallel environment.
However, this problem can be solved by utilizing MPI where the processes are
separate programs communicating with one another using message passing.
Therefore, one MPI process calls only one instance of the ZVODE method.

The ZVODE package was implemented in the Fortran language. The QTM package,
implemented in the C++ language, uses the ZVODE package efficiently. However,
we made some effort to assure that functions offered by the QTM package are
easily accessible as it is in the packages QuTIP and QuantumOptics.jl.
Naturally, programs (delivered as examples or written by the user) utilizing
the QTM have to be compiled. This should not be a problem because of used
makefile mechanism for simplifying this process.

The paper is organized in the following way: in Section 2, we shortly present
selected mathematical features of the QTM. The BDF numerical methods for
solving IVP, used in the presented package are discussed in subsection 2.1.
Whereas in 2.2 the algorithm for the QTM is presented. There are also some
remarks pointing out where the parallel processing techniques and the MPI
technology are utilized.

Section 3 contains selected remarks referring to the implementation of the QTM
package. In this section, the most important data types implemented in package
are presented. We also describe how the ZVODE package is used to solve ODEs
since solving ODEs poses a significant problem in the QTM.

We analyze the efficiency of our implementation in a comparison with other
recently developed packages in section 4. The most important issue presented
in this section is the scalability of the QTM.

A summary of achieved results is discussed in Section 5. There are also
presented further aims which are planned to be realized as the next steps in
the evolution of the demonstrated implementation of the QTM. The article is
ended with acknowledgments and bibliography.

##  2 Quantum Trajectories Method

Performing the QTM requires solving ODEs. The presented solution, as mentioned
above, uses the ZVODE package which allows utilizing some numerical methods
indispensable for the proper functioning of the QTM. The fundamental
information concerning the ZVODE package is presented in subsection 2.1, while
subsection 2.2 contains the description of the QTM, including a method of
calculating a single trajectory [27], [28] with the use of MPI technology. The
properties of the QTM allow us to easily distribute the computations – this
feature helps to accelerate the calculations and makes the process scalable.

Another important element of the solution is a pseudorandom number generator.
In subsection 2.3, we describe the generator which we chose to utilize in the
implementation of the QTM package.

###  2.1 Methods for solving ODEs

The ZVODE package was built on a basis of the VODE package [1]. The VODE
package was based on the LSODE package. The numerical methods, used for the
LSODE (Livermore Solver for Ordinary Differential Equations) package
implementation, are the computational routines based on the group of Linear
Multistep Methods (LMMs). In general, LSODE uses the Adams methods (so-called
predictor-corrector methods) for the non-stiff ODEs solving. In case of stiff
ODEs, the Backward Differentiation Formula (BDF) methods are used. In this
section, the most important mathematical aspects of the aforementioned methods
are briefly presented – the detailed report concerning standard (i.e.
sequential) implementation of LSODE may be found in [19].

The choice of the ZVODE package is determined by the fact that the QTM needs
using complex numbers and the ZDOVE package allows performing calculations on
complex numbers of single and double precision.

For the Initial Value Problem (IVP) in general:

$\begin{array}[]{lcl}y^{\prime}&=&f(x,y),\\\ y(x_{0})&=&y_{0}.\end{array}$ (1)

The LMMs, approximating the problem’s solution, may be described as:

$\sum_{i=0}^{k}\alpha_{(k-i)}y_{(n+1-i)}=h\sum_{i=0}^{k}\beta_{(k-i)}f(x_{(n+1-i)},y_{(n+1-i)}),$
(2)

where $\alpha_{i},\beta_{i}\in\mathbb{R}$ and $\alpha_{k}\neq 0$ and
$|\alpha_{0}|+|\beta_{0}|>0$. The parameter $h$ represents the step width for
integration. Naturally, a value of $h$ is selected during the solver’s work
with the use of adaptive methods.

The LSODE routine uses two methods based on the LMMs: the Adams-Moulton method
and the BDF method.

The implicit Adams-Moulton method, this means the values of $\beta_{i}\neq 0$
(the detailed description of used symbols may be found in [24]), may be
presented as:

$y_{n+1}=y_{n}+h\sum_{i=0}^{k}{\gamma_{i}}\triangledown^{i}f_{n+1}$ (3)

where

${\gamma_{i}}={(-1)}^{i}\int_{0}^{1}{-s+1\choose i}ds$ (4)

The $\triangledown^{i}$ symbol stands for the backward differences. It should
be pointed out that: $\triangledown^{i}f_{n}=f_{n}$ and
$\triangledown^{i+1}f_{n}=\triangledown^{i}f_{n}-\triangledown^{i}f_{n-1}$.
The value $k$ expresses the number of steps, i.e. the values of $y_{i}$ used
in the specified method.

The BDF method is defined as:

$\sum_{i=1}^{k}\frac{1}{i}\triangledown^{i}y_{n+1}=hf(x_{(n+1)},y_{(n+1)})$
(5)

The above notation means that $\beta_{k}\neq 0$ but $\beta_{i}=0$ for
$i=0,1,2,\ldots,k-1$. The values of $\alpha_{i}$ are arbitrarily defined and
may be found in many publications concerning the methods of solving ODEs.

In both cases, for the Adams-Moulton and the BDF method, the problem is to
estimate the value of $f_{n+1}$. It should be stressed that this value is
needed to estimate itself. The Newton method [18] is used in this case to
calculate the estimation. The Newton method is quite rapidly convergent,
though for the system of equations it needs to calculate the Jacobian. This
method for (ZV/LS)ODE may be expressed as:

$g(y^{(1)}_{n+1})=-\left(I-\frac{\partial f}{\partial
f}\right)(y^{(2)}_{n+1}-y^{(1)}_{n+1})$ (6)

where $y^{(2)}_{n+1}$ and $y^{(1)}_{n+1}$ stand for the next approximations of
$y_{n+1}$. The function $g(\cdot)$ represents the function approximating
values of $y_{n+1}$.

Figure 1: A general overview of the ODE solver and the methods used in the
(ZV/LS)ODE implementation

The BDF method is used for stiff ODE problems and the Adams-Moulton method for
non-stiff ODE problems. Using these both approaches together makes the
(ZV/LS)ODE a hybrid method – see Fig. 1. The presented QTM implementation
allows choosing the method, respectively to solving easier (non-stiff ODE) or
more difficult (stiff ODE) QTM problems. Moreover, in ZVODE, for every method
the method’s order can also be controlled automatically. For the Adams-Moulton
method, the orders from the first to the twelfth are available, and for the
BDF method the orders from the first to the sixth.

###  2.2 Quantum Trajectories Method

A description of quantum states’ dynamics may be presented for two fundamental
situations. The first case is an evolution of a closed quantum system, whereas
the second case concerns an evolution of an open quantum system. As an example
of a closed quantum system, we may refer to a model of a quantum circuit. If
we want to consider a quantum system where its dynamics is affected by the
influence of an external environment, we deal with an open quantum system.

In this subsection, we do not aim to describe the mathematical models of
dynamics in open and closed quantum systems – the details concerning this
subject may be found in [29] and [16]. For clarity, it should be recalled that
for closed systems, the evolution is a unitary operation, and it can be
denoted as Schrödinger equation:

$(C1):\;\;i\hbar\frac{\partial}{\partial
t}\Psi=\hat{H}\Psi,\;\;\;(C2):\;\;i\hbar\frac{d}{dt}\left|\psi\right>=H\left|\psi\right>,$
(7)

where (C1) is the form of a partial differential equation and (C2) is the form
which is convenient to use in numerical simulations. In (C2) $H$ is a
Hamiltonian describing system’s dynamics, and $|\psi\rangle$ stands for the
initial system’s state.

The von Neumann equation describes the quantum system’s evolution if an
influence of an external environment has to be considered:

$\dot{\rho}_{\rm tot}(t)=-\frac{i}{\hbar}[H_{\rm tot},\rho_{\rm
tot}(t)],\;\;\;H_{tot}=H_{sys}+H_{env}+H_{int},$ (8)

where $H_{sys}$ denotes the dynamics of a closed/core system, $H_{env}$ stands
for the environment’s dynamics, and $H_{int}$ describes the interaction’s
dynamics between the external environment and the system. The environment’s
influence can be removed from (8) by the partial trace operation. In such a
case, we obtain an equation which describes the dynamics of the core system.
Such a system can be expressed by the Lindblad master equation:

$\dot{\rho}(t)=-\frac{i}{\hbar}[H(t),\rho(t)]+\sum_{n}\frac{1}{2}\left[2C_{n}\rho(t)C_{n}^{+}-\rho(t)C_{n}^{+}C_{n}-C_{n}^{+}C_{n}\rho(t)\right],$
(9)

where $C_{n}$ stands for a set of collapse operators. These operators
represent the influence of an external environment which affects a simulated
system. Naturally, applying a collapse operator causes irreversible
modification of a quantum state. It should be emphasized that simulating the
behavior of a quantum system needs exponentially growing memory capacity
according to the system’s dimension.

The QTM is a method which facilitates reducing the memory requirements during
the simulation. Of course, any simulation of a quantum system’s behavior needs
calculating many single trajectories – the more, the better – because they
will be averaged to one final trajectory, and a greater number of trajectories
ensures improved accuracy. However, every trajectory may be simulated
separately, and this fact provides a natural background to utilize a parallel
approach while implementing the QTM.

If we would like to compare simulating the behavior of a quantum system with
the use of Lindblad master equation and the QTM, we should consider the
requirements of both methods on computational resources. The Lindblad master
equation methods utilize the density matrix formalism and the QTM is based on
a wave function of $n$-dimensional state’s vector (termed as a pure state). A
number of this vector’s entries grows exponentially, but using sparse matrices
facilitates efficient simulation of a quantum system’s behavior. It should be
stressed that a simulation based on the wave function concerns only one state
of a quantum system what seems to be a disadvantage of this solution because
for the Lindblad master equation methods, the density matrix describes many
different states of the same system. Unfortunately, using density matrices in
most cases is not possible because of memory requirements – the size of the
density matrix grows exponentially when the dimension of a simulated system
increases. While the QTM enables monitoring the influence of an external
environment on a quantum state by modifying the state’s vector with the use of
a collapse operator.

For the QTM, an evolution of a quantum system is described by so-called
effective Hamiltonian $H_{\rm eff}$, defined with the use of a set of collapse
operators $C_{n}$, and a system’s Hamiltonian $H_{\rm sys}$:

$H_{\rm eff}=H_{\rm sys}-\frac{i\hbar}{2}\sum_{n}C^{+}_{n}C_{n}.$ (10)

The set of operators $C_{n}$ determines the probability of a quantum jump
occurrence. This phenomenon is caused by a single collapse operator acting on
a current quantum state of the system. A probability of a quantum jump
occurrence is:

$\delta p=\delta t\sum_{n}\left<\psi(t)|C^{\dagger}_{n}C_{n}|\psi(t)\right>.$
(11)

If a quantum jump takes place, the system’s state at the moment of time
$(t+\delta t)$ – that is just after the collapse operation – can be expressed
as:

$\left|\psi(t+\delta
t)\right>=\frac{C_{n}\left|\psi(t)\right>}{\sqrt{\left<\psi(t)|C_{n}^{\dagger}C_{n}|\psi(t)\right>}}.$
(12)

Furthermore, if many collapse operators may be applied in the considered
model, the probability of using $i$-th operator is:

$P_{i}(t)=\frac{\left<\psi(t)|C_{i}^{\dagger}C_{i}|\psi(t)\right>}{\delta p}.$
(13)

Of course, simulating the phenomenon of a system’s collapse needs a random
numbers generator to ensure its probabilistic character (this issue is
discussed in the following subsection).

Let us emphasize that all the considered calculations correspond to operations
implemented on matrices and vectors within the QTM package. the needed
matrices are usually band matrices. Namely, it is convenient, in terms of the
memory consumption and the speed of calculation, to use sparse matrices. In
consequence, the compressed sparse row (CSR) format may be utilized – it also
gives an additional speed-up, especially when many matrix-vector
multiplications have to be realized.

By summarizing the above remarks, we can formulate an algorithm which presents
how a single quantum trajectory is calculated:

######  Algorithm 1.

Computation of a single trajectory
The algorithm computing a single trajectory, according to the QTM, may be
described as four computational steps:

  1. (A)

a value $0\leq r<1$ is computed by a pseudorandom number generator ($r$
denotes the probability of a quantum jump occurrence);

  2. (B)

to get the state’s vector at the moment $t$ the Schrödinger equation is
integrated with the Hamiltonian $H_{\rm eff}$ (providing that the state’s
vector norm has to be equal or greater to $r$:
$\langle\psi(t)|\psi(t)\rangle\geq r$);

  3. (C)

if a quantum jump is realized, then the system’s state projection at the
moment $t$, to one of the states given by Eq (12), is calculated. The operator
$C_{n}$ is selected to meet the following relation for the adequate $n$:
$\sum_{i=1}^{n}P_{i}(t)\geq r$ and $P_{i}(t)$ is given by Eq (13);

  4. (D)

the state obtained by the projection of wave-function in the previous step is
a new initial value corresponding to the moment of time $t$; next, the new
value of $r$ is randomly selected, and the procedure repeats the process of
the quantum trajectory generation, starting from the step (B) – more
precisely: the simulation is performed again, but starts from the previously
given value of $t$.

The presented approach is based on [11], [12], [27], [28], [6], [5].

###  2.3 Pseudorandom number generator

The QTM package described in this work utilizes a pseudorandom number
generator from the Generalized Feedback Shift Register (GFSR) class. More
precisely, we chose the LFSR113 method which is defined by using recurrence
over the field $F_{2}$ consisting of elements ${0,1}$:

$x_{n}=(a_{1}x_{n-1}+a_{2}x_{n-2}+\ldots+a_{k}x_{n-k})\mod 2,$ (14)

where $a_{j_{1}}$ are generator’s parameters ($j_{1}=1,...,k$) and $x_{j_{2}}$
are generator’s seeds ($j_{2}=n-1,...,n-k$). The generator’s period is
$r=2^{k}-1$ if and only if the characteristic polynomial of recurrence:

$P(z)=z^{k}-a_{1}z^{k-1}-\ldots-a_{k},$ (15)

is primitive.

The generated values, for $n\geq 0$, may be expressed as:

$u_{n}=\sum^{L}_{i=1}x_{ns+i-1}2^{-i},$ (16)

where $s$ denotes the step size, and $L$ stands for the number of bits in a
generated word. If $(x_{0},x_{1},\ldots,x_{k-1})\neq 0$, and $s$ is coprime to
$r$, then we obtain a periodic sequence of $u_{n}$ (with a period denoted as
$r$).

Of course, the quality of a generator is determined by the sequence $x_{n}$.
The proper choice of seeds is described in [14] where it is also shown that
four seeds ($k=4$) are sufficient to generate high-quality pseudorandom
numbers. The LFSR113 generator’s realization is very fast because of the
bitwise operations usage – this feature also does not collide with ensuring a
sufficient period for generated numbers: $2^{113}$. However, it should be
emphasized that the selection of generator’s seeds is crucial – there are four
initial values, and they have to be integers greater than 1, 7, 15 and 127,
respectively.

##  3 General implementation remarks

In the QTM, we calculate many independent trajectories. Because there is no
relations between trajectories, it is easy to implement them utilizing
parallel computing. Naturally, parallel computation shortens the time of
calculations in comparison to serial computation. The calculated trajectories
must be averaged to obtain the final trajectory.

Fig. 2 depicts the idea of computing trajectories in many computational nodes.
The exchange of necessary information between nodes is realized with use of
MPI protocol. This protocol offers a scalable solution, what means that our
QTM package works efficiently within a cluster of workstations, connected with
the use of the Ethernet network, and also with one multi-core personal
computer. In section 4, we show the acceleration of computations carried out
with the use of th QTM package – the acceleration is noticeable regardless of
whether we work with a cluster of workstations or with a one multiprocessor
computer.

Figure 2: Computing and averaging trajectories with $n+1$ computational nodes.
We assume that $n$ nodes calculate trajectories, and one node (denoted as node
$n$) plays a master role – it coordinates the process of computations for the
QTM. The figure presents also the general scheme of message passing during the
realization of the QTM. The master node sends tasks, e.g. commands to
calculate trajectories, and computing nodes response with results of their
calculations

The process of calculating the final trajectory is presented in Fig. 2. The
diagram shows the flow of messages during the calculations. The messages are
sent by the master node and other nodes respond to these messages. The main
part of the calculations is preceded by some preliminary activities, e.g.
initialization of pseudorandom numbers generators. Then, the body of
calculations begins, and computational nodes calculate and average
trajectories. The number of trajectories is influenced by the number of
computational nodes and by the number of trajectories planned to be generated
to simulate an analyzed problem. We assume an uniform load for every
computational node. Each computational node calculates a number of
trajectories and averages them. Finally, the averaged trajectories are sent to
the master node where are averaged once again in order to obtain the final
trajectory.

Calculating one trajectory is directly connected with the ZVODE package.
Despite the fact that this package was implemented in Fortran, the functions
from the ZVODE library may be called in a code written in C++, just like other
functions implemented in C++. To achieve that, we have to prepare an
intermediary function based on a template. The exemplary function called
zvode_method_for_mc is presented in Fig.3.

⬇

template <typename TYPE, size_t SIZE, size_t AlphaSize>

int zvode_method_for_mc(TYPE h, TYPE &_Tout_par, int steps,

simpleComplex<TYPE> &_T_par,

uVector< simpleComplex<TYPE>, SIZE > &_Y_par,

int (*fnc)(long int *NEQ, TYPE *T,

dblcmplx *Y,

dblcmplx *YDOT,

dblcmplx *RPAR, long int *IPAR) )

{

int i = 0;

int errLvl = 0;

zvode_((fncFP)fnc,

&neq, _Y_par.m, &_T_par.re, &_Tout_par, &itol,

&rtol, &atol_val, &itask, &istate, &iopt,

zwork, &lzw, rwork, &lrw, iwork, &liw,

(fncFP)dummyjex,

&mf, &rpar, &ipar);

if(istate < 0 ) {

if(verbose_mode == 1)

cerr << ”HALT:␣Error␣in␣zvode,␣value␣ISTATE␣=␣”

<< istate << endl ;

errLvl = -1;

}

return errLvl;

}

Figure 3: A function solving IVP (for a time-independent Hamiltonian) with the
use of the zvode_ method from the ZVODE package. The underscore placed in the
name of the method is a requirement imposed by the Fortran language

The scheme of tasks realized during the calculation of a single trajectory,
according to Algorithm 1, is shown in Fig. 4. The algorithm is implemented in
the C++ language, but it uses the zvode_method_for_mc method to solve a system
of ODEs. The ZVODE package does not offer the reentrant111The function’s
reentrant property allows many threads to utilize the same function. This
property is realized by avoidance of using shared and global variables in
function’s implementation. property, therefore, one MPI process may call only
one instance of any method from ZVODE. However, this feature does not pose a
problem because we may run many MPI processes at the same time.

Figure 4: The general scheme presenting generation of a single quantum
trajectory (PRN – pseudo-random numbers; H – Hamiltonian data;
zvode_method_for_mc – subroutine to ODEs solving; _T_par – time variable;
_Y_par – actual state of the system), fnc – calculates the right side of ODEs

It is a very important assumption for our implementation that we use templates
and static memory allocation (this applies to the code written by the user in
order to, for example, describe some operators) – these techniques enable
implementation of a code with an easier declaration of objects (the code is
more similar to a code of a script language than to a C++ code with dynamic
memory allocation). This makes using the package easier, but equally
efficient. We even expect better efficiency because the static size of data
structures is known during the compilation process, therefore, the compiler is
able to optimize the numerical computation.

The QTM package offers three basic data types. The first one is
`simpleComplex<T>`. It is dedicated to the operations on complex numbers,
where `T` may be a float or a double. If we want to use `simpleComplex<T>`
type together with the functions from the ZVODE library, then `T` is always
double. The second type `uVector<T>(int size)` is dedicated to vectors – the
first parameter describes the type of vector elements, and a constructor
parameter represents the number of entries. Similarly for matrices, a type
`uMatrix<T>(int size)` was introduced, where the first parameter describes the
type of matrix elements, and the constructor parameter stands for the
dimension (only square matrices are utilized in QTM). Package also utilize the
static and dynamic memory allocation. However, dynamic aspects of memory
management were hidden in the implementation layer. Users of the QTM package
is not obliged to create objects in dynamic approach.

An exemplary declaration of four matrices describing collapse operators is:

⬇

uMatrix< simpleComplex<double> > c_ops[ 4 ]={ {WV_LEAD_DIM},

{WV_LEAD_DIM}, {WV_LEAD_DIM}, {WV_LEAD_DIM}};

A value of WV_LEAD_DIM expresses, in this case, the dimension of operators and
matrices.

In the QTM package, a `uCSRMatrix<T>(size, rowptr, colind)` type was also
implemented. This type represents a column-row oriented sparse matrix. The
sparse matrices are often utilized to describe processes taking place in open
quantum systems. Using sparse matrices allows increasing the efficiency of
computation and decreasing the amount of memory needed to hold, for example,
collapse operators or expectation values. Furthermore, sparse matrices in CSR
format offer a shorter time of multiplying these matrices by vectors.

The QTM package utilizes the MPI standard. However, after defining the initial
structures describing a simulated problem, the user is not obliged to deal
with details of the MPI protocol. The whole computational process is realized
by a function `mpi_main`:

⬇

template<size_t N, size_t Ntrj, size_t _WV_LEAD_DIM,

size_t _WV_LEAD_DIM_SQR, size_t _C_OPS_SIZE>

int mpi_main(int argc, char *argv[], int verbose_mode,

double _from_time, double _to_time,

int use_colappse_operator,

int use_expecation_operator,

extra_options opt);

In the above form of `mpi_main` function, the Hamiltonian is time-independent.
If we would like the Hamiltonian to be time-dependent, in the call of
`mpi_main`, we point an additional function which is used during the
calculation of trajectories. The whole process of communication with the use
of MPI is automatically realized within the `mpi_main` function. The
parameters of the final trajectory may be directed to standard output or to a
text file, what is determined by a value of the parameter `opt`. This
parameter also serves to indicate the numerical method for solving ODEs. The
ZVODE package offers two methods: ADAMS and BDF. A method is selected as
below:

⬇

opt.ode_method = METBDF;

or

⬇

opt.ode_method = METADAMS;

The BDF method is dedicated to solving stiff ODEs but the ADAMS method may
also be used in the QTM (sometimes it needs to calculate more trajectories or
to increase the accuracy of the Adams’ ODEs solver).

##  4 Selected problems – implementation and performance

During the realization of the QTM implementation, the essential task was to
improve its performance in comparison to other existing QTM implementations.
The usage of C++ programming language and especially the MPI technology
facilitated us to obtain a stable solution with decent performance, thanks to
the parallel processing of trajectories.

The performance of the presented solution was compared with two recently
developed packages the QuTIP and the QuantumOptics.jl, which also fully
supports the QTM (and the QuTIP also utilizes the ZVODE method). We have
prepared two examples to examine the efficiency of the unitary Hamiltonian and
the trilinear Hamiltonian. We also show computations referring to the Jaynes-
Cummings model. The fourth example presents the results of the experiment also
conducted in [7], and shows the accuracy of the QTM package.

###  4.1 Unitary Hamiltonian

The first example concerns a simulation of a system described by the following
unitary Hamiltonian:

$H_{\mathrm{SYS}}=\frac{2\pi}{10}\sigma_{x},\;\;\;\mbox{and}\;\;\;\sigma_{x}=\left[\begin{array}[]{cc}0&1\\\
1&0\end{array}\right]$ (17)

where $\sigma_{x}$ represents Pauli operator $X$ (also termed as the $NOT$
operator). The initial state of the analyzed system is

$|\psi_{0}\rangle=|0\rangle=\left[\begin{array}[]{c}1\\\ 0\end{array}\right].$
(18)

The collapse operator and the expectation value operator, used in the
simulation, are given below:

$C_{0}=\frac{5}{100}\sigma_{x},\;E_{0}=\sigma_{z},\;\;\;\mathrm{and}\;\;\;\sigma_{z}=\left[\begin{array}[]{cc}1&0\\\
0&-1\end{array}\right]$ (19)

where $\sigma_{z}$ stands for the Pauli operator $Z$.

We ran the experiment on a PC equipped with the Intel i7-4950k 4.0 Ghz
processor under the Ubuntu 16.04 LTS operating system. The utilized version of
the QuTIP package is 4.2 and the QuantumOptics.jl ran on Julia 0.6.4. Although
the dimensions of the above structures are quite small, the simulation of 1000
trajectories with the QuTIP package, when only one core is running, needs
about $\approx 0.90$ seconds. Utilizing e.g. two cores for the calculation
does not reduce the time of process because the QuTIP needs time for
coordinating two threads. Of course, with the greater number of trajectories,
it is easy to observe that the time of simulation is shorter when more cores
are active. If the QTM package is used, calculating 1000 trajectories with the
use of e.g. of eight MPI computational nodes cores takes 2-3 seconds. This
time is mainly consumed by starting the MPI processes. It should be mentioned
that the QuantumOptics.jl package needs more time for calculation because of
the JIT compiler usage. However, for a such small problems the time of
calculation is almost the same, irrespective of the used package.

Fig. 5 shows the evolution of the expectation value with the Pauli Z gate in
time.

The most important pieces of the code are written in a direct form (we define
an effective Hamiltonian form directly filling the matrix $H$ with suitable
entries) because the QTM package allows specifying data structures directly.
Constants WD_LD and WD_LD_SQR respectively stand for the vector state
dimension, and its square and they are equal to two and four; co represents
the collapse operator, eo stands for the expectation operator value and H is
the Hamiltonian:

⬇

int r = 0;

co[0] = msc( 0.0, 0.0 ); co[1] = msc( 0.05, 0.0 );

co[2] = msc( 0.05, 0.0 ); co[3] = msc( 0.0, 0.0 );

eo[0] = msc( 1.0, 0.0); eo[1] = msc( 0.0, 0.0);

eo[2] = msc( 0.0, 0.0); eo[3] = msc(-1.0, 0.0);

alpha[0] = msc( 1.0, 0.0);

alpha[1] = msc( 0.0, 0.0);

H[0] = msc( -0.00125, 0.0);

H[1] = msc( 0.0, -0.62831853);

H[2] = msc( 0.0, -0.62831853);

H[3] = msc( -0.00125, 0.0);

c_ops[0].rows=2; c_ops[0].cols=2;

c_ops[0].m = co;

opt.type_output = OUTPUT_FILE;

opt.only_final_trj = 1;

opt.ode_method = METADAMS;

opt.tolerance = 1e-7;

opt.file_name = strdup(”output-data.txt”);

opt.fnc = &myfex_fnc_f1;

r = mpi_main<N, Ntrj, WV_LD, WV_LD_SQR, 1>(argc,

argv, 1, 0, 10, 1, 1, opt);

Figure 5: The figure showing changes of the unitary Hamiltonian’s expectation
value in the time

A significant part of the source code is a function calculating the right side
of the ODEs. In the case of the example for a unitary Hamiltonian, we may
utilize a direct approach that is multiplying the matrix form of the
Hamiltonian H by the vector Y. The product is assigned to the variable YDOT.
We can describe these calculation in details as follows:

⬇

int myfex_fnc_f1( long int *NEQ,

double *T,

dblcmplx *Y,

dblcmplx *YDOT,

dblcmplx *RPAR,

long int *IPAR)

{

simpleComplex<double> o0, o1, out0, out1;

o0.re=0.0; o0.im=0.0; o1.re=0.0; o1.im=0.0;

out0.re=0.0; out0.im=0; out1.re=0.0; out1.im=0.0;

o0 = Y[0] * H[0]; o1 = Y[1] * H[1];

out0 = o0 + o1;

o0.re=0.0; o0.im=0.0; o1.re=0.0; o1.im=0.0;

o0 = Y[0] * H[2]; o1 = Y[1] * H[3];

out1 = o0 + o1;

YDOT[0] = out0; YDOT[1] = out1;

return 0;

}

The approach presented above naturally allows further optimization of the code
during the operation of multiplying the matrix by the vector.

###  4.2 Trilinear Hamiltonian

In this example, we first assume that the Hamiltonian H is also time-
independent. The simulated process is a time evolution of an optical
parametric amplifier given by the following trilinear Hamiltonian [4]:

$H_{\mathrm{SYS}}=\mathbf{i}K(ab^{\dagger}c^{\dagger}-a^{\dagger}bc).$ (20)

The symbols $a$, $b$ and $c$ stand for the boson annihilation operators
corresponding to the pump, signal and idler fields respectively. The variable
$K$ represents the value of the coupling constant. For the purpose of the
tests, we assumed that $K=1$, and $\mathbf{i}$ represents an imaginary unit.
The initial state of analyzed system is:

$|\psi_{0}\rangle=|\alpha\rangle_{a}|0\rangle_{b}|0\rangle_{c},$ (21)

it is a coherent state for the pump mode ($a$), and the vacuum states for the
signal ($b$) and the idler ($c$) modes.

We also utilize three expectation operators defined as:

$\displaystyle a_{0}=a\otimes I_{b}\otimes I_{c},$ $\displaystyle
a_{1}=I_{a}\otimes b\otimes I_{c},$ $\displaystyle a_{2}=I_{a}\otimes
I_{b}\times c,$ (22) $\displaystyle
n_{0}=a^{\dagger}_{0}a_{0},\;n_{1}=a^{\dagger}_{1}a_{1},\;n_{2}=a^{\dagger}_{2}a_{2}.$

where $a$, $b$, $c$ still represent the boson annihilation operators with the
same or different dimensionality. $I_{a}$, $I_{b}$ and $I_{c}$ are the
identity operators for fields: pump_mode (a), vacuum (b) and idler (c). The
collapse operators are denoted as:

$\displaystyle c_{0}=\sqrt{2\cdot\gamma_{0}}a_{0}$ $\displaystyle
c_{1}=\sqrt{2\cdot\gamma_{1}}a_{1}$ (23) $\displaystyle
c_{2}=\sqrt{2\cdot\gamma_{2}}a_{2}$

where $\gamma_{0}=\gamma_{1}=0.1,\gamma_{2}=0.4$.

The obtained expectation values of the photons’ number, during the experiment
with 1000 trajectories, are presented in Fig. 6.

Figure 6: The expectation values describing the number of photons for an
experiment with generation of 1000 trajectories in three different expectation
operators (n0, n1, n2) defined in Eq (4.2) for a trilinear Hamiltonian

Fig. (7) depicts the simulations’ duration for the QuTIP, the QuantumOptics.jl
and the QTM packages with the different number of trajectories. We also
diversify the number of dimensions for the initial state in our experiments.
The calculation for the MPI protocol was conducted with the use of nine
computers equipped with Intel Core i7-4790K 4.0 GHZ processing units, and
working under the Ubuntu 16.04 LTS operating system. Each processor has four
cores, therefore, we were able to run the 32 MPI processes in eight
computational nodes. The last processing unit serves as a master node.

We can observe a direct profit thanks to dividing tasks between many
computational nodes. We should expect a significant speed-up of calculations
if the number of computational nodes increases, therefore, we consider this to
be a very good result. Especially for high-dimensional systems utilizing
parallel computing and many computational nodes, the execution of the tasks
directly translates into better performance i.e. shorter computation time with
the QTM package. It should be emphasized that in the case of calculating
10.000 – 20.000 trajectories 32 MPI processes still offer sufficient computing
power to increase the size of the problem or the number of trajectories. For
the QuTIP package, this computation was run with the use of one computational
node and the whole computing power was consumed. The same situation can be
observed in the case of the QuantumOptics.jl package. It also utilizes the JIT
system of compilation and offers a very efficient usage of a processing unit.
However, only calculating many trajectories using the MPI technology causes
significant reduction of the calculation time.

The simulation of a trilinear Hamiltonian requires using sparse matrices in
order to maintain both low memory requirements and high performance. The QTM
package offers basic data types to simplify the basic transformations and
realize the calculation (also the definitions of operators may be given
directly, as it was done for the Hamiltonian). A few selected steps connected
with the preparation of the Hamiltonian representations are presented below:

⬇

// Hamiltonian preparations

const int N0=8, N1=8, N2=8;

double K=1.0, gamma0=0.1, gamma1=0.1, gamma2=0.4;

double alpha_triham=sqrt(3);

uMatrix< simpleComplex<double> > d1(N0);

…

uMatrix< simpleComplex<double> > a0(N0*N1*N2), C0(N0*N1*N2), num0(N0*N1*N2);

…

uMatrix< simpleComplex<double> > H(N0*N1*N2);

…

destroy_operator(d1);

eye_of_matrix(d2);

eye_of_matrix(d3);

a0=tensor(d1,d2,d3);

…

H=unity*K*(a0*dagger(a1)*dagger(a2)-dagger(a0)*a1*a2);

Heff=prepare_effective_H(H, C0, C1, C2);

H = convertToCSRMatrix(Heff);

…

vacuum=tensor(basis(N0,0), basis(N1, 0), basis(N2,0));

D=exp_of_matrix(alpha_triham * dagger(a0)

- cojugate(alpha_triham)*a0);

alpha=D*vacuum;

…

r = mpi_main<100, Ntrj, WV_LD, WV_LD_SQR, 3>(argc,

argv, 1, 0, 4, 1, 1, opt);

The changes also apply to the function calculating the right side of the ODEs.
Luckily, only the function realizing multiplication has to be changed to the
one supporting the CSR matrices. Therefore, the function calculating the right
side of the ODEs is:

⬇

int myfex_fnc_f1(long int *NEQ,

double *T,

dblcmplx *Y,

dblcmplx *YDOT,

dblcmplx *RPAR,

long int *IPAR)

{

size_t i, j;

for ( i=0; i < WAVEVECTOR_LEAD_DIM; i++)

{

YDOT[i].re = 0.0;

YDOT[i].im = 0.0;

for ( j=H.row_ptr[i] ; j < H.row_ptr[i+1] ; j++)

{

YDOT[i]= YDOT[i] + (H.values[j] * Y[H.col_ind[j]]);

}

}

return 0;

}

Referring to the function calculating the right side of the ODEs, it should be
emphasized that the access to the current time value (parameter $T$) is
possible, what facilitates using a time-dependent Hamiltonian.

The QTM package is implemented in the C++ language and the presented examples
have to be compiled, therefore, the computation is carried out with a high
efficiency.

Figure 7: The duration of the simulation for the trilinear Hamiltonian. The
experiments were run for the different number of trajectories with use of the
QuTIP ( four CPU cores were used), the QuantumOptics.jl (termed as QuanOpt,
with only one CPU core used during the computations), and the QTM package
(with two different numbers of the computational nodes)

###  4.3 Jaynes-Cummings Model

The third example refers to the Jaynes-Cummings Model (JCM). This problem may
also be simulated in the QTM package. In this case, a single trajectory
corresponds to a computation run with the use of the master equation. Let us
assume that the system’s dimension is $N=40$. The dimensions of the operators
are given in subscripts.

$\displaystyle g=1,\;\Delta=-0.1,\;\alpha=4.0,$ $\displaystyle
a=d^{-}_{N}\otimes I_{2},\;b=I_{N}\otimes\sigma^{-},$ (25) $\displaystyle
H_{\mathrm{SYS}}=\Delta a^{\dagger}a+g(a^{\dagger}b+ab^{\dagger}),$

where $I_{2}$ stands for the identity operator sized $2\times 2$, $\sigma^{-}$
represents the annihilation operator for Pauli spins. The values $\Delta$ and
$g$ represent coupling strength between an atom and a cavity. The initial
state is expressed as:

$\displaystyle|\psi_{0}\rangle=|\chi^{\alpha}_{N}\rangle\otimes|1\rangle$ (26)

where we make tensor product between a coherent state $\chi^{\alpha}_{N}$ and
a single qubit in the state $|1\rangle$.

The results obtained during the simulation of the JCM are presented in Fig. 8.
It should be emphasized that we utilize only one trajectory, and needed
operators were represented by sparse matrices.

Figure 8: The values of the spin excitation in the JCM obtained after
calculating one trajectory without using collapse operators

The most important pieces of the code for the JCM are listed as follows:

⬇

const int N = 40;

double g = 1.0, delta = -0.1, alpha = 4.0

uMatrix< simleComplex<double>, N> d1(N);

uMatrix< simleComplex<double>, 2> eye2(2), sigmam(2);

…

sigmam_operator(sigmam);

destroy_operator(d1); eye_of_matrix(eye2);

…

a = tensor(d1, eye2);

sigmaminus = tensor(eyeN, sigmam);

expect = dagger(sigmamnus)*sigmaminus;

H=delta * dagger(a) * a + g*(dagger(a)*sigmaminus \+

a * dagger(sigmaminus));

Heff=prepare_effective_H(H);

H = convertToCSRMatrix(Heff);

…

r = mpi_main<600, 1, N, N*N, 0>(argc,

argv, 1, 0, 35, 0, 1, opt);

###  4.4 The birth and death of a photon in a cavity

The fourth and the last example refers to the convergency of the simulation of
photon’s birth and death in a cavity, and it is based on paper [7]. Let $N$ be
a number of task’s dimensions, and $N=5$. Then, the Hamiltonian and the
collapse operators are expressed as:

$\displaystyle a=d^{-},\;H_{\mathrm{SYS}}=a^{\dagger}\cdot
a,\;|\psi_{0}\rangle=[0,1,0,0,0]^{\dagger},\;\kappa=1.0/0.129,\;\mathrm{t}=0.063,$
(27) $\displaystyle c_{0}=\sqrt{\kappa\cdot(1+t)}\cdot
a,\;c_{1}=\sqrt{\kappa\cdot t}\cdot a^{\dagger}$

where $d^{-}$ denotes the destroy operator, H – the Hamiltonian and $c_{0}$,
$c_{1}$ represent the collapse operators, t – the temperature of the
environment.

Fig. 9 shows that the number of generated trajectories improves the accuracy
for a solved problem. Naturally, it also confirms the correctness of the
realized QTM implementation.

Figure 9: The probability of the decay for the one-photon state obtained with
the different numbers of trajectories. Increasing the number of trajectories
improves the accuracy of calculation

The source code for this example is very similar to the previously presented
pieces of the code. We utilize dense matrices because the system’s dimension
is low, and it does not influence the performance.

##  5 Conclusions

A package which implements the Quantum Trajectories Method approach was
presented in this article. The package is dedicated to examining the
properties of the quantum open systems. The implementation is based on the MPI
standard. The package was prepared in the C++ programming language, but the
implementation does not require from the final user any advanced programming
techniques. Utilizing the MPI standard allows using the package within systems
realizing high-performance computing, but also small systems like personal
computers because the communication processes introduced by the MPI package do
not virtually increase demand of computing powers.

The current version of our package is the first version that has been made
public. Naturally, further works and development are planned. We would like to
implement a version based on the CUDA/OpenCL technology, which will be able to
utilize computing powers of graphics processing units.

Recently, a new approach to the QTM was presented [15]. In next versions of
the package, this novelty will be considered: the version of used QTM will be
matched in accordance with a given problem in order to obtain a shorter time
of calculations.

## Acknowledgments

We would like to thank for useful discussions with colleagues at the Institute
of Control and Computation Engineering (ISSI) of the University of Zielona
Góra, Poland. We would like also to thank anonymous referees for useful
comments on the preliminary version of this paper. The numerical results were
done using the hardware and software available at the ”FPGA/GPU $\mu$-Lab”
located at the Institute of Control and Computation Engineering of the
University of Zielona Góra, Poland.

## References

  * [1] P. N. Brown, G. D. Byrne, and A. C. Hindmarsh.  Vode, a variable- coefficient ode solver.  SIAM J. Sci. Stat. Comput., 10:1038 – 1051, 1989.
  * [2] Lyndon Clarke, Ian Glendinning, and Rolf Hempel.  The mpi message passing interface standard.  In Karsten M. Decker and René M. Rehmann, editors, Programming Environments for Massively Parallel Distributed Systems, pages 213–218, Basel, 1994. Birkhäuser Basel.
  * [3] A.J. Daley.  Quantum trajectories and open many-body quantum systems.  Advances in Physics, 63(2):77–149, 2014.
  * [4] F. Dell’Anno, S. De Siena, and F. Illuminati.  Multiphoton quantum optics and quantum state engineering.  Physics Reports, 428(2):53 – 168, 2006.
  * [5] R. Dum, R. Zoller, and H. Ritsch.  Monte carlo simulation of the atomic master equation for spontaneous emission.  Phys. Rev. A, 45:4879, 1992.
  * [6] B.M. Garraway and P.L. Knight.  Evolution of quantum superpositions in open environments: Quantum trajectories, jumps, and localization in phase space.  Phys. Rev. A, 50:2548 – 2563, 1994.
  * [7] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deléglise, U. Busk Hoff, M. Brune, J.M. Raimond, and S. Haroche.  Quantum jumps of light recording the birth and death of a photon in a cavity.  Nature, 446:297 – 300, 2007.
  * [8] E. Hairer, S.P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations. I, 2nd ed. Springer, Berlin, Heidelberg, 1993.  Springer Series in Computational Mathematics.
  * [9] E. Hairer, S.P. Nørsett, and G. Wanner.  Solving Ordinary Differential Equations. II, 2nd ed. Springer, Berlin, Heidelberg, 1993.  Springer Series in Computational Mathematics.
  * [10] A.C. Hindmarsh.  LSODE and LSODI, two new initial value ordinary differential equation solvers.  ACM-SIGNUM Newsletter, 15(4):10 – 11, 1980.
  * [11] J.R. Johansson, P.D. Nation, and F. Nori.  Qutip: An open-source python framework for the dynamics of open quantum systems.  Comp. Phys. Comm., 183:1760 – 1772, 2012.
  * [12] J.R. Johansson, P.D. Nation, and F. Nori.  Qutip 2: A python framework for the dynamics of open quantum systems.  Comp. Phys. Comm., 184:1234 – 1240, 2013.
  * [13] S. Krämer, D. Plankensteiner, L. Ostermann, and H. Ritsch.  Quantumoptics.jl: A julia framework for simulating open quantum systems.  Computer Physics Communications, 227:109 – 116, 2018.
  * [14] P. L’Ecuyer.  Tables of maximally-equidistributed combined lfsr generators.  Mathematics of Computation, 68(255):261 – 269, 1999.
  * [15] A. Liniov, V. Volokitin, M. Meyerov, I.and Ivanchenko, and S. Denisov.  Increasing performance of the quantum trajectory method by grouping trajectories.  In V. Voevodin and S. Sobolev, editors, Supercomputing, pages 136–150, Cham, 2017. Springer International Publishing.
  * [16] M.A. Nielsen and I.L. Chuang.  Quantum Computation and Quantum Information: 10th Anniversary Edition.  Cambridge University Press, 2010.
  * [17] V.G. Pimenov and A.S. Hendy.  A numerical solution for a class of time fractional diffusion equations with delay.  International Journal of Applied Mathematics and Computer Science, 27(3):477 – 488, 2017.
  * [18] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery.  Numerical recipes (3rd ed.).  Cambridge University Press, Cambridge/New York, 2007.
  * [19] K. Radhakrishnan and A.C. Hindmarsh.  Description and use of LSODE, the livermore solver for ordinary differential equations.  Technical Report UCRL-ID-113855, Lawrence Livermore National Laboratory, Livermore, CA, 1993.
  * [20] Raimar Sandner and András Vukics.  C++qedv2 milestone 10: A c++/python application-programming framework for simulating open quantum dynamics.  Computer Physics Communications, 185(9):2380 – 2382, 2014.
  * [21] Marek Sawerwain and Joanna Wiśniewska.  Github repository for qtm package, October 2018.  https://github.com/qMSUZ/QTM.
  * [22] Marek Sawerwain and Joanna Wiśniewska.  The programming library implementing Quantum Trajectories Method (QTM package), October 2018.  doi: 10.5281/zenodo.1443581.
  * [23] R. Schacka and T.A. Brun.  A c++ library using quantum trajectories to solve quantum master equations.  Comp. Phys. Comm., 102:210 – 228, 1997.
  * [24] K. Soetaert, J. Cash, and F. Mazzia.  Solving Differential Equations in R.  Springer-Verlag, Berlin Heidelberg, 2012.
  * [25] Sze M Tan.  A computational toolbox for quantum and atomic optics.  Journal of Optics B: Quantum and Semiclassical Optics, 1(4):424, 1999.
  * [26] L. Vasiliu, F. Pop, F. Negru, M. Mocanu, V. Cristea, and J. Kolodziej.  A hybrid scheduler for many task computing in big data systems.  International Journal of Applied Mathematics and Computer Science, 27(2):385 – 399, 2017.
  * [27] J. Wiśniewska and M. Sawerwain.  Gpu: Accelerated computation routines for quantum trajectories method.  In V. Kindratenko, editor, Numerical Computations with GPUs, pages 299 – 318. Springer, 2014.
  * [28] J. Wiśniewska, M. Sawerwain, and W. Leoński.  High performance computing and quantum trajectory method in cpu and gpu systems.  Journal of Physics: Conference Series, 574(1):012127, 2015.
  * [29] R. E. Wyatt.  Quantum Dynamics with Trajectories.  Springer, New York, 2005.
  * [30] K.W. Yip, T. Albash, and D.A. Lidar.  Quantum trajectories for time-dependent adiabatic master equations.  Phys. Rev. A, 97:022116, 2018.