# A solvable class of non-Markovian quantum multipartite dynamics

Adrián A. Budini  Consejo Nacional de Investigaciones Científicas y Técnicas
(CONICET), Centro Atómico Bariloche, Avenida E. Bustillo Km 9.5, (8400)
Bariloche, Argentina, and Universidad Tecnológica Nacional (UTN-FRBA), Fanny
Newbery 111, (8400) Bariloche, Argentina  Juan P. Garrahan  School of Physics
and Astronomy, University of Nottingham, Nottingham, NG7 2RD, UK  Centre for
the Mathematics and Theoretical Physics of Quantum Non-Equilibrium Systems,
University of Nottingham, Nottingham, NG7 2RD, UK

###### Abstract

We study a class of multipartite open quantum dynamics for systems of
arbitrary number of qubits. The non-Markovian quantum master equation can
involve arbitrary single or multipartite and time-dependent dissipative
coupling mechanisms, expressed in terms of strings of Pauli operators. We
formulate the general constraints that guarantee the complete positivity of
this dynamics. We characterize in detail underlying mechanisms that lead to
memory effects, together with properties of the dynamics encoded in the
associated system rates. We specifically derive multipartite “eternal” non-
Markovian master equations that we term hyperbolic and trigonometric due to
the time dependence of their rates. For these models we identify a transition
between positive and periodically divergent rates. We also study non-Markovian
effects through an operational (measurement-based) memory witness approach.

##  I Introduction

In the theory of open quantum systems, the formulation of quantum Markovian
master equations is completely determined by the theory of quantum semigroups
alicki . In contrast, the study of non-Markovian memory effects presents two
problems. The first one is that the most general structure of a quantum master
equation that captures memory effects, and at the same time is consistent with
the completely positive (CP) condition of the solution map breuerbook ; vega ;
wiseman , is not known. The second one is that different inequivalent memory
witnesses can be used to define and measure non-Markovian effects BreuerReview
; plenioReview .

The first problem has been known for many years. In fact, arbitrary non-
Markovian quantum master equations may lead to unphysical solutions wilkie ;
barnett ; budini ; cresserJD  where the average state (the density matrix)
being not positive definite. For tackling this issue a broad class of
phenomenological and theoretical approaches has been formulated vega , dealing
with both time-convoluted and convolutionless master equations LocalNonLocal .
Examples include the dynamics induced by stochastic Hamiltonians defined by
non-white noises GaussianNoise , phenomenological single memory kernels
shabani ; petruccioneLidarEq ; salo ; kossaDariusz , interaction with
incoherent degrees of freedom lindbladrate ; PostMarkovian ; boltzman ; pekola
; megier ; maximal  and arbitrary ancilla systems swf ; hush , related quantum
collisional models embedding ; collisionVacchini ; ciccarello ;
palmaMultipartito ; strunz ; portugal ; brasilCollisional ; brasil , quantum
generalizations of semi-Markov processes Semi ; andrez , and random unitary
dynamics wudarski ; Polonia , together with some exact derivations from
underlying (microscopic or effective) unitary dynamics TwoQubits ;
exactDecayTLS ; DivergingRatesJCModel ; deltaCorrelated ; ferialdi ;
exactChina ; plenio ; additivity ; smirne .

Despite these advances wilkie ; barnett ; budini ; cresserJD ; LocalNonLocal ;
GaussianNoise ; shabani ; petruccioneLidarEq ; salo ; kossaDariusz ;
lindbladrate ; PostMarkovian ; boltzman ; pekola ; megier ; maximal ;
embedding ; collisionVacchini ; ciccarello ; palmaMultipartito ; strunz ;
portugal ; brasilCollisional ; brasil ; Semi ; andrez ; wudarski ; Polonia ;
TwoQubits ; exactDecayTLS ; DivergingRatesJCModel ; deltaCorrelated ; ferialdi
; exactChina ; plenio ; additivity ; smirne ; swf ; hush  most studies of non-
Markovian evolutions are restricted in general to single or bipartite systems.
In fact, in general checking the CP condition of the dynamics is a non trivial
task, whose difficulty in turn increases with the system’s Hilbert space
dimension. However, quantum information intrinsically requires multipartite
processing, and as a consequence the formulation of multipartite non-Markovian
dynamics is of interest from both theoretical and practical points of view.

Our main goal in this paper is to formulate and study a class of solvable
multipartite non-Markovian master equations. The class of systems we consider
are defined in terms an arbitrary number $N$ of qubits, whose interaction with
the environment can be taken into account through arbitrary Pauli channels.
The evolution of the system’s density matrix $\rho_{t}$ is given by the time-
local master equation $(d/dt)\rho_{t}=\mathcal{L}[\rho_{t}],$ where the
generator of the evolution has the general structure

$\displaystyle\ \mathcal{L}[\bullet]\\!$ $\displaystyle=$
$\displaystyle\\!\\!\\!\sum_{\begin{subarray}{c}i=1,\cdots N\\\
\alpha=x,y,z\end{subarray}}\\!\\!\Gamma_{i}^{\alpha}(t)(\sigma_{i}^{\alpha}\bullet\sigma_{i}^{\alpha}-\bullet)\\!$
$\displaystyle\\!\\!\\!\\!+\\!\\!\\!\\!\sum_{\begin{subarray}{c}i=1,\cdots
N\\\
\alpha,\beta=x,y,z\end{subarray}}\\!\\!\Gamma_{i}^{\alpha\beta}(t)(\sigma_{i}^{\alpha}\sigma_{i+1}^{\beta}\bullet\sigma_{i+1}^{\beta}\sigma_{i}^{\alpha}-\bullet)$
$\displaystyle\\!\\!\\!\\!+\\!\\!\\!\\!\sum_{\begin{subarray}{c}i=1,\cdots
N\\\
\alpha,\beta,\gamma=x,y,z\end{subarray}}\\!\\!\\!\\!\Gamma_{i}^{\alpha\beta\gamma}(t)(\sigma_{i}^{\alpha}\sigma_{i+1}^{\beta}\sigma_{i+2}^{\gamma}\bullet\sigma_{i+2}^{\gamma}\sigma_{i+1}^{\beta}\sigma_{i}^{\alpha}-\bullet)$
$\displaystyle\\!\\!\\!\\!+\cdots.$

Here, $\sigma_{i}^{\alpha}$ is the $\alpha$-th Pauli operator $(\alpha=x,y,z)$
acting on qubit $i$, while $\Gamma_{i}^{\alpha\cdots\beta}(t)$ define local
and bipartite time-dependent (coupling) rates. In general, these rate
functions may take both positive and negative values. The problem is to
characterize which constraints must be fulfilled by them in order to obtain
physically valid solutions. Interestingly, the resolution of this issue leads
us to consider all possible multipartite interaction terms, that is,
decoherence channels that involve coupling between an arbitrary number of
qubits. We also explore which rates emerge when the memory effects arise from
different underlying mechanisms based on coupling with incoherent degrees of
freedom megier ; maximal . The explicit formulation of an operational
(measurement based) memory witness modi ; budiniCPF ; BIF  further provides an
alternative characterization of non-Markovian effects.

As a specific example we study a family of “hyperbolic” and “trigonometric”
eternal multipartite non-Markovian master equations where some rates are
negative or develop divergences at all times, respectively. These cases
provide a non-trivial extension and generalization of previous results valid
for single systems canonicalCresser .

The paper is structured as follows. In Sec. II we present the general class of
multipartite dynamics we consider, characterizing solution of the master
equation, resolving in consequence the constraints that guarantee the CP
condition of the map. General properties are derived for this class of models.
In Sec. III the eternal multipartite dynamics are characterized. In Sec. IV we
study memory effects through an operational memory witness. In Sec. V we
provide our Conclusions. The Appendixes give details of derivations and also
obtain the rates associated to different underlying memory mechanisms.

##  II Multipartite dynamics

The system of interest consists of an arbitrary number $N$ of qubits. For
notational convenience we define a set of Pauli strings
$S_{\mathbf{a}}\equiv\sigma_{a_{1}}\otimes\sigma_{a_{2}}\otimes\sigma_{a_{N}},$
each one associated to the vector $\mathbf{a}=(a_{1},a_{2},\cdots,a_{N}).$
Each component $a_{k}$ $(k=1,2,\cdots N)$ assumes the values
$a_{k}=(0,1,2,3)\leftrightarrow(\mathrm{I},\sigma_{x},\sigma_{y},\sigma_{z}),$
each one being associated to the (two-dimensional) identity matrix and the
standard three Pauli matrices.

The evolution of the system’s density matrix $\rho_{t}$ is written in a local-
in-time way. Arbitrary multipartite decoherence channels are considered,

$\frac{d}{dt}\rho_{t}=\mathcal{L}[\rho_{t}]=\sum_{\mathbf{a\neq
0}}\gamma_{t}^{\mathbf{a}}(S_{\mathbf{a}}\rho_{t}S_{\mathbf{a}}-\rho_{t}).$
(2)

The set of functions $\\{\gamma_{t}^{\mathbf{a}}\\}$ define the rates
associated to the multipartite Pauli channel. In general, there are $4^{N}-1$
different rate functions, as the vector $\mathbf{0}=(0,0,\cdots,0)$ is
associated to the identity operator in the full Hilbert space. Our goal is to
characterize the different aspects of this general evolution. A time-
convoluted formulation of the above dynamics is provided in Appendix A.

###  II.1 Subsystem dynamics

Given the evolution above, we ask about the dynamics of any particular
subsystem. Introducing the splitting
$\mathbf{a}=(\mathbf{a}_{\mathbf{s}},\mathbf{a}_{\mathbf{e}}),$ where
$\mathbf{a}_{\mathbf{s}}$ corresponds to the set of local operators that
define the marginal Pauli string of the subsystem of interest, and
$\mathbf{a}_{\mathbf{e}}$ that of the rest of qubits (now considered as part
of the environment), from Eq. (2) the subsystem density matrix
$\rho_{t}^{\mathbf{s}}=\mathrm{Tr}_{\mathbf{e}}[\rho_{t}]$ (where
$\mathrm{Tr}[\bullet]$ is the trace operation) reads

$\frac{d}{dt}\rho_{t}^{\mathbf{s}}=\sum_{\mathbf{a}_{\mathbf{s}}}\gamma_{t}^{\mathbf{a}_{\mathbf{s}}}(S_{\mathbf{a}_{\mathbf{s}}}\rho_{t}^{\mathbf{s}}S_{\mathbf{a}_{\mathbf{s}}}-\rho_{t}^{\mathbf{s}}),\
\ \ \ \ \
\gamma_{t}^{\mathbf{a}_{\mathbf{s}}}\equiv\sum_{\mathbf{a}_{\mathbf{e}}}\gamma_{t}^{\mathbf{a}_{\mathbf{s}},\mathbf{a}_{\mathbf{e}}}.$
(3)

From this equation we conclude that any subsystem, even when in general is
correlated with the complementary part, has an independent self-evolution. In
addition, the structure of this evolution belongs to the same class as that of
the full system [Eq. (2)]. Consequently, the following results can be
particularized for any subsystem of arbitrary size.

###  II.2 Solution map and completely positive condition

We now show that by using the method of damping bases or spectral
decomposition eigen , the solution map $\rho_{0}\rightarrow\rho_{t}$
corresponding to Eq. (2) can be obtained in an exact way. In order the apply
this technique, first we establish a set of relations fulfilled by the (two
dimensional) Pauli operators. Maintaining the notation
$(\sigma_{0},\sigma_{1},\sigma_{2},\sigma_{3})\leftrightarrow(\mathrm{I},\sigma_{x},\sigma_{y},\sigma_{z}),$
it is easy to check that

$\sigma_{a}\mathrm{Tr}[\sigma_{a}\bullet]=\frac{1}{2}\sum_{b}H_{ab}(\sigma_{b}\bullet\sigma_{b}),$
(4)

where the input $[\bullet]$ is an arbitrary two dimensional operator and
$b=0,1,2,3.$ The inverse relation reads

$\sigma_{a}\bullet\sigma_{a}=\frac{1}{2}\sum_{b}H_{ab}\
\sigma_{b}\mathrm{Tr}[\sigma_{b}\bullet].$ (5)

In these expressions, the coefficients $\\{H_{ab}\\}$ define a four
dimensional Hadamard matrix $H,$ which reads

$H\equiv\left(\begin{array}[]{cccc}1&1&1&1\\\ 1&1&-1&-1\\\ 1&-1&1&-1\\\
1&-1&-1&1\end{array}\right).$ (6)

In deriving Eq. (5), we used that its inverse reads $H^{-1}=H/4.$ Also notice
that $H=H^{T}.$

Now, we introduce an extra rate $\gamma_{t}^{\mathbf{0}},$ which is associated
to the identity string in the full Hilbert space

$\gamma_{t}^{\mathbf{0}}\equiv-\sum_{\mathbf{a\neq
0}}\gamma_{t}^{\mathbf{a}}.$ (7)

With this definition, the Lindbladian-like structure of Eq. (2)] can
straightforwardly be written as

$\mathcal{L}[\bullet]=\sum_{\mathbf{a}}\gamma_{t}^{\mathbf{a}}(S_{\mathbf{a}}\bullet
S_{\mathbf{a}}),$ (8)

where the sum now includes the (identity) string $\mathbf{a=0}$. Written in
this way, applying the “vectorial extension” of Eq. (5) to the Hilbert space
of $N$ qubits, it follows that

$\mathcal{L}[\bullet]=\frac{1}{2^{N}}\sum_{\mathbf{a}}S_{\mathbf{a}}\mathrm{Tr}[S_{\mathbf{a}}\bullet]\sum_{\mathbf{b}}H_{\mathbf{ab}}\gamma_{t}^{\mathbf{b}},$
(9)

where $H_{\mathbf{ab}}\equiv H_{a_{1}b_{1}}H_{a_{2}b_{2}}\cdots
H_{a_{N}b_{N}}$ can be read as the matrix elements of the external product of
$N$ single Hadamard matrices, cf. Eq. (6). From this last expression, by using
that
$\mathrm{Tr}[S_{\mathbf{a}}S_{\mathbf{b}}]=2^{N}\delta_{\mathbf{a},\mathbf{b}},$
it is straightforward to determine the eigenvalues and eigenoperators of
$\mathcal{L}[\bullet].$ They read

$\mathcal{L}[S_{\mathbf{a}}]=\mu_{t}^{\mathbf{a}}S_{\mathbf{a}},\ \ \ \ \ \ \
\ \ \
\mu_{t}^{\mathbf{a}}=\sum_{\mathbf{b}}H_{\mathbf{ab}}\gamma_{t}^{\mathbf{b}}.$
(10)

Consequently, any Pauli string $S_{\mathbf{a}}$ is a right eigenoperator with
eigenvalue $\mu_{t}^{\mathbf{a}}.$ Given that $\mathcal{L}[\bullet]$ also
defines the adjoint evolution (as the “jump operators” are Hermitian) eigen ,
$S_{\mathbf{a}}$ is also a left eigenoperator. Notice also that by using the
inverse of the Hadamard matrix, the inverse relation
$\gamma_{t}^{\mathbf{a}}=\sum_{\mathbf{b}}H_{\mathbf{ab}}\mu_{t}^{\mathbf{b}}/4^{N}$
follows.

From the method of damping bases eigen , Eq. (10) allows us to write the
solution of Eq. (2) as

$\rho_{t}=\frac{1}{2^{N}}\sum_{\mathbf{a}}\exp\left[\int_{0}^{t}dt^{\prime}\mu_{t^{\prime}}^{\mathbf{a}}\right]S_{\mathbf{a}}\mathrm{Tr}[S_{\mathbf{a}}\rho_{0}].$
(11)

One can see that the conditions
$\mathrm{Tr}[\rho_{t}]=\mathrm{Tr}[\rho_{0}]=1$ are satisfied after noting
that $\mathrm{Tr}[S_{\mathbf{a}}]=2^{N}\delta_{\mathbf{a},\mathbf{0}}$ and
$\mu_{t^{\prime}}^{\mathbf{0}}=0.$ This last equality follows from Eqs. (7)
and (10) jointly with the property $H_{\mathbf{0b}}=1$ $\forall\mathbf{b.}$ By
using the vectorial extension of Eq. (4), we get the density matrix written in
a Kraus representation alicki ; breuerbook ,

$\rho_{t}=\sum_{\mathbf{a}}p_{t}^{\mathbf{a}}(S_{\mathbf{a}}\rho_{0}S_{\mathbf{a}}).$
(12)

The weights are
$p_{t}^{\mathbf{a}}=4^{-N}\sum_{\mathbf{b}}H_{\mathbf{ab}}\exp[\int_{0}^{t}dt^{\prime}\mu_{t^{\prime}}^{\mathbf{b}}],$
which from Eq. (10) can explicitly be written in terms of the time-dependent
rates as

$p_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\exp\left[\sum_{\mathbf{c}}H_{\mathbf{bc}}\int_{0}^{t}dt^{\prime}\gamma_{t^{\prime}}^{\mathbf{c}}\right].$
(13)

The final expressions (12) and (13) are the main results of this section. They
completely characterize the solution map in terms of the set of rates
$\\{\gamma_{t}^{\mathbf{a}}\\}$ and the initial condition $\rho_{0}.$ In
addition, they naturally provide a constraint that the rates must to fulfill
in order to obtain a CP map, that is, one that gives physical solution. In
fact, the Kraus representation theorem alicki ; breuerbook  implies the
conditions $0\leq p_{t}^{\mathbf{a}}\leq 1,$ which means that
$\\{p_{t}^{\mathbf{a}}\\}$ are a set of normalized probabilities. In the
single qubit case $(N=1),$ previously obtained constraints are recovered
wudarski . In the general case, $4^{N}$ inequalities must be fulfilled. We
notice that a sufficient, but not necessary, condition is
$\int_{0}^{t}dt^{\prime}\gamma_{t^{\prime}}^{\mathbf{a}}\geq 0$
$\forall\mathbf{a\neq 0.}$ In fact, this constraint implies that all
eigenvalues, cf. Eq. (10), satisfy $\mu_{t}^{\mathbf{a}}\leq 0$
$(\mathbf{a}\neq\mathbf{0}).$ Consequently, taking an arbitrary but fixed time
$t,$ the solution (11) of the non-Markovian dynamics, via the association
$\int_{0}^{t}dt^{\prime}\mu_{t^{\prime}}^{\mathbf{a}}=t\mu_{M}^{\mathbf{a}},$
is equivalent to the solution of a (well behaved) Markovian dynamics generated
by a Lindbladian with eigenvalues $\\{\mu_{M}^{\mathbf{a}}\\}.$

###  II.3 Non-Markovianity and time-dependent rates

Different (inequivalent) memory witnesses based only on the system propagator
can be used to define non-Markovianity BreuerReview ; plenioReview  such as
for example the trace distance between two different initial conditions
BreuerFirst  or those based on the $k$-positivity of the solution map
DarioSabrina . Here, as the dynamics is written naturally in a canonical form
canonicalCresser , memory effects can also be defined by the negativity of the
time-dependent rates $\\{\gamma_{t}^{\mathbf{a}}\\}.$ In this way, it is of
interest to determine these elements for any well behaved solution defined by
the probabilities $\\{p_{t}^{\mathbf{a}}\\}$ in Eq. (12).

We can invert Eq. (13),

$\mu_{t}^{\mathbf{a}}=\frac{d}{dt}\ln\left[\sum_{\mathbf{b}}H_{\mathbf{ab}}p_{t}^{\mathbf{b}}\right],$
(14)

and using Eq. (10) we get explicit expressions for the set of rates
$\\{\gamma_{t}^{\mathbf{a}}\\}$ in terms of the normalized time-dependent
weights $0\leq p_{t}^{\mathbf{c}}\leq 1$,

$\gamma_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\frac{d}{dt}\ln\left[\sum_{\mathbf{c}}H_{\mathbf{bc}}p_{t}^{\mathbf{c}}\right].$
(15)

The signs of $\\{\gamma_{t}^{\mathbf{a}}\\}$ can be taken as a signature of
departure from a Markovian regime canonicalCresser . Alternatively, in Sec. V
we study operational measures for non-Markovianity. We notice that Eqs. (13)
and (15) provide a multipartite generalization of the case $N=1$ studied in
Ref. wudarski .

###  II.4 Additivity of non-Markovian master equations

Given two sets of (arbitrary) normalized probabilities
$\\{p_{t}^{\mathbf{a}}\\}$ and $\\{\tilde{p}_{t}^{\mathbf{a}}\\},$ the
relation (15) allows us to obtain the corresponding sets of rates
$\\{\gamma_{t}^{\mathbf{a}}\\}$ and $\\{\tilde{\gamma}_{t}^{\mathbf{a}}\\}.$
From these we can obtain a new master equation defined by Eq. (2) with rates
$\\{\gamma_{t}^{\mathbf{a}}+\tilde{\gamma}_{t}^{\mathbf{a}}\\}$. In fact, it
is always possible to associate a set of probabilities
$\\{q_{t}^{\mathbf{a}}\\}$ to these added rates, that is,

$\\{p_{t}^{\mathbf{a}}\\}\leftrightarrow\\{\gamma_{t}^{\mathbf{a}}\\},\ \ \
\\{\tilde{p}_{t}^{\mathbf{a}}\\}\leftrightarrow\\{\tilde{\gamma}_{t}^{\mathbf{a}}\\},\
\ \Rightarrow\ \
\exists\\{q_{t}^{\mathbf{a}}\\}\leftrightarrow\\{\gamma_{t}^{\mathbf{a}}+\tilde{\gamma}_{t}^{\mathbf{a}}\\}.$
(16)

Consequently, as occurs to Markovian Lindblad equations breuerbook , for our
class of models arbitrary well behaved evolutions (defined by a given set of
rates) can be added in an arbitrary way. The validity of this result follows
from the commutation of two arbitrary propagators, Eq. (12), a property
supported by the relation

$S_{\mathbf{a}}S_{\mathbf{b}}\bullet
S_{\mathbf{b}}S_{\mathbf{a}}=S_{\mathbf{b}}S_{\mathbf{a}}\bullet
S_{\mathbf{a}}S_{\mathbf{b}}=S_{\mathbf{c}}\bullet
S_{\mathbf{c}}^{{\dagger}},$ (17)

which is valid for arbitrary Pauli strings $S_{\mathbf{a}}$ and
$S_{\mathbf{b}},$ and where $S_{\mathbf{c}}=S_{\mathbf{a}}S_{\mathbf{b}}$ or
equivalently $S_{\mathbf{c}}=S_{\mathbf{b}}S_{\mathbf{a}}.$ Eq. (17) can be
straightforwardly demonstrated from Eq. (5).

###  II.5 Coupling with incoherent degrees of freedom

Memory effects are induced whenever extra degrees of freedom are traced out.
Here, we consider a general coupling with incoherent degrees of freedom. Based
on Ref. lindbladrate , the more general case can always be described by
writing the system density matrix $\rho_{t}$ and the probabilities of the
incoherent system $\\{q_{t}^{\mathbf{h}}\\}$ as

$\rho_{t}=\sum_{\mathbf{h}}\rho_{t}^{\mathbf{h}},\ \ \ \ \ \ \ \ \ \ \
q_{t}^{\mathbf{h}}=\mathrm{Tr}[\rho_{t}^{\mathbf{h}}],$ (18)

where the auxiliary states $\\{\rho_{t}^{\mathbf{h}}\\}$ correspond to the
system state given that the extra (hidden) incoherent degrees of freedom are
in the particular state $\mathbf{h.}$ The evolution of the states
$\\{\rho_{t}^{\mathbf{h}}\\}$ may involve coupling between all of them
lindbladrate . Given the structure Eq. (2), each auxiliary state
$\rho_{t}^{\mathbf{h}}$ must to assume the form

$\rho_{t}^{\mathbf{h}}=\sum_{\mathbf{\alpha}}g_{\mathbf{\alpha}}^{\mathbf{h}}(t)(S_{\mathbf{\alpha}}\rho_{0}S_{\mathbf{\alpha}}),$
(19)

where the parameter $\mathbf{\alpha}$ runs over a set of Pauli strings that
depends on each specific problem. The functions
$g_{\mathbf{\alpha}}^{\mathbf{h}}(t)$ in turn obey a classical master equation
whose structure also depends on each specific model. The initial conditions
read $\rho_{0}^{\mathbf{h}}=\rho_{0}q_{0}^{\mathbf{h}},$ where $\rho_{0}$ is
the initial system state and $q_{0}^{\mathbf{h}}$ is the initial probability
of the incoherent degrees of freedom. In fact, at time $t,$
$q_{t}^{\mathbf{h}}=\sum_{\mathbf{\alpha}}g_{\mathbf{\alpha}}^{\mathbf{h}}(t).$
On the other hand, the system density matrix evolution [Eq. (12)] is defined
by the probabilities
$p_{t}^{\mathbf{\alpha}}=\sum_{\mathbf{h}}g_{\mathbf{\alpha}}^{\mathbf{h}}(t).$
A general treatment is not possible. Relevant examples are worked out in
Appendix B such as a mapping with a classical Markovian master equation,
stochastic Hamiltonians, and statistical mixtures of Markovian evolutions. In
all cases, explicit expressions for the rates [Eq. (15)] can be obtained. A
representative class of dynamics is studied in the next section.

##  III Multipartite eternal non-Markovianity

For a single qubit, $N=1$, the system density matrix evolution, Eq. (2), may
involve rates that are negative at all times. This property was called
“eternal non-Markovianity” canonicalCresser ; megier . The results of Appendix
B [see Eqs. (49), (54), and (60)] and Appendix C [see Eqs. (66) and (68)]
guarantee that this property also emerges in multipartite dynamics, $N>1$,
which have $4^{N}-1$ rates.

In order to provide simple (multipartite) examples, here we restrict to the
case where the evolution is

$\displaystyle\mathcal{L}[\bullet]$ $\displaystyle=$
$\displaystyle\Big{\\{}\gamma_{t}^{\underline{\mathbf{a}}}(S_{\underline{\mathbf{a}}}\bullet
S_{\underline{\mathbf{a}}}-\bullet)+\gamma_{t}^{\underline{\mathbf{b}}}(S_{\underline{\mathbf{b}}}\bullet
S_{\underline{\mathbf{b}}}-\bullet)$ (20)
$\displaystyle+\gamma_{t}^{\underline{\mathbf{c}}}(S_{\underline{\mathbf{c}}}\bullet
S_{\underline{\mathbf{c}}}^{\dagger}-\bullet)\Big{\\}},$

where $S_{\underline{\mathbf{a}}}\mathbf{\ }$and $S_{\underline{\mathbf{b}}}$
are two arbitrary multipartite Pauli strings, while
$S_{\underline{\mathbf{c}}}=S_{\underline{\mathbf{a}}}S_{\underline{\mathbf{b}}}.$
Depending on the time-dependence of the rates we define what we term
“hyperbolic” and “trigonometric” cases of eternal non-Markovianity.

###  III.1 Hyperbolic eternal non-Markovianity

The system density matrix is written as the addition of two auxiliary states
$\rho_{t}=\rho_{t}^{(1)}+\rho_{t}^{(2)}$ [Eq. (18)], whose evolution reads

$\displaystyle\frac{d\rho_{t}^{(1)}}{dt}$ $\displaystyle=$
$\displaystyle-\gamma\rho_{t}^{(1)}+\gamma
S_{\underline{\mathbf{a}}}\rho_{t}^{(1)}S_{\underline{\mathbf{a}}},$ (21a)
$\displaystyle\frac{d\rho_{t}^{(2)}}{dt}$ $\displaystyle=$
$\displaystyle-\varphi\rho_{t}^{(2)}+\varphi
S_{\underline{\mathbf{b}}}\rho_{t}^{(2)}S_{\underline{\mathbf{b}}}.$ (21b) The
initial conditions for the auxiliary states are taken to be
$\rho_{0}^{(1)}=\rho_{0}^{(2)}=\rho_{0}/2.$ Given that the auxiliary states do
not couples, the rates of the non-Markovian evolution follow Eq. (15) with
probabilities
$p_{t}^{\mathbf{a}}=[p_{1}^{\mathbf{a}}(t)+p_{2}^{\mathbf{a}}(t)]/2$, with the
sets $\\{p_{1}^{\mathbf{a}}(t)\\}$ and $\\{p_{2}^{\mathbf{a}}(t)\\},$ via Eq.
(13), associated to $\rho_{t}^{(1)}$ and $\rho_{t}^{(2)}$, respectively.
Taking $\varphi=\gamma,$ we get [see also derivation from Eq. (60) in Appendix
C]
$\gamma_{t}^{\underline{\mathbf{a}}}=\gamma_{t}^{\underline{\mathbf{b}}}=\frac{1}{2}\gamma,\
\ \ \ \ \ \ \gamma_{t}^{\underline{\mathbf{c}}}=-\frac{1}{2}\gamma\tanh(\gamma
t).$ (22)

This result provides a multipartite generalization, $(N>1)$, of the single
qubit case $(N=1)$ studied in Ref. canonicalCresser . Similarly to the results
of Ref. megier  we notice that in this particular case alternative dynamics
such as the mapping to a classical master equation [see Eq. (49)] and
stochastic Hamiltonians [see Eq. (54)] also lead to the same rates.

###  III.2 Trigonometric eternal non-Markovianity

Based on Eq. (18), instead of the evolution (21), here we consider

$\displaystyle\frac{d\rho_{t}^{(1)}}{dt}$ $\displaystyle=$
$\displaystyle-\gamma\rho_{t}^{(1)}+\varphi
S_{\underline{\mathbf{b}}}\rho_{t}^{(2)}S_{\underline{\mathbf{b}}},$ (23a)
$\displaystyle\frac{d\rho_{t}^{(2)}}{dt}$ $\displaystyle=$
$\displaystyle-\varphi\rho_{t}^{(2)}+\gamma
S_{\underline{\mathbf{a}}}\rho_{t}^{(1)}S_{\underline{\mathbf{a}}}.$ (23b) The
initial conditions are taken as
$\rho_{0}^{(1)}=[\varphi/(\varphi+\gamma)]\rho_{0}$ and
$\rho_{0}^{(2)}=[\gamma/(\varphi+\gamma)]\rho_{0},$ where $\rho_{0}$ is the
system initial state. Notice that the incoherent transitions
$(1)\leftrightarrow(2)$ imply the system transformations $\rho\rightarrow
S_{\underline{\mathbf{a}}/\underline{\mathbf{b}}}\rho
S_{\underline{\mathbf{a}}/\underline{\mathbf{b}}}$.

Taking into account Eq. (19), in order to solve Eq. (23) each auxiliary state
is written as $(h=1,2)$

$\rho_{t}^{(h)}=g_{\mathbf{0}}^{(h)}\rho_{0}+g_{\underline{\mathbf{a}}}^{(h)}S_{\underline{\mathbf{a}}}\rho_{0}S_{\underline{\mathbf{a}}}+g_{\underline{\mathbf{b}}}^{(h)}S_{\underline{\mathbf{b}}}\rho_{0}S_{\underline{\mathbf{b}}}+g_{\underline{\mathbf{c}}}^{(h)}S_{\underline{\mathbf{c}}}\rho_{0}S_{\underline{\mathbf{c}}}^{\dagger},$
(24)

where as before
$S_{\underline{\mathbf{c}}}=S_{\underline{\mathbf{a}}}S_{\underline{\mathbf{b}}}$,
and $g_{\mathbf{\alpha}}^{(h)}$ are time-dependent functions. Using Eq. (17),
it is possible to derive a classical master equation for the (eight)
$g$-functions which involves coupling between pairs of them. The corresponding
solutions allow to obtain the probabilities
$p_{t}^{\mathbf{a}}=\sum_{h}g_{\mathbf{a}}^{(h)}(t).$ Finally, the rates
associated to the non-Markovian evolution follow from Eq. (15)

$\gamma_{t}^{\underline{\mathbf{a}}}=\gamma_{t}^{\underline{\mathbf{b}}}=\frac{\varphi\gamma(\varphi+\gamma)}{e^{t(\varphi+\gamma)}(\varphi-\gamma)^{2}+4\varphi\gamma}.$
(25)

Furthermore,

$\displaystyle\gamma_{t}^{\underline{\mathbf{c}}}$ $\displaystyle=$
$\displaystyle\varphi\gamma\\{\delta_{+}\Upsilon^{2}(1-e^{t\Upsilon})-\delta_{-}^{2}[\Upsilon(1+e^{t\Upsilon})$
(26)
$\displaystyle+e^{t\delta_{+}}[(\delta_{+}-\Upsilon)-e^{t\Upsilon}(\delta_{+}+\Upsilon)]]\\}$
$\displaystyle\times\\{(e^{t\delta_{+}}\delta_{-}^{2}+4\varphi\gamma)[(1+e^{t\Upsilon})\Upsilon\delta_{+}$
$\displaystyle-(1-e^{t\Upsilon})\delta_{-}^{2}]\\}^{-1},$

where the coefficients are

$\Upsilon\equiv(\varphi^{2}-6\varphi\gamma+\gamma^{2})^{1/2},\ \ \ \ \ \ \
\delta_{\pm}\equiv\varphi\pm\gamma.$ (27)

Depending on the ratio $\varphi/\gamma,$ different characteristic behaviors
are obtained. In Fig. 1 we plot both rates. Consistent with Eq. (25),
$\gamma_{t}^{\underline{\mathbf{a}}}$ and
$\gamma_{t}^{\underline{\mathbf{b}}}$ are always positive functions. However,
this is not the case for $\gamma_{t}^{\underline{\mathbf{c}}},$ Eq. (26),
which depending on $\varphi/\gamma$ develops a transition between positivity
[Figs. 5(a) and 5(d)] and a periodic divergent behavior [Figs. 5(b) and 5(c)].
From Eq. (27) we deduce that this change occurs in the boundaries of the
interval $3-\sqrt{8}<(\varphi/\gamma)<3+\sqrt{8}$, with
$\gamma_{t}^{\underline{\mathbf{c}}}$ developing divergences in this interval,
while being positive outside it.

From the plots it is also evident that $\gamma_{t}^{\underline{\mathbf{a}}}$
and $\gamma_{t}^{\underline{\mathbf{b}}}$ approach a constant when
$\varphi\approx\gamma.$ In fact, when $\varphi=\gamma,$ the previous
expressions reduce to

$\gamma_{t}^{\underline{\mathbf{a}}}=\gamma_{t}^{\underline{\mathbf{b}}}=\frac{1}{2}\gamma,\
\ \ \ \ \ \ \gamma_{t}^{\underline{\mathbf{c}}}=\frac{1}{2}\gamma\tan(\gamma
t).$ (28)

Based on Eq. (22), we name this case as a _trigonometric eternal non-
Markovian_. The probabilities $\\{p_{t}^{\mathbf{a}}\\}$ [Eq. (12)] also
assume a simple form,

$\displaystyle p_{t}^{\mathbf{0}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-\gamma t}[\cosh(\gamma t)+\cos(\gamma t)],$ (29a)
$\displaystyle p_{t}^{\underline{\mathbf{a}}}$ $\displaystyle=$ $\displaystyle
p_{t}^{\underline{\mathbf{b}}}=\frac{1}{4}[1-e^{-2\gamma t}],$ (29b)
$\displaystyle p_{t}^{\underline{\mathbf{c}}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-\gamma t}[\cosh(\gamma t)-\cos(\gamma t)].$ (29c)
These solutions apply to arbitrary multipartite Pauli strings
$\underline{\mathbf{a}}$ and $\underline{\mathbf{b}}.$

Figure 1: Time dependent rates [Eqs. (25) and (26)] corresponding to the
multipartite trigonometric eternal non-Markovian evolution [Eq. (23)] for
different values of the rate ratio $\varphi/\gamma.$

###  III.3 Adding non-Markovian evolutions

Added to the previous examples (see also Appendix C), the possibility of
adding arbitrary (well defined) rates [Eq. (16)] gives us a procedure for
constructing a large family of well behaved dynamics. For example, we write

$\displaystyle\mathcal{L}[\bullet]$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{N}\frac{\gamma_{i}}{2}\Big{\\{}(\sigma_{i}^{x}\sigma_{i+1}^{x}\bullet\sigma_{i+1}^{x}\sigma_{i}^{x}-\bullet)$
$\displaystyle\ \ \ \ \ \ \
+(\sigma_{i}^{y}\sigma_{i+1}^{y}\bullet\sigma_{i+1}^{y}\sigma_{i}^{y}-\bullet)$
$\displaystyle\ \ \ \ \ \ \
+f_{i}(t)(\sigma_{i}^{z}\sigma_{i+1}^{z}\bullet\sigma_{i+1}^{z}\sigma_{i}^{z}-\bullet)\Big{\\}}.$

In this traslational invariant generator (say with periodic boundaries in one
dimension), we may chose $f_{i}(t)=-\tanh(\gamma_{i}t)$ or alternatively
$f_{i}(t)=\tan(\gamma_{i}t)$ [see Eqs. (22) and (28) respectively].

One interesting aspect of using additivity for constructing multipartite
evolutions is that, even when the underlying evolutions have a clear memory
mechanism (see also Appendix B), the resulting dynamics does not necessarily.
For example, while our approach guarantees that Eq. (III.3) leads to a
completely positive dynamics [with solution defined by Eqs. (12) and (13)] it
is not evident which underlying processes may lead to this master equation. In
addition, in general there may be subsystems that are coupled between then,
one part being Markovian and the other non-Markovian. For example, take
$f_{i}(t)=\gamma_{i}/2$ for $i\leq N_{0},$ and $f_{i}(t)=\tan(\gamma_{i}t)$
for $i>N_{0}.$

##  IV Operational memory witness

An alternative and deeper characterization of quantum non-Markovianity can be
obtained by defining memory effects via measurement based approaches modi ;
budiniCPF ; BIF . Here, we study a conditional past-future (CPF) correlation
budiniCPF . This object relies on performing three successive measurement of
arbitrary system observables and calculating the correlation between the last
(future) and first (past) outcomes conditioned to a given intermediate
(present) outcome. For Markovian dynamics it vanishes identically, while
memory effects leads to a non null CPF correlation.

The measurements, denoted in successive order by $\underline{\mathbf{x}},$
$\underline{\mathbf{y}},$ and $\underline{\mathbf{z}},$ correspond to
observations of three Hermitian operators $S_{\underline{\mathbf{m}}}$ with
eigenvectors $\\{|m\rangle\\}$ and eigenvalues $\\{m\\},$

$S_{\underline{\mathbf{m}}}|m\rangle=m|m\rangle,\ \ \ \ \
\underline{\mathbf{m}}=\underline{\mathbf{x}},\underline{\mathbf{y}},\underline{\mathbf{z}}.$
(31)

The CPF correlation then reads budiniCPF

$C_{pf}(t,\tau)|_{y}=\sum_{z,x}zx[P(z,x|y)-P(z|y)P(x|y)],$ (32)

where $\\{x\\},$ $\\{y\\},$ and $\\{z\\}$ denotes the three sets of successive
outcomes (operators eigenvalues), while $t$ and $\tau$ are the (first and
second) time intervals between the successive measurements. With $P(u|v)$ we
denote the conditional probability of $u$ given $v.$

All probabilities appearing in Eq. (32) can be determine from the (outcomes)
joint probability $P(z,y,x)\leftrightarrow P(z,t+\tau,y,t;x,0),$ which in turn
can be calculated after knowing the underlying system-environment dynamics. In
Appendix D we show that $P(z,y,x)$ and $C_{pf}(t,\tau)|_{y}$ can be calculated
exactly assuming that memory effects emerge due to the coupling with
incoherent degrees of freedom [Eqs. (18) and (19)].

Each specific model [see examples (21) and (23)] is completely defined by the
set of functions $\\{g_{\mathbf{\alpha}}^{\mathbf{h}}(t)\\}$ [Eq. (19)]. Given
that they obey a (linear) classical master equation, they can be written as

$g_{\mathbf{\alpha}}^{\mathbf{h}}(t)=\sum_{\mathbf{h}^{\prime}}f_{\mathbf{\alpha}}^{\mathbf{hh}^{\prime}}(t)q_{0}^{\mathbf{h}^{\prime}}\equiv(\mathbf{h}|\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0}),$
(33)

where the set of functions
$\\{f_{\mathbf{\alpha}}^{\mathbf{hh}^{\prime}}(t)\\}$ are independent of the
initial conditions $\\{q_{0}^{\mathbf{h}}\\}.$ Furthermore, for notational
simplicity, we introduced a vectorial orthogonal base $\\{|\mathbf{h})\\}$ for
the incoherent degrees of freedom, such that
$f_{\mathbf{\alpha}}^{\mathbf{hh}^{\prime}}(t)\leftrightarrow(\mathbf{h}|\mathbb{F}_{\mathbf{\alpha}}(t)|\mathbf{h}^{\prime})$
and $q_{0}^{\mathbf{h}}\leftrightarrow(\mathbf{h}|q_{0}).$

The observables $S_{\underline{\mathbf{m}}}$ [Eq. (31)] may in principle be
defined by arbitrary linear combinations of Pauli strings
$\\{S_{\mathbf{a}}\\}.$ Here, for simplicity they are defined by a unique
Pauli string. In this case, the general expression for the CPF correlation
[Eq. (73)] reduces to (see Appendix D)

$C_{pf}(t,\tau)|_{y}=\delta_{\underline{\mathbf{z}},\mathbf{y}}\delta_{\underline{\mathbf{y}},\underline{\mathbf{x}}}\frac{(1-\langle
x\rangle^{2})}{[2^{N}P(y)]^{2}}\sum_{\mathbf{\alpha,\beta}}H_{\underline{\mathbf{y}}\mathbf{\alpha}}H_{\underline{\mathbf{y}}\mathbf{\beta}}[(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)\mathbb{F}_{\mathbf{\beta}}(t)|q_{0})-(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)|q_{t})(1|\mathbb{F}_{\mathbf{\beta}}(t)|q_{0})].$
(34)

In here,
$|q_{t})=\sum_{\mathbf{\alpha}}\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0})$ define
the probabilities of the incoherent degrees of freedom at time $t,$ while
$(1|\equiv\sum_{\mathbf{h}}(\mathbf{h}|.$ Furthermore, $\langle
x\rangle\equiv\sum_{x}xP(x)$ where $P(x)=\langle x|\rho_{0}|x\rangle.$
Finally, $P(y)$ is the probability for the outcomes of the second measurement.
It is

$P(y)=\frac{1}{2^{N}}\Big{[}1+y\langle
x\rangle\delta_{\underline{\mathbf{y}},\underline{\mathbf{x}}}\sum_{\mathbf{\alpha}}H_{\underline{\mathbf{y}}\mathbf{\alpha}}(1|\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0})\Big{]}.$
(35)

The term
$\delta_{\underline{\mathbf{z}},\mathbf{y}}\delta_{\underline{\mathbf{y}},\underline{\mathbf{x}}}$
in Eq. (34) implies that, for observables defined by unique Pauli strings,
memory effects are detected only when the three observables are the same
$S_{\underline{\mathbf{x}}}=S_{\underline{\mathbf{y}}}=S_{\underline{\mathbf{z}}}.$
This constraint does not emerge when the observables correspond to other basis
of operators (see for example Ref. BIF ).

The general solution Eq. (34) can be specified for the trigonometric eternal
model [Eq. (23)]. Stationary initial conditions are assumed,
$|q_{t})=|q_{0}),$ with $q_{0}^{(1)}=\varphi/(\varphi+\gamma)$ and
$q_{0}^{(2)}=\gamma/(\varphi+\gamma).$ For simplicity, first we consider the
case $N=1.$ When the three measurements are performed in direction
$\mathbf{\underline{a}}$ or $\mathbf{\underline{b}}$ we get

$\displaystyle C_{pf}(t,\tau)|_{y}$ $\displaystyle=$
$\displaystyle-\frac{(1-\langle
x\rangle^{2})}{[2^{N}P(y)]^{2}}\exp\left[-(t+\tau)(\gamma+\varphi)/2\right]$
(36)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\times\frac{4^{2}\gamma^{2}\varphi^{2}}{(\gamma+\varphi)^{2}\Upsilon^{2}}\sinh\left(\frac{\Upsilon
t}{2}\right)\sinh\left(\frac{\Upsilon\tau}{2}\right).\ \ \ \ \ $

When the three measurements are performed in direction
$\mathbf{\underline{c},}$ we get

$\displaystyle C_{pf}(t,\tau)|_{y}$ $\displaystyle=$
$\displaystyle\frac{(1-\langle x\rangle^{2})}{[2^{N}P(y)]^{2}}\
\frac{4\gamma\varphi(\gamma-\varphi)^{2}}{(\gamma+\varphi)^{4}}$
$\displaystyle\times[1-e^{-\tau(\gamma+\varphi)}][1-e^{-t(\gamma+\varphi)}].$

These results allow us to analyze the transition to divergent rates [Eq. (26)]
in a complementary way. In Fig. 2 we plot the CPF correlation Eq. (36). We
observe that when the rate $\gamma_{t}^{\underline{\mathbf{c}}}$ does not
develop divergences [Fig. 2(a)], the CPF correlation is negative for any value
of the time intervals $t$ and $\tau.$ On the other hand, in the interval
$3-\sqrt{8}<(\varphi/\gamma)<3+\sqrt{8}$ where the rate
$\gamma_{t}^{\underline{\mathbf{c}}}$ develops divergences [Fig. 2(b)], the
CPF correlation presents oscillations between positive an negative values.

For the model (23), the generalization to $N>1,$ independently of the chosen
observables, always lead to Eq. (36) or Eq. (IV). This results follows by
noting that in Eq. (34) the coefficients $\mathbf{\alpha}$ and
$\mathbf{\beta}$ only assume the four values
$\mathbf{\alpha,\beta}=(\mathbf{0,}\underline{\mathbf{a}},\underline{\mathbf{b}},\underline{\mathbf{c}})$
[see Eq. (24)]. Furthermore, using that
$H_{\underline{\mathbf{y}}\mathbf{\alpha}}H_{\underline{\mathbf{y}}\mathbf{\beta}}=H_{\underline{\mathbf{y}}\mathbf{\gamma}},$
where $\mathbf{\gamma}$ corresponds to the string
$S_{\mathbf{\gamma}}=S_{\mathbf{\alpha}}S_{\mathbf{\beta}},$ for a fixed
$\underline{\mathbf{y}}$ $(\underline{\mathbf{y}}=\underline{\mathbf{a}},$ or
$\underline{\mathbf{y}}=\underline{\mathbf{b}},$ or
$\underline{\mathbf{y}}=\underline{\mathbf{c}})$ the four matrix elements
$H_{\underline{\mathbf{y}}\mathbf{\gamma}},$ similarly to the case $N=1,$ can
only assume the values $(\pm 1),$ which always lead to Eq. (36) or Eq. (IV).
On the other hand, for $N>1$ accidentally it may also happen that the CPF
correlation vanishes. This occur because we assumed that the incoherent
degrees of freedom are stationary, which implies
$\sum_{\mathbf{\alpha,\beta}}[(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)\mathbb{F}_{\mathbf{\beta}}(t)|q_{0})-(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)|q_{t})(1|\mathbb{F}_{\mathbf{\beta}}(t)|q_{0})]=0.$
Thus, when $H_{\underline{\mathbf{y}}\mathbf{\gamma}}=1,$ it follows that
$C_{pf}(t,\tau)|_{y}=0$ [see Eq. (34)]. These accidental cases can always be
surpassed by considering arbitrary measurement operators written as linear
combinations of the Pauli strings.

Figure 2: CPF correlation [Eq. (36)] corresponding to the eternal non-
Markovian trigonometric model [Eq. (23)], for different values of
$\varphi/\gamma$ and measurement time-interval relations $\tau/t.$ In all
cases, the system initial condition is such that $\langle x\rangle=0.$

As an example, we consider a bipartite case where
$S_{\underline{\mathbf{a}}}=\sigma_{1}^{x}\sigma_{2}^{x},$
$S_{\underline{\mathbf{b}}}=\sigma_{1}^{y}\sigma_{2}^{y},$ and
$S_{\underline{\mathbf{c}}}=\sigma_{1}^{z}\sigma_{2}^{z}.$ The CPF correlation
Eq. (36) is obtained when the three measurement are defined by any of the
bipartite operators
$S_{\underline{\mathbf{m}}}=(\sigma_{1}^{x}\sigma_{2}^{z}),$
$(\sigma_{1}^{y}\sigma_{2}^{z}),$ $(\sigma_{1}^{z}\sigma_{2}^{x}),$
$(\sigma_{1}^{z}\sigma_{2}^{y}),$ Eq. (IV) is obtained when
$S_{\underline{\mathbf{m}}}=(\sigma_{1}^{x}\sigma_{2}^{y}),(\sigma_{1}^{y}\sigma_{2}^{x}),$
while $C_{pf}(t,\tau)|_{y}=0$ when
$S_{\underline{\mathbf{m}}}=(\sigma_{1}^{x}\sigma_{2}^{x}),$
$(\sigma_{1}^{y}\sigma_{2}^{y}),$ $(\sigma_{1}^{z}\sigma_{2}^{z}).$

##  V Summary and conclusions

We studied a class of solvable multipartite non-Markovian master equations
where the system consists of an arbitrary number of qubits and whose structure
is written in terms of arbitrary multipartite Pauli coupling terms. Starting
from a local-in-time representation of the evolution, we found the explicit
solution for the system density matrix, which in turn allowed us to formulate
the constraints that time-dependent rates must obey in order to guarantee the
completely positive condition of the solution map.

We also found explicit analytical expressions for the time-dependent rates
associated to a given evolution. Their sign (positive or negative) can be used
as an indicator of non-Markovianity. Memory effects were also characterized by
operational methods, where a CPF correlation defined by a set of three
consecutive system measurements becomes a memory witness. We showed that this
quantity can be obtained in an exact way for arbitrary measurement processes
and arbitrary interaction with incoherent degrees of freedom.

As application of the previous results, we presented simple underlying
dynamics that lead to the phenomenon of eternal non-Markovianity, that is,
multipartite dynamics where some rates depart at all times from that of a
Markovian regime. Both hyperbolic and trigonometric cases were established,
characterized by a rate that is negative at all times or that develops
periodical divergences. Even when these features develop, the CPF correlation
is always a smooth function.

In the Appendices we found the rates associated to different underlying memory
mechanisms such as a mapping with a classical master equation, stochastic
Hamiltonians and statistical superpositions of Markovian dynamics. We showed
that under particular conditions different mechanisms may lead to the same
time-dependent rates. Nevertheless, these accidental degeneracies do not occur
in general. We also found that the phenomenon of eternal non-Markovianity
becomes quite common in multipartite dynamics.

The class of models we studied here provides a useful solvable framework for
studying quantum non-Markovianity in multipartite settings. This allows to
formulate a wide range of well-behaved multipartite non-Markovian master
equations. The study of diverse memory witness can be tackled starting from
here. Our results also lead to interesting questions such as determining which
kind of underlying dynamics can be associated to an arbitrary non-Markovian
multipartite Pauli evolution.

## Acknowledgments

AAB acknowledges support from CONICET, Argentina. JPG acknowledges financial
support from EPSRC Grant no. EP/R04421X/1 and the Leverhulme Trust Grant No.
RPG-2018-181.

##  Appendix A Time-convoluted approach

Instead of the local-in-time formulation defined by Eq. (2), alternatively one
may start with a time convoluted evolution

$\frac{d}{dt}\rho_{t}=\mathcal{L}[\rho_{t}]=\sum_{\begin{subarray}{c}\mathbf{a}\\\
\mathbf{a\neq
0}\end{subarray}}\int_{0}^{t}dt^{\prime}k_{\mathbf{a}}(t-t^{\prime})(S_{\mathbf{a}}\rho_{t^{\prime}}S_{\mathbf{a}}-\rho_{t^{\prime}}),$
(38)

where the set of time-dependent kernels $\\{k_{\mathbf{a}}(t)\\}$ must to be
constrained such that the solution map is CP. Similarly to Sec. II, by
defining the kernel $k_{\mathbf{0}}(t)\equiv-\sum_{\mathbf{a\ (a\neq
0)}}k_{\mathbf{a}}(t),$ here the weights of the solution (12) can be written
as

$p_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\lambda_{\mathbf{b}}(t),$
(39)

where the coefficients $\lambda_{\mathbf{b}}(t)$ obey the evolution

$\frac{d}{dt}\lambda_{\mathbf{b}}(t)=\int_{0}^{t}dt^{\prime}k_{\mathbf{b}}(t-t^{\prime})\lambda_{\mathbf{b}}(t^{\prime}).$
(40)

The inverse relations for determining the kernels $\\{k_{\mathbf{a}}(t)\\}$ a
function of probabilities $\\{p_{\mathbf{a}}(t)\\}$ can be written in a
Laplace domain $[f(z)=\int_{0}^{\infty}dte^{-zt}f(t)]$ as

$k_{\mathbf{a}}(z)=\frac{z\lambda_{\mathbf{a}}(z)-1}{\lambda_{\mathbf{a}}(z)},\
\ \ \ \ \ \
\lambda_{\mathbf{a}}(z)=\sum_{\mathbf{b}}H_{\mathbf{ab}}p_{\mathbf{b}}(z).$
(41)

##  Appendix B Non-Markovian underlying mechanisms

Here, we consider different mechanisms that lead to memory effects. The
present analysis provides nontrivial multipartite extensions of some results
developed in Ref. megier  for the case $N=1.$

###  B.1 Mapping with a classical Markovian master equation

The solution map [Eq. (12)] is defined by a set of normalized probabilities
$\\{p_{t}^{\mathbf{a}}\\}.$ It is possible to formulate an underlying
mechanism such that $\\{p_{t}^{\mathbf{a}}\\}$ correspond to the solution of
an arbitrary Markovian classical master equation with $4^{N}$ different
states.

We assume that the system density matrix interacts with an incoherent system
whose states, in contrast to Eq. (18), can be put in one-to-one correspondence
with the Pauli string vectors $\\{\mathbf{a}\\}.$ Therefore, the system
density matrix $\rho_{t}$ can be written in terms of a set of auxiliary states
$\\{\rho_{t}^{\mathbf{a}}\\}$ lindbladrate  such that

$\rho_{t}=\sum_{\mathbf{a}}\rho_{t}^{\mathbf{a}}.$ (42)

The evolution of the auxiliary states is Markovian and involves coupling
between all of them. We write

$\frac{d}{dt}\rho_{t}^{\mathbf{a}}=-\sum_{\begin{subarray}{c}\mathbf{b}\\\
\mathbf{b}\neq\mathbf{a}\end{subarray}}\phi_{\mathbf{ba}}\rho_{t}^{\mathbf{a}}+\sum_{\begin{subarray}{c}\mathbf{b}\\\
\mathbf{b}\neq\mathbf{a}\end{subarray}}\phi_{\mathbf{ab}}S_{\mathbf{a}}S_{\mathbf{b}}\rho_{t}^{\mathbf{b}}S_{\mathbf{b}}S_{\mathbf{a}}.$
(43)

Here, $\\{\phi_{\mathbf{ba}}\\}$ are arbitrary rates. The stochastic
interpretation of this equation is quite simple. Whenever the incoherent
system undergoes the transition $\mathbf{b}\rightarrow\mathbf{a},$ the quantum
system undergoes the transformation $\rho\rightarrow
S_{\mathbf{a}}S_{\mathbf{b}}\rho S_{\mathbf{b}}S_{\mathbf{a}}.$ Between
transition the system is frozen. The average system dynamics is given by Eq.
(43), where $\rho_{t}^{\mathbf{a}}$ corresponds to the conditional system
state given that the incoherent one is in the state associated to
$\mathbf{a}.$

It is simple to check that the solutions $\\{\rho_{t}^{\mathbf{a}}\\}$ of Eq.
(43) can be written as

$\rho_{t}^{\mathbf{a}}=p_{t}^{\mathbf{a}}(S_{\mathbf{a}}\rho_{0}S_{\mathbf{a}}),$
(44)

where the weights $p_{t}^{\mathbf{a}}$ must to fulfill the classical master
equation

$\frac{d}{dt}p_{t}^{\mathbf{a}}=-\sum_{\begin{subarray}{c}\mathbf{b}\\\
\mathbf{b}\neq\mathbf{a}\end{subarray}}\phi_{\mathbf{ba}}p_{t}^{\mathbf{a}}+\sum_{\begin{subarray}{c}\mathbf{b}\\\
\mathbf{b}\neq\mathbf{a}\end{subarray}}\phi_{\mathbf{ab}}p_{t}^{\mathbf{b}}.$
(45)

Consequently, from Eqs. (42) and (44) we recover the solution Eq. (12)
$[\rho_{t}=\sum_{\mathbf{a}}p_{t}^{\mathbf{a}}(S_{\mathbf{a}}\rho_{0}S_{\mathbf{a}})],$
where the probabilities $\\{p_{t}^{\mathbf{a}}\\}$ fulfill a the classical
master equation (45). For consistence, its initial condition must be
$p_{0}^{\mathbf{a}}=\delta_{\mathbf{a,0}}.$

Particular case: Given that Eq. (45) is arbitrary, it is not possible to find
a general expression for the rates $\\{\gamma_{t}^{\mathbf{a}}\\}$ [Eq. (15)]
in terms of the underlying ones $\\{\phi_{\mathbf{ba}}\\}.$ Nevertheless, this
mapping can be performed, for example, when Eq. (45) assumes the form

$\frac{d}{dt}p_{t}^{\mathbf{0}}=-\phi
p_{t}^{\mathbf{0}}+\varphi\sum_{\begin{subarray}{c}\mathbf{a}\\\ \mathbf{a\neq
0}\end{subarray}}p_{t}^{\mathbf{a}},\ \ \ \ \
\frac{d}{dt}p_{t}^{\mathbf{a}}=-\varphi p_{t}^{\mathbf{a}}+x_{\mathbf{a}}\phi
p_{t}^{\mathbf{0}},$ (46)

where $\phi$ and $\varphi$ are arbitrary rates and the weights
$\\{x_{\mathbf{a}}\\}$ satisfies
$\sum_{\mathbf{a(a}\neq\mathbf{0)}}x_{\mathbf{a}}=1.$ The probabilities, with
initial condition $p_{0}^{\mathbf{a}}=\delta_{\mathbf{a,0}},$ can be written
as

$p_{t}^{\mathbf{a}}=p_{\infty}^{\mathbf{a}}[1-\exp(-\Phi
t)]+\delta_{\mathbf{a,0}}\exp(-\Phi t),$ (47)

where $\Phi\equiv(\phi+\varphi),$ and the stationary values are

$p_{\infty}^{\mathbf{0}}=\frac{\varphi}{\phi+\varphi},\ \ \ \ \ \
p_{\infty}^{\mathbf{a}}=\frac{x_{\mathbf{a}}\phi}{\phi+\varphi}.$ (48)

From the solutions (47), the general expression (15), after some calculations
steps tanh , lead to

$\gamma_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}\frac{\Phi}{2}H_{\mathbf{ab}}\left[\tanh\left(\frac{t\Phi}{2}+\zeta_{\mathbf{b}}\right)-1\right],$
(49)

where the parameters are

$\zeta_{\mathbf{b}}\equiv\frac{1}{2}\ln\left(\frac{h_{\infty}^{\mathbf{b}}}{1-h_{\infty}^{\mathbf{b}}}\right),\
\ \ \ \ \ \ \ \
h_{\infty}^{\mathbf{b}}\equiv\sum_{\mathbf{c}}H_{\mathbf{bc}}p_{\infty}^{\mathbf{c}}.$
(50)

It is simple to check that, due to probability normalization,
$h_{\infty}^{\mathbf{0}}=1.$ Hence, in Eq. (49) the term with $\mathbf{b=0}$
cancels out. Furthermore if $h_{\infty}^{\mathbf{b}}=0,$ it follows
$\tanh(t\Phi/2+\zeta_{\mathbf{b}})\rightarrow-1.$ In general, the time
dependence of the rate $\gamma_{t}^{\mathbf{a}}$ arise from a linear
combination of hyperbolic tangent functions with coefficient that are $\pm 1.$
Thus, in general some rates can be negative at any time.

###  B.2 Stochastic Hamiltonians

We consider a stochastic evolution, where the system wave vector
$|\psi_{t}\rangle$ is driven by a stochastic Hamiltonian,

$\frac{d|\psi_{t}\rangle}{dt}=-iH_{st}|\psi_{t}\rangle=-i\frac{1}{2}\xi_{t}^{\mathbf{\alpha}}S_{\mathbf{\alpha}}|\psi_{t}\rangle.$
(51)

The Hamiltonian $H_{st}$ is characterized by a noise with an arbitrary
statistics but null average
$\langle\langle\xi_{t}^{\mathbf{\alpha}}\rangle\rangle=0.$ The index
$\mathbf{\alpha\leftrightarrow\alpha}_{t}$ run overs all possible Pauli
strings. Its time variation is very slow such that over a single realization
it can be considered as a frozen parameter. Thus, the average state
$\rho_{t}^{\mathbf{\alpha}}=\langle\langle|\psi_{t}\rangle\langle\psi_{t}|\rangle\rangle$
for a given $\mathbf{\alpha}$ reads
$\rho_{t}^{\mathbf{\alpha}}=(1/2)[1+G_{t}^{\mathbf{\alpha}}]\rho_{0}+(1/2)[1-G_{t}^{\mathbf{\alpha}}](S_{\mathbf{\alpha}}\rho_{0}S_{\mathbf{\alpha}}),$
where

$G_{t}^{\mathbf{\alpha}}\equiv\Big{\langle}\Big{\langle}\exp\Big{(}i\int_{0}^{t}dt^{\prime}\xi_{t^{\prime}}^{\mathbf{\alpha}}\Big{)}\Big{\rangle}\Big{\rangle},$
(52)

is the characteristic noise function for a given $\mathbf{\alpha.}$ After
averaging this parameter, the system state can be written as
$\rho_{t}=\sum_{\mathbf{\alpha},(\mathbf{\alpha}\neq\mathbf{0})}x_{\mathbf{\alpha}}\rho_{t}^{\mathbf{\alpha}},$
where
$\sum_{\mathbf{\alpha},(\mathbf{\alpha}\neq\mathbf{0})}x_{\mathbf{\alpha}}=1.$
The parameters $\\{x_{\mathbf{\alpha}}\\}$ correspond to the statistical
weight of each Pauli string during the variation of the coefficient
$\mathbf{\alpha.}$ It is straightforward to check that
$\rho_{t}=\sum_{\mathbf{a}}p_{t}^{\mathbf{a}}(S_{\mathbf{a}}\rho_{0}S_{\mathbf{a}}),$
which recovers Eq. (12) with

$p_{t}^{\mathbf{0}}=\frac{1}{2}(1+\sum_{\begin{subarray}{c}\mathbf{a}\\\
\mathbf{a}\neq\mathbf{0}\end{subarray}}x_{\mathbf{a}}G_{t}^{\mathbf{a}}),\ \ \
\ \ \ p_{t}^{\mathbf{a}}=\frac{x_{\mathbf{a}}}{2}(1-G_{t}^{\mathbf{a}}).$ (53)

Similarly to the previous model, it is not possible to find a general simple
expression for the rates $\gamma_{t}^{\mathbf{a}}$ in terms of these
probabilities. Manageable expressions arise in the following situations.

Particular cases: If the noise is the same for all “directions”
$G_{t}^{\mathbf{a}}=G_{t},$ from Eqs. (15) and (53), after some algebra tanh ,
we get the rates

$\gamma_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}\frac{\dot{g}_{t}}{2}H_{\mathbf{ab}}\left[\tanh\left(\frac{g_{t}}{2}+\zeta_{\mathbf{b}}\right)-1\right],$
(54)

where the scalar functions read

$g_{t}=\ln(1/G_{t}),\ \ \ \ \ \ \
\dot{g}_{t}=-\frac{1}{G_{t}}\frac{dG_{t}}{dt},$ (55)

and where $\zeta_{\mathbf{b}}$ is defined by Eq. (50) with, instead of Eq.
(48), with $p_{\infty}^{\mathbf{0}}=1/2,$ and
$p_{\infty}^{\mathbf{a}}=x_{\mathbf{a}}/2.$

For a stationary Gaussian white noise, where
$\langle\langle\xi_{t}\xi_{t^{\prime}}\rangle\rangle=\Phi\delta(t-t^{\prime}),$
Eq. (52) becomes $G(t)=\exp(-\Phi t).$ It is simple to check that in this
situation Eq. (54) recovers the solution (49) of the previous model with
$\varphi=\phi.$ This results show that there are different underlying models
that may lead to the same system density matrix evolution. This degeneracy is
not universal and clearly depends on the underlying parameters.

For a stationary symmetric dichotomic noise with amplitude $A$ and switching
rate $\eta,$ the characteristic noise function [Eq. (52)] is

$G_{t}=e^{-\eta t}[\cosh(\chi t)+\frac{\eta}{\chi}\sinh(\chi t)],\ \ \ \ \
\chi\equiv\sqrt{\eta^{2}-A^{2}}.$ (56)

In contrast to the previous cases, here the rates defined by Eq. (54) may
develop divergences. In fact, the functions (55) become

$g_{t}=\ln(1/G_{t}),\ \ \ \ \ \ \
\dot{g}_{t}=\frac{A^{2}}{\eta+\chi[1/\tanh(\chi t)]}.$ (57)

Hence, divergent rates are found whenever $\eta<A.$

###  B.3 Statistical mixtures of Markovian evolutions

Departures with respect to a Markovian regime emerge whenever the system
evolution is written as the statistical superposition of different Markovian
propagators. Hence, we write

$p_{t}^{\mathbf{a}}=\sum_{k=1}^{n}q_{k}p_{k}^{\mathbf{a}}(t),$ (58)

where $\\{q_{k}\\}$ are normalized positive weights $(\sum_{k=1}^{n}q_{k}=1),$
and each set of probabilities $\\{p_{k}^{\mathbf{a}}(t)\\}$ is associated to a
Markovian solution of Eq. (2) with time-independent positive rates
$\\{\gamma_{k}^{\mathbf{a}}\\}.$

From Eq. (15), the non-Markovian evolution is characterized by the rates

$\gamma_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\frac{\sum_{k=1}^{n}q_{k}\mu_{k}^{\mathbf{b}}\exp(t\mu_{k}^{\mathbf{b}})}{\sum_{k^{\prime}=1}^{n}q_{k^{\prime}}\exp(t\mu_{k^{\prime}}^{\mathbf{b}})}.$
(59)

where $\mu_{k}^{\mathbf{b}}$ are eigenvalues of the $k$-Markovian dynamics,
$\mu_{k}^{\mathbf{b}}=\sum_{\mathbf{c}}H_{\mathbf{bc}}\gamma_{k}^{\mathbf{c}}.$
The specific properties of these rates strongly depend on the considered
Markovian evolutions and statistic weights.

Particular cases: In the two-state case, $n=2,$ the probabilities are
$p_{t}^{\mathbf{a}}=q_{1}p_{1}^{\mathbf{a}}(t)+q_{2}p_{2}^{\mathbf{a}}(t),$
where each solution is associated to the rates $\gamma_{1}^{\mathbf{a}}$ and
$\gamma_{2}^{\mathbf{a}},$ and $q_{1}+q_{2}=1.$ From Eq. (59), after some
algebra tanh , we get

$\gamma_{t}^{\mathbf{a}}=\frac{1}{2}(\gamma_{1}^{\mathbf{a}}+\gamma_{2}^{\mathbf{a}})+\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\Delta_{\mathbf{b}}\tanh(t\Delta_{\mathbf{b}}+\zeta),$
(60)

where the parameters are

$\Delta_{\mathbf{b}}\equiv\frac{1}{2}\sum_{\mathbf{c}}H_{\mathbf{bc}}(\gamma_{1}^{\mathbf{c}}-\gamma_{2}^{\mathbf{c}}),\
\ \ \ \ \ \ \zeta\equiv\frac{1}{2}\ln(\frac{q_{1}}{q_{2}}).$ (61)

In this case, many rates may also be negative at all times (see next section).

In the other extreme, a continuos-state case can be considered. Thus, Eq. (58)
is rewritten as

$p_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\Big{\langle}\prod_{\mathbf{c}}\exp(tH_{\mathbf{bc}}\gamma^{\mathbf{c}})\Big{\rangle},$
(62)

where we used the explicit expression (13) and the replacement
$\sum_{k=1}^{n}q_{k}\rightarrow\left\langle\cdots\right\rangle.$ The symbol
$\left\langle\cdots\right\rangle$ denotes an average over the set of random
rates $\\{\gamma^{\mathbf{c}}\\},$ each “realization” defining a Markov
solution. Assuming that all rates are independent random variables it follows
that
$\left\langle\cdots\right\rangle\rightarrow\int_{0}^{\infty}d\gamma^{\mathbf{c}}\cdots
P(\gamma^{\mathbf{c}}),$ where $P(\gamma^{\mathbf{c}})$ is the corresponding
probability density. By assuming an exponential probability density
$P(\gamma^{\mathbf{c}})=\tau_{\mathbf{c}}\exp(-\gamma^{\mathbf{c}}\tau_{\mathbf{c}}),$
by using that
$\gamma^{\mathbf{0}}=-\sum_{\mathbf{c}(\mathbf{c}\neq\mathbf{0})}\gamma^{\mathbf{c}},$
[see Eq. (7)] we get

$p_{t}^{\mathbf{a}}=\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\prod_{\begin{subarray}{c}\mathbf{c}\\\
\mathbf{c}\neq\mathbf{0}\end{subarray}}\frac{\tau_{\mathbf{c}}}{\tau_{\mathbf{c}}+(1-H_{\mathbf{bc}})t},$
(63)

where we have used that $H_{\mathbf{b0}}=1.$ From Eq. (15), the corresponding
rates associated to the non-Markovian evolution are

$\gamma_{t}^{\mathbf{a}}=-\frac{1}{4^{N}}\sum_{\mathbf{b}}H_{\mathbf{ab}}\sum_{\begin{subarray}{c}\mathbf{c}\\\
\mathbf{c}\neq\mathbf{0}\end{subarray}}\frac{(1-H_{\mathbf{bc}})}{\tau_{\mathbf{c}}+(1-H_{\mathbf{bc}})t}.$
(64)

We notice that both $\\{p_{t}^{\mathbf{a}}\\}$ and
$\\{\gamma_{t}^{\mathbf{a}}\\}$ develop a power-law behavior. In spite of this
feature the rates are positive at all times, $\gamma_{t}^{\mathbf{a}}>0$
$(\mathbf{a}\neq\mathbf{0}).$ While most of the memory witnesses BreuerReview
; plenioReview  associate this property to a Markovian regime, from
operational approaches it is possible to detect and infer the presence of
memory effects budiniCPF ; BIF .

##  Appendix C Bipartite and tripartite eternal non-Markovian evolutions

Besides the previous examples, the developed approach allow us to show that
master equations characterized by eternal non-Markovian effects are quite
common for multipartite systems. As an example, we consider the statistical
superposition of two different Markovian dynamics characterized by the rates
$\gamma_{1}^{\mathbf{a}}$ and $\gamma_{2}^{\mathbf{a}}$ and equal weights
$[q_{1}=q_{2}$ in Eq. (60)]. Taking
$\gamma_{1}^{\mathbf{a}}=\gamma(\delta_{\mathbf{a},\underline{\mathbf{a}}}-\delta_{\mathbf{a},\mathbf{0}}),$
and
$\gamma_{2}^{\mathbf{a}}=\gamma(\delta_{\mathbf{a},\underline{\mathbf{b}}}-\delta_{\mathbf{a},\mathbf{0}}),$
and using that
$(H_{\mathbf{\alpha}\underline{\mathbf{a}}}-H_{\mathbf{\alpha}\underline{\mathbf{b}}})/2=(\pm
1,0),$ and
$H_{\mathbf{\alpha}\underline{\mathbf{a}}}H_{\mathbf{\alpha}\underline{\mathbf{b}}}=H_{\mathbf{\alpha}\underline{\mathbf{c}}},$
from Eq. (60) we recover the rates defined in Eq. (22). When each (vectorial)
rate involves different Pauli channels more complex expressions are obtained.

As a first example, take a bipartite system $(N=2)$ with

$\displaystyle\gamma_{1}^{\mathbf{a}}$ $\displaystyle=$
$\displaystyle\gamma(\delta_{\mathbf{a},10}+\delta_{\mathbf{a},01}-2\delta_{\mathbf{a},00}),$
(65a) $\displaystyle\gamma_{2}^{\mathbf{a}}$ $\displaystyle=$
$\displaystyle\gamma(\delta_{\mathbf{a},20}+\delta_{\mathbf{a},02}-2\delta_{\mathbf{a},00}).$
(65b) Thus, each dynamics is defined by a local (single) dephasing local
mechanism acting alternatively in $x$\- and $y$-directions. From Eq. (60) we
obtain  $\gamma_{t}^{\mathbf{a}_{0}}=\frac{1}{2}\gamma,\ \ \ \ \ \ \
\gamma_{t}^{\mathbf{a}_{\pm}}=\pm\frac{1}{4}\gamma\tanh(2\gamma t),$ (66)

where $\mathbf{a}_{0}$ and $\mathbf{a}_{\pm}$ correspond to the following
Pauli strings, $\mathbf{a}_{0}=(10),\ (01),$ $(20),$ $(02),$ and
$\mathbf{a}_{+}=(11),\ (22),$ while $\mathbf{a}_{-}=(30),$ $(03),$ $(12),$
$(21).$ Furthermore,

$\gamma_{t}^{33}=-\frac{\gamma}{4}[2\tanh(\gamma t)-\tanh(2\gamma
t)]=-2\gamma\frac{\sinh^{4}(\gamma t)}{\sinh(4\gamma t)},$

while $\gamma_{t}^{\mathbf{a}}=0$ if
$\mathbf{a}\neq(\mathbf{a}_{0},\mathbf{a}_{+},\mathbf{a}_{-}).$ There are
eleven non-null rates out of the fifteen possible ones, five of them being
negative at all times.

As a second example we consider a tripartite system $(N=3),$ where

$\displaystyle\gamma_{1}^{\mathbf{a}}$ $\displaystyle=$
$\displaystyle\gamma(\delta_{\mathbf{a},110}+\delta_{\mathbf{a},101}+\delta_{\mathbf{a},011}-3\delta_{\mathbf{a},000}),$
(67a) $\displaystyle\gamma_{2}^{\mathbf{a}}$ $\displaystyle=$
$\displaystyle\gamma(\delta_{\mathbf{a},220}+\delta_{\mathbf{a},202}+\delta_{\mathbf{a},022}-3\delta_{\mathbf{a},000}).$
(67b) Hence, each Markovian evolution correspond to dephasing in $x$\- and
$y$-directions but now considering all pairs of bipartite dephasing operators.
From Eq. (60) we get  $\displaystyle\gamma_{t}^{\mathbf{a}_{+}}$
$\displaystyle=$ $\displaystyle\frac{1}{4}\gamma[2+\tanh(2\gamma t)],$ (68a)
$\displaystyle\gamma_{t}^{\mathbf{a}_{-}}$ $\displaystyle=$
$\displaystyle-\frac{1}{4}\gamma\tanh(2\gamma t),$ (68b) where
$\mathbf{a}_{\pm}$ correspond to the following Pauli strings,
$\mathbf{a}_{+}=(110),\ (101),$ $(011),\ (220),\ (202),\ (022),$ while
$\mathbf{a}_{-}=(330),$ $(303),$ $(033),$ $(123),$ $(132),$ $(213),$ $(231),$
$(312),$ $(321),$ and $\gamma_{t}^{\mathbf{a}}=0$ if
$\mathbf{a}\neq\mathbf{a}_{+},\mathbf{a}_{-}.$ In this case, out of sixty-
three possible rates, fifteen are non-null, nine of them being negative at all
times.

##  Appendix D CPF correlation calculus

For a system coupled to incoherent degrees of freedom [Eq. (18)], the
(bipartite) system-environment state
$\rho_{t}^{se}=\sum_{\mathbf{h}}\rho_{t}^{\mathbf{h}}|\mathbf{h}),$ from Eqs.
(19) and (33), reads

$\rho_{t}^{se}=\sum_{\mathbf{\alpha}}(S_{\mathbf{\alpha}}\rho_{0}S_{\mathbf{\alpha}})\
\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0}).$ (69)

This evolution defines the system-environment dynamics between measurements.
The measurement of operator $S_{\underline{\mathbf{m}}}$ [Eq. (31)] leads to
the transformation
$\rho^{se}=\sum_{\mathbf{h}}\rho^{\mathbf{h}}|\mathbf{h})\rightarrow|m\rangle\langle
m||q_{m}),$ where $|q_{m})=\sum_{\mathbf{h}}\ \langle
m|\rho_{t}^{\mathbf{h}}|m\rangle/\mathrm{Tr}[\langle
m|\rho_{t}^{\mathbf{h}}|m\rangle]|\mathbf{h}).$ With these ingredients, the
calculation of the joint probability can be performed in a standard way. We
get,

$\frac{P(z,y,x)}{P(x)}=\sum_{\mathbf{\alpha,\beta}}|\langle
z|\sigma_{\mathbf{\alpha}}|y\rangle|^{2}|\langle
y|S_{\mathbf{\beta}}|x\rangle|^{2}(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)\mathbb{F}_{\mathbf{\beta}}(t)|q_{0}),$
(70)

where $P(x)=\langle x|\rho_{0}|x\rangle$ and
$(1|\equiv\sum_{\mathbf{h}}(\mathbf{h}|.$ This result is valid for arbitrary
Hermitian system observables.

Using Bayes rule, the conditional probabilities that define the CPF
correlation [Eq. (32)] can be written as $P(z,x|y)=P(z,y,x)/P(y),$ where
$P(y)=\sum_{z,x}P(z,y,x).$ Furthermore, $P(z|y)=\sum_{x}P(z,x|y),$ and
$P(x|y)=\sum_{z}P(z,x|y).$ From Eq. (70), and using

$\displaystyle\sum_{z}z|\langle z|S_{\mathbf{\alpha}}|y\rangle|^{2}$
$\displaystyle=$ $\displaystyle\langle
y|S_{\mathbf{\alpha}}S_{\underline{\mathbf{z}}}S_{\mathbf{\alpha}}|y\rangle,$
(71a) $\displaystyle\sum_{x}x|\langle y|S_{\mathbf{\beta}}|x\rangle|^{2}P(x)$
$\displaystyle=$ $\displaystyle\langle
y|S_{\mathbf{\beta}}S_{\underline{\mathbf{x}}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\beta}}|y\rangle,$
(71b) $\displaystyle\sum_{x}|\langle y|S_{\mathbf{\beta}}|x\rangle|^{2}P(x)$
$\displaystyle=$ $\displaystyle\langle
y|S_{\mathbf{\beta}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\beta}}|y\rangle,$
(71c) where the system state $\rho_{\underline{\mathbf{x}}}$ is
$\rho_{\underline{\mathbf{x}}}\equiv\sum_{x}P(x)\ |x\rangle\langle
x|=\sum_{x}\langle x|\rho_{0}|x\rangle\ |x\rangle\langle x|,$ (72)

the CPF correlation can be written as

$C_{pf}(t,\tau)|_{y}=\frac{1}{P(y)^{2}}\sum_{\mathbf{\alpha,\beta,\gamma}}\Theta^{\mathbf{\alpha\beta\gamma}}|_{y}\Lambda_{\mathbf{\alpha\beta\gamma}}(t,\tau).$
(73)

The coefficients $\Theta^{\mathbf{\alpha\beta\gamma}}|_{y}$ are

$\Theta^{\mathbf{\alpha\beta\gamma}}|_{y}=\langle
y|S_{\mathbf{\alpha}}S_{\underline{\mathbf{z}}}S_{\mathbf{\alpha}}|y\rangle\langle
y|S_{\mathbf{\beta}}S_{\underline{\mathbf{x}}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\beta}}|y\rangle\langle
y|\sigma_{\mathbf{\gamma}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\gamma}}|y\rangle,$

while the time-dependence follows from

$\displaystyle\Lambda_{\mathbf{\alpha\beta\gamma}}(t,\tau)$ $\displaystyle=$
$\displaystyle+(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)\mathbb{F}_{\mathbf{\beta}}(t)|q_{0})(1|\mathbb{F}_{\mathbf{\gamma}}(t)|q_{0})$
$\displaystyle-(1|\mathbb{F}_{\mathbf{\alpha}}(\tau)\mathbb{F}_{\mathbf{\gamma}}(t)|q_{0})(1|\mathbb{F}_{\mathbf{\beta}}(t)|q_{0}),$

where $|q_{t})=\sum_{\mathbf{\alpha}}\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0}),$
and the probability $P(y)$ is

$P(y)=\sum_{\mathbf{\alpha}}(1|\mathbb{F}_{\mathbf{\alpha}}(t)|q_{0})\langle
y|S_{\mathbf{\alpha}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\alpha}}|y\rangle.$
(74)

The expression (73) is valid for arbitrary observables
$\sigma_{\underline{\mathbf{m}}}$ [Eq. (31)]. In general, they can be written
as linear combinations of Pauli strings $S_{\mathbf{a}}.$ Assuming, for
simplicity, that each $S_{\underline{\mathbf{m}}}$ correspond to a unique
Pauli string operator, from Eq. (5) it follows the relations

$\displaystyle\langle
y|S_{\mathbf{\alpha}}S_{\underline{\mathbf{z}}}\sigma_{\mathbf{\alpha}}|y\rangle$
$\displaystyle=$ $\displaystyle
H_{\mathbf{\alpha}\underline{\mathbf{y}}}\delta_{\underline{\mathbf{z}},\mathbf{y}}a_{y},$
(75a) $\displaystyle\langle
y|S_{\mathbf{\beta}}S_{\underline{\mathbf{x}}}\rho_{\underline{\mathbf{x}}}S_{\mathbf{\beta}}|y\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2^{N}}(H_{\mathbf{\beta}\underline{\mathbf{y}}}\delta_{\underline{\mathbf{y}},\underline{\mathbf{x}}}a_{y}+\langle
x\rangle),$ (75b) $\displaystyle\langle
y|S_{\mathbf{\gamma}}\rho_{\underline{\mathbf{x}}}\sigma_{\mathbf{\gamma}}|y\rangle$
$\displaystyle=$
$\displaystyle\frac{1}{2^{N}}(1+H_{\mathbf{\gamma}\underline{\mathbf{y}}}\delta_{\underline{\mathbf{y}},\underline{\mathbf{x}}}a_{y}\langle
x\rangle),\ \ \ $ (75c) where $\langle
x\rangle\equiv\mathrm{Tr}[S_{\underline{\mathbf{x}}}\rho_{\underline{\mathbf{x}}}].$
By introducing these equalities in Eq. (73), after some algebra we get Eq.
(34). Generalization to arbitrary observables can be worked out in a similar
way from Eq. (73).

## References

  * (1) R. Alicki and K. Lendi, Quantum Dynamical Semigroups and Applications, Lecture Notes in Physics 286 (Springer, Berlin, 1987).
  * (2) H. P. Breuer and F. Petruccione, The theory of open quantum systems, (Oxford University Press, 2002).
  * (3) I. de Vega and D. Alonso, Dynamics of non-Markovian open quantum systems, Rev. Mod. Phys. 89, 015001 (2017).
  * (4) L. Li, M. J. W. Hall, and H. M. Wiseman, Concepts of quantum non-Markovianity: A hierarchy, Phys. Rep. 759, 1 (2018).
  * (5) H. P. Breuer, E. M. Laine, J. Piilo, and V. Vacchini, Colloquium: Non-Markovian dynamics in open quantum systems, Rev. Mod. Phys. 88, 021002 (2016).
  * (6) A. Rivas, S. F. Huelga, and M. B. Plenio, Quantum non-Markovianity: characterization, quantification and detection, Rep. Prog. Phys. 77, 094001 (2014).
  * (7) J. Wilkie, Positivity preserving non-Markovian master equations, Phys. Rev. E 62, 8808 (2000).
  * (8) S. M. Barnett and S. Stenholm, Hazards of reservoir memory, Phys. Rev. A 64, 033808 (2001).
  * (9) A. A. Budini, Stochastic representation of a class of non-Markovian completely positive evolutions, Phys. Rev. A 69, 042107 (2004).
  * (10) S. Daffer, K. Wodkiewicz, J.D. Cresser, and J.K. McIver, Depolarizing channel as a completely positive map with memory, Phys. Rev. A 70, 010304(R) (2004).
  * (11) D. Chruściński and A. Kossakowski, Non-Markovian Quantum Dynamics: Local versus Nonlocal, Phys. Rev. Lett. 104, 070406 (2010).
  * (12) A. A. Budini, Quantum systems subject to the action of classical stochastic fields, Phys. Rev. A 64, 052110 (2001); J. I. Costa-Filho, R. B. B. Lima, R. R. Paiva, P. M. Soares, W. A. M. Morgado, R. Lo Franco, and D. O. Soares-Pinto, Enabling quantum non-Markovian dynamics by injection of classical colored noise, Phys. Rev. A 95, 052126 (2017); 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); A. Kiely, Exact classical noise master equations: Applications and connections, Euro Phys. Lett. 134, 10001 (2021).
  * (13) A. Shabani and D. A. Lidar, Completely positive post-Markovian master equation via a measurement approach, Phys. Rev. A 71, 020101(R) (2005); C. Sutherland, T. A. Brun, and D. A. Lidar, Non-Markovianity of the post-Markovian master equation, Phys. Rev. A 98, 042119 (2018).
  * (14) S. Maniscalco and F. Petruccione, Non-Markovian dynamics of a qubit, Phys. Rev. A 73, 012111 (2006).
  * (15) J. Salo, S. M. Barnett, and S. Stenholm, Non-Markovian thermalization of a two-level system, Op. Comm. 259, 772 (2006).
  * (16) D. Chruściński and A. Kossakowski, Sufficient conditions for a memory-kernel master equation, Phys. Rev. A 94, 020103(R) (2016).
  * (17) A. A. Budini, Lindblad rate equations, Phys. Rev. A 74, 053815 (2006).
  * (18) A. A. Budini, Post-Markovian quantum master equations from classical environment fluctuations, Phys. Rev. E 89, 012147 (2014).
  * (19) B. Vacchini, Non-Markovian dynamics for bipartite systems, Phys. Rev. A 78, 022112 (2008).
  * (20) N. Megier, D. Chruściński, J. Piilo, and W. T. Strunz, Eternal non-Markovianity: from random unitary to Markov chain realisations, Sci. Rep. 7, 6379 (2017).
  * (21) A. A. Budini, Maximally non-Markovian quantum dynamics without environment-to-system backflow of information, Phys. Rev. A 97, 052133 (2018).
  * (22) B. Donvil, P. Muratore-Ginanneschi, and J. P. Pekola, Hybrid master equation for calorimetric measurements, Phys. Rev. A 99, 042127 (2019).
  * (23) H. P. Breuer, B. Kappler, and F. Petruccione, Stochastic wave-function method for non-Markovian quantum master equations, Phys. Rev. A 59, 1633 (1999); H. P. Breuer, Genuine quantum trajectories for non-Markovian processes, Phys. Rev. A 70, 012106 (2004).
  * (24) M. R. Hush, I. Lesanovsky, and J. P. Garrahan, Generic map from non-Lindblad to Lindblad master equations, Phys. Rev. A 91, 032113 (2015).
  * (25) A. A. Budini, Embedding non-Markovian quantum collisional models into bipartite Markovian dynamics, Phys. Rev. A 88, 032115 (2013); A. A. Budini and P. Grigolini, Non-Markovian nonstationary completely positive open-quantum-system dynamics, Phys. Rev. A 80, 022103 (2009).
  * (26) B. Vacchini, Non-Markovian master equations from piecewise dynamics, Phys. Rev. A 87, 030101(R) (2013).
  * (27) V. Giovannetti and G. M. Palma, Master Equations for Correlated Quantum Channels, Phys. Rev. Lett. 108, 040401 (2012).
  * (28) N. K. Bernardes, A. R. R. Carvalho, C. H. Monken, and M. F. Santos, Environmental correlations and Markovian to non-Markovian transitions in collisional models, Phys. Rev. A 90, 032111 (2014).
  * (29) F. Ciccarello, G. M. Palma, and V. Giovannetti, Collision-model-based approach to non-Markovian quantum dynamics, Phys. Rev. A 87, 040103(R) (2013); S. Lorenzo, F. Ciccarello, and G. M. Palma, Class of exact memory-kernel master equations, Phys. Rev. A 93, 052111 (2016); S. Lorenzo, F. Ciccarello, and G. M. Palma, Composite quantum collision models, Phys. Rev. A 96, 032107 (2017).
  * (30) S. Kretschmer, K. Luoma, and W. T. Strunz, Collision model for non-Markovian quantum dynamics, Phys. Rev. A 94, 012106 (2016).
  * (31) B. Çakmak, M. Pezzutto, M. Paternostro, and Ö. E. Müstecaplıoglu, Non-Markovianity, coherence, and system-environment correlations in a long-range collision model, Phys. Rev. A 96, 022109 (2017).
  * (32) R. Ramirez Camasca and G. T. Landi, Memory kernel and divisibility of Gaussian collisional models, Phys. Rev. A 103, 022202 (2021).
  * (33) H. P. Breuer and B. Vacchini, Quantum Semi-Markov Processes, Phys. Rev. Lett. 101, 140402 (2008); H. P. Breuer and B. Vacchini, Structure of completely positive quantum master equations with memory kernel, Phys. Rev. E 79, 041147 (2009); B. Vacchini, Generalized Master Equations Leading to Completely Positive Dynamics, Phys. Rev. Lett. 117, 230401 (2016).
  * (34) D. Chruściński and A. Kossakowski, Generalized semi-Markov quantum evolution, Phys. Rev. A 95, 042131 (2017); D. Chruściński and A. Kossakowski, From Markovian semigroup to non-Markovian quantum evolution, Euro Phys. Lett. 97, 20005 (2012).
  * (35) D. Chruściński and F. A. Wudarski, Non-Markovian random unitary qubit dynamics, Phys. Lett. A 377, 1425 (2013); D. Chruściński and F. A. Wudarski, Non-Markovianity degree for random unitary evolution, Phys. Rev. A 91, 012104 (2015).
  * (36) F. A. Wudarski, P. Nalezyty, G. Sarbicki, and D. Chruściński, Admissible memory kernels for random unitary qubit evolution, Phys. Rev. A 91, 042105 (2015); F. A. Wudarski and D. Chruściński, Markovian semigroup from non-Markovian evolutions, Phys. Rev. A 93, 042120 (2016); D. Chruściński and K. Siudzińska, Generalized Pauli channels and a class of non-Markovian quantum evolution, Phys. Rev. A 94, 022118 (2016); K. Siudzińska and D. Chruściński, Memory kernel approach to generalized Pauli channels: Markovian, semi-Markov, and beyond, Phys. Rev. A 96, 022129 (2017); K. Siudzińska, Markovian semigroup from mixing noninvertible dynamical maps, Phys. Rev. A 103, 022605 (2021).
  * (37) E. Ferraro, M. Scala, R. Migliore, and A. Napoli, Non-Markovian dissipative dynamics of two coupled qubits in independent reservoirs: Comparison between exact solutions and master-equation approaches, Phys. Rev. A 80, 042112 (2009).
  * (38) B. Vacchini and H. P. Breuer, Exact master equations for the non-Markovian decay of a qubit, Phys. Rev. A 81, 042103 (2010).
  * (39) D. Maldonado-Mundo, P. Öhberg, B. W. Lovett, and E. Andersson, Investigating the generality of time-local master equations, Phys. Rev. A 86, 042107 (2012).
  * (40) J. Lankinen, H. Lyyra, B. Sokolov, J. Teittinen, B. Ziaei, and S. Maniscalco, Complete positivity, finite-temperature effects, and additivity of noise for time-local qubit dynamics, Phys. Rev. A 93, 052103 (2016).
  * (41) L. Ferialdi, Exact non-Markovian master equation for the spin-boson and Jaynes-Cummings models, Phys. Rev. A 95, 020101(R) (2017).
  * (42) H. Z. Shen, D. X. Li, Shi-Lei Su, Y. H. Zhou, and X. X. Yi, Exact non-Markovian dynamics of qubits coupled to two interacting environments, Phys. Rev. A 96, 033805 (2017).
  * (43) D. Tamascelli, A. Smirne, S. F. Huelga, and M. B. Plenio, Nonperturbative Treatment of non-Markovian Dynamics of Open Quantum Systems, Phys. Rev. Lett. 120, 030402 (2018).
  * (44) N. Megier, A. Smirne, and B. Vacchini, The interplay between local and non-local master equations: exact and approximated dynamics, New J. Phys. 22, 083011 (2020).
  * (45) D. Burgarth, P. Facchi , M. Ligabò, and D. Lonigro, Hidden non-Markovianity in open quantum systems, Phys. Rev. A 103, 012203 (2021).
  * (46) 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).
  * (47) A. A. Budini, Quantum Non-Markovian Processes Break Conditional Past-Future Independence, Phys. Rev. Lett. 121, 240401 (2018); A. A. Budini, Conditional past-future correlation induced by non-Markovian dephasing reservoirs, Phys. Rev. A 99, 052125 (2019).
  * (48) A. A. Budini, Detection of bidirectional system-environment information exchanges, Phys. Rev. A 103, 012221 (2021).
  * (49) M. J. W. Hall, J. D. Cresser, L. Li, and E. Andersson, Canonical form of master equations and characterization of non-Markovianity, Phys. Rev. A 89, 042120 (2014).
  * (50) H. J. Briegel and B. G. Englert, Quantum optical master equations: The use of damping bases, Phys. Ref. A 47, 3311 (1993); S. M. Barnetts and S. Stenholm, Spectral decomposition of the Lindblad operator, J. Mod. Optics 47, 2869 (2000).
  * (51) 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).
  * (52) D. Chruściński and S. Maniscalco, Degree of Non-Markovianity of Quantum Evolution, Phys. Rev. Lett. 112, 120404 (2014).
  * (53) We used the equality $(d/dt)\ln[pe^{at}+qe^{bt}]=(a+b)/2+\Delta\tanh[t\Delta+\zeta],$ where $\Delta=(a-b)/2$ and $\zeta=(1/2)\ln(p/q).$