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.

1. 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.

Figure 1

Standard collision model.

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 past-future 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, , for the state of d-dimensional ancillas. This includes, for instance, and 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 time-bin 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 complex parameters to be specified, not 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.

2. 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 of a finite dimension d. The state is fully defined by complex numbers in the decomposition where is an orthonormal basis in . A collection of complex numbers 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 connected in a line; see Figure 2.

Figure 2

Tensor diagram for a matrix product state.

and are rank-2 tensors with elements and , respectively, whereas for all the tensor has rank 3 and is composed of elements . On the other hand, if the physical index is fixed, then can be viewed as a matrix with elements . Similarly, if and are fixed, then and can be viewed as a row matrix and a column matrix with matrix elements and , respectively. Arrows in Figure 2 also indicate the order for multiplication of matrices. The contraction yields which explains the MPS name. A number of the values that the virtual index can take is not related to the physical dimension d of the particles. We will refer to the maximal number 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 , . The resulting diagram is depicted in Figure 3.

Figure 3

Equivalent diagram to that in Figure 2. ⊤ denotes transposition with respect to virtual indices.

Nonetheless, if indices are fixed, then the c-number does not change because

3. 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 , , and the environment consists of n uncorrelated ancillas in a pure state , . Each elementary collision is described by a unitary operator , 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.

Figure 4

Matrix product state of the system and ancillas induced by collisions in the standard collision model.

As a result of n collisions, we get a correlated state of 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, . Therefore, we can associate each virtual index with a vector , so that a collection of vectors for a fixed k forms an orthonormal basis in . The very diagram in Figure 4 serves as the proof for the following result.

Let the system and n ancillas be initially in the pure states where

The result of Proposition 1 explains the previously known observations of Ref. [8] that the partial swap interactions () generate the W-type of entanglement, whereas the controlled unitary interactions () 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].

Let the system and ancillas be qubits. The system is initially in the excited state Note that the matrix

The explicit relation between the unitary operator U and tensors , 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 to get rid of its degrees of freedom. The resulting state of n particles is , where 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 . Clearly, each collision could be described by its own unitary operator. Then, one should replace in the formula for in Proposition 1. Numerical tools for optimization over many unitary operators are presented, e.g., in Refs. [81,82].

Let us consider entanglement of the state with respect to a bipartition into ancillas on one side and ancillas and the system on the other side, i.e., the left-right bipartition with the boundary in between the ancillas k and . 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, , where . The reduced density operator for k ancillas is presented in the form of a tensor diagram in Figure 5a.

Figure 5

(a) Tensor diagram for the reduced density operator . Overline denotes complex conjugation. (b) Simplified tensor diagram for due to the right normalization condition.

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 conjugation. In fact, if , then because by Proposition 1. If , then Proposition 1 implies

If , then we deal with dummy indices and . Hence, we have proved the following result.

MPS

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. Indeed, Equation (5) is equivalent to , 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 .

Entanglement entropy of the state with respect to the cut in between the ancillas k and reads

Note that because is an upper bound for the Schmidt rank of .

4. 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 . The initial state of n ancillas is given by a factorized density operator . Collisional dynamics with the elementary unitary operator U drives the system and ancillas to the state where the subscript in the notation means that U nontrivially acts on the system and the k-th ancilla. A tensor diagram for Equation (8) is presented in Figure 6a.

Figure 6

(a) Tensor diagram for Equation (8). (b) Matrix product density operator. (c) Tensor decomposition guaranteeing positive semidefiniteness of the matrix product density operator.

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,b are different. The operator in Figure 6b reads

Here, and are rank-4 tensors, whereas is a rank-6 tensor for all , …, n. If indices and are fixed, then we treat as a row matrix with elements . Similarly, if indices and are fixed, then we treat as a column matrix with elements . If and indices and are fixed, then we treat as a matrix with elements , i.e., is a row multiindex and is a column multiindex. The explicit expressions for M-tensors follow from Figure 6a,b and read

Let be the spectral decomposition for the system initial state. Let 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 , the decomposition takes place, with and if , and if , and if . Tensor diagrams in 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 one-dimensional 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 for k ancillas from future collisions at times . The results of this section are summarized as follows.

The standard collision model with initially mixed states of the system (

The main benefit of the constructed MPDO representation is that it exploits only parameters instead of 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].

Let the system and ancillas be qubits. The system is initially in the excited state where where

5. 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 be a right-canonical MPDO for n ancillas shown in Figure 7a. Here, we have added a formal density operator (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, is the trivial 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 with respect to horizontal virtual indices.

Figure 7

(a) Matrix product density operator representation for a generally correlated state of ancillas, where we have formally added an auxiliary density matrix for the bond degrees of freedom to redefine the first MPDO tensor (dotted region). (b) Recurrence relation of density operators for bond degrees of freedom in the interaction-free evolution. (c) Tensor diagram for the system density operator after k collisions with correlated ancillas.

Figure 7b illustrates that the tensors 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., and ) provided is also a density matrix. This follows from the right normalization condition (16).

Figure 7c depicts the system density operator after k collisions. The partial trace over ancillas corresponds to a vertical connecting line for the bond degrees of freedom (blue arrows in Figure 7c). The partial trace over ancillas 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 where ∘ denotes a map concatenation. Each map is completely positive and trace preserving because it adopts the diagonal sum representation The trace preserving property 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 qubit-ancilla 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 environment or , where either and , or and ; , , are nonintersecting subsequences of collision numbers. The latter environment reproduces an arbitrary Pauli dynamical map [51].

Consider the GHZ state of 3-dimensional ancillas The qubit system is initially in the state where

Substituting Equations ( where We portray the typical behaviour of

6. Master Equation

Equation (18) defines the discrete dynamical map that transforms the initial system density operator to the system density operator after k collisions, i.e., . If , then we interpret as a continuous time t. A time-local master equation can be derived if is invertible; namely, . 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 be a collection of the density operators for the bond degrees of freedom generated by Equation (17). For each k, define the following map acting on both the system and the bond degrees of freedom:

Since , we have , so is a projection. A pictorial representation of the projection is given in Figure 9a, which shows that breaks the left-right correlations between bunches of ancillas and . Clearly,

Figure 9

(a) Projection . (b) Reduced density operator for the -th ancilla in the initial correlated state of ancillas. (c) Reduced density operator for the -th ancilla and -th ancilla in the initial correlated state of ancillas.

A complementary projection is defined through where is the identity transformation for the system and bond degrees of freedom. is also a projection because . We can also rewrite in the form where 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 , we readily get

To simplify the notation, let us introduce the system-bond density operator after the k-th collision, . Then, for all k. Applying 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 to the both sides of equation , we obtain

Recalling the relation and taking partial trace over the bond indices in Equation (34), we get

Subtracting 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 relates the density matrix increment (in between the times and ) with the past density operator at time . If , then we have a time-local term giving the density operator increment caused by the latest collision (among those that have already happened): with being a reduced density operator for the -th ancilla in the initial state; see Figure 9b. If , then describes a nontrivial effect of preceding collisions on the system evolution and reads

If there were no correlations in the environment, then would be the only contribution to the kernel because would vanish for all . Indeed, the MPDO rank equals 1 for a factorized environment state, so for all k. Each is unambiguously defined because is the trivial identity matrix in this case, and for any matrix F. If the environment is correlated, then the memory contribution in general.

Consider the GHZ state of 3-dimensional ancillas

The memory kernel

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.

7. Effect of Two-Point Correlations

An important simplification comes from a series expansion for with respect to the interaction strength between the system and an individual environment particle. Let be the system-particle interaction Hamiltonian during the m-th collision, where ℏ is the reduced Planck constant, g has the physical dimension of frequency, and is a dimensionless Hermitian operator with the operator norm . Then, the elementary unitary interaction in the m-th collision is . The map in Equation (19) has a contribution of both U and , so we have where and denote the commutator and the anticommutator, respectively. Substituting Equation (41) into Equation (38), we keep track of the leading terms in the memory kernel, namely,

The term vanishes because see Equation (30). Physically, the 0-th order of involves no system-environment interaction and, consequently, no contribution to the memory kernel.

To calculate the term , we should fix for all but one of . If , then we have a zero contribution to because of Equation (44). Suppose ; then, , and because is a trace preserving map due to the right-normalization condition, whereas nullifies the trace of any operator (see Equation (30)). Therefore, the term also vanishes.

Similar considerations for the term lead to a conclusion that may be nonzero only if we fix for and , whereas for all we fix . This results in where the coefficient reads

Recalling the definition (30), we get

Tensor representation in Figure 9c justifies that i.e., we get the reduced density operator for the -th ancilla and -th ancilla in the initial correlated state of ancillas. Similarly, i.e., we get a tensor product of individual reduced density operators for the -th ancilla and -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 by a conventional formula we readily see that , where and . 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.

8. 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 coincide for all , and MPDO is fully described by the density matrix and the tensor M. In this case, all local density operators for individual ancillas coincide ; however, the two-ancilla density operator . If ancillas are initially in such a homogeneous state, we can expect that the kernel 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 in continuous time , namely, [66]. The correlations among ancillas are irrelevant in this scenario because in the considered limit.

To reveal a nonunitary system dynamic at a long timescale, one should consider a different limit , [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 explodes in the master equation because ; however, this problem disappears in a proper interaction picture [66].

In the stroboscopic limit, one cannot simply replace in the left hand side of Equation (36) by if the term does not vanish, because , and the second summand cannot be neglected. However, if the expression vanishes, then the characteristic frequency of system dynamics is so that . The second summand vanishes in the (first-order) stroboscopic limit because . The time-local memory kernel component should also be considered in this limit, i.e.,

The higher-order contributions to the memory kernel (43) vanish only if the correlation length (in the chain of ancillas) is finite. If this is the case, then for .

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 being defined by the second-largest eigenvalue of the transfer matrix (in absolute values). If the correlation length is finite, then represents a sum of exponentially decaying terms, where are eigenvalues of the transfer matrix T such that and are the corresponding maps.

To explicitly find the kernel 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 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 (), then the system “feels” not only the individual ancillas (which results in the local term ) but also a somewhat averaged correlated state (which results in the nonlocal term ). We summarize these results as follows.

Let the system collisionally interact with an array of ancillas in the homogeneous MPDO with a finite correlation length. If the expression

Consider an infinite chain of spin-1 particles (ancillas) in the AKLT state [ Each individual ancilla has a reduced density operator

At time Averaging over single-ancilla degrees of freedom yields

The local term is given by formula (

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 Note that and has eigenvalues 1 (of multiplicity 1) and where

Substituting Equation ( On the other hand,

Finally, Equation (

The reader may notice the formally negative rate

9. Effect of Multipoint Correlations in the Higher-Order Stroboscopic Limit

We start with an example stimulating the discussion of the higher-order stroboscopic limit.

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 reads

For such an interaction,

The exact treatment of the problem via Equation ( with the depolarization function

A feature of the depolarizing dynamical map ( i.e., the image of the system Bloch ball shrinks after the first collision and then expands after the second collision. If

If

We will refer to the limit , as the n-th order stroboscopic limit. Surely, the expressions g, , …, explode in this limit; however, if their contribution to the system dynamics vanishes, then the limit is well defined. The higher-order contributions , , …vanish in the n-th order stroboscopic limit. In Equation (55), for , 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 , which describes -point correlations among ancillas. Recalling the notation in Equation (41), we get, for instance, the following expression for the third-order memory-kernel: where is a three-point correlation function, which is determined by tensor diagrams similar to those in Figure 9c and reads with being the third-order Waldenfelds cumulant [113,114]. In the second-order stroboscopic limit, the expression reduces to the time-local generator , which contributes to the final GKSL master equation . 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, , …, 2 determine the system dynamics in the -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.

10. 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.