# Tightly Integrated Motion Classification and State Estimation in Foot-
Mounted Navigation Systems

Isaac Skog  Dept. of Electrical Engineering
Uppsala University
Uppsala, Sweden
<EMAIL_ADDRESS>Gustaf Hendeby  Dept. of Electrical Engineering
Linköping University
Linköping, Sweden
<EMAIL_ADDRESS>Manon Kok  Delft Center for Systems and Control
Delft University of Technology
Delft, The Netherlands
<EMAIL_ADDRESS>
###### Abstract

A framework for tightly integrated motion mode classification and state
estimation in motion-constrained inertial navigation systems is presented. The
framework uses a jump Markov model to describe the navigation system’s motion
mode and navigation state dynamics with a single model. A bank of Kalman
filters is then used for joint inference of the navigation state and the
motion mode. A method for learning unknown parameters in the jump Markov
model, such as the motion mode transition probabilities, is also presented.
The application of the proposed framework is illustrated via two examples. The
first example is a foot-mounted navigation system that adapts its behavior to
different gait speeds. The second example is a foot-mounted navigation system
that detects when the user walks on flat ground and locks the vertical
position estimate accordingly. Both examples show that the proposed framework
provides significantly better position accuracy than a standard zero-velocity
aided inertial navigation system. More importantly, the examples show that the
proposed framework provides a theoretically well-grounded approach for
developing new motion-constrained inertial navigation systems that can learn
different motion patterns.

###### Index Terms:

Inertial navigation, Zero-velocity detection, Constant height detection,
Filter bank, Motion-constraints.

##  I Introduction

Current state-of-the-art technology for zero-velocity aided inertial
navigation systems is based upon a strategy of loose integration between the
zero-velocity detector and the inertial navigation filter [1]. That is, the
zero-velocity detector and the inertial navigation filter are treated as
separate functions, with a one-directional flow of information from the former
to the latter. From an information theoretical perspective this is suboptimal
because: (a) the estimated navigation states carry information about the
system’s motion mode that is not used in the zero-velocity detector, and (b)
the test statistic used in the zero-velocity detector is quantized before it
is used in the navigation filter. Thus, information that could be used both to
improve the detection of zero-velocity events and to control how the zero-
velocity updates are performed, is lost. The same argumentation also applies
to other motion constraints commonly applied to foot-mounted inertial
navigation systems and where an external motion classifier is used to
determine when to apply these constraints. A few examples of commonly used
motion constraints are that the system keeps a constant height [2], constant
heading [3], or constant speed [4].

The fact that the estimated navigation state carries information that can be
used to improve the zero-velocity detection process is utilized in [5], where
the velocity estimates from the inertial navigation filter are used as a prior
for a Bayesian zero-velocity detector. And in [6], the test statistics from
the zero-velocity detector are used to control the magnitude of the
measurement covariance in the inertial navigation filter. Thereby,
quantization of the zero-velocity test statistics is avoided and the zero-
velocity update process can adapt to different gait conditions. Still, both
[5] and [6] employ a loose integration strategy where the zero-velocity
detector and the inertial navigation filter are treated as separate functions.

This paper instead proposes a framework for tightly integrated motion mode
classification and navigation state estimation, of which tightly integrated
zero-velocity aided inertial navigation is a special case. The core of the
framework is a jump Markov model that includes both the navigation state and
the motion mode. That is, a single model is used to describe both the
kinematics of the system under various motion modes and the probability of
transitioning between the motion modes. A bank of Kalman filters is then used
to jointly estimate the navigation state and the motion mode. A method to
automatically learn unknown parameters in the jump Markov model is also
presented. The proposed framework is evaluated on two data sets consisting of
various gait conditions and motion modes.

Reproducible research: The data and code used to produce the presented results
can be downloaded at:
https://gitlab.liu.se/open-shoe/filterbanks.

##  II Signal Model

Let the system state $x_{k}$ and input vector $u_{k}$ at time instant $k$ be
defined as

$x_{k}=\begin{bmatrix}r_{k}\\\ v_{k}\\\ q_{k}\\\
\xi_{k}\end{bmatrix}\quad\text{and}\quad u_{k}=\begin{bmatrix}s_{k}\\\
\omega_{k}\end{bmatrix},$ (1)

respectively. Here $r_{k}$, $v_{k}$, and $q_{k}$ denote the position,
velocity, and attitude quaternion, respectively. Further, $\xi_{k}$ denotes
potential auxiliary states, such as sensor biases, needed to describe the
behavior of the system. Moreover, $s_{k}$ and $\omega_{k}$ denote the specific
force and angular velocity, respectively.

Next, define the discrete state

$\delta_{k}\in\\{1,\ldots,L\\}$ (2)

that indicates the current motion mode of the navigation system, e.g., the
system is stationary, keeping a constant height, or moving at constant
velocity. A jump Markov nonlinear model that can be used to the describe the
dynamics of the inertial navigation system is then given by [7]

$\displaystyle x_{k+1}=f(x_{k},\delta_{k},u_{k},\eta_{k})$ (3a) $\displaystyle
y_{k}=h(x_{k},\delta_{k},u_{k})+e_{k}$ (3b)
$\displaystyle\eta_{k}\overset{\text{\tiny
i.i.d.}}{\sim}\mathcal{N}\bigl{(}n_{k};0,Q_{k}\bigr{)}$ (3c) $\displaystyle
e_{k}\overset{\text{\tiny
i.i.d.}}{\sim}\mathcal{N}\bigl{(}e_{k};0,R(\delta_{k})\bigr{)}$ (3d)
$\displaystyle p(\delta_{k+1}|\delta_{k})=\Pi_{\delta_{k+1},\delta_{k}}.$ (3e)

Here $y_{k}\triangleq 0$ is used as a pseudo-observation to impose a set of
motion mode dependent “stochastic constraints”, such as zero-velocity
constraints, on the navigation state $x_{k}$. The exact form of these
constraints, defined by $h(\cdot)$, depends on the assumed motion modes and
will be discussed later. The “hardness” of the imposed constraints is
controlled via the motion mode dependent covariance $R(\delta_{k})$. Further,
$\mathcal{N}\bigl{(}\,\cdot\,;\mu,\Sigma\bigr{)}$ denotes a multivariate
normal distribution parameterized by the mean vector $\mu$ and covariance
matrix $\Sigma$. The function $f(\cdot)$ describes the system dynamics and is
given by the inertial navigation equations and the dynamics of the auxiliary
states $\xi_{k}$. Moreover, $Q_{k}$ denotes the covariance of the process
noise, and $\Pi_{i,j}$ denotes the $i$:th row and $j$:th column entry of the
mode transition probability matrix. Finally, $p(a|b)$ denotes the probability
density (mass) function of $a$ given $b$. Note that it is straightforward to
extend the model (3) to also include real observations from various sensors,
but it is out of the scope of this paper, as the focus is on motion-
constrained inertial navigation.

##  III State Estimation Using a Filter Bank

A variety of inference techniques for jump Markov models exist. Here a
commonly used filter banks solution to the inference problem will be
recapitulated; for details the reader is referred to [8] and [9]. Thereafter,
necessary approximations needed to adapt the filter bank solution to the
considered navigation problem will be presented.

###  III-A Linear model with normal distributed noise

The goal of the inference process is to estimate the a posteriori distribution
$p(x_{k}|y_{1:k})$ of the state $x_{k}$ given all the observations up until
time $k$, denoted as $y_{1:k}$. If the motion mode sequence $\delta_{1:k}$ is
known, the state-space model in the Markov jump system is linear, and the
process and measurement noise are normally distributed, then the a posteriori
distribution can be estimated with the Kalman filter. That is,

$p(x_{k}|\delta_{1:k},y_{1:k})=\mathcal{N}\bigl{(}x_{k};\hat{x}^{\delta_{1:k}}_{k|k},P^{\delta_{1:k}}_{k|k}\bigr{)},$
(4)

where $\hat{x}_{k|k}^{\delta_{1:k}}$ and $P_{k|k}^{\delta_{1:k}}$ denote the
Kalman filter state estimate and state covariance given the mode sequence
$\delta_{1:k}$, respectively. In reality, the motion mode sequence is unknown
and must be estimated from the measurements. Let $p(\delta_{1:k}|y_{1:k})$
denote the a posteriori distribution of the motion mode sequences
$\delta_{1:k}$ given the measurements $y_{1:k}$. The sought-after a posteriori
distribution of the state $x_{k}$ can then be found as

$\begin{split}p(x_{k}|y_{1:k})&=\sum_{i=1}^{L^{k}}p(x_{k}|\delta^{i}_{1:k},y_{1:k})p(\delta^{i}_{1:k}|y_{1:k})\\\
&=\sum_{i=1}^{L^{k}}w^{i}_{k}\mathcal{N}\bigl{(}x_{k};\hat{x}^{\delta^{i}_{1:k}}_{k|k},P^{\delta^{i}_{1:k}}_{k|k}\bigr{)},\end{split}$
(5)

where $w^{i}_{k}\triangleq p(\delta^{i}_{1:k}|y_{1:k})$ is probability of the
$i$:th motion mode sequence. Hence, the posteriori distribution of the state
$x_{k}$ is given by a mixture of $L^{k}$ weighted normal distributions where
the expected value and covariance of each mixture component are calculated via
a Kalman filter. That is, the posteriori distribution can be calculated via a
bank of Kalman filters where the size of the filter bank grows according to a
tree structure. The weights $w^{i}_{k}$ in the mixture, i.e., the probability
for each branch in the tree, can be recursively calculated as

$w^{i}_{k}\propto
p(y_{k}|y_{1:k-1},\delta^{i}_{1:k})p(\delta^{i}_{k}|\delta^{i}_{k-1})w^{i}_{k-1},$
(6a) where $\sum_{i=1}^{L^{k}}w^{i}_{k}=1,$ (6b) and
$p(y_{k}|y_{1:k-1},\delta^{i}_{1:k})=\mathcal{N}\bigl{(}y_{k};\hat{y}_{k|k-1}^{\delta^{i}_{1:k}},S_{k}^{\delta^{i}_{1:k}}\bigr{)}.$
(6c)

Here $\hat{y}_{k|k-1}^{\delta_{1:k}}$ and $S_{k}^{\delta_{1:k}}$ denote the
Kalman filter measurement prediction and innovation covariance, respectively.

From the posteriori distribution, the minimum variance estimate of the state
$x_{k}$ can be calculated as

$\hat{x}^{\text{\scriptsize
mv}}_{k}=\sum_{i=1}^{L^{k}}w^{i}_{k}\hat{x}_{k|k}^{\delta^{i}_{1:k}}$ (7a) and
the associated conditional covariance matrix as $P^{\text{\scriptsize
mv}}_{k}=\sum_{i=1}^{L^{k}}w^{i}_{k}\bigl{(}P_{k|k}^{\delta^{i}_{1:k}}+(\hat{x}_{k|k}^{\delta^{i}_{1:k}}-\hat{x}^{\text{\scriptsize
mv}}_{k})(\hat{x}_{k|k}^{\delta^{i}_{1:k}})-\hat{x}^{\text{\scriptsize
mv}}_{k})^{\top}\bigr{)}.$ (7b)

###  III-B Approximative solution

The outlined filter bank solution to the inference problem can in general not
be used without modifications due to the growing number of Kalman filters
needed. Therefore, many strategies have been developed for pruning and merging
branches in the growing filter bank tree so that a fixed complexity is
achieved [10].

Beyond the challenges with the exponentially increasing complexity, the
outlined general solution cannot still be applied to the Markov jump system
model in (3). The system is nonlinear and the attitude states belong to a
manifold of Euclidean space. This implies that $p(x_{k}|\delta_{1:k},y_{1:k})$
in (4) is generally not a normal distribution, neither can the minimum
variance estimate of the attitude $q_{k}$ be calculated as in (7a). A common
way to handle that the system is nonlinear and the attitude not belonging to
Euclidian space in the filtering process is to use an error state Kalman
filter [11]. The a posteriori distribution $p(x_{k}|\delta_{1:k},y_{1:k})$ is
then, at every time instant, approximated as normal distributed with the mean
and covariance given by the error state Kalman filter [12].

One way to calculate a point estimate of the attitude in terms of the Euler
angles $\phi$ is via

$\hat{\phi}^{\text{\scriptsize
mv}}_{k}=\operatorname*{arg\,min}_{\phi\in\Omega_{\phi}}\sum_{i}^{L^{k}}w^{i}_{k}\|C(\hat{q}^{\delta^{i}_{1:k}}_{k|k})-C(\phi)\|^{2}_{F},$
(8)

where $C(q)$ denotes the rotation matrix that transforms a vector from the
body frame to the navigation frame [13]; to simplify the notation, a slight
misusage of notation is admitted, and the rotation matrix is here
interchangeably parameterized by the attitude quaternion $q$ and the
corresponding Euler angles $\phi$. Further, the set
$\Omega_{\phi}=[0,2\pi)\times[-\pi/2,\pi/2)\times[0,2\pi)$. The covariance of
$\hat{\phi}^{\text{\scriptsize mv}}_{k}$ can be approximately calculated as

$\text{Cov}(\hat{\phi}^{\text{\scriptsize
mv}}_{k})\approx\sum_{i}^{L^{k}}w^{i}_{k}\bigl{(}\Sigma^{\delta^{i}_{1:k}}_{k|k}+\Delta\phi^{\delta^{i}_{1:k}}_{k}(\Delta\phi^{\delta^{i}_{1:k}}_{k})^{\top}\bigr{)}.$
(9)

Here $\Sigma^{\delta^{i}_{1:k}}_{k|k}$ denotes the covariance of the attitude
(in terms of Euler angles) calculated by the error state Kalman filter for the
branch corresponding to the motion sequence $\delta^{i}_{1:k}$. Further, the
difference between the attitude estimate calculated by the error state Kalman
filter and the minimum variance attitude estimate in (8) is given by

$\Delta\phi^{\delta^{i}_{1:k}}_{k|k}=\\{\phi\in\Omega_{\phi}~{}\text{s.t.}~{}~{}C(\hat{\phi}^{\text{\scriptsize
mv}}_{k})=C(\phi)C(\hat{q}^{\delta^{i}_{1:k}}_{k|k})\\}.$ (10)

Note that the optimization problem in (8) can be efficiently solved using a
polar decomposition [13].

##  IV Learning of Unknown Model Parameters

The signal model in (3) may have unknown parameters, such as the transition
probability matrix $\Pi$, whose values can be hard to specify accurately.
Instead of resorting to cumbersome hand-tuning of these parameters, they may
be learned from data. A variety of methods to learn model parameters in jump
Markov models has been suggested [14, 15]. Here a maximum likelihood method,
similar to that presented in [12], for learning the parameters from the data
will be outlined.

Let $\theta$ denote the unknown parameters in the model. The maximum
likelihood estimate of these parameters is then given by

$\hat{\theta}=\operatorname*{arg\,max}_{\theta}p(y_{1:k};\theta),$ (11a) where
$p(y_{1:k};\theta)=\prod_{n=2}^{k}p(y_{n}|y_{1:n-1})p(y_{1}),$ (11b) and
$p(y_{n}|y_{1:n-1})=\sum_{i=1}^{L^{n}}\mathcal{N}\bigl{(}y_{n};\hat{y}_{n|n-1}^{\delta^{i}_{1:n}},S_{n}^{\delta^{i}_{1:n}}\bigr{)}p(\delta^{i}_{n}|\delta^{i}_{n-1})w^{i}_{n-1}.$
(11c)

Here $p(y_{1})$ denotes the a priori probability of the observation $y_{1}$.
Note that the number of modes in the full mixture distribution, which
subsequently must be summed to marginalize away the mode dependence, increases
exponentially in $n$. However, in practice, the actual number of modes to
consider is reduced using pruning or merging in the filter bank. Hence, the
number of modes to consider when evaluating (11c) can be kept tractable.

##  V Application Examples

Next, the proposed inference framework is used to realize two foot-mounted
inertial navigation systems that incorporate different motion modes. The first
system adaptively selects the covariance matrix used in the zero-velocity
update, and the second system automatically detects if the system returns to
the same height after a step. For both systems, the motion mode state
transition matrix is learned from training data using the method outlined in
Sec. IV. The performances of both systems are compared to the OpenShoe system
presented in [16].

Unconstrainedmotion $\delta=1$Almost
stationary$\delta=2$Stationary$\delta=3$$\Pi_{1,1}$$\Pi_{2,2}$$\Pi_{3,3}$$\Pi_{2,1}$$\Pi_{3,2}$$\Pi_{2,3}$$\Pi_{1,2}$

(a) Varying gait speed example.

Unconstrainedmotion $\delta=1$Stationary &new height$\delta=2$Stationary &same
height$\delta=3$$\Pi_{1,1}$$\Pi_{3,1}$$\Pi_{3,3}$$\Pi_{2,1}$$1$$\Pi_{1,3}$

(b) Return to same height example.

Figure 1: Mode transition diagrams for the application examples.

###  V-A Example: Varying gait speed

The challenge of designing and tuning a foot-mounted zero-velocity aided
inertial navigation system so that it works well for multiple gait speeds is
well-known [1]. To use the proposed framework to design a system that
adaptively selects the detection threshold, as well as the covariance matrix
used in the zero-velocity updates, we define the following three motion modes:
the unconstrained motion mode ($\delta_{k}=1$), the almost stationary motion
mode ($\delta_{k}=2$), and the stationary motion mode ($\delta_{k}=3$). A
motion mode transition diagram illustrating how the system is assumed to
transition between these motion modes is shown in Fig. 1(a). The associated
motion mode transition probabilities $\Pi_{i,j}$ are also shown.

Next, introduce the following system dynamics

$f(x_{k},\delta_{k},u_{k},\eta_{k})=\begin{bmatrix}r_{k}+\Delta tv_{k}\\\
v_{k}+\Delta t\bigl{(}C(q_{k})(s_{k}+\eta^{s}_{k})+g\bigr{)}\\\
\Omega\bigl{(}\Delta t\,(\omega_{k}+\eta^{\omega}_{k})\bigr{)}q_{k}\\\
\end{bmatrix}.$ (12a) Here $\Delta t$ and $g$ denote the sampling period and
gravity vector, respectively. Further, $\Omega(\cdot)$ denotes the quaternion
update matrix (see [17] for details). Moreover, $\eta^{s}_{k}$ and
$\eta^{\omega}_{k}$ denote the process noise associated with the
accelerometers and gyroscopes, respectively.

Finally, define the stochastic motion constraints imposed during the different
motion modes as follows

$h(x_{k},u_{k},\delta_{k})=\begin{cases}\hfil~{}~{}~{}0\hfil,&\delta_{k}=1\\\
h_{0}(x_{k},u_{k}),&\delta_{k}=\\{2,3\\}\end{cases},$ (12b)

where

$h_{0}(x_{k},u_{k})=\begin{bmatrix}v_{k}\\\ \omega_{k}\\\
C(q_{k})s_{k}+g\end{bmatrix}.$ (12c)

The “hardness” of the constraints is controlled by the covariance matrix

$R(\delta_{k})=\begin{cases}\sigma^{2}_{nc}I,&\delta_{k}=1\\\
R_{0}(\delta_{k}),&\delta_{k}=\\{2,3\\}\end{cases},$ (12d)

where

$R_{0}(\delta_{k})=\sigma^{2}_{v}({\delta_{k}})I\oplus\sigma^{2}_{\omega}(\delta_{k})I\oplus\sigma^{2}_{s}(\delta_{k})I.$
(12e)

Here $\oplus$ denotes the direct sum matrix operator and $I$ denotes an
identity matrix of appropriate size. Further, $\sigma^{2}_{v}({\delta_{k}})$,
$\sigma^{2}_{\omega}(\delta_{k})$, and $\sigma^{2}_{s}(\delta_{k})$ denote the
mode-dependent variance associated with the stochastic constraints on the
velocity, angular velocity, and acceleration, respectively. Moreover,
$\sigma^{2}_{nc}$ is a design parameter that controls the measurement
likelihood associated with no constraints.

A filter bank with the jump Markov model defined by (12) was designed and two
data sets were recorded with a sensor array consisting of 32 InvenSense
MPU9150 IMUs while a person walked and ran back and forth along a straight
line. From the recorded data, two data sets with 50 measurement sequences each
were created by repeatedly drawing four random IMUs and averaging their
measurements. These measurement sequences were then processed using the
OpenShoe system algorithm with the zero-velocity detector threshold tuned to
generate a detection at least once per gait cycle. The measurements were also
processed using the designed filter bank. Pruning was used to keep the tree
size to a maximum of nine leaves, and a pruning strategy where the most
probable leaves were retained was used.

The measurement sequences in the first data set were used to learn the mode
transition matrix $\Pi$. The initial value and end result of the learning
process were

$\widehat{\Pi}^{\text{\tiny init}}=\begin{bmatrix}1/2&1/3&0\\\ 1/2&1/3&1/2\\\
0&1/3&1/2\end{bmatrix}~{}\text{and}~{}~{}\widehat{\Pi}^{\text{\tiny
end}}=\begin{bmatrix}0.993&0.073&0\\\ 0.007&0.893&0.005\\\ 0&0.034&0.995\\\
\end{bmatrix},$

respectively, and the optimization converged in less than 10 iterations. The
matrix $\widehat{\Pi}^{\text{\tiny end}}$ corresponds to a mode transition
system with low-pass characteristics. That is, the probability of staying in a
mode is much higher than transitioning to another mode. On average, the
percentage of time spent in the motion modes one, two, and three is 57%, 5%,
and 38%, respectively.

The learned mode transition matrix was then used in the filter bank when
processing the measurements in the second data set. The result is shown in
Fig. 2. In Fig. 2(a) the speed and most likely mode sequence
$\hat{\delta}_{1:k}$ estimated by the filter bank from one of the measurement
sequences are shown. And in Fig. 2(b) the horizontal plane errors at the end
of the trajectory are shown.

From the figures, the following things can be observed. The filter bank
adaptively selects different motion modes depending on the gait speed. That
is, when walking, the filter bank frequently selects the stationary motion
mode as the most likely mode, whereas when running, the filter bank never
selects the stationary mode. Comparing the horizontal positioning error of the
OpenShoe system and the filter bank the systems have about the same cross-
track error, whereas the filter bank has a significantly smaller along-track
error. The OpenShoe system has an along-track bias error of about minus one
meter. This is likely due to the high detection threshold used in the zero-
velocity detector, causing zero-velocity updates to be applied even when the
system is moving. This, in turn, causes part of each step to be cut away,
especially when the user walks and the foot’s transition from stationary to
moving is less distinct. Thanks to the ability of the filter bank to
adaptively select the detection threshold and the covariance used in the zero-
velocity update, this effect is reduced. However, looking at the vertical root
mean square (RMS) error, both systems have an error of several meters. How to
reduce this error will be illustrated next.

(a) Estimated speed and most likely motion mode versus time.

(b) Horizontal position error at end of the trajectory. Also shown are the
mean error (blue dot) and 95% confidence interval (blue circle) calculated
from the filter covariance, assuming the error to be normally distributed.

Figure 2: Results from varying gait speed application example.

(a) Estimated height and most likely motion mode versus time.

(b) Height versus horizontal position error at end of the trajectory. Also
shown are the mean error (blue dot) and 95% confidence interval (blue circle)
calculated from the filter covariance, assuming the error to be normally
distributed.

Figure 3: Results from return to same height application example.

###  V-B Example: Return to same height

The vertical position error often grows faster than the horizontal position
error in foot-mounted zero-velocity aided inertial navigation systems. This is
because the beginning and end of each step, when the foot is mainly moving in
the vertical direction, are cut away by the zero-velocity updates. However,
the vertical error growth can be significantly reduced by noting that most of
the time a person is walking on flat ground and the foot should return to the
same height [3]. To design a system using the proposed framework that
incorporates this information, we introduce the following three motion modes:
the unconstrained motion mode ($\delta_{k}=1$), the stationary at new height
motion mode ($\delta_{k}=2$), and the stationary at the same height motion
mode ($\delta_{k}=3$). A motion mode transition diagram illustrating how the
system is assumed to transition between these motion modes is shown in Fig.
1(b). Here it has been assumed that if the system becomes stationary at a new
height, the system will be stationary for at least two samples. Hence the
probability of transitioning from mode $\delta_{k}=2$ to mode $\delta_{k}=3$
is one.

Next, we introduce the following extended system dynamics

$f_{\text{\tiny
ext}}(x_{k},\delta_{k},u_{k},\eta_{k})=\\!\begin{bmatrix}f(x_{k},\delta_{k},u_{k},\eta_{k})\\\
\mathds{1}(\delta_{k}=1)\xi_{k}+\mathds{1}(\delta_{k}\neq 1)[p_{k}]_{3}\\\
\end{bmatrix},$ (13a) where $\mathds{1}(\cdot)$ denotes the indicator function
and $[a]_{j}$ denotes $j$:th element of the vector $a$. Furthermore, the
auxiliary state $\xi_{k}$ stores the height from the time when the system was
last stationary.

The stochastic motion constraints are as follows

$h_{\text{\tiny
ext}}(x_{k},u_{k},\delta_{k})=\begin{cases}\hfil~{}~{}~{}0\hfil,&\delta_{k}=1\\\
~{}~{}h_{0}(x_{k},u_{k}),&\delta_{k}=2\\\ \begin{bmatrix}h_{0}(x_{k},u_{k})\\\
[p_{k}]_{3}-\xi_{k}\end{bmatrix},&\delta_{k}=3\\\ \end{cases}.$ (13b)

The “hardness” of these constraints is controlled by the covariance matrices

$R_{\text{\tiny
ext}}(\delta_{k})=\begin{cases}\sigma^{2}_{nc}I,&\delta_{k}=1\\\
R_{0}(\delta_{k}),&\delta_{k}=2\\\
R_{0}(\delta_{k})\oplus\sigma^{2}_{h}(\delta_{k})I,&\delta_{k}=3\\\
\end{cases},$ (13c)

where $\sigma^{2}_{h}(\delta_{k})$ denotes the variance associated with the
stochastic constraints on height.

The jump Markov model defined by (13a), (13b), and (13c) was used to design a
filter bank, and two data recordings were collected while a person walked on
flat ground and then climbed and descended a stair. Following the same
procedure as in the previous example, the data recordings were used to create
two data sets with 50 measurement sequences each. The measurement sequences in
the first data set were used to learn the mode transition matrix $\Pi$. The
initial value and end result of the learning process were

$\widehat{\Pi}^{\text{\tiny init}}=\begin{bmatrix}1/3&0&1/2\\\ 1/3&0&0\\\
1/3&1&1/2\end{bmatrix}~{}\text{and}~{}~{}\widehat{\Pi}^{\text{\tiny
end}}=\begin{bmatrix}0.976&0&0.031\\\ 0.003&0&0\\\
0.021&1&0.969\end{bmatrix},$

respectively, and the optimization converged in less than 10 iterations. Once
again the matrix corresponds to a mode transition system with low-pass
characteristics, where the average percentages of time spent in motion modes
one, two, and three is 56.2%, 0.2%, and 43.6%, respectively.

The learned mode transition matrix was then used in the filter bank when
processing the measurements in the second data set. Once again, pruning was
used to keep the tree size to a maximum of nine leaves and a pruning strategy
where the most leaves were retained was used. The result is shown in Fig. 3.
In Fig. 3(a) the height and the most likely mode sequence $\hat{\delta}_{1:k}$
estimated by the filter bank is shown. And in Fig. 3(b) the height versus
horizontal plane errors at the end of the trajectory are shown.

From the figures, the following things can be observed. The filter bank can
detect whether the user walks on flat ground or not, i.e., whether the foot
returns to the same height as in the previous step or not. Therefore it
effectively reduces the height error from several meters to a few decimeters.
A minor increase in the horizontal position error is observed with the filter
bank.

##  VI Discussion and Conclusions

A framework for tightly integrated motion classification and state estimation
in motion-constrained inertial navigation systems has been presented. The
framework provides a structured and theoretically sound way to design motion-
constrained inertial navigation systems that can learn different motion
patterns. The application of the framework has been illustrated via two
examples of foot-mounted zero-velocity aided inertial navigation systems. The
examples show that a significant performance gain can be achieved compared to
a standard foot-mounted inertial navigation system. However, the price paid is
an increased computational complexity (proportional to the number of filters
in the filter bank) and an increased number of system parameters to tune. The
challenges with tuning the system parameters can to some extent be alleviated
via the proposed parameter learning method. Still, tuning the system
parameters can be challenging, especially if many motion modes are included in
the system model.

##  VII Outlook and Future Research

A basic version of the filter bank framework has been presented and tested
with models that have a few motion modes. Many possible improvements and open
research questions exist. First and foremost, the framework should be compared
against other adaptive zero-velocity detector frameworks. More complex models
that include a variety of motion modes should also be explored. Related to
that, hieratical model structures constructed of several linked small hidden
Markov models for the motion states are of special interest to obtain a
computationally attractive algorithm. Further, to obtain a smoother transition
between motion modes and a more robust algorithm framework, the feasibility of
substituting the normal distribution in (6c) with a more heavy-tailed
distribution, such as the Student-t distribution, should be explored.
Moreover, since the probability of transitioning between motion modes varies
both with time and the motion dynamics, the use of a constant mode transition
probability matrix $\Pi$ is clearly suboptimal. Therefore, the feasibility of
learning a time-varying and motion-dynamic dependent transition matrix should
be explored.

## Acknowledgment

This work has been partially funded by the Swedish Research Council project
2020-04253 _Tensor-field based localization_ and the Dutch Research Council
(NWO) research program Veni project 18213 _Sensor Fusion For Indoor
Localisation Using The Magnetic Field_.

## References

  * [1] J. Wahlström and I. Skog, “Fifteen years of progress at zero velocity: A review,” _IEEE Sensors J._ , vol. 21, no. 2, pp. 1139–1151, Jan. 2021\.
  * [2] W. Zhang, D. Wei _et al._ , “The improved constraint methods for foot-mounted PDR system,” _IEEE Access_ , vol. 8, pp. 31 764–31 779, Feb. 2020.
  * [3] K. Abdulrahim, C. Hide _et al._ , “Using constraints for shoe mounted indoor pedestrian navigation,” _The Journal of Navigation_ , vol. 65, no. 1, p. 15–28, 2012.
  * [4] N. Kronenwett, S. Qian _et al._ , “Elevator and escalator classification for precise indoor localization,” in _Int. Conf. on Indoor Positioning and Indoor Navigation (IPIN)_ , Nantes, France, Sep. 2018, pp. 1–8.
  * [5] J. Wahlström, I. Skog _et al._ , “Zero-velocity detection – A Bayesian approach to adaptive thresholding,” _IEEE Sensors L._ , vol. 3, no. 6, pp. 1–4, Jun. 2019.
  * [6] C.-S. Jao and A. M. Shkel, “ZUPT-aided INS bypassing stance phase detection by using foot-instability-based adaptive covariance,” _IEEE Sensors J._ , vol. 21, no. 21, pp. 24 338–24 348, Nov. 2021.
  * [7] E. Mazor, A. Averbuch _et al._ , “Interacting multiple model methods in target tracking: a survey,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 34, no. 1, pp. 103–123, Jan. 1998.
  * [8] H. Blom and Y. Bar-Shalom, “The interacting multiple model algorithm for systems with Markovian switching coefficients,” _IEEE Trans. Autom. Control_ , vol. 33, no. 8, pp. 780–783, Aug. 1988.
  * [9] X. Rong Li and V. Jilkov, “Survey of maneuvering target tracking. part v. multiple-model methods,” _IEEE Trans. Aerosp. Electron. Syst._ , vol. 41, no. 4, pp. 1255–1321, Oct. 2005.
  * [10] X.-R. Li and Y. Bar-Shalom, “Multiple-model estimation with variable structure,” _IEEE Transactions on Automatic Control_ , vol. 41, no. 4, pp. 478–493, Apr. 1996.
  * [11] J. Sola, “Quaternion kinematics for the error-state kalman filter,” _arXiv preprint arXiv:1711.02508_ , 2017.
  * [12] M. Kok, J. D. Hol _et al._ , “Using inertial sensors for position and orientation estimation,” _Foundations and Trends in Signal Processing_ , vol. 11, no. 1, 2017.
  * [13] M. Moakher, “Means and averaging in the group of rotations,” _SIAM J. Matrix Anal. Appl._ , vol. 24, no. 1, p. 1–16, Jan. 2002.
  * [14] V. Jilkov and X. Li, “Online Bayesian estimation of transition probabilities for Markovian jump systems,” _IEEE Trans. Signal Process._ , vol. 52, no. 6, pp. 1620–1630, Jun. 2004.
  * [15] E. Özkan, F. Lindsten _et al._ , “Recursive maximum likelihood identification of jump Markov nonlinear systems,” _IEEE Trans. Signal Process._ , vol. 63, no. 3, pp. 754–765, Feb. 2015.
  * [16] J.-O. Nilsson, A. K. Gupta _et al._ , “Foot-mounted inertial navigation made easy,” in _Int. Conf. on Indoor Positioning and Indoor Navigation (IPIN)_ , Busan, South Korea, Oct. 2014, pp. 24–29.
  * [17] J. A. Farrell and M. Barth, _The Global Positioning System and Inertial Navigation_. McGraw-Hill, 1998.


# Transfer-tensor description of memory effects in open-system dynamics and
multi-time statistics

Stefano Gherardini<EMAIL_ADDRESS>Dipartimento di Fisica e
Astronomia, Università degli Studi di Firenze, via G. Sansone 1, I-50019 Sesto
Fiorentino, Italy  LENS, CNR-INO, and QSTAR, via N. Carrara 1, I-50019 Sesto
Fiorentino, Italy  INFN Sezione di Firenze, via G. Sansone 1, I-50019 Sesto
Fiorentino, Italy  Andrea Smirne<EMAIL_ADDRESS>Institute of
Theoretical Physics, and IQST, Universität Ulm, Albert-Einstein-Allee 11,
89069 Ulm, Germany  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  Susana
Huelga<EMAIL_ADDRESS>Institute of Theoretical Physics, and IQST,
Universität Ulm, Albert-Einstein-Allee 11, 89069 Ulm, Germany  Filippo Caruso
<EMAIL_ADDRESS>Dipartimento di Fisica e Astronomia, Università degli
Studi di Firenze, via G. Sansone 1, I-50019 Sesto Fiorentino, Italy  LENS,
CNR-INO, and QSTAR, via N. Carrara 1, I-50019 Sesto Fiorentino, Italy

###### Abstract

The non-Markovianity of an arbitrary open quantum system is analyzed in
reference to the multi-time statistics given by its monitoring at discrete
times. On the one hand, we exploit the hierarchy of inhomogeneous transfer
tensors, which provides us with relevant information about the role of
correlations between the system and the environment in the dynamics. The
connection between the transfer-tensor hierarchy and the CP-divisibility
property is then investigated, by showing to what extent quantum Markovianity
can be linked to a description of the open-system dynamics by means of the
composition of 1-step transfer tensors only. On the other hand, we introduce
the set of stochastic transfer tensor transformations associated with local
measurements on the open system at different times and conditioned on the
measurement outcomes. The use of the transfer-tensor formalism accounts for
different kinds of memory effects in the multi-time statistics and allows us
to compare them on a similar footing with the memory effects present in non-
monitored non-Markovian dynamics, as we illustrate on a spin-boson case study.

##  I Introduction

In a dynamical quantum system the interaction with the external environment
typically leads to non-Markovian behaviours, which, broadly speaking, are
associated with the presence of memory effects in the system evolution [1, 2,
3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22].
Although by now many distinct definitions of quantum non-Markovianity have
been introduced, two main pathways can be identified [23].

In one case, the focus is on the evolution of the open-system state at
different times, as fixed, e.g., by the dynamical maps or the master equations
that characterize the open-system dynamics [24, 25], while in the other case
the focus is on the statistics associated with a sequence of local
measurements (or other active interventions) performed on the open system at
subsequent times. These two approaches to non-Markovianity are indeed
inherently different, as the former addresses the predictions related with
observables at a single time, while the latter concerns the multi-time
statistics. In addition to such a difference, which already appears at the
classical level [26, 27], the quantum nature of the system at hand implies a
further, fundamental distinction between the single-time and the multi-time
notions of non-Markovianity. Contrary to what happens classically, in quantum
systems it is generally not possible, even in principle, to access the
information associated with sequential measurements without disturbing the
system that is being measured. Indeed, quantum measurements alter the current
state of the system, as well as its correlations with the environment, which
will later result in a modified evolution [7, 9, 29, 32, 28, 30, 34, 33, 31,
35]. As a consequence, the presence of intermediate measurements can bring
along specific forms of memory, which combine in a non-trivial way with the
memory strictly due to the system-environment interaction.

In this paper, we investigate the relation between the memory effects
appearing in an open quantum system dynamics and those associated with the
multi-time statistics due to sequential measurements, by means of the transfer
tensor (TT) method [36, 37, 38, 39]. We will show that the latter, which was
introduced to treat efficiently the long-time dynamics of open quantum
systems, also allows one to treat memory effects on the dynamics and the
multi-time statistics on a similar footing. First, we prove to what extent the
divisibility property of the dynamical maps, which is the defining property of
quantum non-Markovianity according to [7, 4], is linked to the hierarchy of
TTs also in the inhomogeneous case, i.e., beyond the time-translational
invariant regime explored in the original paper [36]. Then, we extend the
definition of the TTs to the situation where the open system is measured at
subsequent times, via the corresponding conditional stochastic dynamics. This
allows us to identify different forms of memory in the multi-time statistics,
related with the interplay between the correlations and environmental-state
transformations induced by, respectively, the interaction between the open
system and the environment and the sequential measurements. After introducing
proper quantifiers related with the $L_{2}$-norm of the multi-step (i.e., step
greater than $1$) TTs, we compare in a case study the relevance of such memory
effects for different kinds of measurements, including the case without any
intermediate measurement.

The rest of the paper is organized as follows. In Sec.II, we recall the main
features of the TTs formalism that are relevant for our analysis, focusing in
particular on the recursive construction of the TTs. In Sec.III, we prove that
if only multi-step TTs are different from zero the open-system dynamics is
divisible, while the converse statement does not hold, as we show by means of
an explicit example. In Sec.IV, we first generalize the TTs formalism to the
case of multi-time measurements, by introducing a family of TTs conditioned on
the sequence of measurement outcomes, and then we exploit it to characterize
the different kinds of memory effects associated with the multi-time
statistics. Finally, the conclusions and outlooks of our analysis are
discussed in Sec.V

##  II Inhomogeneous transfer tensors hierarchy

We begin by briefly recalling the TT formalism for the dynamics of an open
quantum system [36] and the microscopic characterization of the TTs recently
derived in [38], which allows to build up the complete TTs hierarchy and to
connect it with the system-environment correlations. Such a construction will
be extended to an evolution conditioned on the occurrence of sequential
measurements in the following of the paper.

###  II.1 General microscopic definition of TTs

Let us consider a quantum system $\mathbb{S}$ in interaction with an
environment $\mathbb{E}$. The dynamics of the composite system
$\mathbb{S}+\mathbb{E}$ is governed by a possibly time-dependent Hamiltonian
$H(t)=H_{0}(t)+H_{\rm int}$, where $H_{\rm int}$ takes into account the
interaction of the system with the environment, while
$H_{0}(t)=H_{\mathbb{S}}(t)+H_{\mathbb{E}}$ concerns the uncoupled time
evolution of $\mathbb{S}$ and $\mathbb{E}$. The interaction Hamiltonian
$H_{\rm int}$ is fixed but unknown, and $H_{\mathbb{S}}(t)$ can be time-
dependent. Assuming an initial product state between the open system and the
environment and a fixed initial state of $\mathbb{E}$, i.e.,
$\rho_{t_{0}}=\rho_{\mathbb{S},t_{0}}\otimes\rho_{\mathbb{E},t_{0}}$, the
dynamics of $\mathbb{S}$ can be described by means of the formalism of quantum
dynamical maps [24, 25, 8]. In fact, one can define a one parameter family of
completely positive, trace-preserving (CPTP) dynamical maps
$\left\\{\Phi(t,t_{0})\right\\}_{t\geq t_{0}}$, with
$\Phi(t,t_{0}):\mathcal{S}(\mathcal{H}_{\mathbb{S}})\rightarrow\mathcal{S}(\mathcal{H}_{\mathbb{S}})$,
where $\mathcal{S}(\mathcal{H}_{\mathbb{S}})$ denotes the sets of density
operators (non-negative operators with unit trace) acting on
$\mathcal{H}_{\mathbb{S}}$ and $\Phi(t,t_{0})$ is CPTP at any time $t$.
Hereafter, the action of the dynamical maps of the system will be taken into
account at the discrete time instants (not necessarily equally spaced)
$t_{k}$, $k=1,\ldots,m$, and we will use the notation
$\Phi(t_{k},t_{0})\equiv\Phi_{k}$ and

$\rho_{\mathbb{S},k}=\rho_{\mathbb{S}}(t_{k})=\Phi_{k}[\rho_{\mathbb{S},0}].$
(1)

This will make the comparison with the sequential-measurement scenario more
transparent and will further allow us to directly apply the TT method [36].
The latter relies in fact on a family of maps, the TTs $T_{k,j}$, which are
defined by the relation

$T_{k,0}=\Phi_{k}-\sum_{j=1}^{k-1}T_{k,j}\Phi_{j},$ (2)

equivalently expressed by

$\Phi_{k}=\sum_{j=0}^{k-1}T_{k,j}\Phi_{j}.$ (3)

In other terms, the state at the time instant $t_{k}$ is obtained by
propagating the states at $t_{j}$ according to the equation
$\rho_{\mathbb{S},k}=\sum_{j=0}^{k-1}T_{k,j}\rho_{\mathbb{S},j}$. We stress
that, at variance with the original formulation, we are not assuming time-
invariance; thus our analysis will also apply to time-dependent Hamiltonians
and non-stationary initial environmental states [36].

A useful expression for generic, possibly inhomogeneous TTs was recently
derived in [38]. Given the global unitary dynamics of the composite system
$\mathbb{S}+\mathbb{E}$ from $t_{s}$ to $t_{k}$, i.e.,
$\mathcal{U}_{t_{s}:t_{k}}[\rho]\equiv
U_{t_{s}:t_{k}}\rho\,U_{t_{s}:t_{k}}^{\dagger}$ with
$U_{t_{s}:t_{k}}\equiv\mathcal{T}\exp(-i\int_{t_{s}}^{t_{k}}H(\tau)d\tau)$
($\mathcal{T}$ is the time ordering operator), let us introduce the CPTP maps
$\Gamma_{k|k-n}:\mathcal{S}(\mathcal{H}_{\mathbb{S}})\rightarrow\mathcal{S}(\mathcal{H}_{\mathbb{S}})$
defined as

$\Gamma_{k|k-n}[\sigma_{\mathbb{S}}]\equiv{\rm
Tr}_{\mathbb{E}}\left[\mathcal{U}_{t_{k-n}:t_{k}}[\sigma_{\mathbb{S}}\otimes\sigma_{\mathbb{E},k-n}]\right],$
(4)

where ${\rm Tr}_{\mathbb{S}/\mathbb{E}}[\cdot]$ denotes the partial trace
w.r.t. the Hilbert space of $\mathbb{S}$ or $\mathbb{E}$,
$\sigma_{\mathbb{S}}$ denotes a density operator on $\mathcal{H}_{\mathbb{S}}$
and $\sigma_{\mathbb{E},k-n}$ is a time-dependent density operator on
$\mathcal{H}_{\mathbb{E}}$. Here, we will focus on the case where
$\sigma_{\mathbb{E},k-n}$ is the state of the environment at $t_{k-n}$. It can
then be proven that the $n-$step TT $T_{k,k-n}$ is given by the recursive
relation [38]

$T_{k,k-n}=\Gamma_{k|k-n}-\sum_{j=1}^{n-1}T_{k,k-j}\Gamma_{k-j|k-n}$ (5)

that allows one to reconstruct the whole hierarchy of TTs from the CPTP maps
$\Gamma_{k|k-n}$, as shown explicitly in [38]. In the next paragraph, we
repeat the construction for the lowest orders of the hierarchy, both for the
sake of illustration and since it will be useful in the following when we will
move to the stochastic case.

###  II.2 Reconstruction of the TTs hierarchy

First, let us combine together the definition of the dynamical maps of the
system at the discrete time instants $t_{k}$, $k=1,\ldots,n$, i.e., Eq. (1),
and the recursive general expression of the TT transformation given by Eq.(3).
In this way, the latter can be written via the following recursive relation:

$\Phi_{k}=\underline{T}_{k}\underline{\Xi}_{k-1}[\rho_{\mathbb{S},0}],$ (6)

where

$\underline{T}_{k}\equiv\left(T_{k,0},T_{k,1},\ldots,T_{k,k-1}\right)$

and

$\underline{\Xi}_{k-1}[\rho_{\mathbb{S},0}]\equiv\begin{pmatrix}\mathbbm{1}_{\mathbb{S}}[\rho_{\mathbb{S},0}]\\\
\Phi_{1}=\underline{T}_{1}\underline{\Xi}_{0}[\rho_{\mathbb{S},0}]\\\
\vdots\\\
\Phi_{k-1}=\underline{T}_{k-1}\underline{\Xi}_{k-2}[\rho_{\mathbb{S},0}]\end{pmatrix},$

which are valid for $n\geq 1$
($\underline{T}_{0}=\underline{\Xi}_{-1}=\mathbbm{1}_{\mathbb{S}}[\rho_{\mathbb{S},0}]$,
identity operator on $\mathbb{S}$). For the sake of clarity, we show the first
terms of the recursive expansion of Eq. (6):

$\begin{cases}\Phi_{0}=\mathbbm{1}_{\mathbb{S}}[\rho_{\mathbb{S},0}]\\\
\Phi_{1}=T_{1,0}\Phi_{0}=T_{1,0}[\rho_{\mathbb{S},0}]\\\
\Phi_{2}=T_{2,0}\Phi_{0}+T_{2,1}\Phi_{1}=\left(T_{2,0}+T_{2,1}T_{1,0}\right)[\rho_{\mathbb{S},0}]\\\
\Phi_{3}=T_{3,0}\Phi_{0}+T_{3,1}\Phi_{1}+T_{3,2}\Phi_{2}\\\
\,\,\,\,\,\,\,\,\,=\left(T_{3,0}+T_{3,1}T_{1,0}+T_{3,2}T_{2,0}+T_{3,2}T_{2,1}T_{1,0}\right)[\rho_{\mathbb{S},0}].\\\
\end{cases}$ (7)

By combining Eqs. (1) and (6), we immediately find that

$T_{1,0}[\rho_{\mathbb{S},0}]=\Phi_{1}[\rho_{\mathbb{S},0}],$ (8)

so that $\rho_{\mathbb{S},1}=T_{1,0}[\rho_{\mathbb{S},0}]$. Then, by
continuing the recursion, one has

$\rho_{\mathbb{S},2}=T_{2,0}[\rho_{\mathbb{S},0}]+T_{2,1}[T_{1,0}[\rho_{\mathbb{S},0}]].$

Here, it is worth observing that $T_{2,1}[T_{1,0}[\rho_{\mathbb{S},0}]]$
returns the component of $\rho_{\mathbb{S},2}$ conditioned on finding the
system in the states $\rho_{\mathbb{S},0}$ and $T_{1,0}[\rho_{\mathbb{S},0}]$
at the time instants $t_{0}$ and $t_{1}$ and no system-environment
correlations were present. In fact, given the maps $\Gamma_{k|k-n}$ defined in
Eq. (4), the 1-step TTs $T_{k,k-1}$ satisfy

$T_{k,k-1}=\Gamma_{k|k-1}$ (9)

as a direct consequence of Eq. (5) (for the trivial case $n=1$). In other
terms, we have

$T_{2,0}[\rho_{\mathbb{S},0}]=\rho_{\mathbb{S},2}-\Gamma_{2|1}[\rho_{\mathbb{S},1}],$
(10)

i.e.,

$T_{2,0}=\Phi_{2}-\Gamma_{2|1}\Phi_{1}.$ (11)

For the sake of brevity, from here on we will remove the subscript
$\mathbb{S}$ from the notation; it will used (together with $\mathbb{E}$) only
if necessary.
Moving on to the level $k=3$ of the hierarchy, by applying Eq. (7) one has
that the reduced density operator of $\mathbb{S}$ at $t_{3}$ is equal to

$\displaystyle\rho_{3}$ $\displaystyle=$ $\displaystyle
T_{3,0}[\rho_{0}]+T_{3,1}[T_{1,0}[\rho_{0}]]+T_{3,2}[T_{2,0}[\rho_{0}]]$
$\displaystyle+$ $\displaystyle T_{3,2}[T_{2,1}[T_{1,0}[\rho_{0}]]].$

Using the identity in Eq. (9) for $T_{3,2}$ and $T_{2,1}$, as well as Eq.
(10), we have

$\begin{cases}T_{3,2}[T_{2,0}[\rho_{0}]]=\Gamma_{3|2}[\rho_{2}-\Gamma_{2|1}[\rho_{1}]]\\\
T_{3,2}[T_{2,1}[T_{1,0}[\rho_{0}]]]=\Gamma_{3|2}\Gamma_{2|1}[\rho_{1}].\end{cases}$
(12)

Moreover, $T_{3,1}$ can be expressed by exploiting Eq. (5) (for $n=2$) that
along with Eq. (9) gives

$T_{3,1}\equiv\Gamma_{3|1}-\Gamma_{3|2}\Gamma_{2|1}.$ (13)

Thus, $T_{3,1}$ can be written as a function of a term involving only one-step
$\Gamma$s, i.e., $\Gamma_{3|2}\Gamma_{2|1}$, and of the $2-$step
$\Gamma_{3|1}$. By combining together all terms of Eqs. (12) and (13), one has

$\rho_{3}=T_{3,0}[\rho_{0}]+\Gamma_{3|1}[\rho_{1}]+\Gamma_{3|2}[\rho_{2}]-\Gamma_{3|2}\Gamma_{2|1}[\rho_{1}],$

i.e.,

$T_{3,0}=\Phi_{3}-\Gamma_{3|2}\Phi_{2}+(\Gamma_{3|2}\Gamma_{2|1}-\Gamma_{3|1})\Phi_{1}.$
(14)

As mentioned, the procedure can be generalized to any $k$ via Eq. (5), so that
the whole hierarchy of TTs can be expressed in terms of the maps
$\Gamma_{k|k-n}$ [38] (indeed, $\Phi_{k}=\Gamma_{k|0}$). Importantly, this
construction reveals the influence of the correlations between the system and
the environment into the dynamics of the open quantum system. In fact, the
CPTP maps (4) generate a hierarchy of $1-$, $2-$,…,$n-$step TTs conditioned to
the fact that the system passes through product states at the different steps
taken into account. As we will show in the next section, this is naturally
linked to the presence of memory effects in the open-system dynamics.

##  III Transfer tensors and Markovianity of the dynamics

Here, we derive a definite connection between a property of the TTs – the
possibility to build up the whole hierarchy only in terms of $1-$step TTs –
and the non-Markovianity of the dynamics. In particular, we rely on the
definition introduced in [4], which identifies Markovian dynamics as those
described by a family of CP-divisible dynamical maps. Adapting the original
definition in [4] to the case of a discrete set of times, then the following
definition can be stated:
The (discrete) dynamics $\left\\{\Phi_{k}\right\\}_{k=1,\ldots,m}$ is
Markovian when for any $k\geq j\geq 0$ there is a CPTP map $\mathcal{E}_{k,j}$
such that $\Phi_{k}=\mathcal{E}_{k,j}\Phi_{j}$.

Hence, if we assume that $T_{k,k-n}=0$ for $n\geq 2$, then by applying
recursively Eq. (3) we have

$\Phi_{k}=T_{k,k-1}T_{k-1,k-2}\cdots T_{j+1,j}\Phi_{j}$ (15)

for any $k\geq j\geq 0$. But then Eq. (9) implies

$\Phi_{k}=\Gamma_{k|k-1}\Gamma_{k-1|k-2}...\Gamma_{j+1|j}\Phi_{j}.$ (16)

Any conditional map $\Gamma_{k|j}$ is CPTP by construction, see Eq.(4), so
that we can conclude that any map $\Phi_{k}$ can be decomposed as
$\Phi_{k}=\mathcal{E}_{k,j}\Phi_{j}$, with $\mathcal{E}_{k,j}$ CPTP map:
_If only one-step TTs are different from zero, the resulting (discrete-time)
evolution is CP-divisible_.
We further note that if we restrict to the case of equally-spaced time
instants, $t_{k}=k\Delta$, and translational invariant TTs (i.e.,
$T_{k,k-n}=T_{n,0}$), then $T_{k,k-n}=0$ for $n\geq 2$ implies that the open-
system dynamics is not only Markovian, but also a semigroup, i.e.,
$\Phi_{k}=(\Phi_{1})^{k}$, thus recovering what shown in [36].

The previous result, besides linking a property of the hierarchy of TTs to the
Markovianity of the corresponding dynamics, provides us with an explicit
illustration of one of the physical meanings of such property. In fact, from
Eq.(16) we see that when one-step TTs are the only non-zero TTs, the dynamical
maps can be obtained by using the maps $\Gamma$ only. System-environment
correlations, despite being present as a consequence of the interaction term
within the unitary operators, will not affect the reduced dynamics of the open
system, so that at any time $t_{k-n}$ the actual global state can be
effectively replaced by a product state
$\rho_{\mathbb{S},k-n}\otimes\sigma_{\mathbb{E},k-n}$, see Eq.(4).
Experimentally, Eq. (16) can be validated by independently measuring the maps
$\\{\Phi_{k}\\}$ and the 1-step $\\{\Gamma_{k,k-1}\\}$, obtained by breaking
the system-environment correlations at the time instants $t_{k}$. In
particular, the former can be measured by means of quantum tomography
processes, while the latter, e.g., by preparing a fresh copy of the system in
the state $\sigma_{\mathbb{S},k}$ [13], which is known by virtue of the
previously reconstructed map $\Phi_{k}$.

The fact that only one-step TTs are non-zero is in general a stronger
requirement than the CP-divisibility of the dynamics, as follows from the
analysis of [23] (see also [40]) and the explicit example provided below.
Indeed, having only one-step TTs different from zero is only a _sufficient
criterion_ (and thus not necessary) for CP-divisible discrete-time evolutions.
In this regard, let us consider the simplest case, namely $m=2$, whereby we
recall that $m$ is the last (and greater) value of the index $k$ in the
sequence of time instants $t_{k}$. Hence, the only non-trivial requirement for
CP-divisibility is the existence of a CPTP map $\mathcal{E}_{2,1}$ obeying the
relation

$\mathcal{E}_{2,0}=\Phi_{2}=\mathcal{E}_{2,1}\Phi_{1}=\mathcal{E}_{2,1}T_{1,0}$

since by definition $T_{1,0}\equiv\Phi_{1}$. Now, according to the definition
of the TT transformation, $\Phi_{2}$ is also equal to

$\Phi_{2}=T_{2,0}+T_{2,1}T_{1,0}.$

Therefore, assuming that the inverse of $T_{1,0}$ (i.e., $(T_{1,0})^{-1}$)
does exist, one has

$\mathcal{E}_{2,1}=T_{2,0}\left(T_{1,0}\right)^{-1}+T_{2,1}.$ (17)

What we will show is that $\mathcal{E}_{2,1}$ can be CPTP also if $T_{2,0}\neq
0$ and even if $\left(T_{1,0}\right)^{-1}$ is not be generally CPTP. In fact,
consider the reduced dynamics of a two-level system, fixed by the unitary
dynamics given by the global Hamiltonian
$H=\sigma_{z}+H_{\mathbb{E}}+\sigma_{z}\otimes B$, where the free
environmental Hamiltonian $H_{\mathbb{E}}$ and the interaction term commute:
$[H_{\mathbb{E}},B]=0$. Now, also assume that the global state at the initial
time $t_{0}=0$ is $\rho(0)=\rho_{\mathbb{S}}(0)\otimes\rho_{\mathbb{E}}(0)$,
where the environmental state is stationary with respect to the free dynamics
such that $[H_{\mathbb{E}},\rho_{\mathbb{E}}]=0$. As a consequence, on the one
hand, the reduced dynamics of the open system is a pure dephasing, namely
$\rho_{11}(t)=\rho_{11}(0)$ and $\rho_{10}(t)=k(t)\rho_{10}(0)$, where
$\rho_{ij}(t)=\langle i|\rho_{\mathbb{S}}(t)|j\rangle$ ($\left\\{\lvert
1\rangle,\lvert 0\rangle\right\\}$ are the eigenvectors of $\sigma_{z}$) and
$k(t)=e^{i\omega t}\mbox{tr}_{\mathbb{E}}[e^{2iBt}\rho_{\mathbb{E}}]$. On the
other hand, the state of the environment will not evolve in time, also when
interacting with the system, i.e.,
$\rho_{\mathbb{E}}(t)=\rho_{\mathbb{E}}(0)$. For any two times $t_{2}\geq
t_{1}\geq 0$, we will then have

$T_{2,1}=\Gamma_{2|1}=\Phi_{1},$

where we used Eq.(9) and Eq.(4), respectively. As a result (see also Eq.(11)),

$T_{2,0}=\Phi_{2}-\Phi^{2}_{1}.$

Specifically, using the matrix representation in the basis $\left\\{\lvert
1\rangle,\lvert 0\rangle\right\\}$, for a generic state $\rho_{\mathbb{S}}$ we
have

$\mathcal{E}_{2,1}[\rho_{\mathbb{S}}]=\begin{pmatrix}\rho_{11}&\frac{k(t_{2})}{k(t_{1})}\rho_{10}\\\
\frac{k^{*}(t_{2})}{k^{*}(t_{1})}\rho_{01}&\rho_{00}\end{pmatrix},$

which is positive definite for any $\rho_{\mathbb{S}}$ if and only if
$|k(t_{2})|\leq|k(t_{1})|$. Note that this implies that $\mathcal{E}_{2,1}$ is
also CPTP under the same condition, since for pure dephasing positivity and
complete positivity coincide. On the other hand,

$T_{2,0}[\rho_{\mathbb{S}}]=\begin{pmatrix}0&(k(t_{2})-k(t_{1})^{2})\rho_{10}\\\
(k^{*}(t_{2})-k^{*}(t_{1})^{2})\rho_{01}&0\end{pmatrix}$

that is equal to 0 for any $\rho_{\mathbb{S}}$ if and only if
$k(t_{2})=k(t_{1})^{2}$. The latter condition implies but is not implied by
the former (indeed, $|k(t)|\leq 1$), which proves our claim:
$\mathcal{E}_{2,1}$ can be CPTP, and hence the dynamics CP-divisible, even
though the $2$-step transfer tensor $T_{2,0}$ is different from zero.

##  IV Transfer tensors and multi-time statistics

In this paragraph, we move on to the second part of our investigation, where
we consider an open system that, besides interacting with the environment, is
measured at some discrete instants of time. We will first introduce a proper
counterpart of the transfer tensors in such dynamical regime, and will then
show how it can be used to account for the memory present in the multi-time
statistics. Before that, let us specify the framework we refer to.

Assume that the open quantum system $\mathbb{S}$ is monitored at the same $m$
instants of time where the discrete dynamics has been evaluated so far. In
particular, we take a sequence of quantum measurements locally performed on
$\mathbb{S}$ according to the observables $\mathcal{O}_{k}\equiv
F_{\theta_{k}}\otimes I_{\mathbb{E}}$; $\\{\theta_{k}\\}$ is the set of the
possible measurement outcomes, and $\\{F_{\theta_{k}}\\}$ denotes the set of
positive semi-definite operators on $\mathcal{H}_{\mathbb{S}}$ satisfying the
relation $\sum_{\theta_{k}}F_{\theta_{k}}=\mathbbm{1}_{\mathbb{S}}$ $\forall
k$. The probability that the outcome $\theta_{k}$ associated with the
measurement operator $F_{\theta_{k}}$ occurs is equal to ${\rm
Tr}[\rho_{\mathbb{S},k}F_{\theta_{k}}]$, while the post-measurement state of
$\mathbb{S}$ equals to
$\widehat{\rho}_{\mathbb{S},k}\equiv\mathcal{M}_{\theta_{k}}[\rho_{\mathbb{S},k}]/{\rm
Tr}[\mathcal{M}_{\theta_{k}}[\rho_{\mathbb{S},k}]]$, where we introduced the
super-operator $\mathcal{M}_{\theta_{k}}[\rho_{\mathbb{S},k}]\equiv
M_{\theta_{k}}\rho_{\mathbb{S},k}M^{\dagger}_{\theta_{k}}$ and
$M_{\theta_{k}}$ fulfills the identity
$F_{\theta_{k}}=M^{\dagger}_{\theta_{k}}M_{\theta_{k}}$. Note that
$\rho_{\mathbb{S},k}$ ($\widehat{\rho}_{\mathbb{S},k}$) denotes the state of
the system before (after) the measurement at the time step $k$.

Importantly, a measurement of the open system affects not only its current
state, but also its future evolution as a consequence of the change in the
correlations between the system and the environment due to the measurement
itself [9]. The very notion of dynamics of the open system becomes more
subtle, since it cannot be clearly separated from the results of the
sequential measurements. It is thus useful to introduce the notion of
conditional dynamics, whereby the system admits a different dynamics for each
sequence of measurement outcomes. Explicitly, the conditional dynamics of the
system (for a given global evolution $\mathcal{U}_{t_{s}:t_{k}}$ and initial
environmental state $\rho_{\mathbb{E},0}$) is fixed by the CP map

$\displaystyle\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}[\rho_{\mathbb{S},0}]$
$\displaystyle\equiv$ $\displaystyle{\rm
Tr}_{\mathbb{E}}\left[\mathcal{U}_{t_{k-1}:t_{k}}(\mathcal{M}_{\theta_{k-1}}\otimes\mathbbm{1}_{\mathbb{E}})\cdots\right.$
(18) $\displaystyle\cdots$
$\displaystyle\left.(\mathcal{M}_{\theta_{1}}\otimes\mathbbm{1}_{\mathbb{E}})\,\mathcal{U}_{t_{0}:t_{1}}[\rho_{\mathbb{S},0}\otimes\rho_{\mathbb{E},0}]\right],$

depending on the time instants $\underline{t}\equiv(t_{1},\ldots,t_{k})$ as
well as on the measurement outcomes
$\underline{\theta}\equiv(\theta_{1},\ldots,\theta_{k-1})$. In other terms,
$\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}$ has to be understood
as a stochastic map, so that we could effectively define different
trajectories in the set of CP maps, each of them associated with a different
sequence of measurement outcomes. In general
$\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}$ is not trace
preserving.

The joint probability distributions to get the measurement outcomes
$\theta_{1},\theta_{2},\ldots,\theta_{k}$ at the time instants
$t_{1},t_{2},\ldots,t_{k}$ is directly linked to the stochastic map
$\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}$ by the following
relation:

$q_{k}\equiv{\rm Prob}(\theta_{k},t_{k};\ldots;\theta_{1},t_{1})={\rm
Tr}[\mathcal{M}_{\theta_{k}}\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}[\rho_{\mathbb{S},0}]].$
(19)

These quantities define the multi-time statistics associated to sequential
measurements at different times and, as recalled in the introduction, suitable
notions of Markovianity can be attributed to them [29, 23, 15, 13, 14, 28, 30,
34, 31].

###  IV.1 Conditional transfer tensors

As anticipated, in order to characterize the conditional dynamics originated
by different sequences of outcomes, we can introduce a stochastic version of
the TTs. Let us denote them as
$\widetilde{T}^{\underline{\theta},\underline{t}}_{k,j}$, with
$k,j=1,\ldots,m$, explicitly pointing out their dependence on the instants and
outcomes of the repeated measurements.

The basic idea is to express the definition of the conditional dynamical map
in a recursive fashion. As first step, note that the global state after the
first measurement is proportional (via the normalization factor) to

$\widehat{\rho}_{1}\sim(M_{\theta_{1}}\otimes
I_{\mathbb{E}})\mathcal{U}_{t_{0}:t_{1}}[\rho_{0}](M_{\theta_{1}}^{\dagger}\otimes
I_{\mathbb{E}})\equiv\mathcal{V}_{t_{0}:t_{1}}[\rho_{0}]$ (20)

that also provides the formal definition of the map $\mathcal{V}[\rho]$.
Accordingly, the stochastic quantum map of $\mathbb{S}$ at time $t_{1}$, i.e.,
$\widetilde{\Phi}_{1}$, can be implicitly linked to
$\mathcal{V}_{t_{0}:t_{1}}[\rho_{0}]$ through the relation

$\mathcal{M}_{\theta_{1}}\widetilde{\Phi}_{1}[\rho_{\mathbb{S},0}]\equiv{\rm
Tr}_{\mathbb{E}}(\mathcal{V}_{t_{0}:t_{1}}[\rho_{0}]).$ (21)

Analogously, the $k-$th dynamical map of $\mathbb{S}$ after a sequence of
quantum measurements is defined by

$\displaystyle\widehat{\rho}_{k}$ $\displaystyle\sim$
$\displaystyle(M_{\theta_{k}}\otimes
I_{\mathbb{E}})\mathcal{U}_{t_{k-1}:t_{k}}[\widehat{\rho}_{k-1}](M_{\theta_{k}}^{\dagger}\otimes
I_{\mathbb{E}})$ (22) $\displaystyle\equiv$
$\displaystyle\mathcal{V}_{t_{k-1}:t_{k}}[\widehat{\rho}_{k-1}]$

with

$\mathcal{M}_{\theta_{k}}\widetilde{\Phi}_{k}[\rho_{\mathbb{S},0}]={\rm
Tr}_{\mathbb{E}}(\mathcal{V}_{t_{k-1}:t_{k}}[\widehat{\rho}_{k-1}])$ (23)

and $\widehat{\rho}_{0}\equiv\rho_{0}$. For the sake of clarity, we recall
that, by definition, the CP maps $\widetilde{\Phi}_{k}$ are obtained by
tracing the state of the composite system $\mathbb{S}+\mathbb{E}$ (w.r.t. the
environment $\mathbb{E}$) just after applying the evolution superoperator
$\mathcal{U}_{t_{k-1}:t_{k}}$. As done with the definition of transfer tensors
linking together the quantum maps $\Phi_{k}$ of the system at time instants
$t_{k}$, $k=1,\ldots,m$, we can thus introduce the stochastic transformations
$\widetilde{T}^{\underline{\theta},\underline{t}}_{k,j}$ that relate the
conditional dynamics of $\mathbb{S}$ after each measurement of the sequence.
In particular, Eq. (23) has the same structure as the not monitored dynamics,
with $\Phi_{k}$ replaced by $\mathcal{M}_{\theta_{k}}\widetilde{\Phi}_{k}$.
Hence, we define the stochastic TTs via

$\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}=\sum_{j=0}^{k-1}\widetilde{T}^{\underline{\theta},\underline{t}}_{k,j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{j}\,.$
(24)

Note that if there are no measurements, $\mathcal{M}_{\theta_{k}}=\mathbbm{1}$
and $\widetilde{\Phi}^{\underline{\theta},\underline{t}}_{k}=\Phi_{k}$ (see
Eq. (18)) with the result that the TTs
$\widetilde{T}^{\underline{\theta},\underline{t}}_{k,j}$ are no longer
conditional objects and can be identified with the $T_{k,j}$ since Eq. (24)
reduces to Eq. (3). From now on, we will drop the label
$(\cdot)^{\underline{\theta},\underline{t}}$, which is implied in all the
expressions with a tilde.

As before, by expanding Eq. (24) we can write the expression for the $k-$th
dynamical map of $\mathbb{S}$ as a function of $\rho_{\mathbb{S},0}$, i.e.,

$\widetilde{\Phi}_{k}=\underline{\widetilde{T}}_{k}\underline{\mathcal{M}}_{\theta_{k-1}}\underline{\Upsilon}_{k-1}[\rho_{\mathbb{S},0}],$
(25)

where

$\displaystyle\underline{\widetilde{T}}_{k}$ $\displaystyle\equiv$
$\displaystyle\left(\widetilde{T}_{k,0},\widetilde{T}_{k,1},\ldots,\widetilde{T}_{k,k-1}\right)$
$\displaystyle\underline{\mathcal{M}}_{\theta_{k-1}}$ $\displaystyle\equiv$
$\displaystyle{\rm
diag}\left(\mathcal{M}_{\theta_{0}},\mathcal{M}_{\theta_{1}},\ldots,\mathcal{M}_{\theta_{k-1}}\right)$
(26)

and

$\underline{\Upsilon}_{k-1}[\rho_{\mathbb{S},0}]\equiv\begin{pmatrix}\mathbbm{1}_{\mathbb{S}}[\rho_{\mathbb{S},0}]\\\
\widetilde{\Phi}_{1}=\underline{\widetilde{T}}_{1}\underline{\Upsilon}_{0}[\rho_{\mathbb{S},0}]\\\
\widetilde{\Phi}_{2}=\underline{\widetilde{T}}_{2}\underline{\mathcal{M}}_{\theta_{1}}\underline{\Upsilon}_{1}[\rho_{\mathbb{S},0}]\\\
\vdots\\\
\widetilde{\Phi}_{k-1}=\underline{\widetilde{T}}_{k-1}\underline{\mathcal{M}}_{\theta_{k-2}}\underline{\Upsilon}_{k-2}[\rho_{\mathbb{S},0}]\end{pmatrix}$
(27)

with $\underline{\Upsilon}_{-1}$ and $\mathcal{M}_{\theta_{0}}$ equal to the
identity operator $\mathbbm{1}_{\mathbb{S}}$ (the first measurement of the
sequence, indeed, is not performed at time $t_{0}$, but at $t_{1}$) and
$\underline{\Upsilon}_{0}=\underline{\Xi}_{0}=\Phi_{0}$.

For the sake of clarity, we show also in this case the first terms of the
recursive expansion of Eq. (25):

$\begin{cases}\widetilde{\Phi}_{0}=\Phi_{0}=\mathbbm{1}_{\mathbb{S}}[\rho_{\mathbb{S},0}]\\\
\widetilde{\Phi}_{1}=\widetilde{T}_{1,0}\widetilde{\Phi}_{0}=\widetilde{T}_{1,0}[\rho_{\mathbb{S},0}]\\\
\widetilde{\Phi}_{2}=\widetilde{T}_{2,0}\widetilde{\Phi}_{0}+\widetilde{T}_{2,1}\mathcal{M}_{\theta_{1}}\widetilde{\Phi}_{1}=\left(\widetilde{T}_{2,0}+\widetilde{T}_{2,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}\right)[\rho_{\mathbb{S},0}]\\\
\widetilde{\Phi}_{3}=\widetilde{T}_{3,0}\widetilde{\Phi}_{0}+\widetilde{T}_{3,1}\mathcal{M}_{\theta_{1}}\widetilde{\Phi}_{1}+\widetilde{T}_{3,2}\mathcal{M}_{\theta_{2}}\widetilde{\Phi}_{2}\\\
\,\,\,\,\,\,\,\,\,=\left(\widetilde{T}_{3,0}+\widetilde{T}_{3,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}+\widetilde{T}_{3,2}\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,0}\right.\\\
\,\,\,\,\,\,\,\,\left.+\,\widetilde{T}_{3,2}\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}\right)[\rho_{\mathbb{S},0}].\end{cases}$

In addition, let us observe that also the multi-time statistics for
$\mathbb{S}$, given by the joint probability distribution in Eq.(19), can be
expressed by means of these recursive relations. Indeed, by using again Eq.
(25), one has that

$q_{k}={\rm
Tr}\left(\mathcal{M}_{\theta_{k}}\underline{\widetilde{T}}_{k}\underline{\mathcal{M}}_{\theta_{k-1}}\underline{\Upsilon}_{k-1}[\rho_{\mathbb{S},0}]\right),$
(28)

$\forall k=1,\ldots,m$, where the first terms of the recursion for $q_{k}$ are
given by

$\begin{cases}q_{0}={\rm Tr}\left(\rho_{\mathbb{S},0}\right)=1\\\ q_{1}={\rm
Tr}\left(\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}[\rho_{\mathbb{S},0}]\right)\\\
q_{2}={\rm
Tr}\left(\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,0}[\rho_{\mathbb{S},0}]\right)+{\rm
Tr}\left(\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}[\rho_{\mathbb{S},0}]\right)\\\
q_{3}={\rm
Tr}\left(\mathcal{M}_{\theta_{3}}\widetilde{T}_{3,0}[\rho_{\mathbb{S},0}]\right)+{\rm
Tr}\left(\mathcal{M}_{\theta_{3}}\widetilde{T}_{3,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}[\rho_{\mathbb{S},0}]\right)\\\
\,\,\,\,\,\,\,+\,\,{\rm
Tr}\left(\mathcal{M}_{\theta_{3}}\widetilde{T}_{3,2}\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,0}[\rho_{\mathbb{S},0}]\right)\\\
\,\,\,\,\,\,\,+\,\,{\rm
Tr}\left(\mathcal{M}_{\theta_{3}}\widetilde{T}_{3,2}\mathcal{M}_{\theta_{2}}\widetilde{T}_{2,1}\mathcal{M}_{\theta_{1}}\widetilde{T}_{1,0}[\rho_{\mathbb{S},0}]\right).\end{cases}$

Finally, also the set $\\{\widetilde{T}\\}$ of stochastic TTs admits a
hierarchic structure, very similar to that of Eq. (5), based on a
generalization of the CPTP maps $\Gamma$s. The main difference is that now we
need a map acting on the reduced dynamics of $\mathbb{S}$ at each $t_{k}$
after that a quantum measurement has been performed on the system. Thus, such
map will depend on the resulting post-measurement environmental state. In
particular, following the same construction as in [38], we get

$\widetilde{T}_{k,k-n}=\widetilde{\Gamma}_{k|k-n}-\sum_{j=1}^{n-1}\widetilde{T}_{k,k-j}\mathcal{M}_{\theta_{k-j}}\widetilde{\Gamma}_{k-j|k-n},$
(29)

where the CPTP maps

$\widetilde{\Gamma}_{k|k-n}[\sigma_{\mathbb{S}}]={\rm
Tr}_{\mathbb{E}}[\mathcal{U}_{t_{k-n}:t_{k}}[\sigma_{\mathbb{S}}\otimes\widetilde{\sigma}_{\mathbb{E},k-n}]]$
(30)

also include a dependence on the whole sequence of outcomes, fixing the
environmental state $\widetilde{\sigma}_{\mathbb{E},k-n}$ (which motivates the
use of the tilde).

###  IV.2 Conditional transfer tensors and memory effects

The relation in Eq. (29) allows us to bring forward the connection between TTs
and CP-divisibility to the stochastic level that describes the effects
sequential measurements.

In particular, let us consider the situation where, for any sequence of
outcomes, all the stochastic TTs $\widetilde{T}_{k,k-n}$ are equal to 0 for
$n\geq 2$. The stochastic TTs hierarchy of Eq. (29) implies that each $1$-step
transfer tensor $\widetilde{T}_{k,k-1}$ is equal to the corresponding $1$-step
$\widetilde{\Gamma}$, i.e.,
$\widetilde{T}_{k,k-1}=\widetilde{\Gamma}_{k|k-1}$, which by definition is
completely positive. In addition, recursively applying
$\widetilde{\Phi}_{k}=\sum_{j=0}^{k-1}\widetilde{T}_{k,j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}_{j}$,
we get

$\widetilde{\Phi}_{k}=\widetilde{T}_{k|k-1}\mathcal{M}_{\theta_{k-1}}\widetilde{T}_{k-1|k-2}\mathcal{M}_{\theta_{k-2}}\cdots\widetilde{T}_{j+1|j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}_{j}$
(31)

for any $k\geq j\geq 0$. Therefore, we have

$\widetilde{\Phi}_{k}=\widetilde{\Gamma}_{k|k-1}\mathcal{M}_{\theta_{k-1}}\widetilde{\Gamma}_{k-1|k-2}\mathcal{M}_{\theta_{k-2}}\cdots\widetilde{\Gamma}_{j+1|j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}_{j},$
(32)

which implies that for any $k\geq j\geq 0$ there is a (conditional) CP map
$\widetilde{\mathcal{E}}_{k,j}$ such that

$\widetilde{\Phi}_{k}=\widetilde{\mathcal{E}}_{k,j}\widetilde{\Phi}_{j}$

i.e., the family of maps $\widetilde{\Phi}_{k}$ is CP-divisibile [41]. We
conclude that, if all the stochastic TTs $\widetilde{T}_{k,k-n}$ are equal to
0 for $n\geq 2$, the corresponding conditional dynamics
$\\{\widetilde{\Phi}_{k}\\}_{k=1,\ldots,m}$ is CP-divisible for any fixed
sequence of outcomes.

As one can directly check from the definition in Eq. (18), the condition in
Eq. (32) is satisfied whenever the global state after the measurements is a
product state, e.g., if the set $\\{\mathcal{M}_{\theta_{k}}\\}$ describe
projective measurements of a non-degenerate system’s observables. This means
that in these situations the conditional dynamics $\widetilde{\Phi}_{k}$ will
be automatically CP-divisible, simply due to the kind of employed measurement.
As we will see in the next paragraph by means of an explicit example, CP-
divisibility of the conditional dynamics, and in particular the validity of
Eq. (32), does not imply that only one-step conditional TTs are different from
zero.

The stochastic TTs $\widetilde{T}_{k,k-n}$ with $n\geq 2$ enclose the
influence of the possible correlations resulting both from the system-
environment interaction and the measurement on the open system up to a time
step $j$ on the conditional dynamics at the later time step $k$. In general,
such correlations depend on the whole sequence of measurement outcomes
performed up to the $j-$th time instant, so that their influence on the
subsequent conditional dynamics can be read as a signature of memory in the
multi-time statistics defined by Eq. (19). However, even in the case where Eq.
(32) holds, some memory can be present in the conditional dynamics, and hence
in the multi-time statistics. In fact, the post-measurement environmental
state $\widetilde{\sigma}_{\mathbb{E}}$, which defines the conditional maps
$\widetilde{\Gamma}$ in Eq. (30), might (and in general will) depend on the
sequence of outcomes. As we will see in the next paragraph, the proper
description of such memory effects can require non-zero TTs
$\widetilde{T}_{k,k-n}$ with $n\geq 2$.

Finally, it is worth observing that, if in addition to the validity of Eq.
(32) the dependence of $\widetilde{\sigma}_{\mathbb{E}}$ on the previous
outcomes can be neglected (i.e., no memory of the sequence of measurements is
left), Eq. (32) automatically implies the quantum regression theorem [42, 43,
44, 45, 46]. Essentially, the latter means that the whole multi-time
statistics can be reconstructed only from the initial open-system state and
open-system dynamical maps that are not conditioned on the sequence of
outcomes. Similar considerations have been recently discussed also in [15].

###  IV.3 Spin-boson case study

In this paragraph, we show an application of the previously introduced
stochastic transfer tensors, with the aim to illustrate the capability of the
method in giving a quantitative account of the memory effects both in the
conditional dynamics and in the multi-time statistics. In addition, the use of
TTs also allows to compare on a similar footing the memory effects in the
presence of different kinds of measurements, or even no intermediate
measurements at all.

We consider a single spin in interaction with $5$ quantum harmonic oscillators
and monitored by a sequence of $M=5$ quantum measurements at regular intervals
$\Delta$ with $t_{k}=k\Delta$ and $k=1,\ldots,M$. The Hamiltonian $H$ of the
composite system is given by

$H=\frac{1}{2}\sigma^{z}_{\mathbb{S}}+\sum_{k=1}^{5}\omega_{\mathbb{E},k}b_{\mathbb{E},k}^{\dagger}b_{\mathbb{E},k}+\sum_{k=1}^{5}g_{k}(\sigma_{\mathbb{S}}^{-}b_{\mathbb{E},k}^{\dagger}+\sigma_{\mathbb{S}}^{+}b_{\mathbb{E},k})$
(33)

where $\sigma^{z}$ is the Pauli matrix in the $z$ direction,
$b_{\mathbb{E}}^{\dagger}$ and $b_{\mathbb{E}}$ denote respectively the
creation and annihilation operators associated to each harmonic oscillator of
$\mathbb{E}$, and $\sigma_{\mathbb{S}}^{+}$ and $\sigma_{\mathbb{S}}^{-}$ are
the raising and lowering operators of the spin $\mathbb{S}$. Moreover,
regarding the frequencies $\omega_{\mathbb{E},k}$ of the oscillators and the
interaction couplings $g_{k}$, for the numerical simulations here performed we
have chosen the values
$\\{\omega_{\mathbb{E},k}\\}_{k=1}^{5}\approx(1.99,0.73,0.89,2.04,1.58)$ and
$\\{g_{k}\\}_{k=1}^{5}\approx(1.67,1.32,2.15,2.70,1.07)$. All these values are
expressed in units ensuring that $\hbar=1$. The former have been uniformly
sampled from the interval $[0,5]$, while the latter are sampled from the
probability distribution $\alpha g\exp(-g/\beta)$, with $\alpha=1$ and
$\beta=2$. The parameters $\alpha$ and $\beta$ tune the intensity,
respectively, of the interaction between $\mathbb{S}$ and $\mathbb{E}$ and of
the value of $g$ that corresponds to the peak of the distribution. Note that
this model accounts for both decoherence and dissipation in the dynamics of
$\mathbb{S}$ due to the interaction with the external environment, yet with
the form of the interaction preserving the total number of excitations. On the
other hand, the excitations are generally not preserved by the action of the
measurements. In our simulations, each measurement at $t=t_{k}$ is provided by
the following POVM operators [47]:

$F_{\pm}\equiv(1-\lambda)|\pm\rangle\\!\langle\pm|+\frac{\lambda}{2}\mathbbm{1}_{\mathbb{S}}$
(34)

where $|\pm\rangle\equiv(|0\rangle\pm|1\rangle)/\sqrt{2}$ denote the
eigenstates of the Pauli matrix $\sigma^{x}$. This choice is motivated by the
fact that, by changing the parameter $\lambda$, it is possible to analyze the
case without measurements ($\lambda=1$), the one with projective measurements
($\lambda=0$) and all the intermediate cases ($0<\lambda<1$) corresponding to
a partial collapse of the spin wave-function.

The evolution of the composite system has been evaluated starting from an
initial product state of the ground states of both $\mathbb{S}$ and
$\mathbb{E}$, and taking the time interval $\Delta$ between measurements equal
to 1, 2 and 3 (in natural units allowing for $\hbar=1$). Conversely, the
hierarchy of transfer tensors has been derived by numerically implementing Eq.
(29). The latter is based on the computation of the maps
$\widetilde{\Gamma}_{k|k-n}$, which are all obtained by propagating an initial
product state composed by the reduced states of $\mathbb{S}$ and $\mathbb{E}$
after each possible quantum measurement at the time instants $t_{k}$.

Figure 1: Spin-boson dynamics: Ground state population
$\rho_{\mathbb{S}}^{(1,1)}$ (blues solid lines) and coherence terms
$|\rho_{\mathbb{S}}^{(1,2)}|$ (red dash-dotted lines) of the single spin
$\mathbb{S}$ as a function of time and $4$ distinct $\lambda$ values for
$\Delta=1$ in natural units. Depending on the value of $\lambda$ that fixes
the intermediate measurement operator, the coherence and population of the
spin has a different behaviour. In particular, revivals occur for $\lambda=1$,
while a more and more pronounced damping is observed for decreasing values of
$\lambda$. Indeed, sudden changes in the population and coherence terms denote
the action of the measurements, and they are more evident in the case of
projective measurements ($\lambda=0$). Figure 2: Comparison between the
$L_{2}$-norm of the stochastic TTs $\widetilde{T}_{k,0}$ (coloured dots) and
the figure of merit $\mathcal{D}_{k}$ (circles, squares, diamonds and x-marks)
for $k=2,\ldots,5$ (corresponding to the discrete time instants
$t_{2},\ldots,t_{5}$) and $\lambda=0,0.25,0.5,1$ (identified, respectively, by
the colour black, red, blue and magenta). The former quantifier takes into
account the influence of the $k$-step TTs all referring to $t=t_{0}$, while
the latter sums up the impact of all the involved $k-\ell$ steps TTs with
$\ell=0,\ldots,k-2$ and they both account for the non-Markovianity of the
implemented quantum dynamics. If all the stochastic TTs
$\widetilde{T}_{k,k-n}$ are equal to 0 for $n\geq 2$, then both quantifiers
are equal to zero (the opposite does not generally hold).

In Fig. 1 we plot both the population and coherence terms of the single spin
$\mathbb{S}$ as a function of time and $\lambda$ values
($\lambda\in\\{0,0.25,0.5,1\\}$) for $\Delta=1$ (all in natural units). As one
can observe from the figure, depending on the value of $\lambda$ that
corresponds to different intermediate measurement operators, the spin quantum
coherence tends to decrease for smaller values of $\lambda$ or to be
maintained oscillating over time otherwise. In the latter case, repeated
“coherence revivals” are originated by the interaction of $\mathbb{S}$ with
only the environment $\mathbb{E}$ composed of 5 quantum harmonic oscillators.
Instead, in the former case also the external measurement apparatus plays a
crucial role. Hence, from Fig. 1 one can clearly appreciate that monitoring
the open quantum system $\mathbb{S}$ through a sequence of quantum
measurements radically changes its reduced dynamics, as well as the
corresponding non-Markovian behaviors. In particular, if the quantum coherence
of $\mathbb{S}$ is damping over time, then also high-order step memory effects
can be considered negligible. Otherwise, each state of $\mathbb{S}$ also
depends on past contributions, even up to the initial instant in the extreme
case.

Let us now analyze more quantitatively these considerations by introducing
proper quantifiers of non-Markovianity. In this regard, in Fig. 2 we plot two
quantifiers of the impact of the stochastic $n$-step TTs, with $n\geq 2$, on
the conditional dynamics of the system. In each panel, we consider different
POVMs, defined as in Eq. (34) for $\lambda=0,0.25,0.5,1$ and identified
respectively by the colours black, red, blue and magenta, as a function of the
discrete time instants $t_{k}$ with $k=1,\ldots,5$. The three panels
correspond to three different values of the time interval $\Delta$. The two
quantifiers are the $L_{2}$-norm $\|\cdot\|_{2}$ of the stochastic transfer
tensors $\widetilde{T}_{k,0}$ (coloured dots in the figure) and the figure of
merit (circles, squares, diamonds and x-marks)

$\mathcal{D}_{k}\equiv\frac{1}{k-2}\sum_{\ell=0}^{k-2}\|\widetilde{T}_{k,\ell}\mathcal{M}_{\theta_{\ell}}\widetilde{\Phi}_{\ell}\|_{2}.$
(35)

Starting from the decomposition of the conditional map $\widetilde{\Phi}_{k}$
in correspondence of $t=t_{k}$ as given in Eq. (24), one can observe that the
former quantifier focuses on the role of the $k$-step TT for $k\geq 2$ (note,
indeed, that $\widetilde{\Phi}_{0}=\mathbbm{1}$, i.e., no measurement is
performed at the initial time instant), while the latter sums up the influence
of all the involved $k-\ell$ steps TTs with $k\geq 2$ and $\ell=0,\ldots,k-2$
(thus, the $1$-step TT is not considered). If all the stochastic TTs
$\widetilde{T}_{k,k-n}$ are equal to 0 for $n\geq 2$, then both quantifiers
are equal to zero; the opposite does not generally hold.

Overall, from Fig. 2 we can observe that the quantities
$\|\widetilde{T}_{k,0}\|_{2}$ and $\mathcal{D}_{k}$ are comparable; namely, in
the considered parametric regime the $k$-th step TT provides with similar
information on the memory effects influencing $\widetilde{\Phi}_{k}$, compared
to the sum of all the contributions from the TTs with step larger than 1. They
both indicate that non-Markovian behaviours, originated from the interaction
of the spin both with the environment and the observer, become relevant for
$t/\Delta\geq 3$ with all the considered values of $\Delta$. Moreover, both
the non-Markovianity quantifiers reach their maximum values in the case where
no intermediate measurement is performed (magenta dots and x-marks),
corresponding in Fig. 1 to a nearly periodic evolution of the coherence term
$|\rho_{\mathbb{S}}^{(1,2)}|$ over time. Thus, comparing Figs. 1 and 2, we can
deduce that such periodic behaviours are originated by high-order step memory
effects, presumably even from $t_{0}$. While without applying any measurement
both quantifiers typically increase with time, a more pronounced non-monotonic
behavior is observed in the presence of intermediate measurements, which
reduce the memory effects caused by the interaction between $\mathbb{S}$ and
$\mathbb{E}$. Although a decreased impact of multi-step TTs in the presence of
intermediate measurements (corresponding to a decreased impact of system-
environment correlations on the conditional dynamics) might have been
expected, the quantifiers do not necessarily reach their minimum values in the
case of projective measurements ($\lambda=0$). Rather, the minimum values are
reached more often for $\lambda=0.25$ or $\lambda=0.5$ (the behavior for these
two POVMs is quite similar). However, as anticipated, $n$-step TTs for $n\geq
2$ can be non-zero also in the case of projective measurements, i.e., when Eq.
(32) holds. Such TTs account for memory effects due to the dependence of the
post-measurement environmental state on the previous outcomes. Even more, the
latter can exceed the memory effects – due to both system-environment
correlations and changes in the environmental states – occurring in the case
of non-projective measurements, as quantified by $\|\widetilde{T}_{k,0}\|_{2}$
and $\mathcal{D}_{k}$.

Figure 3: $L_{2}$-norm of the difference between the left- and the right-and-
side of equation (32), namely
$\|\widetilde{\Phi}_{k}-\widetilde{\Gamma}_{k|k-1}\mathcal{M}_{\theta_{k-1}}\widetilde{\Gamma}_{k-1|k-2}\mathcal{M}_{\theta_{k-2}}\cdots\widetilde{\Gamma}_{j+1|j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}_{j}\|_{2}$,
computed for $\Delta=2$ and $\lambda=0,0.25,0.5,1$ (magenta dots, blue
x-marks, red diamonds, black squares, respectively) in correspondence of the
discrete time-instants $t_{k}=k\Delta$ with $k\in\\{2,3,4,5\\}$. This
represents the $L_{2}$-norm of the difference between the actual stochastic
map $\widetilde{\Phi}_{k}$ and the CP-divisible map provided by the product
states of the dynamical evolution of the spin after each measurement. During
the transient of the spin dynamics, Eq. (32) is violated as the value of
$\lambda$ increases (thus, for less invasive intermediate measurements);
however, such tendency can be reversed in case the spin approaches an
equilibrium asymptotic state or a stable periodicity phase.

In order to point out the memory effects originated by system-environment
correlations, we consider the $L_{2}$-norm of the difference between the left-
and the right-and-side of Eq. (32), namely between the actual stochastic map
$\widetilde{\Phi}_{k}$ (in correspondence of the $k$-th time instant $t_{k}$)
and the CP-divisible map
$\widetilde{\Gamma}_{k|k-1}\mathcal{M}_{\theta_{k-1}}\widetilde{\Gamma}_{k-1|k-2}\mathcal{M}_{\theta_{k-2}}\cdots\widetilde{\Gamma}_{j+1|j}\mathcal{M}_{\theta_{j}}\widetilde{\Phi}_{j}$
associated with the product states after each measurement. The results are
plotted in Fig. 3. It is worth noting that Eq. (32) is satisfied only when
projective measurements are performed, so that the composite system is
actually in a product state after each measurement and the conditional
dynamics of $\mathbb{S}$ can be fully described by means of the composition of
$1$-step maps $\widetilde{\Gamma}$. In this case, the computation of the
operators $\widetilde{\Gamma}_{k|k-1}$ is enough by itself. However, as said,
some $n$-step stochastic TTs with $n\geq 2$ can be still different from zero,
despite overall their combination is canceled to ensure the validity of Eq.
(32).

As final remark, from Fig. 3 we can also observe that for $k=2,3,4$ the
violation of Eq. (32) increases with the value of $\lambda$, i.e., when moving
toward less invasive measurements; indeed, the largest violation is reached
for the case without intermediate measurements. Moreover, quite interestingly,
the situation is reversed in the correspondence of the last observed value of
time, $t_{k}/\Delta=5$, for which $\lambda=0.25$ induces the largest
violation, while the unconditional dynamics the smallest (apart, of course,
from the case of projective measurements). The latter behaviour of the
unconditional dynamics is usually originated by a tendency towards an
equilibrium asymptotic state or to a stable periodicity phase. In our case
study, being the environment composed by a finite number (five) of quantum
harmonic oscillators, the simulated unconditional dynamics is practically
periodic, with period approximately of $t=10$ (in natural units).

##  V Conclusion and outlook

In this paper, we have studied the non-Markovian dynamics of open quantum
systems and their multi-time statistics by means of the transfer-tensor
approach, which allowed us to treat the memory effects in the two different
scenarios on a similar footing. After showing the connection between the
hierarchy of the TTs and the divisibility of the dynamics, we have extended
the definition of TTs to the case where the open system of interest undergoes
quantum measurements at subsequent instants of time. We have thus introduced a
stochastic family of TTs, depending on the sequence of measurement outcomes,
and a related hierarchy that captures how the multi-time statistics is
influenced by both the system-environment correlations and the dependence of
the current environmental state on the previous outcomes. Finally, we have
defined two quantifiers of multi-time memory effects, which rely directly on
the various contributions to the hierarchy of TTs, and we have investigated
their behavior in a paradigmatic case study, comparing different kinds of
intermediate measurements, as well as the case with no monitoring at all.

The precise relation between the memory effects present, respectively, in non-
Markovian quantum dynamics and in their multi-time statistics remains still to
be addressed. In particular, more quantitative and general results
characterizing different kinds of multi-time statistics, such as those obeying
the quantum-regression theorem, seem to be needed in order to distinguish the
different sources of memory and to compare them with those leading to non-
Markovian quantum dynamics. Hopefully, our results will provide useful
insights to look for such rigorous connections. More in general, our analysis
suggests the usefulness of the TT approach also to deal with multi-time
statistics; an example in this direction might be to use the hierarchy of
stochastic TTs to evaluate the non-Markovianity along the so-called most
probable trajectory of the system, namely the trajectory originated with
higher probability by a sequence of quantum measurements [48, 49]. In
addition, the TT formalism addressed in this paper is also expected to provide
a powerful tool to evaluate experimentally the degree of non-Markovianity [50]
in monitored quantum dynamics [51, 52], as well as to test the validity of the
quantum regression theorem or possible generalizations [53, 54, 11]. One
further advantage of the presented approach would be to track memory effects
also in the reduced state evolution of the environment $\mathbb{E}$. This
feature might be useful in those quantum technology devices (e.g., even the
commercial ones now available online as IBM or Rigetti [55, 56, 57]) in which
the number of the environmental degrees of freedom is not much larger than the
ones of the system $\mathbb{S}$ under analysis. Accordingly, in such a case,
also the state of $\mathbb{E}$ is subjected to frequent and relevant changes
due to the presence of $\mathbb{S}$, and the bi-directional exchange of
information between the system and the environment likely leads to departures
from Markovian dynamics and the quantum regression theorem. Finally, it is
worth also mentioning how stochastic transfer tensors may find application in
modeling the interaction between a quantum system $\mathbb{S}$ and one or more
external thermal baths, giving rise to thermalization processes in the large-
time limit. During the transient of the dynamics, indeed, $\mathbb{S}$ may
exhibit non-Markovian behaviours and non-trivial memory effects. In this
regard, if for example we resort to quantum collision models [58, 59], a
thermal baths is assumed of being composed by a large number of small
subsystems, such that, at discrete time instants, the dynamics between the
quantum system and the bath take place through successive “collisions”
provided by pairwise short interactions. Similar thermalization effects are
also obtained by monitoring $\mathbb{S}$ through a sequence of projective
measurements, in the limit of many measurements [60]. Accordingly, in such a
contexts, the TT formalism is expected to shed light in quantitative terms on
the effective influence of transient non-Markovian dynamics, due to the
interaction with thermal baths, within the quantum system under analysis in
allowing for specific thermalization quantum processes.

#### Note added:

During the completion of this work, the related paper [61] appeared, in which
a generalization of the transfer-tensor formalism to multi-time measurements
is introduced by means of the process-tensor description of the open-system
evolution.

###### Acknowledgements.

S.G. and F.C. were financially supported from PATHOS EU H2020 FET-OPEN grant
No. 828946, the Fondazione CR Firenze through the project Quantum-AI and UNIFI
grant Q-CODYCES.

## References

  * [1] M.M. Wolf, and J.I. Cirac. Dividing Quantum Channels. Commun. Math. Phys. 279, 147-168 (2008).
  * [2] 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).
  * [3] E.-M. Laine, J. Piilo, and H.-P. Breuer. Measure for the non-Markovianity of quantum processes. Phys. Rev. A 81, 062115 (2010).
  * [4] A. Rivas, S.F. Huelga, and M.B. Plenio. Entanglement and non-Markovianity of quantum evolutions. Phys. Rev. Lett. 105, 050403 (2010).
  * [5] B.-H. Liu, L. Li, Y.-F. Huang, C.-F. Li, G.-C. Guo, E.-M. Laine, H.-P. Breuer, and J. Piilo. Experimental control of the transition from Markovian to non-Markovian dynamics of open quantum systems. Nat. Phys. 7, 931-934 (2011).
  * [6] D. Chruscinski, and S. Maniscalco. Degree of Non-Markovianity of Quantum Evolution. Phys. Rev. Lett. 112, 120404 (2014).
  * [7] A. Rivas, S.F. Huelga, and M.B. Plenio. Quantum non-Markovianity: characterization, quantification and detection. Rep. Prog. Phys. 77, 094001 (2014).
  * [8] F. Caruso, V. Giovannetti, C. Lupo, and S. Mancini. Quantum channels and memory effects. Rev. Mod. Phys. 86, 1203 (2014).
  * [9] 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).
  * [10] M.M. Müller, S. Gherardini, and F. Caruso. Stochastic quantum Zeno-based detection of noise correlations. Scientific Reports 6, 38650 (2016).
  * [11] I. de Vega, and D. Alonso. Dynamics of non-Markovian open quantum systems. Rev. Mod. Phys. 89, 015001 (2017).
  * [12] S. Cialdi, M.A.C. Rossi, C. Benedetti, B. Vacchini, D. Tamascelli, S. Olivares, and M.G.A. Paris. All-optical quantum simulator of qubit noisy channels. Appl. Phys. Lett. 110, 081107 (2017).
  * [13] F.A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi. Operational Markov condition for quantum processes. Phys. Rev. Lett. 120, 040405 (2018).
  * [14] F.A. Pollock, C. Rodríguez-Rosario, T. Frauenheim, M. Paternostro, and K. Modi, Non-Markovian quantum processes: Complete framework and efficient characterization. Phys. Rev. A 97, 012127 (2018).
  * [15] A.A. Budini. Quantum Non-Markovian Processes Break Conditional Past-Future Independence. Phys. Rev. Lett. 121, 240401 (2018).
  * [16] S. Campbell, F. Ciccarello, G.M. Palma, and B. Vacchini. System-environment correlations and Markovian embedding of quantum non-Markovian dynamics. Phys. Rev. A 98, 012142 (2018).
  * [17] P. Taranto, F.A. Pollock, S. Milz, M. Tomamichel, and K. Modi. Quantum Markov Order. Phys. Rev. Lett. 122, 140401 (2019).
  * [18] H.-V. Do, C. Lovecchio, I. Mastroserio, N. Fabbri, F.S. Cataliotti, S. Gherardini, M.M. Müller, N. Dalla Pozza, and F. Caruso. Experimental proof of Quantum Zeno-assisted Noise Sensing. New J. Phys. 21, 113056 (2019).
  * [19] S. Campbell, M. Popovic, D. Tamascelli, and B. Vacchini. Precursors of non-Markovianity. New J. Phys. 21, 053036 (2019).
  * [20] S. Cialdi, C. Benedetti, D. Tamascelli, S. Olivares, M. G. A. Paris, and B. Vacchini. Experimental investigation of the effect of classical noise on quantum non-Markovian dynamics. Phys. Rev. A 100, 052104 (2019).
  * [21] S. Gherardini, S. Marcantoni, and F. Caruso. Irreversibility mitigation in unital non-Markovian quantum evolutions. Phys. Rev. Res. 2, 033250 (2020).
  * [22] S. Martina, S. Gherardini, and F. Caruso. Machine learning approach for quantum non-Markovian noise classification. arXiv:2101.03221.
  * [23] L. Li, M.J. W. Hall, H.M. Wiseman. Concepts of quantum non-Markovianity: a hierarchy. Phys. Rep. 759, 1-51 (2018).
  * [24] H. Breuer, and F. Petruccione. The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2002).
  * [25] Á. Rivas, and S.F. Huelga. Open Quantum Systems: An Introduction (Springer, Berlin, 2012).
  * [26] B. Vacchini, A. Smirne, E.-M. Laine, J. Piilo, and H.-P. Breuer. Markovian and non-Markovian dynamics in quantum and classical systems. New J. Phys. 13, 093004 (2011).
  * [27] A. Smirne, A. Stabile, and B. Vacchini. Signatures of non-Markovianity in classical single-time probability distributions. Phys. Scr. T153, 014057 (2013).
  * [28] S. Milz, M. S. Kim, F. A. Pollock, K. Modi. Completely Positive Divisibility Does Not Mean Markovianity. Phys. Rev. Lett. 123, 040401 (2019).
  * [29] A. Smirne, D. Egloff, M.G. Díaz, M.B. Plenio, and S.F. Huelga. Coherence and non-classicality of quantum Markov processes. Quantum Sci. Technol. 4, 01LT01 (2019).
  * [30] P. Strasberg and M.G. Díaz. Classical quantum stochastic processes. Phys. Rev. A. 100, 022120 (2019).
  * [31] S. Milz, F. Sakuldee, F.A. Pollock, and K. Modi. Kolmogorov extension theorem for (quantum) causal modelling and general probabilistic theories. Quantum 4, 255 (2020).
  * [32] A. Smirne, T. Nitsche, D. Egloff, S. Barkhofen, S. De, I. Dhand, C. Silberhorn, S.F. Huelga, and M.B. Plenio. Experimental Control of the Degree of Non-Classicality via Quantum Coherence. Quantum Sci. Technol. 5 04LT01 (2020).
  * [33] M.G. Díaz, B. Desef, M. Rosati, D. Egloff, J. Calsamiglia, A. Smirne, M. Skotiniotis, and S. F. Huelga. Accessible coherence in open quantum system dynamics. Quantum 4, 249 (2020).
  * [34] S. Milz, D. Egloff, P. Taranto, T. Theurer, M.B. Plenio, A. Smirne, S.F. Huelga. When is a non-Markovian quantum process classical? Phys. Rev. X 10, 041049 (2020).
  * [35] S. Gherardini, A. Smirne, M.M. Müller, and F. Caruso. Advances in sequential measurements and control of open quantum systems. Proceedings 12 (1), 11 (2019).
  * [36] J. Cerrillo, and J. Cao. Non-Markovian dynamical maps: Numerical processing of open quantum trajectories. Phys. Rev. Lett. 112, 110401 (2014).
  * [37] R. Rosenbach, J. Cerrillo, S.F. Huelga, J. Cao, and M.B. Plenio. Efficient simulation of non-Markovian system-environment interaction. New J. Phys. 18, 023035 (2016).
  * [38] F.A. Pollock, and K. Modi. Tomographically reconstructed master equations for any open quantum dynamics. Quantum 2, 76 (2018).
  * [39] Y.-Q. Chen, K.-L. Ma, Y.-C. Zheng, J. Allcock, S. Zhang, and C.-Y. Hsieh. Non-Markovian Noise Characterization with the Transfer Tensor Method. Phys. Rev. Applied 13, 034045 (2020).
  * [40] In fact, Eq.(16) implies the property named No Information Backflow in [23], which is sufficient but not necessary for CP-divisibility, so that the latter cannot imply $T_{k,k-n}=0$ for $n\geq 2$.
  * [41] We use here the term CP-divisible referred to a family of maps that are not necessarily trace-preserving. Hence, we also relax the requirement that the intermediate maps $\widetilde{\mathcal{E}}_{k,j}$ are trace preserving.
  * [42] M. Lax. Quantum Noise. XI. Multitime Correspondence between Quantum and Classical Stochastic Processes. Phys. Rev. 172, 350 (1968).
  * [43] N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Amsterdam,1992).
  * [44] H. Carmichael, An Open Systems Approach to Quantum Optics (Springer-Verlag, Berlin, 1993).
  * [45] C. W. Gardiner and P. Zoller, Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics (Springer, Berlin, 2004).
  * [46] G. Guarnieri, A. Smirne, and B. Vacchini. Quantum regression theorem and non-Markovianity of quantum dynamics. Phys. Rev. A 90, 022110 (2014).
  * [47] C. Addis, T. Heinosaari, J. Kiukas, E.-M. Laine, and S. Maniscalco. Dynamics of incompatibility of quantum measurements in open systems. Phys. Rev. A 93, 022114 (2016).
  * [48] S.J. Weber, A. Chantasri, J. Dressel, A.N. Jordan, K.W. Murch, and I. Siddiqi. Mapping the optimal route between two quantum states. Nature (London) 511, 570 (2014).
  * [49] S. Gherardini. Exact nonequilibrium quantum observable statistics: A large-deviation approach. Phys. Rev. A 99 (6), 062105 (2019).
  * [50] C.-F. Li, G.-C. Guo, J. Piilo. Non-Markovian quantum dynamics: What is it good for? EPL 128, 30001 (2019).
  * [51] G.A.L. White, C.D. Hill, F.A. Pollock, L.C.L. Hollenberg, and K. Modi. Demonstration of non-Markovian process characterisation and control on a quantum processor. Nat. Comm. 11, 6301 (2020).
  * [52] L. Xiang, Z. Zong, Z. Zhan, Y. Fei, C. Run, Y. Wu, W. Jin, C. Xiao, Z. Jia, P. Duan, J. Wu, Y. Yin, G. Guo. Quantify the non-Markovian process with intermediate projections in a superconducting processor. arXiv:2105.03333.
  * [53] H.-S. Goan, P.-W. Chen, and C.-C. Jian. Non-Markovian finite-temperature two-time correlation functions of system operators: Beyond the quantum regression theorem. J. Chem. Phys. 12, 124112 (2011).
  * [54] D.P.S. McCutcheon. Optical signatures of non-Markovian behavior in open quantum systems. Phys. Rev. A 93, 22119 (2016).
  * [55] J. Morris, F.A. Pollock, and K. Modi. Non-Markovian memory in IBMQX4. arXiv:1902.07980.
  * [56] G. García-Pérez, M.A.C. Rossi, and S. Maniscalco IBM Q Experience as a versatile experimental testbed for simulating open quantum systems. npj Quantum Inf. 6, 1 (2020).
  * [57] S. Martina, et al. In preparation.
  * [58] F. Ciccarello, S. Lorenzo, V. Giovannetti, and G.M. Palma. Quantum collision models: open system dynamics from repeated interactions. arXiv:2106.11974.
  * [59] O. Arisoy, S. Campbell, and Ö.E. Müstecaplioğlu. Thermalization of finite many-body systems by a collision model. Entropy 21, 1182 (2019).
  * [60] S. Gherardini, G. Giachetti, S. Ruffo, and A. Trombettoni. Thermalization processes induced by quantum monitoring in multi-level systems. arXiv:2012.15216.
  * [61] M. R. Jørgensen and F. A. Pollock. A discrete memory-kernel for multi-time correlations in non-Markovian quantum processes. Phys. Rev. A 102, 052206 (2020).