Multipartite Correlations in Quantum Collision Models

Quantum collision models have proved to be useful for a clear and concise description of many physical phenomena in the field of open quantum systems: thermalization, decoherence, homogenization, nonequilibrium steady state, entanglement generation, simulation of many-body dynamics, and quantum thermometry. A challenge in the standard collision model, where the system and many ancillas are all initially uncorrelated, is how to describe quantum correlations among ancillas induced by successive system-ancilla interactions. Another challenge is how to deal with initially correlated ancillas. Here we develop a tensor network formalism to address both challenges. We show that the induced correlations in the standard collision model are well captured by a matrix product state (a matrix product density operator) if the colliding particles are in pure (mixed) states. In the case of the initially correlated ancillas, we construct a general tensor diagram for the system dynamics and derive a memory-kernel master equation. Analyzing the perturbation series for the memory kernel, we go beyond the recent results concerning the leading role of two-point correlations and consider multipoint correlations (Waldenfelds cumulants) that become relevant in the higher-order stroboscopic limits. These results open an avenue for the further analysis of memory effects in collisional quantum dynamics.


Introduction
The standard collision model, introduced as early as in 1963 in Ref. [1], considers a quantum system that sequentially interacts with identical uncorrelated ancillary particles or oscillator modes. Each elementary system-particle interaction lasts for a finite period of time τ and is described by an elementary unitary evolution operator U; see Figure 1. However simple this model may look like, it (i) naturally describes the system dynamics induced by repeated interactions, e.g., in the micromaser theory [2]; (ii) gives an intuitively clear picture of various phenomena, such as thermalization [3,4], decoherence [4][5][6], homogenization [7,8], nonequilibrium steady state [9][10][11], and entanglement generation [11][12][13]; and (iii) is amenable to analytical treatment, which makes it possible to derive time-continuous master equations in appropriate limits on the system-environment interaction strength and the collision duration [14][15][16] (in the standard collision model, the system dynamics are Markovian and completely positive divisible due to a pastfuture independence of ancillary particles [17]). Ideas of repeated interactions underlie the discrete-time open quantum walks and their continuous-time limit [18][19][20][21][22][23]. Hence, it is no wonder that quantum collision models are becoming increasingly popular in quantum information, quantum technology, and mathematical physics communities. Mysteriously, the quantum physics community and mathematical physics community do not know much about each other and sometimes conduct rather isolated research on highly interrelated topics. Mathematical physicists usually refer to the standard quantum collision model as the repeated interaction model and treat it as a particular model of non-equilibrium quantum statistical mechanics [24]. In addition to the derivation of the master equation, the interest of mathematical physicists is also focused on the asymptotic state in the limit of large times [25,26] and the study of random repeated interactions [27,28]. On the other hand, quantum physicists find new applications of quantum collision models in simulations of open quantum many-body dynamics [29,30] (including simulations on noisy intermediate-scale quantum processors [31]), relaxation processes caused by the dilute gas environment [32], quantum thermodynamics [33], and quantum thermometry [34,35]. The collisional picture of repeated interactions also takes place in quantum optics and waveguide quantum electrodynamics, where the electromagnetic field is represented in the form of discrete time-bin modes interacting with a quantum emitter [36][37][38][39][40][41][42][43][44][45][46][47][48][49]; however, the time-bin modes constituting the radiation field can be correlated so that the system dynamics becomes non-Markovian and exhibits memory effects in general. Besides the initially correlated state of ancillary particles or modes [46][47][48][49][50][51], memory effects in quantum collision models also appear as a result of two-ancilla collisions in between the system-ancilla collisions, where the latest involved ancilla interacts with the one that would interact with the system during the next collision [52][53][54][55]. An alternative scenario for non-Markovian dynamics (e.g., due to random telegraph noise) assumes that a system is composed of the very open system under study and an auxiliary sybsystem, which alternately interacts with a fresh reservoir ancilla and the system under study [56]. Another approach considers repeated interactions of the system with the particles it has already collided (including many-body collisions) [57][58][59]. Quantum channels with memory can also be viewed in terms of quantum collision models [60][61][62][63][64][65][66][67]. The presented list of possible modifications for quantum collision models is far from being complete; in this regard, we refer the interested reader to the recent review papers [68,69]. Nonetheless, the reader can see a great flexibility of quantum collision models to describe a variety of physical situations in a rather simple way.
One of the current challenges in the standard collision model is related to quantum correlations among ancillas that are induced by successive system-ancilla interactions. These correlations lead to an advantage in the collisional quantum thermometry [34]. However, a direct numerical simulation of the output ancillas' state is possible for a relatively small number n of ancillas because of an exponentially growing dimension, d n , for the state of d-dimensional ancillas. This includes, for instance, d = 2 and n ≤ 12 in Ref. [34]. Another challenge appears if the ancillas are initially correlated. This scenario takes place, e.g., when the second system starts interacting with an array of ancillas that were originally uncorrelated but previously interacted with the first system in the standard collision model [65,66]. Alternatively, the ancillas can represent timebin correlated modes in the structured electromagnetic radiation [36][37][38][39][40][41][42][43][44][45][46][47][48][49] or particles in a correlated spin chain, e.g., spin-1 particles in the ground state of the Affleck-Kennedy-Lieb-Tesaki (AKLT) antiferromagnetic Hamiltonian [70]. Ref. [71] reports that the correlations can break convergence of the system state to the same state of all locally identical ancillas (such a convergence-known as homogenization-would have taken place under appropriate conditions, were the ancillas uncorrelated). Again, the exponential increase in Hilbert-space dimension limits the numerical study in Ref. [71] to 16 ancillas. Therefore, we face a general problem of how to deal with correlations among ancillas (either induced by the system or initially present).
The first goal of this paper is to represent the system-induced correlations among ancillas (in the standard collision model) by developing the tensor network formalism applied recently in Ref. [72]. The main idea behind the tensor network representation (in the form of the matrix product state [73][74][75][76]) is that many n-partite states of d-dimensional ancillas require only about ndr 2 complex parameters to be specified, not d n parameters. As we show in this paper, r equals the system dimension in the standard collision model. Our second goal is to develop the ideas of Ref. [72] and derive a more general master equation for the system dynamics in the nonstandard collision model with an initially correlated environment. The point of Ref. [72] is that two-point correlations among ancillas play a leading role in the system dynamics if each elementary unitary evolution operator slightly deviates from the identity operator. However, it may happen that the leading contribution vanishes for a specific interaction, and we demonstrate such an example in this paper. Therefore, one needs to consider higher-order correlations among ancillas and their effect on the system dynamics. We close this gap and provide a recipe for how to derive a master equation valid in the corresponding perturbation order for the elementary unitary evolution operator.

Tensor Network Notation
Tensor network representation of quantum states is reviewed in a number of papers [73][74][75][76][77][78] and a book [79]. Consider a pure state |ψ of n particles, where each particle is associated with a Hilbert space H of a finite dimension d. The state is fully defined by d n complex numbers C i 1 i 2 ...i n in the decomposition where {|i k } i k =1,...,d is an orthonormal basis in H. A collection of d n complex numbers {C i 1 i 2 ...i n } can be viewed as a rank-n tensor C with a picture representation involving a letter "C" with n legs. To distinguish the ket-vector |ψ from the bra-vector ψ|, we add arrows to the legs; namely, we associate outcoming arrows with ket-vectors and incoming legs with bra-vectors. A tensor diagram concisely depicts a contraction of tensors: the connected lines are summed over. The tensor diagram for an n-partite matrix product state (MPS) with open boundary conditions contains n tensors A [1] , . . . , A [n] connected in a line; see Figure 2. A [1] and A [n] are rank-2 tensors with elements A [1],i 1 a 1 and A [n],i n a n−1 , respectively, whereas for all k = 2, . . . , n − 1 the tensor A [k] has rank 3 and is composed of elements A [n],i n a n−1 ,1 , respectively. Arrows in Figure 2 also indicate the order for multiplication of matrices. The contraction yields which explains the MPS name. A number |{a k }| of the values that the virtual index a k can take is not related to the physical dimension d of the particles. We will refer to the maximal number max k=1,...,n−1 |{a k }| as the bond dimension. Clearly, the MPS representation for a given state |ψ is not unique in general; however, the less the bond dimension, the easier the calculations and the analysis. In view of this, the minimal bond dimension among all possible MPS representations is called the MPS rank and denoted by r. The greater r, the more entangled the state |ψ can be with respect to the left-right bipartitions [80].
Arrows in tensor diagrams simplify their interpretation. For instance, changing the direction of arrows from left to right in the connecting lines in Figure 2, we get the transposed matrices (A [k],i k ) , k = 1, . . . , n. The resulting diagram is depicted in Figure 3. Nonetheless, if indices i 1 , . . . , i n are fixed, then the c-number C i 1 i 2 ...i n does not change because

Matrix Product State Correlations in the Standard Collision Model
We begin with the simplest scenario, in which the system is initially in a pure state |ϕ ∈ H S , dimH S = d S , and the environment consists of n uncorrelated ancillas in a pure state |ψ 1 ⊗ |ψ 2 ⊗ . . . ⊗ |ψ n ∈ H ⊗n , dimH = d. Each elementary collision is described by a unitary operator U : H S ⊗ H → H S ⊗ H, which is viewed as a 4-rank tensor. After n collisions, the system and ancillas get entangled, and their composite state is given by a tensor diagram in Figure 4. As a result of n collisions, we get a correlated state of n + 1 particles: n ancillas and one system particle. Dotted lines in Figure 4 denote tensors that should be contracted to get the matrix product state structure. The rightmost dotted region depicts an identity operator I. Clearly, the bond dimension equals the number of the system degrees of freedom, d S . Therefore, we can associate each virtual index a k with a vector |a k ∈ H S , so that a collection of vectors {|a k } for a fixed k forms an orthonormal basis in H S . The very diagram in Figure 4 serves as the proof for the following result. Proposition 1. Let the system and n ancillas be initially in the pure states |ϕ , |ψ 1 , . . . , |ψ n . Then the output state |Ψ ∈ H ⊗n ⊗ H S of the system and ancillas in the standard collision model with the elementary unitary operator U adopts an MPS representation where A [1],i 1 1,a 1 = a 1 | ⊗ i 1 | U |ϕ ⊗ |ψ 1 , A [k],i k a k−1 ,a k = a k | ⊗ i k | U |a k−1 ⊗ |ψ k for all k = 2, . . . , n, and A [n+1],i n+1 a n ,1 = δ a n ,i n+1 .
The result of Proposition 1 explains the previously known observations of Ref. [8] that the partial swap interactions (U = exp[−igτ ∑ i,j |ij ji|]) generate the W-type of entanglement, whereas the controlled unitary interactions (U = ∑ i U i ⊗ |i i|) generate entanglement of the Greenberger-Horne-Zeilinger (GHZ) type. In fact, both W and GHZ states of many qubits are particular forms of the matrix product states with the bond dimension 2 [73][74][75][76][77][78]. Example 1. Let the system and ancillas be qubits. The system is initially in the excited state |ϕ = |↑ . Each ancilla is initially in the ground state, i.e., |ψ k = |↓ for all k.
The explicit relation between the unitary operator U and tensors A [k] , which we establish in Proposition 1, enables one to approach the quantum engineering problem too. Suppose one wants to create an entangled state |Ψ of n particles that adopts a matrix product representation with the bond dimension r. Then one needs to take an r-dimensional quantum system and let it sequentially interact with the initially uncorrelated particles. Finally, one performs a projective measurement on the system in the basis {|i n+1 } to get rid of its degrees of freedom. The resulting state of n particles is i n+1 |Ψ , where i n+1 is the measurement outcome and |Ψ is given by Equation (4). The unitary operator U should be optimized in such a way as to maximize the overlap |( Ψ | ⊗ i n+1 |) |Ψ | 2 . Clearly, each collision could be described by its own unitary operator. Then, one should replace U → U k in the formula for A [k] in Proposition 1. Numerical tools for optimization over many unitary operators {U k } n k=1 are presented, e.g., in Refs. [81,82]. Let us consider entanglement of the state |Ψ with respect to a bipartition into ancillas 1, . . . , k on one side and ancillas k + 1, . . . , n and the system on the other side, i.e., the left-right bipartition with the boundary in between the ancillas k and k + 1. Entanglement of a pure state with respect to a bipartition is quantified by the entanglement entropy that equals the von Neumann entropy of either reduced density operator, The reduced density operator 1...k = tr k+1,...,n+1 |Ψ Ψ| for k ancillas is presented in the form of a tensor diagram in Figure 5a. The tensor diagram in Figure 5a gets simpler if we take into account the following important property (referred to as the right-normalization condition [75]): Here, the overline denotes the complex conjugation and † denotes the Hermitian con- 1,a n = δ a n a n because A [n+1],i n+1 a n ,1 = δ a n ,i n+1 by Proposition 1. If m = 2, . . . , n, then Proposition 1 implies If m = 1, then we deal with dummy indices a 0 = a 0 = 1 and ∑ i 1 ,a 1 A Hence, we have proved the following result. (4) satisfies the right-normalization condition (5).

Proposition 2. MPS |Ψ in Equation
An MPS satisfying the right normalization condition is also called right-canonical [75]. The advantage of the right-canonical form is that the partial trace over the rightmost particles corresponds to a single connecting line in the tensor diagram; see Figure 5b.
,i m = I, which is exactly the vertical connecting line in Figure 5b. Physically, the reduced density operator for k ancillas does not depend on future system collisions with other ancillas that happen after time kτ.
Entanglement entropy E(Ψ) of the state |Ψ with respect to the cut in between the ancillas k and k + 1 reads Note that E(Ψ) ≤ log d S because d S is an upper bound for the Schmidt rank of |Ψ .

Generalization to Mixed States of the System and Ancillas
Let us consider the standard collision model, where the system and ancillas are generally mixed. This scenario is especially relevant to the task of quantum thermometry [34,35]. The initial state of the system is given by the density operator S . The initial state of n ancillas is given by a factorized density operator n k=1 k . Collisional dynamics with the elementary unitary operator U drives the system and ancillas to the state where the subscript Sk in the notation U Sk means that U nontrivially acts on the system and the k-th ancilla. A tensor diagram for Equation (8) is presented in Figure 6a. Dotted regions in Figure 6a show tensor contractions or tensor combinations that effectively lead to the equivalent tensor diagram depicted in Figure 6b. Note, however, that the arrows in the upper horizontal line in Figure 6a Let S = ∑ l λ l S |ϕ l S ϕ l S | be the spectral decomposition for the system initial state. Let k = ∑ m λ m k |ψ m k ψ m k | be the spectral decomposition for the initial state of the k-th ancilla. Then one readily obtains the representation A tensor diagram for Equation (14) is depicted in Figure 6c. We see that for any Figure 6b,c define the so-called matrix product density operator (MPDO) [83,84], which is automatically Hermitian and positive semidefinite. MPDOs are successfully used to study the dissipative dynamics and the Gibbs states of onedimensional quantum chains [83][84][85][86]. Among other questions, Ref. [86] addresses an important question how to prepare MPDO states experimentally. Our results show one more method to prepare an MPDO state via the standard collision model. In our construction, the MPDO is right canonical, i.e., it additionally satisfies the right-normalization condition Equation (16) mathematically shows the independence of the reduced density operator 1...k (kτ) for k ancillas from future collisions at times t > kτ. The results of this section are summarized as follows.
The main benefit of the constructed MPDO representation is that it exploits only parameters needed for a description of a general state of the system and n ancillas. In other words, computational resources scale linearly (not exponentially) with the number of ancillas if one uses the MPDO representation. This fact opens an avenue for further numerical studies in the collisional quantum thermometry [34,35]. If the system interacts with a thermal reservoir in between the collisions with ancillas, one can readily include such a system-reservoir interaction in the tensor network representation in the form of a quantum channel [87]. Example 2. Let the system and ancillas be qubits. The system is initially in the excited state |ϕ = |↑ . Each ancilla is initially in the Gibbs state where k B is the Boltzmann constant, T is the temperature, and E ↑ and E ↓ are the energy levels for the ancilla states |↑ and |↓ , respectively. Consider the energy exchange unitary U = exp[gτ(|↓↑ ↑↓| − |↑↓ ↓↑|)]. After n collisions, the mixed state of the system and ancillas is fully described by a right-canonical MPDO with D = 2. The explicit form for this MPDO is given by Proposition 3 and reads cos gτ 0 0 1 , where k = 2, . . . , n.

Collision Model with a Generally Correlated State of Ancillas
Let us consider a more complicated collision model, in which ancillas are initially correlated. Surprisingly enough, any pure state of n ancillas adopts an MPS representation [73][74][75][76][77][78]. However, the MPS rank for a generally correlated state grows exponentially with n. On the other hand, many important states of correlated ancillas such as few-photon wavepackets [46][47][48][49], artificial photonic tensor network states [37,[88][89][90][91][92][93], and ground states of gapped one-dimensional local Hamiltonians for the spin chains [94] are described by MPSs with a low MPS rank. As the states of ancillas are mixed in general, we exploit the MPDO formalism. We pay little attention to the rank of decomposition as our further goal is to reveal the effect of ancillas' correlations on the system dynamics. Note that the correlations can be either quantum (genuinely entangled ancillas) or classical (fully separable state of ancillas); however, both types strongly affect the system dynamics (see an example in Ref. [50]).
Let the intial state 1...n be a right-canonical MPDO for n ancillas shown in Figure 7a. Here, we have added a formal density operator χ 0 (i.e., a positive semidefinite operator with unit trace) for the bond degrees of freedom (blue arrows in Figure 7a). In the conventional MPDO notation, χ 0 is the trivial 1 × 1 identity matrix for dummy indices; however, in our construction, it can be an arbitrary density matrix such that the tensor contraction is well defined. Note that we changed the direction of arrows in the upper line in Figure 7a. This implies transposition of matrices B [k],i k b with respect to horizontal virtual indices. Figure 7b illustrates that the tensors {B [k] } k define a "free evolution" for the bond degrees of freedom if ancillas do not interact with the system; namely, the matrix is a valid density matrix (i.e., χ † k+1 = χ k+1 ≥ 0 and tr[χ k ] = 1) provided χ k is also a density matrix. This follows from the right normalization condition (16). Figure 7c depicts the system density operator S (kτ) after k collisions. The partial trace over ancillas k + 1, . . . , n corresponds to a vertical connecting line for the bond degrees of freedom (blue arrows in Figure 7c). The partial trace over ancillas 1, . . . , k corresponds to vertical connecting lines for the ancillary degrees of freedom (green arrows in Figure 7c). Ref. [72] discusses the natural Markovian embedding for the system dynamics that follows from the diagram in Figure 7c. In our case, we have The trace preserving property ∑ j m b K † j m b K j m b = I follows from the right normalization condition (16) and unitarity of U.
The tensor diagram in Figure 7c is a particular form of the process tensor-a recently developed approach to an operational description of non-Markovian quantum dynamics [59,[95][96][97][98]. The complexity of the non-Markovian dynamics simulation depends on the dimension of the effective reservoir in the Markovian embedding [99,100]: the less the dimension of the Markovian embedding, the simpler the simulation. In our model, the role of the effective reservoir is played by the bond degrees of freedom that specify correlations among the ancillas.
The emergence of non-Markovian dynamics in the case of correlated ancillas was demonstrated in Ref. [50], where an exemplary indecomposable qubit channel was realized as a result of a qubit's collisional interactions with many qutrit ancillas in the GHZ state. The analytical treatment in Ref. [50] was only possible due to a peculiar controlled-unitary qubitancilla interaction. Were the qubit-ancilla interaction different from the controlled-unitary type, the methods of Ref. [50] would not provide any analytical expression for the qubit system dynamics (nor would it be possible to study its non-Markovianity). As we show in the example below, the developed tensor network formalism enables us to resolve that difficulty and analytically derive the qubit dynamics even for non-controlled unitary collisions. Since any environment state adopts an MPDO form, our results generalize those of Ref. [51], where non-Markovian qubit dynamics are induced by a specific correlated environ- where either i k = 0 and i k = 1, or i k = 1 and i k = 0; {k x }, {k y }, {k z } are nonintersecting subsequences of collision numbers. The latter environment reproduces an arbitrary Pauli dynamical map [51].
The qubit system is initially in the state S ≡ S (0). Each qubit-qutrit collision is described by the unitary operator where (σ 1 , σ 2 , σ 3 ) ≡ (σ x , σ y , σ z ) is the conventional set of Pauli operators, and (J 1 , J 2 , J 3 ) ≡ (J x , J y , J z ) is a set of SU(2) generators for a qutrit (spin-1 particle). In the conventional orthonormal basis (|1 , |2 , |3 ), the corresponding matrices are Substituting Equations (20) and (21) into Equation (19), we get the map E [m] ≡ E , which does not depend on the collision number m. Then Equation (18) results in the following qubit system density operator after k collisions: where λ(k) is a scaling coefficient for the x and y components of the qubit Bloch vector, and λ z (k) is a scaling coefficient for the z component of the qubit Bloch vector. The center of the Bloch ball is a steady point under such dynamics. The explicit formulas for the scaling coefficients are We portray the typical behaviour of λ(k) and λ z (k) in Figure 8. Whenever |λ(k)| increases with increasing k, we observe a positive indivisible dynamic, which is often treated as an indication of essential non-Markovianity [101]. We refer the interested reader to Ref.

Master Equation
Equation (18) defines the discrete dynamical map Υ kτ that transforms the initial system density operator S (0) ≡ S to the system density operator S (kτ) after k collisions, i.e., Υ kτ [ S (0)] = S (kτ). If τ → 0, then we interpret kτ as a continuous time t. A time-local master equation Although such a master equation correctly describes the system evolution, it conceals the major role of correlations among ancillas. To reveal the physics of how these correlations affect the system dynamics, we resort to the conventional projection operator techniques [103] and derive the Nakajima-Zwanzig memory kernel equation [104,105] for our model.
Let {χ k } k be a collection of the density operators for the bond degrees of freedom generated by Equation (17). For each k, define the following map P k acting on both the system and the bond degrees of freedom: Since tr[χ k ] = 1, we have P 2 k = P k , so P k is a projection. A pictorial representation of the projection P k is given in Figure 9a, which shows that P k breaks the left-right correlations between bunches of ancillas (1, . . . , k) and (k + 1, . . . , n). Clearly, A complementary projection Q k is defined through where Id is the identity transformation for the system and bond degrees of freedom. Q k is also a projection because Q 2 k = Q k . We can also rewrite Q k in the form where Q bond#k is a projection for the k-th bond degrees of freedom that acts on an operator F for the bond degrees of freedom as follows: Since Q bond#0 [χ 0 ] = 0, we readily get To simplify the notation, let us introduce the system-bond density operator after the k-th collision, Applying Q k to the both sides of the latter equation, we get The recurrent Equation (32) with the initial condition (31) has the following formal solution: If we apply P k+1 to the both sides of equation Recalling the relation P k [R(kτ)] = S (kτ) ⊗ χ k and taking partial trace over the bond indices in Equation (34), we get Subtracting S (kτ) from both sides of Equation (35) and dividing the result by the collision time τ, we get a discrete-time version of the celebrated Nakajima-Zwanzig master equation, namely, where the memory kernel K km relates the density matrix increment (in between the times kτ and (k + 1)τ) with the past density operator at time (k − m)τ. If m = 0, then we have a time-local term K k0 giving the density operator increment caused by the latest collision (among those that have already happened): with k+1 being a reduced density operator for the (k + 1)-th ancilla in the initial state; see Figure 9b. If m ≥ 1, then K km describes a nontrivial effect of preceding collisions on the system evolution and reads If there were no correlations in the environment, then K k0 would be the only contribution to the kernel because K km would vanish for all m ≥ 1. Indeed, the MPDO rank equals 1 for a factorized environment state, so dimH bond#k = 1 for all k. Each χ k is unambiguously defined because χ k is the trivial 1 × 1 identity matrix in this case, and Q bond#k [F] = 0 for any 1 × 1 matrix F. If the environment is correlated, then the memory contribution K km [ S ] = 0 in general.

Example 4.
Consider the GHZ state of 3-dimensional ancillas 1...n = |GHZ GHZ|, |GHZ = 1 √ 3 ∑ j=1,2,3 |j ⊗n and the controlled unitary system-ancilla interaction U = ∑ j=1,2,3 e −igτσ j ⊗ |j j|, where gτ quantifies the dimensionless system-ancilla interaction strength and (σ 1 , σ 2 , σ 3 ) ≡ (σ x , σ y , σ z ) is the conventional set of Pauli operators. This is a scenario considered also in Ref. [50]. A direct calculation yields The memory kernel K km does not decay with the increase of m due to the infinite correlation length in the GHZ state. In view of this, even if the interaction strength gτ 1, one cannot truncate a series expansion for K km with respect to a small parameter gτ. Instead, all orders of gτ are significant for reproducing the system dynamics.
If the correlation length is finite, then it is possible to derive a continuous-time master equation in the appropriate limit for τ and g. This is discussed in what follows.

Effect of Two-Point Correlations
An important simplification comes from a series expansion for U Sm with respect to the interaction strength gτ between the system and an individual environment particle. Let ghH m be the system-particle interaction Hamiltonian during the m-th collision, wherē h is the reduced Planck constant, g has the physical dimension of frequency, and H m is a dimensionless Hermitian operator with the operator norm H m ≤ 1. Then, the elementary unitary interaction in the m-th collision is U = exp(−igτH m ). The map E [m] in Equation (19) has a contribution of both U and U † , so we have The term K see Equation (30). Physically, the 0-th order of Φ [k] i k ,i k involves no system-environment interaction and, consequently, no contribution to the memory kernel.
To calculate the term K Id S , then we have a zero because is a trace preserving map due to the right-normalization condition, whereas Q bond#k nullifies the trace of any operator (see Equation (30)). Therefore, the term K (1) km also vanishes.
Similar considerations for the term K km lead to a conclusion that K km may be nonzero only if we fix Φ where the coefficient C Recalling the definition (30), we get Tensor representation in Figure 9c justifies that i.e., we get the reduced density operator for the (k − m + 1)-th ancilla and (k + 1)-th ancilla in the initial correlated state of ancillas. Similarly, i.e., we get a tensor product of individual reduced density operators for the (k − m + 1)-th ancilla and (k + 1)-th ancilla in the initial correlated state of ancillas. Combining (47)-(50), we obtain a surprisingly simple though exact result, namely, Introducing the environment two-point correlation function for operators O and O by a conventional formula we readily see that C Combining all the findings of this section, we get Equation (54) provides an important physical link between the two-point correlation function of ancillas and the memory kernel.

Stroboscopic Limit
To simplify the analysis, let us assume that the correlated ancillas are initially in the homogeneous right-canonical MPDO, i.e., the tensors B [k] coincide for all k = 1, . . . , n, and MPDO is fully described by the density matrix χ 0 and the tensor M. In this case, all local density operators for individual ancillas coincide 1 = . . . = n ; however, the two-ancilla density operator 12 = 1 ⊗ 2 . If ancillas are initially in such a homogeneous state, we can expect that the kernel K km depends on m only and does not depend on k.
Suppose the collision duration τ tends to zero while the coupling strength g remains constant. Then we get the Hamiltonian dynamics for the system S (t) in continuous time [66]. The correlations among ancillas are irrelevant in this scenario because g 2 τ → 0 in the considered limit.
To reveal a nonunitary system dynamic at a long timescale, one should consider a different limit gτ → 0, g 2 τ = const [56,65,66], which we refer to as the (first-order) stroboscopic limit that is also used in the analysis of dynamics induced by indirect repeated measurements [106,107]. The Hamiltonian part −ig[ H anc , S (t)] explodes in the master equation because g → ∞; however, this problem disappears in a proper interaction picture [66].
In the stroboscopic limit, one cannot simply replace 1 , and the second summand cannot be neglected. However, if the expression −ig[ H anc , S (t)] vanishes, then the characteristic frequency of system dynamics is g 2 τ so that 1 . The second summand vanishes in the (first-order) stroboscopic limit because g 4 τ 3 = (g 2 τ) 2 τ → 0. The time-local memory kernel component should also be considered in this limit, i.e., The higher-order contributions K km , K km , . . . to the memory kernel (43) vanish only if the correlation length l corr (in the chain of ancillas) is finite. If this is the case, then l corr g N+1 τ N → 0 for N = 3, 4, . . .. Therefore, in the first-order stroboscopic limit, we get the time-continuous master equation The correlations are known to decay exponentially in an MPS and an MPDO [73][74][75][76][77], with the correlation length l corr being defined by the second-largest eigenvalue of the transfer matrix T = ∑ i M ii (in absolute values). If the correlation length is finite, then K m represents a sum of exponentially decaying terms, where {λ j } j are eigenvalues of the transfer matrix T such that |λ j | < 1 and {L (j) nonlocal } are the corresponding maps.
To explicitly find the kernel K(t ) in Equation (56), we resort to the Laplace transform (which is often used for the memory kernel master equations [108][109][110]); namely, The result does not depend on s, which means the kernel K(t ) becomes local in the stroboscopic limit and the final master equation takes the form Physically, Equation (61) shows that if the system quickly interacts with ancillas (gτ 1), then the system "feels" not only the individual ancillas (which results in the local term L local ) but also a somewhat averaged correlated state (which results in the nonlocal nonlocal ). We summarize these results as follows.
Proposition 4. Let the system collisionally interact with an array of ancillas in the homogeneous MPDO with a finite correlation length. If the expression −ig[ H anc , S (t)] vanishes, then in the first-order stroboscopic limit gτ → 0, g 2 τ = const, the system dynamics are governed by the master equation (61), where the local and nonlocal contributions to the generator are defined by Equations (55), (58) and (59).

Example 5.
Consider an infinite chain of spin-1 particles (ancillas) in the AKLT state [70] that adopts the following homogenous right canonical MPS representation with the MPS rank 2 [75]: Each individual ancilla has a reduced density operator 1 = 1 3 I; however, the global state is correlated. At time t = 0, the qubit system collides with one of the intermediate ancillas, then collides with its right neighbor and so on. Each collision lasts time τ. The system-particle interaction Hamiltonian ish Averaging over single-ancilla degrees of freedom yields H anc = 1 3 ∑ j=1,2,3 σ j = 0. To use Proposition 4, we set S (0) = 1 2 (I + 1 √ 3 ∑ j=1,2,3 σ j ), so that −ig[ H anc , S (0)] = 0. As we will see later, the latter commutation relation remains valid for all times t, i.e., −ig[ H anc , S (t)] = 0, and the use of Proposition 4 is justified.
The local term is given by formula (55) and reads To find the nonlocal term, we should take correlations into account. Since the system interacts with a part of the infinite spin chain, the state 1...∞ of ancillas (spin-1 particles) is mixed and described by a right-canonical homogeneous MPDO with The transfer matrix reads and has eigenvalues 1 (of multiplicity 1) and − 1 3 (of multiplicity 3). The two-spin reduced density matrix reads where J α is an operator for the spin projection (in units ofh) in the α direction, α = x, y, z; see Equation (22). Substituting Equation (67) in Equation (58), we get On the other hand, i.e., in our case we have a single contribution with λ = − 1 3 and Finally, Equation (61) gives the explicit master equation in the stroboscopic limit The reader may notice the formally negative rate − g 2 τ 3 in front of the second dissipator term; however, the total generator does have the Gorini-Kossakowski-Sudarshan-Lindblad form [111,112] because the Kossakowski matrix is positive semidefinite. Therefore, the actual relaxation rates are positive. The effect of the formally negative rate − g 2 τ 3 in front of the second dissipator term in Equation (70) is that correlations among ancillas slow down the relaxation as compared to the case of uncorrelated ancillas (equation Figure 10 illustrates this phenomenon. Figure 10 also shows a good agreement between the exact system dynamics and the system dynamics in the stroboscopic limit and emphasizes the role of correlations. Disregard of correlations among ancillas leads to a wrong result (see the dashed line in Figure 10).

Effect of Multipoint Correlations in the Higher-Order Stroboscopic Limit
We start with an example stimulating the discussion of the higher-order stroboscopic limit. Example 6. Consider a qubit system interacting with an infinite chain of spin-1 particles (ancillas) in the AKLT state as in Example 5, with the only difference that the system-ancilla interaction Hamiltonian now readsh For such an interaction, H anc = 0, so the use of Proposition 4 is justified. Following the lines of Example 5, we similarly calculate L local and L nonlocal in the first-order stroboscopic limit; however, these contributions cancel each other so that the right hand side of Equation (61) vanishes in the first-order stroboscopic limit, and we get d S (t) dt = 0. The exact treatment of the problem via Equation (18) yields the depolarizing system dynamics with the depolarization function x = 2 + 7 cos 3gτ 2 , y = 7 + 2 cos 3gτ 2 , z = 2 y 2 + 27 sin 2 3gτ 2 .
A feature of the depolarizing dynamical map (72) is that it is neither completely positive divisible nor positive divisible for all gτ ≤ 2 3 arccos(− 11 16 ) because q(τ) ≥ 0 and q(2τ) − q(τ) = 2 5 y 3 6 sin 2 3gτ i.e., the image of the system Bloch ball shrinks after the first collision and then expands after the second collision. If gτ = 4πm/3, m ∈ N, then the system experiences no evolution, i.e., S (t) = S (0). If gτ = 2π/3, then the Bloch ball experiences partial inversion with the scaling parameter − 5 27 after each collision. The latter dynamics are, however, completely positive divisible. If gτ 1, then q(t) ≈ exp(− 1 8 g 4 τ 3 t), so the characteristic frequency of the system dynamics is g 4 τ 3 . The first-order stroboscopic limit is unable to reproduce such a behaviour because it is only sensitive to rates ∼ g 2 τ. This example stimulates us to develop (to some extent) the theory of the higher-order stroboscopic limit.
We will refer to the limit gτ → 0, g n+1 τ n = const as the n-th order stroboscopic limit. Surely, the expressions g, g 2 τ, . . . , g n τ n−1 explode in this limit; however, if their contribution to the system dynamics vanishes, then the limit is well defined. The higherorder contributions g n+2 τ n+1 , g n+3 τ n+2 , . . . vanish in the n-th order stroboscopic limit. In Equation (55), for L local , we should throw away the exploding terms (as they will cancel other exploding terms from the memory kernel) and vanishing terms. In the memory kernel series expansion (43), we should keep the term K (n+1) , which describes (n + 1)point correlations among ancillas. Recalling the notation in Equation (41), we get, for instance, the following expression for the third-order memory-kernel: where C with C(O, O , O ) being the third-order Waldenfelds cumulant [113,114]. In the secondorder stroboscopic limit, the expression g 3 τ 2 K km reduces to the time-local generator g 3 τ 2 L nonlocal , which contributes to the final GKSL master equation d S (t) dt = L local [ S (t)] + g 3 τ 2 L nonlocal [ S (t)]. The higher-order Waldenfelds cumulants are expressed through the lower-order ones [114], thus enabling one to achieve a desired stroboscopic order. To correctly describe evolution in Example 6, one needs to consider the third-order stroboscopic limit.
To give a broader view on the achieved result, the language of tensor networks enabled us to relate the memory kernel components with the multipoint correlation functions of the special form (the Waldenfelds cumulants). Multipoint correlations of orders n, n − 1, . . . , 2 determine the system dynamics in the (n − 1)-th order stroboscopic limit. Although multitime correlation functions have been used in the theory of open quantum systems (see, e.g., [115][116][117]), here we have explicitly demonstrated their origin in the collision model. We believe that the tensor network representation opens an avenue for a further analysis of the effect of multipoint correlations on the collisional dynamics, e.g., Wick's theorem for matrix product states [118] can be of great use.

Conclusions
We presented a tensor network approach to challenges in both the standard collision model and the collision model with correlated ancillas. We showed that the system-ancilla interactions in the standard collision model induce a correlated state of the system and ancillas that is naturally described by a right-canonical MPS (if the system and ancillas are initially in pure states) and a right-canonical MPDO (if the system and ancillas are initially in mixed states). Since the description of MPS and MPDO requires many fewer parameters as compared to a general multipartite state, we believe that the revealed representation can find applications in many practically relevant problems, e.g., this representation can allow one to go well beyond 12 collisions in the numerical study of quantum thermometry [34]. As far as initially correlated ancillas are concerned, we reviewed the recently proposed approach to the tensor network description of the system dynamics (with the emphasis on the two-point correlations) and generalized it to the case of multipartite correlations among ancillas. We showed conditions under which the higher-order stroboscopic limit is to be considered and how the Waldenfelds cumulants contribute to the memory-kernel master equation in this case.

Data Availability Statement:
The data presented in this study are available in article.