# Entropic and trace distance based measures of non-Markovianity

Federico Settimo  Dipartimento di Fisica “Aldo Pontremoli”, Università degli
Studi di Milano, via Celoria 16, 20133 Milan, Italy  Heinz-Peter Breuer
<EMAIL_ADDRESS>Physikalisches Institut, Universität Freiburg,
Hermann-Herder-Straße 3, D-79104 Freiburg, Germany  EUCOR Centre for Quantum
Science and Quantum Computing, Universität Freiburg, Hermann-Herder-Straße 3,
D-79104 Freiburg, Germany  Bassano Vacchini<EMAIL_ADDRESS>Dipartimento di Fisica “Aldo Pontremoli”, Università degli Studi di Milano,
via Celoria 16, 20133 Milan, Italy  Istituto Nazionale di Fisica Nucleare,
Sezione di Milano, via Celoria 16, 20133 Milan, Italy

###### Abstract

We analyze and compare different measures for the degree of non-Markovianity
in the dynamics of open quantum systems. These measures are based on the
distinguishability of quantum states which is quantified, on the one hand, by
the trace distance or, more generally, by the trace norm of the Helstrom
matrix, and, on the other hand, by entropic quantifiers: the Jensen-Shannon
divergence, the Holevo or the quantum skew divergence. We explicitly construct
a qubit dynamics for which the trace norm based non-Markovianity measure is
nonzero, while all the entropic measures turn out to be zero. This leads to
the surprising conclusion that the non-Markovianity measure which employs the
trace norm of the Helstrom matrix is strictly stronger than all entropic non-
Markovianity measures.

##  I Introduction

The study of quantum non-Markovian dynamics involves the investigation of the
very notion of stochastic process in the quantum realm, as well as the
characterization of memory effects in open quantum system dynamics [1, 2, 3,
4]. Memory effects in the dynamics of a quantum system interacting with an
external environment can be uniquely traced back to local retrieval of
exchanged information in the approach to non-Markovianity based on the non
monotonic behavior in time of distinguishability of quantum states. This
strategy was introduced in [5] and validated for different distinguishability
quantifiers of quantum states. In particular, while the original approach was
focused on the trace distance, it was later put into evidence that invariance
under translations of this quantifier led to failure in assessing memory
features in certain dynamics [6]. To avoid this difficulty, the trace norm of
the Helstrom matrix was used as a generalized trace distance also sensitive to
translations [7, 8]. A crucial feature associated to the trace norm of the
Helstrom matrix is the fact that its non-monotonicity in time is equivalent to
lack of P divisibility of the considered dynamics, provided the evolution is
invertible as a linear transformation. In such a way a direct relation could
be established between a divisibility and a distinguishability criterion.

More recently, entropic distinguishability quantifiers have also been
introduced and directly connected to the notion of non-Markovianity as due to
information backflow [9]. To this aim, suitable regularizations of the quantum
relative entropy have been considered, which, at variance with the quantum
relative entropy, remain finite for any pair of states, and allow to introduce
triangle-like inequalities which connect revivals of the quantifier to
information backflow, even in the absence of a true triangle inequality as for
distances. Furthermore, these entropic quantifiers are also sensitive to
translations. In this framework, a special role is played by the Jensen-
Shannon divergence, whose square root is a true distance [10, 11, 12].

Given that entropic distinguishability quantifiers are contractions under
positive trace preserving maps which are not necessarily completely positive,
as it happens for the trace distance and the trace norm of the Helstrom
matrix, a natural question is the role of P divisibility in this context.
Importantly, we show by means of an example that the Jensen-Shannon
divergence, as well as the other entropic quantifiers, might fail in detecting
breaking of P divisibility.

Our results imply that the non-Markovianity measure employing the trace norm
of the Helstrom matrix is strictly stronger than all of the entropic non-
Markovianity measures, leading to a nonzero value even for dynamics for which
the entropic measures are zero, while the opposite cannot happen.

The paper is organized as follows. In Sec. II we introduce and exemplify the
general framework for the treatment of non-Markovianity based on
distinguishability quantifiers, together with the associated measures. In Sec.
III we outline the connection between non-Markovianity and divisibility of the
dynamics, and explore this relationship in its dependence on the considered
distinguishability quantifier. In particular, we construct an example of non-P
divisible evolution whose non-Markovianity measure is zero according to
entropic quantifiers. We summarize and discuss the conclusions of our work in
Sec. IV.

##  II Entropic and trace distance based distinguishability quantifiers

Let us begin by introducing the general framework of non-Markovianity for the
dynamics of open quantum systems. The main aim is to compare the well-known
measure of memory effects based on the trace distance with other measures
using alternative different distinguishability quantifiers between quantum
states, in particular those related to the quantum relative entropy.

###  II.1 Trace distance and Helstrom matrix

In the framework of quantum information and statistics there are many
different quantifiers of distinguishability between two quantum states $\rho$
and $\sigma$. A very important one is given by the trace distance (TD) [13]

$D(\rho,\sigma)=\frac{1}{2}\lVert\rho-\sigma\rVert,$ (1)

where the trace norm of any trace-class operator $A$ is defined as $\lVert
A\rVert=\operatorname{tr}\sqrt{A^{\dagger}A}$. The TD is bounded, $0\leq
D(\rho,\sigma)\leq 1$, with $D(\rho,\sigma)=0$ if and only if $\rho=\sigma$,
and $D(\rho,\sigma)=1$ if and only if $\rho\perp\sigma$. Additionally, the TD
obeys the triangle inequality

$D(\rho,\sigma)\leqslant D(\rho,\tau)+D(\tau,\sigma),$ (2)

is contractive under the action of any completely positive trace preserving
(CPTP) map $\Lambda$, as well as of any positive trace preserving map [14]

$D(\Lambda\rho,\Lambda\sigma)\leq D(\rho,\sigma),$ (3)

and it is invariant under unitary and anti-unitary transformations [15]. It is
also invariant under translations, in the sense that

$D(\rho+A,\sigma+A)=D(\rho,\sigma)$ (4)

for any operator $A$. This follows directly from the fact that the TD depends
on the difference between its two arguments.

It is possible to give the TD an interpretation as the bias in favour of a
correct identification between two quantum states, upon performing a single
measurement. Let us suppose that Alice prepares the state $\rho$ or $\sigma$,
each with probability $\frac{1}{2}$, and sends it to Bob; the TD is linked to
Bob’s maximal probability of correctly distinguishing between the two as [16]

$P_{\text{dist}}(\rho,\sigma)=\frac{1}{2}(1+D(\rho,\sigma)).$ (5)

This feature, combined with the contractivity of the TD under CPTP maps (3),
tells us that CPTP maps cannot increase the probability of distinguishing
between quantum states.

The idea of using the trace norm $\lVert\cdot\rVert$ to quantify the bias in
favour of a correct identification can be generalised also to the case in
which the two states $\rho$ and $\sigma$ are not prepared with the same
a-priori probability. In fact, if one supposes that Alice prepares $\rho$ with
probability $p$ and $\sigma$ with probability $1-p$, then Bob’s maximal
probability of distinguishing between the two is given by [17]

$P_{\text{dist}}(\rho,\sigma)=\frac{1}{2}(1+\lVert\Delta\rVert),$ (6)

where

$\Delta=p\rho-(1-p)\sigma$ (7)

is known as the Helstrom matrix [18]. The trace norm of $\Delta$ represents
the bias in favour of a correct identification and Eq. (6) reduces to (5) in
the unbiased case $p=\tfrac{1}{2}$. The Helstrom matrix can be seen as a
generalisation of the TD to generic ensembles $\set{(p,\rho),(1-p,\sigma)}$
and it inherits properties such as boundedness and contractivity from the TD.

###  II.2 Jensen-Shannon and skew divergences

The TD is not the only possible quantifier of distinguishability between
quantum states. A particularly interesting distinguishability quantifier is
the relative entropy

$S(\rho,\sigma)=\begin{cases}\operatorname{tr}[\rho\log\rho-\rho\log\sigma]\quad&\text{if
}\operatorname{supp}\rho\subseteq\operatorname{supp}\sigma\\\
\infty\quad&\text{otherwise}\end{cases},$ (8)

where we take the logarithm in base 2. The relative entropy, just like the TD,
is contractive under both CPTP maps and positive trace preserving maps [19].
However, as it is evident from the definition, it is not bounded.

The relative entropy can also be naturally associated to a distinguishability
task. In particular, let us suppose to be able to prepare and measure the
states an arbitrarily large number $N$ of times. The relative entropy
$S(\rho,\sigma)$ represents the maximal asymptotic rate at which the
probability of erroneously concluding that the state is $\rho$, when it is
actually $\sigma$, decays with the size $N$ of the sample over which a
measurement is performed, so that, for large enough $N$ the probability of
correctly identifying the state is [20, 21, 22]

$P_{N,\text{dist}}(\rho,\sigma)=1-e^{-NS(\rho,\sigma)}.$ (9)

Unboundedness here arises naturally: whenever
$\operatorname{supp}\rho\not\subseteq\operatorname{supp}\sigma$, it is
possible with certainty to distinguish $\rho$ from $\sigma$ with only a finite
number of measurements and hence the rate is infinite.

It is possible to define a smoothed version of the relative entropy, namely
the Jensen-Shannon divergence (JSD), according to [23]

$\begin{split}J(\rho,\sigma)=&\frac{1}{2}S\left(\rho,\frac{\rho+\sigma}{2}\right)+\frac{1}{2}S\left(\sigma,\frac{\rho+\sigma}{2}\right)\\\
=&H\left(\frac{\rho+\sigma}{2}\right)-\frac{1}{2}H(\rho)-\frac{1}{2}H(\sigma),\end{split}$
(10)

where $H$ denotes the von Neumann entropy
$H(\rho)=-\operatorname{tr}\rho\log\rho$. This definition ensures that the JSD
inherits the contractivity under CPTP maps from the relative entropy and,
additionally, it is bounded according to $0\leq J(\rho,\sigma)\leq 1$, with
$J(\rho,\sigma)=0$ if and only if $\rho=\sigma$, while $J(\rho,\sigma)=1$ if
and only if $\rho\perp\sigma$. In particular, it can be bounded by monotonic
functions of the TD as [24, 25]

$\frac{1}{2}D(\rho,\sigma)^{2}\leq J(\rho,\sigma)\leq D(\rho,\sigma).$ (11)

The lower bound directly follows from the Pinsker inequality [21]. Fig. 1
shows these bounds together with the value of the TD and the JSD for randomly
chosen pairs of states. The JSD is not invariant under translations in the
sense of (4), since, unlike the TD, it does not depend solely on the
difference $\rho-\sigma$. This fact is visualized in Fig. 2.

Figure 1: Plot of the JSD and the TD for $10^{5}$ randomly generated pairs of
qubits. The red lines are the upper and lower bounds of Eq. (11), which are
valid also for arbitrarily dimensional Hilbert spaces. Figure 2: Plots of the
TD (left) and of the JSD (right) for qubit states represented by Bloch vectors
of the form $\boldsymbol{r}_{\rho}=(x_{1},0,0)^{\top}$,
$\boldsymbol{r}_{\sigma}=(x_{2},0,0)^{\top}$. The translational invariance of
the TD is reflected by the fact that the plot on the left only depends on the
difference $x_{1}-x_{2}$. Note further the different sensitivity in the
central and corner regions.

The JSD, unlike the TD, is not a distance since it does not obey the triangle
inequality. However, it has been proven that its square root
($\sqrt{\text{JSD}}$) does obey this inequality and is indeed a distance [10,
11, 12]. Even if it does not obeys the triangle inequality, the JSD obeys a
triangle-like inequality

$\begin{split}J(\rho,\sigma)-J(\rho,\tau)&\leq\frac{1+D(\sigma,\tau)}{2}-J(1,D(\sigma,\tau))\\\
&\leq\sqrt[4]{2}\sqrt[4]{J(\sigma,\tau)},\end{split}$ (12)

which follows from the triangle-like inequalities presented in [24].

It is possible to generalise the JSD to generic ensembles
$\set{(\mu,\rho),(1-\mu,\sigma)}$, however, unlike for the TD, such
generalisation is not unique. As suggested in [26] we point to two distinct
generalizations based on a skewed version of the relative entropy, also called
telescopic relative entropy in the quantum setting [27, 28]. We therefore
introduce the Holevo skew divergence

$K_{\mu}(\rho,\sigma)=\frac{\chi_{\mu}(\rho,\sigma)}{h(\mu)},$ (13)

where

$h(p)=-p\log p-(1-p)\log(1-p)$ (14)

is the binary entropy for the distribution $\set{p,1-p}$ and

$\chi_{\mu}(\rho,\sigma)=H(\mu\rho+(1-\mu)\sigma)-\mu
H(\rho)-(1-\mu)H(\sigma)$ (15)

is the Holevo $\chi$ quantity [29] for the considered ensemble, as well as the
quantum skew divergence

$\begin{split}S_{\mu}(\rho,\sigma)=&\frac{\mu}{\log(1/\mu)}S(\rho,\mu\rho+(1-\mu)\sigma)\\\
&+\frac{1-\mu}{\log(1/(1-\mu))}S(\sigma,(1-\mu)\sigma+\mu\rho).\end{split}$
(16)

Both quantities are bounded and they reduce to the JSD in the unbiased case
$\mu=\frac{1}{2}$. Furthermore, they both obey triangle-like inequalities
similar to the ones that hold for the unbiased case (12), namely [9, 26]

$\displaystyle
S_{\mu}(\rho,\sigma)-S_{\mu}(\rho,\tau)\leq\eta_{\mu}^{S}\sqrt[4]{S_{\mu}(\sigma,\tau)}$
(17) $\displaystyle
K_{\mu}(\rho,\sigma)-K_{\mu}(\rho,\tau)\leq\eta_{\mu}^{K}\sqrt[4]{K_{\mu}(\sigma,\tau)},$
(18)

with

$\displaystyle\eta_{\mu}^{S}=\log\left(\frac{1}{\mu(1-\mu)}\right)\sqrt[4]{\frac{\mu(1-\mu)}{2\,h(\mu)\log^{3}(\mu)\log^{3}(1-\mu)}},$
$\displaystyle\eta_{\mu}^{K}=\sqrt[4]{\frac{8\mu(1-\mu)}{h(\mu)^{3}}}.$ (19)

###  II.3 Non-Markovianity measures from distinguishability quantifiers

The unavoidable interaction between a quantum system and its surroundings
leads to system-environment correlations and non-unitary time evolution of the
state. Assuming that at the initial time $t=0$ the global system-environment
state is factorised, the dynamics is described by a one-parameter family of
CPTP dynamical maps $\Phi=\set{\Phi_{t}}{0\leq t\leq T,\Phi_{0}=\mathbbm{1}}$
such that $\rho(t)=\Phi_{t}\rho(0)$. Assuming that $\Phi_{t}^{-1}$ exists at
all times $t\geq 0$, it is possible to define a two-parameter family of maps
as

$\Phi_{t,s}=\Phi_{t}\Phi_{s}^{-1},\quad t\geq s\geq 0,$ (20)

such that $\Phi_{t,0}=\Phi_{t}$, describing the evolution of the state from
time $s$ to time $t$. The dynamics is said to be (C)P divisible if
$\Phi_{t,s}$ is (completely) positive for all times $t\geq s\geq 0$.

The interaction between the system and the environment can lead to memory
effects during the dynamics of the state. If this happens, the dynamics is
said to be non-Markovian. Following [5, 1] we define a family of non-
Markovianity measures, based on some distinguishability quantifier $d$, as

$\mathscr{N}^{d}(\Phi)=\sup\int\limits_{\sigma_{d}(t)\char 62\relax
0}dt\,\sigma_{d}(t),$ (21)

where

$\sigma_{d}(t)=\frac{d}{dt}\,d(\rho_{1}(t),\rho_{2}(t)),$ (22)

and the maximisation is performed over all possible pairs of initial states
$\rho_{1,2}(0)$ and any eventual parameter defining the distinguishability
quantifier $d$, such as the skewing parameter $\mu$ defining $S_{\mu}$ or
$K_{\mu}$. A certain pair of initial states $\rho_{1,2}(0)$ is said to be
optimal, if the maximum of equation (21) is attained on this pair. Thus, a
dynamical map $\Phi$ is Markovian according to the quantifier $d$ if and only
if $\mathscr{N}^{d}(\Phi)=0$ or, equivalently, if $d(\rho_{1}(t),\rho_{2}(t))$
is a monotonic function of time for any initial pair of states
$\rho_{1,2}(0)$. Alternative approaches are indeed possible, such as
violations of divisibility of the dynamical map [1, 2, 3, 30, 31, 32].

Following [26], in order to have a well-defined measure of non-Markovianity we
ask the quantifier $d$ to obey three properties:

  1. 1.

Boundedness and indistinguishability of identical states:

$0\leq d(\rho,\sigma)\leq 1,$ (23)

with $d(\rho,\sigma)=1$ if and only if $\rho\perp\sigma$, and
$d(\rho,\sigma)=0$ if and only if $\rho=\sigma$. Considering bounded
distinguishability quantifiers allows to perform the maximation in Eq. (21)
thus warranting that the measure of non-Markovianity is well-defined.

  2. 2.

Contractivity under CPTP maps:

$d(\Lambda\rho,\Lambda\sigma)\leq d(\rho,\sigma)$ (24)

for any CPTP map $\Lambda$. This property is crucial so that any revival in
$d$ must necessarily correspond to violations of divisibility of the dynamical
map. In fact, if $\Phi$ is CP divisible, $d$ must be monotonically decreasing,
since the map $\Phi_{t,s}$ describing the evolution from $s$ to $t\char
62\relax s$ is always CPTP. Therefore, a revival in $d$ is possible only if
$\Phi$ violates the divisibility.

  3. 3.

Triangle-like inequalities:

$\displaystyle d(\rho,\sigma)-d(\rho,\tau)\leq\phi(d(\sigma,\tau)),$ (25)
$\displaystyle d(\rho,\sigma)-d(\tau,\sigma)\leq\phi(d(\rho,\tau)),$ (26)

where $\phi(x)$ is a strictly positive concave function for $x\char 62\relax
0$, and with $\phi(0)=0$. This property allows for a microscopic
interpretation of the revivals of $d$ as a twofold exchange of information,
which is at first stored in external degrees of freedom and later retrieved in
the open system.

TD, JSD, $\sqrt{\text{JSD}}$, and their generalisations all obey properties
1-3, and hence lead to a well-defined measure of non-Markovianity. For the TD
and the $\sqrt{\text{JSD}}$, which are actually distances, the function $\phi$
is given by the identity, while, for the JSD and the other entropic
quantities, the function $\phi$ is proportional to the fourth root as follows
from Eqs. (12), as well as (17) and (18).

Given two distinguishability quantifiers $d_{1}$ and $d_{2}$ satisfying 1-3,
we say that $\mathscr{N}^{d_{1}}$ is stronger than $\mathscr{N}^{d_{2}}$ if,
for any dynamical map $\Phi$ such that $\mathscr{N}^{d_{2}}(\Phi)\char
62\relax 0$, then $\mathscr{N}^{d_{1}}(\Phi)\char 62\relax 0$. Furthermore,
$\mathscr{N}^{d_{1}}$ is strictly stronger than $\mathscr{N}^{d_{2}}$ if it is
stronger and there exists $\Phi$ such that $\mathscr{N}^{d_{2}}(\Phi)=0$ and
$\mathscr{N}^{d_{1}}(\Phi)\char 62\relax 0$. Viceversa, $\mathscr{N}^{d_{1}}$
is (strictly) weaker than $\mathscr{N}^{d_{2}}$ if $\mathscr{N}^{d_{2}}$ is
(strictly) stronger than $\mathscr{N}^{d_{1}}$. Two measures are said to be
equivalent if $\mathscr{N}^{d_{1}}$ is both stronger and weaker than
$\mathscr{N}^{d_{2}}$.

An important distinguishability quantifier obeying the abovementioned three
properties is the TD. Optimal pairs for this measure must always be orthogonal
and therefore on the border of the set of states [33]. Additionally, the
triangle inequality (2) allows to upper bound the revival of the TD from $s$
to a later time $t\char 62\relax s\char 62\relax 0$ as [34, 35, 36]

$\begin{split}\Delta
D(t,s)=&D(\rho_{S}^{1}(t),\rho^{2}_{S}(t))-D(\rho_{S}^{1}(s),\rho^{2}_{S}(s))\\\
\leq&D(\rho^{1}_{SE}(s),\rho^{1}_{S}(s)\otimes\rho^{1}_{E}(s))\\\
&+D(\rho^{2}_{SE}(s),\rho^{2}_{S}(s)\otimes\rho^{2}_{E}(s))\\\
&+D(\rho^{1}_{E}(s),\rho^{2}_{E}(s)),\end{split}$ (27)

where $\rho^{i}_{SE}(s)$, for $i=1,2$, is the global system-environment state,
and $\rho^{i}_{S}(s)=\operatorname{tr}_{E}\rho^{i}_{SE}(s)$ and
$\rho^{i}_{E}(s)=\operatorname{tr}_{S}\rho^{i}_{SE}(s)$ are, respectively, the
reduced system and environmental states at time $s$. This allows for a
microscopic interpretation of the measure of non-Markovianity [34, 35]: a
revival in the TD is possible only if at time $s$ the two environments are
different or if correlations have built up during the dynamics. Therefore,
information is stored as correlations or as difference between the
environmental states and can later flow back into the open system. A similar
interpretation also holds for the measure of non-Markovianity arising from the
Helstrom matrix [8], as well as for entropic distinguishability quantifiers
[9, 26]. In particular, the non-Markovianity measure obtained according to Eq.
(21) when the quantifier $d$ is the trace norm of the Helstrom matrix Eq. (7),
which we denote as $\mathscr{N}^{\Delta}(\Phi)$, is positive if and only if
$\Phi$ is not P divisible as has been shown in [7, 8] building on results in
[37, 38]. The measure based on the TD, instead, is strictly weaker than
$\mathscr{N}^{\Delta}$, since it can equal zero even for non-P divisible
dynamics, due to its translational invariance [6]. Given that both the TD and
the quantum relative entropy are contractive under positive trace preserving
maps, and the equivalence between a non-Markovianity measure and a
divisibility property was obtained considering positivity, from now on we will
concentrate our attention simply on positivity.

Let us now focus our attention to the entropic distinguishability quantifiers
of Sec. II.1. Except for the relative entropy, which is unbounded, all the
other quantifiers obey properties 1-3 and hence can be used to define a
measure of non-Markovianity. In particular, we want to investigate whether
these measures of non-Markovianity are equivalent to $\mathscr{N}^{\Delta}$.
Namely, we want to know whether the equivalence between positivity of
$\mathscr{N}^{d}(\Phi)$ and lack of P divisibility of $\Phi$ also holds when
choosing $d$ as one of the previously introduced entropic quantifiers. We will
show in Sec. III.3 that this is not the case: we will use a counterexample to
point out that $\mathscr{N}^{\Delta}$ is strictly stronger. Let us focus in
particular on the JSD, since the $\sqrt{\text{JSD}}$ is just a monotonic
function of it and hence $\mathscr{N}^{J}$ and $\mathscr{N}^{\sqrt{J}}$ are
equivalent.

###  II.4 Behavior on unital models

Let us now focus our discussion on qubits, since it suffices considering the
simplest non-trivial case to prove that $\mathscr{N}^{\Delta}$ is strictly
stronger than $\mathscr{N}^{J}$. For qubits, a generic state $\rho$ can be
represented by means of a real three-dimensional Bloch vector
$\boldsymbol{r}_{\rho}$ with $|\boldsymbol{r}_{\rho}|\leq 1$ in the form

$\rho=\frac{1}{2}(\mathbbm{1}+\boldsymbol{r}_{\rho}\cdot\boldsymbol{\sigma}),$
(28)

where $\boldsymbol{\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z})^{\top}$ is the
vector of the Pauli matrices. Under the action of a generic trace and
Hermiticity preserving linear map, the Bloch vector associated to the state
transforms according to

$\boldsymbol{r}\mapsto\boldsymbol{r}_{t}=D(t)\boldsymbol{r}+\boldsymbol{\kappa}(t),$
(29)

where
$D(t)=\operatorname{diag}\\{\lambda_{1}(t),\lambda_{2}(t),\lambda_{3}(t)\\}$
is a real diagonal $3\times 3$ matrix and $\boldsymbol{\kappa}(t)$ is a real
three-dimensional vector [39]. The representation Eq. (29) is valid up to
orthogonal transformations which in the present context plays no role, given
that both TD and JSD are invariant under unitary transformations, and the
measure of non-Markovianity in Eq. (21) is obtained by maximizing over the
possible initial states. Let us focus in particular on unital maps, which are
the maps that preserve the maximally mixed state at any time $t\geq 0$:
$\Phi_{t}\left[\frac{\mathbbm{1}}{2}\right]=\frac{\mathbbm{1}}{2}$.
Alternatively, employing the representation (29), they are the maps such that
$\boldsymbol{\kappa}(t)=0$ at all times $t\geq 0$. A dynamics of this kind is
not P divisible if and only if at least one of the functions $\lambda_{i}(t)$
does not decrease monotonically.

An important feature of unital dynamics is that $\mathscr{N}^{D}(\Phi)\char
62\relax 0$ if and only if $\Phi$ violates P divisibility: any backflow of
information, corresponding to a violation of P divisibility, is witnessed by
the TD, without the need to generalise it to the Helstrom matrix.
Interestingly, this feature also holds for the JSD. Let $\rho_{1,2}(0)$ be the
optimal pair for the TD. Since they must be pure and orthogonal, we have
$\rho_{1}(0)+\rho_{2}(0)=\mathbbm{1}$ and the Bloch vectors representing the
states obey $\boldsymbol{r}_{1}(0)=-\boldsymbol{r}_{2}(0)$. Thanks to
unitality the transformed average state $(\rho_{1}(t)+\rho_{2}(t))/2$ remains
the maximally mixed state, so that
$\boldsymbol{r}_{1}(t)=-\boldsymbol{r}_{2}(t)$ holds at all times. Thus the TD
between the two states reads $D(\rho_{1}(t),\rho_{2}(t))=r(t)$ and both
evolved states have the same von Neumann entropy

$H(\rho_{1}(t))=H(\rho_{2}(t))=h\left(\frac{1-r(t)}{2}\right),$ (30)

where $h$ is the binary entropy introduced in Eq. (14). It is therefore
possible to rewrite the JSD using Eq. (10) as

$J(\rho_{1}(t),\rho_{2}(t))=1-h\left(\frac{1-D(\rho_{1}(t),\rho_{2}(t))}{2}\right).$
(31)

This expression is a monotonic function of the TD and thus a revival in the
JSD is witnessed if and only if it is witnessed by the TD. Therefore, as it
happens for the TD, $\mathscr{N}^{J}(\Phi)\char 62\relax 0$ if and only if
$\Phi$ violates P divisibility. Unlike the TD, the characterisation of optimal
pairs for the JSD is still an open problem. For unital maps acting on qubits,
numerical evidence suggests that they must be pure and orthogonal, just like
for the TD. This feature allows for an interesting interpretation of
$\mathscr{N}^{J}$ in terms of the von Neumann entropy. By employing the first
line of Eq. (31), holding for any pair of pure and orthogonal initial states,
it is therefore possible to rewrite the measure of non-Markovianity as

$\mathscr{N}^{J}(\Phi)=\max_{\rho(0)}\int_{\Gamma_{\rho(0)}}dt\frac{d}{dt}[-H(\rho(t))],$
(32)

where $\Gamma_{\rho(0)}=\set{t\in\mathbb{R}}{\frac{d}{dt}H(\rho(t))<0}$. The
measure of non-Markovianity for the JSD in the case of unital dynamics is
given by the total decrease of entropy for a single state, maximised over all
possible initial states. This is clearly linked to violations of P
divisibility of $\Phi$, since any unital positive map acting on qubits
increases the entropy of the state [21], so that any revival in the entropy
must necessarily correspond to a violation of P divisibility. Unfortunately,
this feature is only true for qubits, since in higher dimensions orthogonal
states do not need to have the same eigenvalues. Additionally, no similar
interpretation holds for $\sqrt{\text{JSD}}$ or for the two generalisations to
ensembles $K_{\mu}$ and $S_{\mu}$.

###  II.5 Robustness of optimal pairs

Let us now study the robustness of optimal pairs, i.e. how the measure of non-
Markovianity changes when moving away from the optimal pair, for the different
distinguishability quantifiers. We will illustrate this by considering a
simple but yet paradigmatic model: the dephasing model. This model consists in
a modification of the coherences without a corresponding change in the
populations:

$\rho(t)=\begin{pmatrix}\rho_{00}&\rho_{01}\gamma(t)e^{-i\omega_{S}t}\\\
\rho_{10}\gamma^{*}(t)e^{i\omega_{S}t}&\rho_{11}\end{pmatrix},$ (33)

where $\gamma$ is called the decoherence function. For this model, non-
Markovianity corresponds to a non-monotonic behaviour of $\lvert\gamma\rvert$.
It is worth stressing that considerations similar to the ones for this model,
are also valid for other models such as, for example, the phase covariant
model which will be introduced in Sec. III.2.

Figure 3: Measure of non-Markovianity for pure and orthogonal states for TD
(left), JSD (middle), and $\sqrt{\text{JSD}}$ (right). The color provide the
value of the measure of non-Markovianity, rescaled according to the maximum
value reached by the measure for the considered distinguishability quantifier
$d$, obtained considering as initial states the corresponding point on the
surface of the Bloch sphere and its antipodal point. Brighter colors
correspond to higher values of the revivals. The reference values are taken to
be $1\times 10^{-2}$ for $\mathscr{N}^{D}(\Phi)$ as well as
$\mathscr{N}^{\sqrt{J}}(\Phi)$, and $1\times 10^{-3}$ for
$\mathscr{N}^{J}(\Phi)$, reflecting the different scales of the revivals. The
two distances (TD and $\sqrt{\text{JSD}}$) behave very similarly. For the JSD,
on the other hand, the value of the measure of non-Markovianity decreases more
slowly when moving away from the optimal pair, i.e. all the pairs of pure and
orthogonal vectors on the equator of the Bloch sphere. The dynamics $\Phi$ is
given by the dephasing model of Eq. (33) with decoherence function $\gamma(t)$
corresponding to a bosonic bath described by a spectral density with an
exponential cut-off of the form
$J(\omega)=\lambda({\omega^{s}}/{\Omega^{s-1}})\exp({-{\omega}/{\Omega}})$ as
considered in [40, 41], with $s=3$, $\lambda=3$ and $\Omega=1$ in inverse
units of time.

Optimal pairs are all the pairs of pure and orthogonal states corresponding to
antipodal vectors on the equator of the Bloch sphere, since the $x$ and $y$
direction are the only ones in which the dynamics is not trivial. In order to
evaluate the robustness of the optimal pairs, Fig. 3 shows the behavior of the
measure of non-Markovianity when moving away from the equatorial plane, but
still considering pure and orthogonal states. It is possible to notice that
the TD and the $\sqrt{\text{JSD}}$, being both distances, behave very
similarly, with the maximum of the measure of non-Markovianity on the
equatorial plane, and quickly decreasing when moving towards the poles of the
Bloch sphere. For the JSD, on the other hand, the situation is qualitatively
different, with a broader region around the equator with a measured value of
non-Markovianity similar to the maximal one, which is about one order of
magnitude smaller than the value obtained for the other quantifiers.

##  III Non-Markovianity and divisibility

In the definition of a measure of non-Markovianity for a generic
distinguishability quantifier $d$, the condition of contractivity (24) implies
that every P divisible dynamics leads to a zero measure of non-Markovianity.
On the other hand, using the Helstrom matrix, as soon as the dynamics violates
P divisibility one has a non-zero measure of non-Markovianity,
$\mathscr{N}^{\Delta}(\Phi)\char 62\relax 0$. In other words,
$\mathscr{N}^{\Delta}$ is stronger than $\mathscr{N}^{d}$ for any other
quantifier $d$. We now want to study whether there exists other quantifiers
$d$ leading to measures that are equivalent to the one arising from the
Helstrom matrix, with particular focus on the entropic quantifiers.

We already know that the properties 1-3 are not sufficient in order to have a
measure of non-Markovianity equivalent to $\mathscr{N}^{\Delta}$, since it is
strictly stronger than $\mathscr{N}^{D}$. In Sec. III.3 we will show with a
counterexample that the same also holds for the JSD and its generalisations.

###  III.1 Positivity and non-contractivity domain

Let us first tackle the question of the behavior of the JSD under non positive
maps. We already know that the JSD is contractive under any positive map. We
now want to investigate if the reverse is also true, namely we want to clarify
whether, for any non positive map $\Lambda$, there exists a pair of states for
which the JSD is strictly non contractive. Non-positivity of $\Lambda$ implies
that there exists some state $\rho$ which is mapped to a non-positive operator
$\Lambda\rho$. However, the JSD, unlike the TD, cannot be extended to non-
positive operators, since it involves the logarithm of the eigenvalues.
Therefore, the search for a non-contractive pair for $\Lambda$ must be
restricted to the set of states that are mapped to states after the action of
the map, i.e. to the positivity domain

$\mathcal{PD}_{\Lambda}=\set{\rho\in\mathcal{S}(\mathscr{H})}{\Lambda\rho\in\mathcal{S}(\mathscr{H})},$
(34)

where $\mathcal{S}(\mathscr{H})$ is the set of quantum states on a Hilbert
space $\mathscr{H}$. In the following, we will only consider qubits
$\mathscr{H}=\mathbb{C}^{2}$, since this will turn out to be sufficient in
order to show that $\mathscr{N}^{J}$ is strictly weaker than
$\mathscr{N}^{\Delta}$. We denote the set of all qubit states, i.e. the Bloch
sphere, as $\mathcal{S}(\mathbb{C}^{2})=\mathcal{S}$.

Considering unital non-positive maps $\Lambda$, it is easy to show that there
always exists a non-contractive pair inside $\mathcal{PD}_{\Lambda}$. Such
maps act on Bloch vectors according to Eq. (29) with $\boldsymbol{\kappa}=0$
and non-positivity implies that some $\lambda_{i}\char 62\relax 1$, which we
take to be $\lambda_{1}$, without loss of generality. The non-contractive pair
is the one represented by the Bloch vectors
$\boldsymbol{r}_{\rho}=(\lambda_{1}^{-1},0,0)^{\top}=-\boldsymbol{r}_{\sigma}$.
In fact, by direct calculation it is easy to show that
$J(\rho,\sigma)<J(\Lambda\rho,\Lambda\sigma)=1$. In the general case, an
analytic proof for the existence of a non-contractive pair is missing.
However, by parameterising the non positive map $\Lambda$ as in Eq. (29) and
performing a sample on all the parameters, we observed numerically that for
any such map it is always possible to find a pair of states
$\rho,\sigma\in\mathcal{PD}_{\Lambda}$ such that
$J(\Lambda\rho,\Lambda\sigma)\char 62\relax J(\rho,\sigma)$.

Turning back to the dynamical point of view, however, the search for the non-
contractive pair might not be extended to all $\mathcal{PD}_{\Lambda}$. In
fact, not all the domain of positivity of $\Lambda=\Phi_{t,s}$ is available,
but only the image at time $s$ of the Bloch sphere $\Phi_{s}(\mathcal{S})$ is.
We stress that $\Phi_{s}(\mathcal{S})$ is in general only a subset of
$\mathcal{PD}_{\Lambda}$. Thus, in order to have $\mathscr{N}^{J}(\Phi)\char
62\relax 0$ for all non-P divisible processes, we would need to be able to
find a non-contractive pair for the JSD inside $\Phi_{s}(\mathcal{S})$. Let us
now define the set of states in which it is possible to find a non-contractive
pair as the non-contractivity domain

$\begin{split}\mathcal{NCD}_{\Lambda,J}=\\{\rho\in&\mathcal{PD}_{\Lambda}\text{
}|\text{ }\exists\sigma\in\mathcal{PD}_{\Lambda},\\\
&J(\Lambda\rho,\Lambda\sigma)\char 62\relax J(\rho,\sigma)\\}.\end{split}$
(35)

Therefore, in order to have non-Markovianity for all non-P divisible dynamical
maps, we would need to have
$\mathcal{NCD}_{\Lambda,J}=\mathcal{PD}_{\Lambda}$: for any state in
$\Phi_{s}(\mathcal{S})\cap\mathcal{PD}_{\Lambda}$ it is always possible to
find another state such that non-contractivity holds. However, this is not the
case, as it is clear from the example of Fig. 4. There, in fact,
$\mathcal{NCD}_{\Lambda,J}$ is a proper subset of $\mathcal{PD}_{\Lambda}$.
Therefore, if we were able to construct a dynamics $\Phi$ such that for times
$t\char 62\relax s\char 62\relax 0$ it acts as this non-positive map
($\Phi_{t,s}=\Lambda$), with a dynamics prior to time $s$ that is P divisible
and with the Bloch sphere that is mapped at time $s$ inside
$\mathcal{PD}_{\Lambda}$ but outside $\mathcal{NCD}_{\Lambda,J}$, i.e.
$\Phi_{s}(\mathcal{S})\subset\mathcal{PD}_{\Lambda}\setminus\mathcal{NCD}_{\Lambda,J}$,
we would construct a non-P divisible dynamics but with
$\mathscr{N}^{J}(\Phi)=0$. This is indeed feasible as we will show in Sec.
III.3 providing explicitly a model which is similar in spirit to the one just
described.

Figure 4: Section at $y=0$ of the Bloch sphere (light blue) for an example of
a non-positive map $\Lambda$ for which the non-contractivity domain
$\mathcal{NCD}_{\Lambda,J}$ (green) is strictly included in the positive
domain $\mathcal{PD}_{\Lambda}$ (blue). This map acts on a Bloch vector
$\boldsymbol{r}=(x,y,z)^{\top}$ as
$\boldsymbol{r}\mapsto(\lambda_{x}x,\lambda_{y}y,\lambda_{z}z)^{\top}$, with
$\lambda_{x}=\lambda_{y}=1.1$, and $\lambda_{z}=0.1$.

###  III.2 Phase covariant dynamics

In order to construct the counterexample of Sec. III.3, let us first set the
theoretical background of the considered dynamics, namely phase covariant
ones. They contain a broad class of dynamics and they involve maps $\Phi$ that
satisfy covariance with respect to phase transformations, namely [42]

$e^{-i\sigma_{z}\theta}\Phi_{t}[\rho]e^{i\sigma_{z}\theta}=\Phi_{t}\left[e^{-i\sigma_{z}\theta}\rho
e^{i\sigma_{z}\theta}\right]$ (36)

for all real $\theta$ and for all states $\rho\in\mathcal{S}(\mathbb{C}^{2})$.
Phase covariant dynamics are in the form [43, 44]

$\Phi_{t}\rho=\frac{1}{2}\big{[}\mathbbm{1}+\eta_{\perp}(t)(\nu_{x}\sigma_{x}+\nu_{y}\sigma_{y})+\eta_{\parallel}(t)\nu_{z}\sigma_{z}+\kappa_{z}(t)\sigma_{z}\big{]},$
(37)

where $\nu_{i}=\operatorname{tr}[\rho\sigma_{i}]$, for $i=x,y,z$. The complete
positivity conditions reads

$\eta_{\parallel}\pm\kappa_{z}\leq 1,\qquad
1+\eta_{\parallel}\geq\sqrt{4\eta_{\perp}^{2}+\kappa_{z}^{2}}.$ (38)

The dynamics can be reformulated in terms of a master equation of the form

$\displaystyle\frac{d\rho}{dt}=$
$\displaystyle\gamma_{+}(t)\left(\sigma_{+}\rho\sigma_{-}-\frac{1}{2}\\{\rho,\sigma_{-}\sigma_{+}\\}\right)$
(39)
$\displaystyle+\gamma_{-}(t)\left(\sigma_{-}\rho\sigma_{+}-\frac{1}{2}\\{\rho,\sigma_{+}\sigma_{-}\\}\right)+\gamma_{z}(t)(\sigma_{z}\rho\sigma_{z}-\rho),$

where

$\gamma_{\pm}(t)=\frac{\eta_{\parallel}(t)}{2}\frac{d}{dt}\left(\frac{1\pm\kappa_{z}(t)}{\eta_{\parallel}(t)}\right),\quad\gamma_{z}(t)=\frac{1}{4}\frac{d}{dt}\ln\frac{\eta_{\parallel}(t)}{\eta_{\perp}^{2}(t)}.$
(40)

The dynamics is CP divisible if and only if $\gamma_{\pm}(t)\geq 0$ and
$\gamma_{z}(t)\geq 0$. P divisibility, instead, is satisfied whenever [42]

$\gamma_{\pm}(t)\geq
0\quad\text{and}\quad\sqrt{\gamma_{+}(t)\gamma_{-}(t)}+2\gamma_{z}(t)\char
62\relax 0.$ (41)

The composition of two phase covariant dynamics is again phase covariant. If
we suppose that the system undergoes a first phase covariant dynamics
$\Phi^{1}$ from $t=0$ to $t=t_{1}$, and later it evolves following $\Phi^{2}$,
then the total dynamics $\Phi=\Phi^{2}\circ\Phi^{1}$, defined as

$\Phi_{t}=\begin{cases}\Phi^{1}_{t},\quad&\text{if }t\leq t_{1}\\\
\Phi^{2}_{t-t_{1}}\Phi^{1}_{t_{1}},\quad&\text{if }t\char 62\relax
t_{1}\end{cases},$ (42)

is again phase covariant, described by the functions

$\displaystyle\eta_{\parallel,\perp}(t)=\begin{cases}\eta_{\parallel,\perp}^{1}(t)\quad&\text{if
}t\leq t_{1}\\\
\eta_{\parallel,\perp}^{2}(t-t_{1})\eta_{\parallel,\perp}^{1}(t_{1})\quad&\text{if
}t\char 62\relax t_{1}\end{cases},$ (43)
$\displaystyle\kappa_{z}(t)=\begin{cases}\kappa^{1}_{z}(t)\quad&\text{if
}t\leq t_{1}\\\
\kappa^{2}_{z}(t-t_{1})+\eta_{\parallel}^{2}(t-t_{1})\kappa^{1}_{z}(t_{1})\quad&\text{if
}t\char 62\relax t_{1}\end{cases},$ (44)

where the superscripts 1 or 2 label the functions defining respectively
$\Phi^{1}$ or $\Phi^{2}$. The composition of two phase covariant dynamics is
not commutative, since in general
$\Phi^{1}\circ\Phi^{2}\neq\Phi^{2}\circ\Phi^{1}$ as it is evident from (43)
and (44). Furthermore, the family of all phase covariant dynamics do not form
a group, since, in general, the inverse of a dynamics $\Phi^{-1}$ is not
positive.

###  III.3 Example showing that $\mathscr{N}^{J}$ is strictly weaker than
$\mathscr{N}^{\Delta}$

Let us now employ the previously introduced phase covariant model in order to
build a counterexample of a dynamics which is not P divisible but yet leading
to a zero measure of non-Markovianity for the JSD. It follows from this
counterexample that $\mathscr{N}^{J}$ is a strictly weaker measure of non-
Markovianity than $\mathscr{N}^{\Delta}$, since $\mathscr{N}^{\Delta}\char
62\relax 0$ for all non-P divisible dynamics. We will actually consider a
dynamics for which the memory effects are already detected by the TD, without
the need to generalise to the Helstrom matrix.

Figure 5: Schematic representation of the dynamics of the counterexample,
obtained by visualizing the Bloch sphere (blue) and the evolved ellipsoid
(red) at different time steps. After time $t=0$ and till time $T_{nM}$ the
image of the Bloch sphere is a uniformly contracted sphere translated in the
$z$ direction, without violating P divisibility. After $T_{nM}$, the image of
the sphere is further contracted and translated in the $z$ direction, but also
expanded in the $x$ and $y$ directions, thus violating P divisibility. In the
third panel, the dashed line represents the border of the non-contractivity
domain for the dynamical map from $T_{nM}$ onwards, corresponding to the green
area in Fig. 4. It thus clearly appears that the Bloch sphere at $T_{nM}$ is
mapped outside the non-contractivity domain. This way, even if the second part
of the dynamics is not P divisible, there is no pair of states available for
the JSD to witness a revival.

Figure 6: Left panel: The functions of Eqs. (45) and (46). Non-Markovianity is
present after $T_{nM}\approx 2.2T$ and corresponds to a non-monotonic behavior
of $\eta_{\perp}$. Right panel: Rates $\gamma_{\pm}$ and $\gamma_{z}$ for the
corresponding dynamics. The red line corresponds to
$\sqrt{\gamma_{+}\gamma_{-}}+2\gamma_{z}$, whose positivity is necessary for
the dynamics to be P divisible, together with the conditions
$\gamma_{\pm}\char 62\relax 0$. Non-Markovianity for $t\char 62\relax T_{nM}$
corresponds to a violation of divisibility. We considered the choice of the
parameters: $\mu_{1}=5$, $\mu_{2}=4$, $A_{\parallel}=0.01$, $A_{\perp}=1.01$,
$A_{\kappa}=0.45$ and $\alpha=5$.

The dynamical map $\Phi$ of such counterexample is described by the functions

$\displaystyle\begin{split}\eta_{\parallel,\perp}(\tau)=&e^{-\mu_{1}\tau}\sigma(1-\tau)\\\
&+e^{-\mu_{1}}e^{-\mu_{2}(\tau-1)}\sigma(\tau-1)\sigma(2-\tau)\\\
&+e^{-\mu_{1}-\mu_{2}}\left[(3-\tau)+A_{\parallel,\perp}(\tau-2)\right]\sigma(\tau-2),\end{split}$
(45)
$\displaystyle\begin{split}\kappa_{z}(\tau)=&A_{\kappa}\tau\sigma(2-\tau)\\\
&+2A_{\kappa}\left[(3-\tau)+A_{\parallel}(\tau-2)\right]\sigma(\tau-2),\end{split}$
(46)

where $\sigma$ is the sigmoid function

$\sigma(\tau)=\frac{1}{1+e^{-\alpha\tau}},$

which is a smooth version of the Heaviside theta function and $\tau=t/T$ is a
dimensionless time parameter, where $T$ is a reference time determining the
duration of the different stages depicted in Fig. 5. These functions together
with the corresponding rates $\gamma_{\pm}$ and $\gamma_{z}$ obtained from Eq.
(40) are shown in Fig. 6.

The idea of the counterexample follows from the considerations of Sec. III.1:
the dynamics in the time interval in which the memory effects arise will
consist of a non-positive map $\Phi_{t,s}$ similar to the one described in
Fig. 4, for which the non-contractivity domain (35) is strictly smaller than
the positivity domain (34). Prior to this time interval, the dynamics is P
divisible and such that it maps the whole Bloch sphere inside
$\mathcal{PD}_{\Phi_{t,s}}$ but outside $\mathcal{NCD}_{\Phi_{t,s},J}$, so
that there is no pair of states available for the JSD to witness the
violations of P divisibility. A schematic representation of such dynamics is
shown in Fig. 5.

The violation of P divisibility takes place for $t\char 62\relax
T_{\text{nM}}$, with $T_{\text{nM}}\approx 2.2T$ as (41) is violated, and is
due to positivity of the time derivative of $\eta_{\perp}$, in turn leading to
negativity of the second condition appearing in (41). This behavior is shown
in Fig. 6. Such violation corresponds to a revival in the coherences, without
a corresponding revival in the population, thus building on a genuine quantum
effect. The fact that memory effects are due to a unital feature of the map,
and not to the translation $\kappa_{z}$, implies that
$\mathscr{N}^{D}(\Phi)\char 62\relax 0$, which, in turn, leads to
$\mathscr{N}^{\Delta}(\Phi)\char 62\relax 0$.

On the other hand, we evaluated numerically that $\mathscr{N}^{J}(\Phi)=0$, as
it can be seen in Fig. 7. The numerical analysis has been performed
considering all the possible initial pairs of states on the Bloch sphere,
studying their time evolution and evaluating any eventual revival of the JSD,
but none has been found. Therefore, we can conclude that the measure of non-
Markovianity arising from the JSD is strictly weaker than the one arising from
the Helstrom matrix. Or, in other words, there exists non-P divisible dynamics
leading to a zero measure of non-Markovianity.

Nevertheless, this fact is also true for the TD. In order to be able to
capture any violation of P divisibility as a revival of some quantifier, one
has to generalise the TD to ensembles, introducing a bias parameter. One might
wonder if something similar also happens for the JSD: if we generalise it to
ensembles as the Holevo skew divergence (13) or as the quantum skew divergence
(16), could we be able to witness all the violations of P divisibility?
Unfortunately, the answer to this question is no. In fact, considering the
same counterexample, one has again
$\mathscr{N}^{K_{\mu}}(\Phi)=\mathscr{N}^{S_{\mu}}(\Phi)=0$. This fact
actually comes unsurprisingly: for the TD, the generalisation to ensembles
breaks the translational symmetry and makes us able to detect violations of P
divisibility due to the translational components of the dynamics; for the JSD,
on the other hand, there is no symmetry to break and thus generalising it to
ensembles does not lead to any qualitative difference. In particular, one can
consider a time evolution $\tilde{\Phi}$ such that
$\mathscr{N}^{J}(\tilde{\Phi})\char 62\relax 0$, while
$\mathscr{N}^{D}(\tilde{\Phi})=0$, as shown in [9] considering a different
phase covariant model.

A natural question is what additional constraints $d$ has to obey in order to
have a measure of non-Markovianity equivalent to $\mathscr{N}^{\Delta}$.
Building on our counterexample, a necessary condition that $d$ must obey is
naturally the existence of a strictly non-contractive pair of states for any
non-positive map $\Lambda$. A second condition, crucial for entropic
distinguishability quantifiers, is the relation between the non-contractivity
domain $\mathcal{NCD}_{\Lambda,d}$, depending on both the map and the
quantifier, and the positivity domain $\mathcal{PD}_{\Lambda}$ of the map. As
we have shown, if $\mathcal{NCD}_{\Lambda,d}$ is a proper subset of
$\mathcal{PD}_{\Lambda}$, detection of violation of divisibility after a time
$s$ can fail for maps whose image at time $s$, by necessity within the
positivity domain, is strictly outside the non-contractivity domain. This is
exactly the feature we exploited to provide the counterexample.

Figure 7: Non-Markovianity measure for the JSD (solid line) and for the TD
(dashed line) for the considered counterexample as a function of time.
Clearly, $\mathscr{N}^{D}(\Phi)\char 62\relax 0$, since a revival in the TD is
witnessed. The JSD, on the other hand is always a monotonic function and hence
$\mathscr{N}^{J}(\Phi)=0$.

##  IV Conclusions and outlook

In this work, we have compared different measures for the degree of non-
Markovianity in the dynamics of open systems based on distinguishability
quantifiers between quantum states. In particular, we have provided evidence
that the measure based on the trace norm of the Helstrom matrix is strictly
stronger than all of the measures of non-Markovianity based on entropic
distinguishability quantifiers, as well as stronger than the measure based on
the trace distance, which is neither stronger nor weaker than the entropic
ones. This is our central result. It means that the value of the measure based
on the trace norm of the Helstrom matrix associated to a dynamical map $\Phi$,
namely $\mathscr{N}^{\Delta}(\Phi)$, is greater than zero whenever this
happens for the measures associated to entropic distinguishability quantifiers
or to the trace distance, while the reverse is not true as we have
demonstrated here. Thus, we can conclude that the different distinguishability
measures exhibit a quite different performance in the detection of non-
positive maps, which is surprising in view of similar physical interpretations
outlined in Secs. II.1 and II.2.

This result has been obtained considering the explicit expression of a qubit
dynamics which is not P divisible, so that it is non-Markovian according to
$\mathscr{N}^{\Delta}$, though strong numerical evidence shows that the
associated entropic non-Markovianity measure $\mathscr{N}^{J}$, based on the
Jensen-Shannon divergence as distinguishability quantifier, is equal to zero.
A purely analytic proof of this property, in the present or in a different
counterexample, would provide further insights on the relationship between
positivity of a map and its contractivity property with respect to entropic
distinguishability quantifiers obtained from the quantum relative entropy,
such as the Jensen-Shannon divergence.

Moreover, it might be very interesting to clarify if the Helstrom based
quantifier is unique, or whether other distinguishability quantifiers are
equivalent to it as measures for quantum non-Markovianity.

###### Acknowledgements.

The authors thank Andrea Smirne for many fruitful discussions and careful
reading of the manuscript.

## References

  * Breuer  _et al._ [2016] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Colloquium : Non-markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * Rivas  _et al._ [2014] Á. Rivas, S. F. Huelga, and M. B. Plenio, Quantum non-markovianity: characterization, quantification and detection, Reports on Progress in Physics 77, 094001 (2014).
  * Li  _et al._ [2018] L. Li, M. J. Hall, and H. M. Wiseman, Concepts of quantum non-markovianity: A hierarchy, Physics Reports 759, 1 (2018).
  * Milz and Modi [2021] S. Milz and K. Modi, Quantum stochastic processes and quantum non-markovian phenomena, PRX Quantum 2, 030201 (2021).
  * Breuer  _et al._ [2009] H.-P. Breuer, E.-M. Laine, and J. Piilo, Measure for the degree of non-markovian behavior of quantum processes in open systems, Phys. Rev. Lett. 103, 210401 (2009).
  * Liu  _et al._ [2013] J. Liu, X.-M. Lu, and X. Wang, Nonunital non-markovianity of quantum dynamics, Phys. Rev. A 87, 042103 (2013).
  * Chruściński  _et al._ [2011] D. Chruściński, A. Kossakowski, and Á. Rivas, Measures of non-markovianity: Divisibility versus backflow of information, Physical Review A 83, 052128 (2011).
  * Wißmann  _et al._ [2015] S. Wißmann, H.-P. Breuer, and B. Vacchini, Generalized trace-distance measure connecting quantum and classical non-markovianity, Phys. Rev. A 92, 042108 (2015).
  * Megier  _et al._ [2021] N. Megier, A. Smirne, and B. Vacchini, Entropic bounds on information backflow, Phys. Rev. Lett. 127, 030401 (2021).
  * Briët and Harremoës [2009] J. Briët and P. Harremoës, Properties of classical and quantum jensen-shannon divergence, Phys. Rev. A 79, 052311 (2009).
  * Sra [2021] S. Sra, Metrics induced by jensen-shannon and related divergences on positive definite matrices, Linear Algebra and its Applications 616, 125 (2021).
  * Virosztek [2021] D. Virosztek, The metric property of the quantum jensen-shannon divergence, Advances in Mathematics 380, 107595 (2021).
  * Heinosaari and Ziman [2011] T. Heinosaari and M. Ziman,  _The Mathematical Language of Quantum Theory: From Uncertainty to Entanglement_ (Cambridge University Press, 2011).
  * Ruskai [1994] M. B. Ruskai, Beyond strong subadditivity? improved bounds on the contraction of generalized relative entropy, Reviews in Mathematical Physics 6, 1147 (1994).
  * Molnár and Timmermann [2002] L. Molnár and W. Timmermann, Isometries of quantum states, Journal of Physics A: Mathematical and General 36, 267 (2002).
  * Fuchs and van de Graaf [1999] C. A. Fuchs and J. van de Graaf, Cryptographic distinguishability measures for quantum-mechanical states, IEEE Transactions on Information Theory 45, 1216 (1999).
  * Helstrom [1976] C. W. Helstrom,  _Quantum detection and estimation theory_ (Academic Press New York, 1976).
  * Helstrom [1967] C. W. Helstrom, Detection theory and quantum mechanics, Information and Control 10, 254 (1967).
  * Müller-Hermes and Reeb [2017] A. Müller-Hermes and D. Reeb, Monotonicity of the quantum relative entropy under positive maps, Annales Henri Poincaré 18, 1777 (2017).
  * Audenaert [2014] K. M. Audenaert, Quantum skew divergence, Journal of Mathematical Physics 55, 112202 (2014).
  * Bengtsson and Życzkowski [2017] I. Bengtsson and K. Życzkowski,  _Geometry of quantum states: an introduction to quantum entanglement_ (Cambridge university press, 2017).
  * Hayashi [2006] M. Hayashi,  _Quantum Information_ (Springer-Verlag, Berlin, 2006).
  * Majtey  _et al._ [2005] A. P. Majtey, P. W. Lamberti, and D. P. Prato, Jensen-shannon divergence as a measure of distinguishability between mixed quantum states, Phys. Rev. A 72, 052310 (2005).
  * Audenaert [2011a] K. M. R. Audenaert, Telescopic relative entropy–ii triangle inequalities, e-print arXiv:1102.3041 (2011a).
  * Pinsker and Feinstein [1964] M. Pinsker and A. Feinstein,  _Information and Information Stability of Random Variables and Processes_ (Holden-Day, 1964).
  * Smirne  _et al._ [2022] A. Smirne, N. Megier, and B. Vacchini, Holevo skew divergence for the characterization of information backflow, Phys. Rev. A 106, 012205 (2022).
  * Audenaert [2011b] K. M. Audenaert, Telescopic relative entropy, in  _Conference on Quantum Computation, Communication, and Cryptography_ (Springer, 2011) pp. 39–52.
  * Lee [1999] L. Lee, Measures of distributional similarity, in  _Proceedings of the 37th Annual Meeting of the Association for Computational Linguistics_ (Association for Computational Linguistics, College Park, Maryland, USA, 1999) pp. 25–32.
  * Holevo [1973] A. S. Holevo, Bounds for the quantity of information transmitted by a quantum communication channel, Problemy Peredachi Informatsii 9, 3 (1973).
  * Pollock  _et al._ [2018] F. A. Pollock, C. Rodriguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Operational markov condition for quantum processes, Phys. Rev. Lett. 120, 040405 (2018).
  * Budini [2018] A. A. Budini, Quantum non-markovian processes break conditional past-future independence, Phys. Rev. Lett. 121, 240401 (2018).
  * Budini [2022] A. A. Budini, Quantum non-markovian environment-to-system backflows of information: Nonoperational vs. operational approaches, Entropy 24, 649 (2022).
  * Wißmann  _et al._ [2012] S. Wißmann, A. Karlsson, E.-M. Laine, J. Piilo, and H.-P. Breuer, Optimal state pairs for non-markovian quantum dynamics, Phys. Rev. A 86, 062108 (2012).
  * Laine  _et al._ [2010] E.-M. Laine, J. Piilo, and H.-P. Breuer, Witness for initial system-environment correlations in open-system dynamics, EPL (Europhysics Letters) 92, 60010 (2010).
  * Amato  _et al._ [2018] G. Amato, H.-P. Breuer, and B. Vacchini, Generalized trace distance approach to quantum non-markovianity and detection of initial correlations, Phys. Rev. A 98, 012120 (2018).
  * Campbell  _et al._ [2019] S. Campbell, M. Popovic, D. Tamascelli, and B. Vacchini, Precursors of non-markovianity, New Journal of Physics 21, 053036 (2019).
  * Kossakowski [1972a] A. Kossakowski, On necessary and sufficient conditions for a generator of a quantum dynamical semi-group, Bull. Acad. Pol. Sci., Sér. Sci. Math. Astron. Phys. 20, 1021 (1972a).
  * Kossakowski [1972b] A. Kossakowski, On quantum statistical mechanics of non-hamiltonian systems, Reports on Mathematical Physics 3, 247 (1972b).
  * King and Ruskai [2001] C. King and M. Ruskai, Minimal entropy of states emerging from noisy quantum channels, IEEE Transactions on Information Theory 47, 192 (2001).
  * Addis  _et al._ [2014] C. Addis, B. Bylicka, D. Chruściński, and S. Maniscalco, Comparative study of non-markovianity measures in exactly solvable one- and two-qubit models, Phys. Rev. A 90, 052103 (2014).
  * Guarnieri  _et al._ [2014] G. Guarnieri, A. Smirne, and B. Vacchini, Quantum regression theorem and non-Markovianity of quantum dynamics, Physical Review A - Atomic, Molecular, and Optical Physics 90, 1 (2014).
  * Filippov  _et al._ [2020] S. N. Filippov, A. N. Glinov, and L. Leppäjärvi, Phase covariant qubit dynamics and divisibility, Lobachevskii J. Math. 41, 617 (2020).
  * Haase  _et al._ [2018] J. F. Haase, A. Smirne, J. Kołodyński, R. Demkowicz-Dobrzański, and S. F. Huelga, Fundamental limits to frequency estimation: a comprehensive microscopic perspective, New Journal of Physics 20, 053009 (2018).
  * Smirne  _et al._ [2016] A. Smirne, J. Kołodyński, S. F. Huelga, and R. Demkowicz-Dobrzański, Ultimate Precision Limits for Noisy Frequency Estimation, Phys. Rev. Lett. 116, 1 (2016).


# Multi-Objective and Model-Predictive Tree Search for Spatiotemporal
Informative Planning

Weizhe Chen and Lantao Liu  Indiana University, Bloomington, IN, USA, 47408
Emails: {chenweiz<EMAIL_ADDRESS>
###### Abstract

Adaptive sampling and planning in robotic environmental monitoring are
challenging when the target environmental process varies over space and time.
The underlying environmental dynamics require the planning module to integrate
future environmental changes so that action decisions made earlier do not
quickly become outdated. We propose a Monte Carlo tree search method which not
only well balances the environment exploration and exploitation in space, but
also catches up to the temporal environmental dynamics. This is achieved by
incorporating multi-objective optimization and a look-ahead model-predictive
rewarding mechanism. We show that by allowing the robot to leverage the
simulated and predicted spatiotemporal environmental process, the proposed
informative planning approach achieves a superior performance after comparing
with other baseline methods in terms of the root mean square error of the
environment model and the distance to the ground truth.

##  I Introduction and Related Work

Environmental monitoring allows humans to obtain a good understanding of the
underlying environmental state which oftentimes is spatially and temporally
dynamic. Capturing and utilizing the environmental dynamics are vital. For
example, when an oil spill happens in oceans, the search-and-rescue crew need
to estimate the oil diffusion and movement by integrating environmental
information such as ocean currents, surface winds, and temperature/weather to
best approximate the spill spreading behavior. Autonomous robotic systems are
increasingly being utilized as scientific data gathering tools, because they
are efficient, economical, and able to monitor the environment that is too
dangerous, or impossible, for humans to tackle [10]. In robotic environmental
monitoring, the mobile robot needs to build a “distribution map” for the
environmental attribute of interest, based on which it is able to plan the
“best” sampling trajectory. Most existing literature uses Gaussian processes
(GPs) to construct the distribution map and uses adaptive sampling and
planning methods to navigate the robot [15, 22]. GPs not only provide an
estimate of the environmental attribute but also quantify the estimation
uncertainty. Given the estimation uncertainty, many information metrics,
including but not limited to entropy, variance reduction, and mutual
information can be derived. Such information-theoretic planning framework is
also called informative planning [19, 12, 17, 7], and it can effectively guide
the mobile robot to the areas with the largest amount of information (or
uncertainty reduction), leading to a good spatial coverage even within a short
time window.

However, in the aforementioned environmental monitoring scenarios, good
spatial coverage is insufficient to achieve satisfying modeling performance.
The robot needs to explore the environment to discover the areas with high
concentration (namely, hotspots) first, and then visit different areas with
differing efforts to obtain the most valuable samples while catching up with
the spatiotemporal environmental dynamics.

Figure 1: Illustration of current estimation based Monte Carlo tree search
(MCTS) and model-predictive tree search in an oil spill monitoring task. A
large volume of spilled oil usually ends into multiple area patches on the
surface which keep moving with eddies, see the oil spill trajectory from NOAA
for example. The first column of this picture shows the unknown ground truth,
where the arrows represent a vector field of the ocean currents, and the
contour lines indicate a hotspot (patch) of spilled oil. In the second column,
the bottom left image shows the robot’s current estimation of the environment
and the MCTS. The model-predictive tree search is carried out in the robot’s
prediction (upper layer) which comes from its current estimation of the oil
density and the (imperfect) ocean current forecast. Model-predictive tree
search encourages the robot to visit the area where the oil is likely to exist
while current estimation based MCTS guides the robot to where the oil has ever
been.

Recently, Monte Carlo tree search (MCTS) has gained great popularity as an
online planning method. It has been applied to spatiotemporal monitoring [18],
planetary exploration mission [1], environmental exploration [8], and active
object recognition [4, 20]. MCTS reveals a significant advantage for improving
the performance of game-playing agents, especially in the game of Go [21],
based on perfect simulating environments with well-defined rules. However,
existing MCTS planners in robotics applications are typically based on the
robot’s current estimation of the environment without considering much on the
(possibly highly) dynamic nature (see Fig. 1). This motivates us to develop an
informative planning method that incorporates environmental prediction model
where environmental dynamics are integrated with the MCTS working process.

Predicting the environmental process over space and time is a challenging
task. Nevertheless, there is a type of spatiotemporal environmental process
whose dynamics are mainly affected by the motion of the environmental medium
(e.g., fluid, air) in which it is located. Fortunately, for these types of
dynamic processes, simulation modeling is able to roughly estimate the
possible trajectory of the feature of interest. To name a few, such phenomena
include spilled oil driven by ocean currents, wildfire affected by wind, and
air pollution propagated by the prevailing wind. Thanks to the efforts in
these domains, we normally have access to some simulators that provide an
(imperfect) estimation of the environmental process given current wind field
or ocean currents data. For example, General NOAA Operational Modeling
Environment (GNOME) [3] is the modeling tool that predicts how wind, currents,
and other processes might move and spread the spilled oil. Fire Area Simulator
(FARSITE) [11] computes wildfire growth and behavior using wind, fuel
moistures, and topography data. An operational algal bloom forecast system is
able to forecast the expansion of harmful algal blooms using a combination of
satellite imaging and wind prediction [23].

This paper presents a multi-objective variant of MCTS with a simulator that
enables the robot to incorporate its prediction of environmental state based
on prior or related information. The proposed method also allows the robot to
handle potentially competing objectives such as exploration versus
exploitation. In addition, the algorithm can return a planning solution even
if it is interrupted prematurely, making it suitable for online planning and
in-sit decision-making.

##  II PRELIMINARY

###  II-A Gaussian Process Regression in Spatial Modeling

In the spatial modeling task, given a training set $\mathcal{D}$ of $n$
observations, $\mathcal{D}=\\{(\bm{x}_{i},y_{i})|i=1,\dots,n\\}$, where
$\bm{x}$ denotes an input location of dimension $2$ (or $3$) and $y$ denotes a
scalar output or target, we want to infer the relationship between input and
output. Samples from the real world are assumed to be noisy, i.e.
$y=f(\bm{x})+\epsilon,$ where $f(\cdot)$ is the unknown function of interest
and $\epsilon$ is the i.i.d. Gaussian noise with zero mean and variance
$\sigma_{n}^{2}$, $\epsilon\sim\mathcal{N}(0,\sigma_{n}^{2})$.

We place a Gaussian process (GP) prior over the unknown function,
$f(\bm{x})\sim\mathcal{GP}\left(\mu(\bm{x}),k(\bm{x},\bm{x}^{\prime}\right)$,
where $\mu(\bm{x})$ is the mean function and $k(\bm{x},\bm{x}^{\prime})$ is
the covariance function or kernel function. The mean function is usually
specified as zero and the data will be standardized when using GP. Therefore,
the estimation of GP is determined by the kernel function and the training
data. Let $X$ be a design matrix containing $n$ training inputs, $\bm{y}$ be
the aggregated targets, and $X_{*}$ be the testing points corresponding to the
function values $\bm{f}_{*}$. The predictive distribution of GP regression is
$\bm{f}_{*}|X,\bm{y},X_{*}\sim\mathcal{N}(\bar{\bm{f}}_{*}\text{cov}(\bm{f}_{*}))$,
where

$\displaystyle\bar{\bm{f}}_{*}$ $\displaystyle=K_{X_{*},X}K_{X}^{-1}\bm{y},$
(1) $\displaystyle\text{cov}(\bm{f}_{*})$
$\displaystyle=K_{X_{*},X_{*}}-K_{X_{*},X}K_{X}^{-1}K_{X,X_{*}}.$ (2)

Here we used the shorthand $K_{X}=[K_{X,X}+\sigma_{n}^{2}I]$.
$K_{X_{*},X},K_{X_{*},X_{*}},\text{and}K_{X,X}$ are the covariance matrices.
The element $(i,j)$ in each matrix is calculated in a component-wise fashion:
${K_{X,X^{\prime}}}{(i,j)}=k(x_{i},x^{\prime}_{j})$.

The covariance function usually contains some adjustable parameters, which are
called hyper-parameters. A widely adopted covariance function is the squared
exponential kernel,

$k_{SE}(\bm{x},\bm{x^{\prime}})=\sigma_{f}^{2}\exp{\left(-\frac{1}{2l^{2}}||\bm{x}-\bm{x^{\prime}}||^{2}\right)},$

which has hyper-parameters $\sigma_{f}^{2}$ and $l$. The optimal set of
parameters $\bm{\theta}^{*}=\\{\sigma_{f}^{2},l,\sigma_{n}^{2}\\}$ is found by
maximizing the log marginal likelihood $\log p(\bm{y}|X,\bm{\theta})$.

Based on the posterior covariance (Eq. (2)), the information metrics, such as
entropy, mutual information, and variance reduction, can be calculated.
Recently, most planning algorithms for robotic environmental monitoring belong
to informative planning which aims to find the trajectory maximizing one of
these information gain metrics.

##  III Problem Formulation

Maximizing information gain tends to encourage exploration, which is not
enough for catching up with the environmental dynamics. As the environment is
constantly changing, if the robot does not exploit the discovered features
(e.g. hotspots) in time, the information contained in the collected data will
“expire” because the environment may have changed significantly. Since
exploration and exploitation are two incompatible objectives, the problem
falls within the spectrum of multi-objective optimization.

In the general case, multi-objective informative planning requires solving the
following maximization problem:

$\displaystyle\bm{a}^{*}=\operatorname*{arg\,max}_{\bm{a}\in\mathcal{A}}\left\\{f_{1}(\bm{a}),f_{2}(\bm{a}),\dots,f_{D}(\bm{a})\right\\}\quad\text{
s.t. }C_{\bm{a}}\leq B$

where $\bm{a}$ is a sequence of actions, $\mathcal{A}$ is the space of
possible action sequences, $B$ is a budget (e.g. time, energy, memory, or
number of iterations), $C_{\bm{a}}$ is the cost of actions, and
$f_{d}(\bm{a}),d=1,\dots,D$ are the objective functions defining which types
of behaviors should the robot exhibit. At least one of the $f_{d}(\bm{a})$
should be the information gain objective.

##  IV METHODOLOGY

Figure 2: Block diagram of the proposed system.

Our objective is to monitor a time-varying phenomenon using an autonomous
vehicle. The robot needs to build a “distribution map” based on the data
collected for the environmental attribute of interest and plan the best
sampling path for the next step. A big challenge in spatiotemporal monitoring
is that, if the sampling path is planned based on the robot’s estimation of
the current environmental state, it becomes obsolete and thus sub-optimal when
the environment changes. To address this, we proposed to unify the online
informative planning (tree searching) mechanism with a forecast of the
environmental state following estimated environmental dynamics.

The system, shown in Fig. 2, consists of three main steps:

  * •

Modeling. As the robot travels along a path, it continuously collects data of
the environment attribute that we are interested in. These newly acquired
sample data are used to refine its estimation of the environment. Modeling is
not the focus of this work, thus we employ the commonly adopted Gaussian
process regression with the squared exponential kernel as described in Section
II-A.

  * •

Prediction. Combining the estimation of the environment with the available
ocean current or wind data, the robot predicts the possible moving path of the
target by performing a forecast in its simulator. Details are presented in
Section IV-A.

  * •

Planning. Given the Gaussian process and the prediction component, the robot
searches a path to follow using Pareto MCTS. This is discussed in Section
IV-D.

###  IV-A Environment Prediction

The prediction essentially refers to the simulation of possible environmental
states in the future based on meteorological data, such as ocean currents and
wind field, as well as the current estimation of the environmental state.
Incorporating the prediction component is achievable because the robot can
request the nowcast and forecast of the environmental conditions (e.g.,
oceanographic information from Operational Forecast System and wind data from
Global Forecast System). Also, some off-the-shelf simulators are available for
modeling the motion of the target in the ocean currents or wind fields, such
as the trajectory of the oil spill [3], propagation of wildfire [11], and
expansion of harmful algal blooms [23].

Formally, let $X$ be the matrix containing all possible sampling locations and
$\bm{v}_{i}^{t}(X),i=1,\dots,V$ be $V$ vector fields representing different
kinds of meteorological data at time $t$. If we denote the current time as
$t_{0}$ and the current estimation from GP (Eq. (1)) as
$\bm{f}_{\text{predict}}^{t_{0}}(X)$, then the environment simulator is a
function that takes $\bm{f}_{\text{predict}}^{t_{0}}(X)$ and
$\bm{v}_{i}^{t_{0}}(X),i=1,\dots,V$ as inputs, and returns the prediction of
the environment for the next time step $\bm{f}_{\text{predict}}^{t_{1}}(X)$.
Repeated use of the simulator allows us to obtain a series of predictions of
the environment
$\bm{f}_{\text{predict}}^{t_{0}}(X),\bm{f}_{\text{predict}}^{t_{1}}(X),\dots$.
These environment predictions will be part of the objective function for the
planning module, guiding the robot to where hotspot may appear. We will
discuss the objective functions in section IV-C.

###  IV-B Primitive Paths

Figure 3: Illustration of robot’s primitive paths. The arrows indicate the
sampling positions and orientations of the robot as it moves along the path.

In this work, the robot is equipped with a pre-computed path generator which
provides a set of primitive paths for every possible position. Each path
consists of a series of sampling points (see Fig. 3). Note that we will use
the term action and primitive path interchangeably, but they mean the same
thing. The primitive paths are generated as follows:

  * •

With some pre-defined radius around the robot, we obtain a given number of
points (called endpoints) that evenly divide the circle. (Points that directly
behind the robot might be discarded.)

  * •

Connect the robot’s current position with the endpoints using Dubin’s model
[9]. We assume the vehicle moves at constant speed in the plane and can
execute turns with a minimum turning radius. Given the position and heading of
the vehicle, Dubin’s model gives the feasible shortest path joining these two
oriented points.

  * •

Environmental sampling data is then collected along a Dubin’s path once it is
chosen.

###  IV-C Objective Functions

Although our proposed approach is applicable to the multi-objective problem,
we will focus on the bi-objective case for an illustrative purpose. Take the
hotspot monitoring task as an example, two behaviors of the robot are expected
in this task. First, the robot needs to collect as much information about the
environment as possible to determine the high-value areas (i.e., hotspots).
Second, it also needs to visit these important areas more frequently in order
to get better modeling performance. To encourage these two behaviors, we need
to define two objective functions (or reward functions).

Given the number of look-ahead steps $\tau$, the hotspot-seeking reward
function at $t_{0}$ is defined as the average of look-ahead prediction of the
sampling points along a path at time $t_{\tau}$. Suppose each primitive path
$a$ has $m$ sampling points $X_{a}$, then the reward of taking action $a$ at
time $t_{0}$ is given by

$f_{\text{hotspot}}^{t_{0}}(a)=\frac{1}{m}\sum\bm{f}_{\text{predict}}^{t_{\tau}}(X_{a}).$
(3)

This objective function is a bridge between environment prediction and
planning algorithms, encouraging the robot to visit where the target value
will be higher in the future.

For the information gain reward function, we use the average reduction in
variance (ARV) [5], because our method does not require the objective function
to have special properties, such as submodularity [22], and it is handy to
calculate. Given the prior covariance matrix $K_{\text{prior}}$ and the
posterior covariance matrix $K_{\text{post}}$ for $m$ sampling points of the
GP, the ARV objective function is defined as

$f_{ARV}(a)=\frac{1}{m}[\text{Trace}(K_{\text{prior}})-\text{Trace}(K_{\text{post}})].$
(4)

Reward function Eq. (3) encourages the robot to the high-value areas and Eq.
(4) guides the robot to its unknown or uncertain regions. The planning
algorithm aims to choose a sequence of actions $\bm{a}$ that maximize these
two reward functions,

$\displaystyle\bm{a}^{*}=\operatorname*{arg\,max}_{\bm{a}\in\mathcal{A}}\left\\{f_{\text{hotspot}}(\bm{a}),f_{ARV}(\bm{a})\right\\}\text{
s.t. }C_{\bm{a}}\leq B$ (5)

Since the GP and the environment prediction are updated as the robot collects
data, the reward functions are constantly changing. Hence, we need an online
planner that can adapt to the changing reward. In addition, the planner should
also consider the impact of the action on future results when choosing the
current best action (namely, nonmyopic planning [22]).

###  IV-D Pareto Monte Carlo Tree Search

(a) Symmetric Tree

(b) Asymmetric Tree

Figure 4: Illustration of a symmetric search tree and an asymmetric search
tree with the same number of tree nodes. The robot is represented as the
arrow. Red and white correspond to low and high reward, respectively. Pareto
Monte Carlo tree search estimates the expected reward of each child node
through Monte Carlo sampling, and puts more search efforts on the parts that
are more likely to yield the best solution. Figure 5: Illustration of Pareto
Monte Carlo Tree Search main steps.

One straightforward solution to Eq. (5) is to build a fixed-depth search tree
that evaluates all possible combinations of actions, then execute the first
best action, and then replan, i.e., receding-horizon tree search. However, the
computational complexity of this method is prohibitive. The Pareto Monte Carlo
tree search (Pareto MCTS) alleviates this burden by building an asymmetric
search tree which allocates its search time to the most important parts of the
tree, as shown in Fig. 4.

The tree is built incrementally through four steps: selection, expansion,
simulation, and backpropagation [6]. The main steps of Pareto MCTS are
depicted in Fig. 5. During each iteration, the algorithm descends through the
tree (red arrows in the selection step) and a new node is attached to the tree
(bold node in the expansion step), where each node represents a sampling
location and each edge denotes an action leading to a subsequent location.
Starting from the location represented by the newly added node, we take some
random actions and calculate the cumulative reward vector of taking these
actions using the reward functions. The cumulative reward vector will be
utilized to update the estimation of expected reward stored in each node along
the backpropagation path (red arrows).

The selection step chooses an expandable node in the tree. An expandable node
has at least one action that has not yet been tried during the search. To find
an expandable node in a best-first search manner, we start from the root node
and recursively apply a child node selection policy until an expandable node
is reached. We proposed an extension of the upper confidence bound (UCB)
policy [2, 13] for selecting the best next child in the multi-objective case
(described below). In the expansion step, we randomly choose an untried action
from the expandable node, calculate the robot’s location after executing this
primitive path, and build a new child node according to this action and
location.

During the simulation phase, the reward of the recently attached node is
estimated by taking random actions and evaluating the cumulative reward of
these actions until a pre-specified number of iterations is reached. The
estimated reward will be saved as a statistic in the node in order to bias
future tree searches. For our problem, the reward functions for evaluating the
actions are defined as Eq. (3) and (4). This also implies that the reward here
is not a scalar but a vector. The standard MCTS assumes that the reward is a
scalar, so special care should be taken during the selection, simulation, and
backpropagation steps (reward is not used in the expansion step). When
arithmetic operations on the reward are needed, we perform element-wise
operations. For example, in the simulation and backpropagation steps, we need
to calculate the cumulative reward. This is done by adding up the
corresponding elements in the reward vector. Since the modification of the
node selection policy is non-trivial, we will discuss it separately later.

The last step is the backpropagation where the operation of updating the
statistics of all nodes along the path from the newly attached node to the
root of the tree is performed. Specifically, we increase the number of visits
of all nodes along the path by one and add the reward evaluated in the
simulation step to the cumulative reward of each node.

We repeat these four steps until the maximum number of iterations is exceeded
or a given time limit is exhausted. Finally, the algorithm returns the best
action corresponding to the optimal child of the root node in terms of
expected reward. Pareto MCTS is an anytime algorithm which is able to return a
valid action at any moment during the searching process. A better action can
be found if it is allowed to run for a longer time.

For the critical node selection step in Pareto MCTS, one needs to balance
exploitation of the currently most promising child nodes and exploration of
alternatives which may turn out to be a superior choice at a later time. We
develop a selection policy that is able to handle the vector-form reward and
balances exploration and exploitation following the optimism in the face of
uncertainty principle [2].

Pareto optimality: In order to choose the best child node, we need to define a
metric first. In the case of scalar reward, the most promising node is simply
the one with the highest expected reward. However, given two reward vectors,
it is not obvious to tell which one is superior to the other. This is where
the Pareto optimal set comes into play.

In a multi-objective optimization problem, there does not typically exist a
solution that maximizes all objective functions simultaneously. Therefore, the
goal is to find the Pareto optimal set instead, or the so-called Pareto front.
The solutions in the Pareto front cannot be improved for any objective without
hurting other objectives. Formally, let $\bar{r}_{d}^{k}$ be the expected
reward corresponding to the $d$-th objective function of child node $k$. A
node $i$ is said to (Pareto) dominate another node $j$ if
$\bar{r}_{d}^{i}\geq\bar{r}_{d}^{j}$ for all $d=1,\dots,D$, and
$\bar{r}_{d}^{i}>\bar{r}_{d}^{j}$ for at least one $d\in\\{1,2,\dots,D\\}$. A
node is called Pareto optimal node if there does not exist another node that
dominates it. The set of Pareto optimal nodes is called Pareto front or Pareto
optimal set. It is worth mentioning that any node in this set is considered
equally best.

Node selection policy: In the node selection step, we evaluate a score vector
for each child and then choose the best child in terms of Pareto optimality.

First, we calculate the following upper confidence bound vector for each child
node [7]:

$\bm{U}(k)=\frac{\bm{r}^{k}}{n_{k}}+\sqrt{\frac{4\ln{n}+\ln{D}}{2n_{k}}},$ (6)

where $\bm{r}^{k}$ is the cumulative reward vector stored in node $k$, $n_{k}$
is the number of times child $k$ has been visited, $D$ is the number of
dimensions of the reward, and $n=\sum_{k=1}^{K}n_{k}$ is the number of times
the parent node has been visited. Eq. (6) is a combination of the empirical
expected reward and a confidence interval. The first term prefers the recently
discovered most promising child whereas the second term emphasizes exploration
of alternatives. When the number of visits of the $k$-th node remains the same
while that of other nodes increases, the second term will become larger,
making node $k$ more likely to be selected.

Second, we build a Pareto optimal set based on the resulting
$\bm{U}(k),k=1,\dots,K$. Alg. 1 outlines the main steps for computing the
Pareto optimal set for a bi-objective case. See [14] for more general multi-
objective cases.

Algorithm 1 Computing Bi-Objective Pareto Front.

1:procedure $\mathtt{compute\\_pareto\\_front}$($\mathbf{U}$)

2: $\triangleright$ $\mathbf{U}$ is a $K\times{2}$ matrix consists of upper
confidence vectors

3: $\mathbf{U}_{\text{sort}}\leftarrow$ sort $\mathbf{U}$ w.r.t. the first
column in descending order

4:
$\mathbf{u}^{\star}\leftarrow\mathbf{U}_{\text{sorted}}[0,:]$$\triangleright$
get the first row vector

5: $\mathbb{P}\leftarrow\\{\text{node corresponds to
}\mathbf{u}^{\star}\\}$$\triangleright$ add to Pareto optimal set

6: for $i=1,\dots,K-1$ do

7: if $\mathbf{u}^{\star}[1]\leq\mathbf{U}_{\text{sort}}[i,1]$
then$\triangleright$ Compare the second dim.

8: $\mathbf{u}^{\star}\leftarrow\mathbf{U}_{\text{sort}}[i,1]$

9: $\mathbb{P}\leftarrow\mathbb{P}\cup\\{\text{node corresponds to
}\mathbf{u}^{\star}\\}$

10: return $\mathbb{P}$

Finally, the best child is chosen from the Pareto optimal set uniformly at
random. The reason for choosing the best child randomly is that the Pareto
optimal nodes are considered equally optimal if no preference information is
specified.

##  V EVALUATION RESULTS

In order to verify the advantages of incorporating multiple objectives and
environment prediction into planning, we compared the UCT [13], Pareto MCTS,
and the proposed look-ahead based multi-objective tree search in the task of
monitoring a spatially and temporally dynamic hotspot.

###  V-A Experiment Setup

The task of the robot is to monitor a hotspot that passively moved by its
fluid medium. One can think about this as monitoring the oil spill or harmful
algal bloom in the ocean. The robot is equipped with $11$ primitive paths,
each of which is $2.5$ m long. It collects data in a two-dimensional
$150\text{ m }\times 150\text{ m }$ environment with range ($x_{\min}$: $-75$,
$x_{\max}$: $75$, $y_{\min}$: $-75$, $y_{\max}$: $75$). The hotspot is
specified by a Gaussian distribution with initial mean
$\bm{m}=[0,-50]^{\top}$. The robot starts from $[-20,-40]$ with no prior
knowledge of the environment, so it needs to explore the environment first to
find the hotspot.

The input of GP only includes the spatial dimension and the spatiotemporal
correlation is not modeled in this work, thus we assume that the environment
is static in a short time period. The GP will be reset and the robot is
allowed to take $5$ actions at each time step. For the $s$-th ($s=0,\dots,4$)
action $a_{s}$ at time $t_{0}$, the hotspot-seeking reward function is defined
as the average s-step look-ahead prediction

$f_{\text{hotspot}}^{t_{0}}(a_{s})=\frac{1}{m}\sum
f_{\text{predict}}^{t_{s}}(X_{*}).$ (7)

We also assume that a forecast of ocean currents with near future is available
although noisy [16]. At each location $(x,y)$ and time $t$, the predicted
ocean current vector is a noisy version of the true ocean current:

$\bm{v}^{t}_{\text{predict}}=\bm{v}^{t}_{\text{true}}+\bm{\epsilon}=\begin{pmatrix}5\cos{\left(2\pi\frac{t\text{
mod }60}{60}\right)}+\epsilon_{x}\\\ 5\sin{\left(2\pi\frac{t\text{ mod
}60}{60}\right)}+\epsilon_{y}\end{pmatrix},$ (8)

where $\epsilon_{x},\epsilon_{y}$ are independent Gaussian noises drawn from
$\mathcal{N}(0,0.2)$.

The Gaussian distribution is linearly moved (disturbed) by the ocean current
vector:

$\bm{m}_{t+1}=\bm{m}_{t}+\bm{v}^{t}_{\text{true}}$ (9)

In the environment prediction component, we first find out the locations where
the GP estimations are greater than the median of all estimations. Then, these
locations are moved by the predicted ocean current vector
$\bm{v}^{t}_{\text{predict}}$ as in Eq. (9).

In the three compared planning algorithms, UCT only uses the ARV objective
function and selects the child node maximizing the following upper confidence
bound in the selection step:

$U(k)=\frac{r^{k}}{n_{k}}+\sqrt{\frac{2\ln{n}}{n_{k}}},$ (10)

where $r^{k}$ is the scalar cumulative reward stored in node $k$ and other
variables have the same meaning as the variables in Eq. (6). Note that when
$D=1$, Eq. (6) degenerates into Eq. (10). Our proposed Eq. (6) is a
multidimensional extension of Eq. (10). Pareto MCTS differs from the proposed
method in that it only uses the current estimation of the environment in its
hotspot-seeking reward function, and does not involve any environmental model
prediction. All three versions of MCTS use the same parameters: $200$ maximum
number of tree search iterations and $2$ random actions in the simulation
step. The hyper-parameters of GP are updated by maximizing the log marginal
likelihood in an online manner.

###  V-B Results

Figure 6: Snapshots of the robot’s behaviors in the hotspot monitoring task
using UCT (first row), Pareto MCTS (second row) and model-predictive tree
search (third row) at time $15,30,45$, and $59$. The arrows indicate the
positions and orientations of the robot. Blue lines are the historical
trajectories. For each time step, the left image is the ground truth of the
environment, and the right one is robot’s current estimation of the
environment.

Fig. 6 presents the snapshots of the robot’s behaviors when using UCT, Pareto
MCTS, and model-predictive tree search. At the beginning of the mission, all
the three algorithms guide the robot to explore the environment, looking for
the hotspot. This can be seen from the historical trajectories of the robot,
because the trajectory will appear to be stretched when the robot is
exploring, and it will be circling during exploitation. When the hotspot is
found, the three algorithms show different characteristics.

UCT with only information gain objective function ignores the hotspot and
continues exploring (See the first row of Fig. 6), because less information
will be gained from a sampling point once it has been observed. Without
exploitation at a finer level, we cannot identify the hotspot feature in the
mean of GP. This explains why it is insufficient to use only the information
gain objective in spatiotemporal monitoring.

Pareto MCTS tends to visit the high-value area more frequently once it has
been discovered (second row of Fig. 6). However, the robot will easily lose
the target and needs to search for it again (second row of time $30$ and
$45$). As a comparison, the proposed method not only better tracks the
hotspot, but also guides the robot to the area that the hotspot may exist in
the next time (third row of time $30$).

The left pictures of the fourth column show the complete historical paths. The
path planned by the proposed method is the most consistent with the ground-
truth underlying hotspot moving path (a circle). It can be seen that the robot
moderately wandered around the hotspot to collect more valuable samples.

These behavioral characteristics directly affect the final modeling
performance. From the right image of each column in Fig. 6, we can see that
UCT fails to model the hotspot, Pareto MCTS overestimates the range of the
hotspot, and the model-predictive tree search provides a good estimate for the
underlying hotspot. We also quantitatively compare the root mean square error
(RMSE) in Fig. 7 and distance to the center of the hotspot in Fig. 8. When the
robot does not collect any sample in the hotspot, the GP predicts the value of
all locations as $0$, resulting in $0.000255$ RMSE. The robot using UCT did
not find the hotspot most of the time, but occasionally encountered hotspot
during the exploration process (Fig. 8). As can be seen from Fig. 7 and Fig.
8, Pareto MCTS lost the hotspot during the time period $[30,40]$, but most of
time it was able to collect data around hotspot. Since it is often unable to
catch up with the dynamics of the environment, the collected samples do not
represent the underlying hotspot well, leading to unsatisfactory RMSE. The
RMSE of the model-predictive tree search drops drastically and fluctuates at
low values, which reveals that the proposed method can quickly find the
hotspot, then track it and exploit it.

Figure 7: The root mean square error (RMSE) of the three version of MCTS.
Figure 8: Distance from the robot’s position to the center of the hotspot.

##  VI CONCLUSION and FUTURE WORK

This paper presents an informative planning framework for robotic monitoring
of spatiotemporal environments. According to the nowcast and forecast of the
environmental model, the robot is able to infer the future changes of the
target and plan its sampling actions accordingly. Within this framework, we
propose a sampling-based multi-objective tree search planning algorithm that
allows the robot to handle multiple potentially conflicting objectives, such
as exploration and exploitation. The results of a hotspot evolution monitoring
experiment show that our proposed method enables the robot to adapt to the
time-varying environment and balance exploration and exploitation behaviors,
achieving smaller modeling error and closer distance to the hotspot center.

## References

  * Arora et al. [2017] A. Arora, P. M. Furlong, R. Fitch, S. Sukkarieh, and T. Fong.  Multi-Modal Active Perception for Information Gathering in Science Missions.  _Autonomous Robots (AuRo)_ , pages 1–27, 2017.
  * Auer et al. [2002] P. Auer, N. Cesa-Bianchi, and P. Fischer.  Finite-time analysis of the multiarmed bandit problem.  _Machine learning_ , 47(2-3):235–256, 2002.
  * Beegle-Krause [2001] J. Beegle-Krause.  General noaa oil modeling environment (gnome): A new spill trajectory model.  In _International Oil Spill Conference_ , pages 865–871, 2001.
  * Best et al. [2018] G. Best, O. M. Cliff, T. Patten, R. R. Mettu, and R. Fitch.  Dec-MCTS: Decentralized Planning for Multi-Robot Active Perception.  _The International Journal of Robotics Research (IJRR)_ , 2018.
  * Binney et al. [2013] J. Binney, A. Krause, and G. S. Sukhatme.  Optimizing waypoints for monitoring spatiotemporal phenomena.  _The International Journal of Robotics Research (IJRR)_ , 32(8):873–888, 2013.
  * Browne et al. [2012] C. B. Browne, E. Powley, D. Whitehouse, S. M. Lucas, P. I. Cowling, P. Rohlfshagen, S. Tavener, D. Perez, S. Samothrakis, and S. Colton.  A survey of monte carlo tree search methods.  _IEEE Transactions on Computational Intelligence and AI in games_ , 4(1):1–43, 2012.
  * Chen and Liu [2019] W. Chen and L. Liu.  Pareto monte carlo tree search for multi-objective informative planning.  In _Proceedings of Robotics: Science and Systems (RSS)_ , 2019.
  * Corah and Michael [2017] M. Corah and N. Michael.  Efficient Online Multi-robot Exploration via Distributed Sequential Greedy Assignment.  In _Robotics: Science and Systems (RSS)_ , 2017.
  * Dubins [1957] L. E. Dubins.  On curves of minimal length with a constraint on average curvature, and with prescribed initial and terminal positions and tangents.  _American Journal of mathematics_ , 79(3):497–516, 1957.
  * Dunbabin and Marques [2012] M. Dunbabin and L. Marques.  Robots for Environmental Monitoring: Significant Advancements and Applications.  _IEEE Robotics & Automation Magazine(RAM)_, 19(1):24–39, 2012.
  * Finney [1998] M. A. Finney.  Farsite: Fire area simulator-model development and evaluation.  _US Department of Agriculture_ , 4, 1998.
  * Hollinger and Sukhatme [2014] G. A. Hollinger and G. S. Sukhatme.  Sampling-based robotic information gathering algorithms.  _The International Journal of Robotics Research (IJRR)_ , 33(9):1271–1287, 2014.
  * Kocsis and Szepesvári [2006] L. Kocsis and C. Szepesvári.  Bandit based monte-carlo planning.  In _European conference on machine learning (ECML)_ , pages 282–293, 2006.
  * Kung et al. [1975] H.-T. Kung, F. Luccio, and F. P. Preparata.  On finding the maxima of a set of vectors.  _Journal of the ACM (JACM)_ , 22(4):469–476, 1975\.
  * Lim et al. [2016] Z. W. Lim, D. Hsu, and W. S. Lee.  Adaptive informative path planning in metric spaces.  _The International Journal of Robotics Research (IJRR)_ , 35(5):585–598, 2016.
  * Liu and Sukhatme [2018] L. Liu and G. S. Sukhatme.  A Solution to Time-Varying Markov Decision Processes.  _IEEE Robotics and Automation Letters (RA-L)_ , 3(3):1631–1638, 2018.
  * Ma et al. [2018] K.-C. Ma, L. Liu, H. K. Heidarsson, and G. S. Sukhatme.  Data-driven learning and planning for environmental sampling.  _Journal of Field Robotics (JFR)_ , 35(5):643–661, 2018.
  * Marchant et al. [2014] R. Marchant, F. Ramos, S. Sanner, et al.  Sequential Bayesian optimisation for Spatial-Temporal Monitoring.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence (UAI)_ , pages 553–562, 2014.
  * Meliou et al. [2007] A. Meliou, A. Krause, C. Guestrin, and J. M. Hellerstein.  Nonmyopic informative path planning in spatio-temporal models.  In _The Association for the Advancement of Artificial Intelligence (AAAI)_ , volume 10, pages 16–7, 2007.
  * Patten et al. [2018] T. Patten, W. Martens, and R. Fitch.  Monte Carlo Planning for Active Object Classification.  _Autonomous Robots (AuRo)_ , pages 391–421, 2018.
  * Silver et al. [2016] D. Silver, A. Huang, C. J. Maddison, A. Guez, L. Sifre, G. Van Den Driessche, J. Schrittwieser, I. Antonoglou, V. Panneershelvam, M. Lanctot, et al.  Mastering the Game of Go with Deep Neural Networks and Tree Search.  _Nature_ , 529(7587):484, 2016.
  * Singh et al. [2009] A. Singh, A. Krause, and W. J. Kaiser.  Nonmyopic adaptive informative path planning for multiple robots.  In _Twenty-First International Joint Conference on Artificial Intelligence (IJCAI)_ , 2009.
  * Stumpf et al. [2009] R. P. Stumpf, M. C. Tomlinson, J. A. Calkins, B. Kirkpatrick, K. Fisher, K. Nierenberg, R. Currier, and T. T. Wynne.  Skill assessment for an operational algal bloom forecast system.  _Journal of Marine Systems_ , 76(1-2):151–161, 2009.