A solvable class of non-Markovian quantum multipartite dynamics

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 [1]. 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 [2][3][4], is not known. The second one is that different inequivalent memory witnesses can be used to define and measure non-Markovian effects [5,6].
Despite these advances  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 ρ t is given by the time-local master equation (d/dt)ρ t = L[ρ t ], where the generator of the evolution has the general structure Here, σ α i is the α-th Pauli operator (α = x, y, z) acting on qubit i, while Γ α···β i (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 [20,21]. The explicit for-mulation of an operational (measurement based) memory witness [46][47][48] 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 [49].
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 a ≡ σ a1 ⊗ σ a2 ⊗ σ aN , each one associated to the vector a = (a 1 , a 2 , · · · , a N ). Each component a k (k = 1, 2, · · · N ) assumes the values a k = (0, 1, 2, 3) ↔ (I, σ x , σ y , σ 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 ρ t is written in a local-in-time way. Arbitrary multipartite decoherence channels are considered, The set of functions {γ a t } define the rates associated to the multipartite Pauli channel. In general, there are 4 N − 1 different rate functions, as the vector 0 = (0, 0, · · · , 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.

A. Subsystem dynamics
Given the evolution above, we ask about the dynamics of any particular subsystem. Introducing the splitting a = (a s , a e ), where a s corresponds to the set of local operators that define the marginal Pauli string of the subsystem of interest, and a e that of the rest of qubits (now considered as part of the environment), from Eq. (2) the subsystem density matrix ρ s t = Tr e [ρ t ] (where Tr [•] is the trace operation) reads 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.

B. Solution map and completely positive condition
We now show that by using the method of damping bases or spectral decomposition [50], the solution map ρ 0 → ρ 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 (σ 0 , σ 1 , σ 2 , σ 3 ) ↔ (I, σ x , σ y , σ z ), it is easy to check that where the input [•] is an arbitrary two dimensional operator and b = 0, 1, 2, 3. The inverse relation reads In these expressions, the coefficients {H ab } define a four dimensional Hadamard matrix H, which reads 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 γ 0 t , which is associated to the identity string in the full Hilbert space With this definition, the Lindbladian-like structure of Eq. (2)] can straightforwardly be written as where the sum now includes the (identity) string a = 0. Written in this way, applying the "vectorial extension" of Eq. (5) to the Hilbert space of N qubits, it follows that where H ab ≡ H a1b1 H a2b2 · · · H aN bN 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 Tr[S a S b ] = 2 N δ a,b , it is straightforward to determine the eigenvalues and eigenoperators of L [•]. They read Consequently, any Pauli string S a is a right eigenoperator with eigenvalue µ a t . Given that L[•] also defines the adjoint evolution (as the "jump operators" are Hermitian) [50], S a is also a left eigenoperator. Notice also that by using the inverse of the Hadamard matrix, the inverse relation γ a t = b H ab µ b t /4 N follows. From the method of damping bases [50], Eq. (10) allows us to write the solution of Eq. (2) as One can see that the conditions Tr[ρ t ] = Tr[ρ 0 ] = 1 are satisfied after noting that Tr[S a ] = 2 N δ a,0 and µ 0 t ′ = 0. This last equality follows from Eqs. (7) and (10) jointly with the property H 0b = 1 ∀b. By using the vectorial extension of Eq. (4), we get the density matrix written in a Kraus representation [1,2], The weights are p a , which from Eq. (10) can explicitly be written in terms of the time-dependent rates as 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 {γ a t } and the initial condition ρ 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 [1,2] implies the conditions 0 ≤ p a t ≤ 1, which means that {p a t } are a set of normalized probabilities. In the single qubit case (N = 1), previously obtained constraints are recovered [35]. In the general case, 4 N inequalities must be fulfilled. We notice that a sufficient, but not necessary, condition is In fact, this constraint implies that all eigenvalues, cf. Eq. (10), satisfy µ a t ≤ 0 (a = 0). Consequently, taking an arbitrary but fixed time t, the solution (11) of the non-Markovian dynamics, via the association t 0 dt ′ µ a t ′ = tµ a M , is equivalent to the solution of a (well behaved) Markovian dynamics generated by a Lindbladian with eigenvalues {µ a M }.
C. Non-Markovianity and time-dependent rates Different (inequivalent) memory witnesses based only on the system propagator can be used to define non-Markovianity [5,6] such as for example the trace distance between two different initial conditions [51] or those based on the k-positivity of the solution map [52]. Here, as the dynamics is written naturally in a canonical form [49], memory effects can also be defined by the negativity of the time-dependent rates {γ a t }. In this way, it is of interest to determine these elements for any well behaved solution defined by the probabilities {p a t } in Eq. (12). We can invert Eq. (13), and using Eq. (10) we get explicit expressions for the set of rates {γ a t } in terms of the normalized time-dependent weights 0 ≤ p c t ≤ 1, The signs of {γ a t } can be taken as a signature of departure from a Markovian regime [49]. 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. [35].

D. Additivity of non-Markovian master equations
Given two sets of (arbitrary) normalized probabilities {p a t } and {p a t }, the relation (15) allows us to obtain the corresponding sets of rates {γ a t } and {γ a t }. From these we can obtain a new master equation defined by Eq. (2) with rates {γ a t +γ a t }. In fact, it is always possible to associate a set of probabilities {q a t } to these added rates, that is, Consequently, as occurs to Markovian Lindblad equations [2], 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 which is valid for arbitrary Pauli strings S a and S b , and where S c = S a S b or equivalently S c = S b S a . Eq. (17) can be straightforwardly demonstrated from Eq. (5).

E. 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. [17], the more general case can always be described by writing the system density matrix ρ t and the probabilities of the incoherent system {q h t } as where the auxiliary states {ρ h t } correspond to the system state given that the extra (hidden) incoherent degrees of freedom are in the particular state h. The evolution of the states {ρ h t } may involve coupling between all of them [17]. Given the structure Eq. (2), each auxiliary state ρ h t must to assume the form where the parameter α runs over a set of Pauli strings that depends on each specific problem. The functions g h α (t) in turn obey a classical master equation whose structure also depends on each specific model. The initial conditions read ρ h 0 = ρ 0 q h 0 , where ρ 0 is the initial system state and q h 0 is the initial probability of the incoherent degrees of freedom. In fact, at time t, q h t = α g h α (t). On the other hand, the system density matrix evolution [Eq. (12)] is defined by the probabilities

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" [20,49] In order to provide simple (multipartite) examples, here we restrict to the case where the evolution is where S a and S b are two arbitrary multipartite Pauli strings, while S c = S a S b . Depending on the time-dependence of the rates we define what we term "hyperbolic" and "trigonometric" cases of eternal non-Markovianity.

A. Hyperbolic eternal non-Markovianity
The system density matrix is written as the addition of two auxiliary states ρ t = ρ The initial conditions for the auxiliary states are taken to be ρ This result provides a multipartite generalization, (N > 1), of the single qubit case (N = 1) studied in Ref. [49].
Similarly to the results of Ref. [20] we notice that in this particular case alternative dynamics such as the mapping to a classical master equation [see Eq. (B8)] and stochastic Hamiltonians [see Eq. (B13)] also lead to the same rates.

B. Trigonometric eternal non-Markovianity
Based on Eq. (18), instead of the evolution (21), here we consider The initial conditions are taken as ρ Taking into account Eq. (19), in order to solve Eq. (23) each auxiliary state is written as (h = 1, 2) where as before S c = S a S b , and g involves coupling between pairs of them. The corresponding solutions allow to obtain the probabilities a (t). Finally, the rates associated to the non-Markovian evolution follow from Eq. (15) Furthermore, where the coefficients are Depending on the ratio ϕ/γ, different characteristic behaviors are obtained. In Fig. 1 we plot both rates. Consistent with Eq. (25), γ From the plots it is also evident that γ Based on Eq. (22), we name this case as a trigonometric eternal non-Markovian. The probabilities {p a t } [Eq. (12)] also assume a simple form, These solutions apply to arbitrary multipartite Pauli strings a and b. family of well behaved dynamics. For example, we write

C. Adding non-Markovian evolutions
In this traslational invariant generator (say with periodic boundaries in one dimension), we may chose f i (t) = − tanh(γ i t) or alternatively f i (t) = tan(γ 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. (30) 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) = γ i /2 for i ≤ N 0 , and f i (t) = tan(γ 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 [46][47][48].
Here, we study a conditional past-future (CPF) correlation [47]. 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 x, y, and z, correspond to observations of three Hermitian operators S m with eigenvectors {|m } and eigenvalues {m}, The CPF correlation then reads [47] C pf (t, τ )| y = 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 τ 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) ↔ P (z, t + τ, 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, τ )| 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 h α (t)} [Eq. (19)]. Given that they obey a (linear) classical master equation, they can be written as where the set of functions {f hh ′ α (t)} are independent of the initial conditions {q h 0 }. Furthermore, for notational simplicity, we introduced a vectorial orthogonal base {|h)} for the incoherent degrees of freedom, such that f hh ′ α (t) ↔ (h|F α (t)|h ′ ) and q h 0 ↔ (h|q 0 ). The observables S m [Eq. (31)] may in principle be defined by arbitrary linear combinations of Pauli strings {S a }. Here, for simplicity they are defined by a unique Pauli string. In this case, the general expression for the CPF correlation [Eq. (D9)] reduces to (see Appendix D) In here, |q t ) = α F α (t)|q 0 ) define the probabilities of the incoherent degrees of freedom at time t, while (1| ≡ h (h|. Furthermore, x ≡ x xP (x) where P (x) = x|ρ 0 |x . Finally, P (y) is the probability for the outcomes of the second measurement. It is The term δ z,y δ y,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 x = S y = S z . This constraint does not emerge when the observables correspond to other basis of operators (see for example Ref. [48]).
When the three measurements are performed in direction c, we get 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 γ c t does not develop divergences [ Fig. 2(a)], the CPF correlation is negative for any value of the time intervals t and τ. On the other hand, in the interval 3− √ 8 < (ϕ/γ) < 3+ √ 8 where the rate γ c t 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 . 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 (34)]. These accidental cases can always be surpassed by considering arbitrary measurement operators written as linear combinations of the Pauli strings.

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.
Instead of the local-in-time formulation defined by Eq. (2), alternatively one may start with a time convoluted evolution where the set of time-dependent kernels {k a (t)} must to be constrained such that the solution map is CP. Similarly to Sec. II, by defining the kernel k 0 (t) ≡ − a (a =0) k a (t), here the weights of the solution (12) can be written as where the coefficients λ b (t) obey the evolution The inverse relations for determining the kernels {k a (t)} a function of probabilities {p a (t)} can be written in a 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. [20] for the case N = 1.

Mapping with a classical Markovian master equation
The solution map [Eq. (12)] is defined by a set of normalized probabilities {p a t }. It is possible to formulate an underlying mechanism such that {p a t } 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 {a}. Therefore, the system density matrix ρ t can be written in terms of a set of auxiliary states {ρ a t } [17] such that

Stochastic Hamiltonians
We consider a stochastic evolution, where the system wave vector |ψ t is driven by a stochastic Hamiltonian, The Hamiltonian H st is characterized by a noise with an arbitrary statistics but null average ξ α t = 0. The index α ↔ α 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 ρ α t = |ψ t ψ t | for a given α reads is the characteristic noise function for a given α. After averaging this parameter, the system state can be written The parameters {x α } correspond to the statistical weight of each Pauli string during the variation of the coefficient α. It is straightforward to check that ρ t = a p a t (S a ρ 0 S a ), which recovers Eq. (12) with Similarly to the previous model, it is not possible to find a general simple expression for the rates γ a t in terms of these probabilities. Manageable expressions arise in the following situations.
Particular cases: If the noise is the same for all "directions" G a t = G t , from Eqs. (15) and (B12), after some algebra [53], we get the rates where the scalar functions read and where ζ b is defined by Eq. (B9) with, instead of Eq. (B7), with p 0 ∞ = 1/2, and p a ∞ = x a /2. For a stationary Gaussian white noise, where ξ t ξ t ′ = Φδ(t − t ′ ), Eq. (B11) becomes G(t) = exp(−Φt). It is simple to check that in this situation Eq. (B13) recovers the solution (B8) of the previous model with ϕ = φ. 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 η, the characteristic noise (B15) In contrast to the previous cases, here the rates defined by Eq. (B13) may develop divergences. In fact, the functions (B14) become . (B16) Hence, divergent rates are found whenever η < A.

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 where {q k } are normalized positive weights ( n k=1 q k = 1), and each set of probabilities {p a k (t)} is associated to a Markovian solution of Eq. (2) with time-independent positive rates {γ a k }. From Eq. (15), the non-Markovian evolution is characterized by the rates where µ b k are eigenvalues of the k-Markovian dynamics, µ b k = c H bc γ c k . 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 a t = q 1 p a 1 (t) + q 2 p a 2 (t), where each solution is associated to the rates γ a 1 and γ a 2 , and q 1 + q 2 = 1. From Eq. (B18), after some algebra [53], we get where the parameters are 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. (B17) is rewritten as where we used the explicit expression (13) and the replacement n k=1 q k → · · · . The symbol · · · denotes an average over the set of random rates {γ c }, each "realization" defining a Markov solution. Assuming that all rates are independent random variables it follows that · · · → and P (x|y) = z P (z, x|y). From Eq. (D6), and using z z| z|S α |y | 2 = y|S α S z S α |y , (D7a) x x| y|S β |x | 2 P (x) = y|S β S x ρ x S β |y , (D7b) x | y|S β |x | 2 P (x) = y|S β ρ x S β |y , where the system state ρ x is the CPF correlation can be written as C pf (t, τ )| y = 1 P (y) 2 α,β,γ Θ αβγ | y Λ αβγ (t, τ ). (D9) The coefficients Θ αβγ | y are Θ αβγ | y = y|S α S z S α |y y|S β S x ρ x S β |y y|σ γ ρ x S γ |y , while the time-dependence follows from where |q t ) = α F α (t)|q 0 ), and the probability P (y) is The expression (D9) is valid for arbitrary observables σ m [Eq. (31)]. In general, they can be written as linear combinations of Pauli strings S a . Assuming, for simplicity, that each S m correspond to a unique Pauli string operator, from Eq. (5) it follows the relations y|S α S z σ α |y = H αy δ z,y a y , y|S β S x ρ x S β |y = 1 2 N (H βy δ y,x a y + x ), (D11b) y|S γ ρ x σ γ |y = 1 2 N (1 + H γy δ y,x a y x ), (D11c) where x ≡ Tr[S x ρ x ]. By introducing these equalities in Eq. (D9), after some algebra we get Eq. (34). Generalization to arbitrary observables can be worked out in a similar way from Eq. (D9).