{'fraction_non_alphanumeric': 0.039250892065728765, 'fraction_numerical': 0.010682060955930817, 'mean_word_length': 4.473317576813036, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A brief journey through collision models for multipartite open quantum dynamics

The quantum collision models are a useful method to describe the dynamics of an open quantum system by means of repeated interactions between the system and some particles of the environment, which are usually termed"ancillas". In this paper, we review the main collision models for the dynamics of multipartite open quantum systems, which are composed of several subsystems. In particular, we are interested in models that are based on elementary collisions between the subsystems and the ancillas, and that simulate global and/or local Markovian master equations in the limit of infinitesimal timestep. After discussing the mathematical details of the derivation of a generic collision-based master equation, we provide the general ideas at the basis of the collision models for multipartite systems, we discuss their strengths and limitations, and we show how they may be simulated on a quantum computer. Moreover, we analyze some properties of a collision model based on entangled ancillas, derive the master equation it generates for small timesteps, and prove that the coefficients of this master equation are subject to a constraint that limits their generality. Finally, we present an example of this collision model with two bosonic ancillas entangled in a two-mode squeezed thermal state.

The quantum collision models are a useful method to describe the dynamics of an open quantum system by means of repeated interactions between the system and some particles of the environment, which are usually termed "ancillas". In this paper, we review the main collision models for the dynamics of multipartite open quantum systems, which are composed of several subsystems. In particular, we are interested in models that are based on elementary collisions between the subsystems and the ancillas, and that simulate global and/or local Markovian master equations in the limit of infinitesimal timestep. After discussing the mathematical details of the derivation of a generic collision-based master equation, we provide the general ideas at the basis of the collision models for multipartite systems, we discuss their strengths and limitations, and we show how they may be simulated on a quantum computer. Moreover, we analyze some properties of a collision model based on entangled ancillas, derive the master equation it generates for small timesteps, and prove that the coefficients of this master equation are subject to a constraint that limits their generality. Finally, we present an example of this collision model with two bosonic ancillas entangled in a two-mode squeezed thermal state.

I. INTRODUCTION
The simplest example of an open quantum system is a single quantum particle (e.g., a qubit) interacting with an external environment, such as a thermal bath. If the system-environment coupling is weak and the autocorrelation functions of the environment decay sufficiently fast, it can be shown that the dynamics of the state of the system is described by a Markovian master equation [1]. The complete characterization of the structure of this master equation is arguably the most celebrated result of the theory of open quantum systems. Specifically, it is known that the most general Markovian master equation of an open quantum system is the Gorini-Kossakowski-Sudarshan-Lindblad (GKLS) master equation [2][3][4]. The non-diagonal GKLS master equation for an open system with Hilbert space H S is written as (note that throughout the paper we will assume = 1): with D S ≡ dim(H S ) and where ρ S (t) is the density matrix describing the state of the system at time t. H S = H † S is a self-adjoint operator on H S , driving the free system dynamics.
with Tr[F j F † k ] = δ jk are some orthonormal traceless operators that together with the identity form a basis in the space of the bounded operators in H S . We refer to them as the Gorini-Kossakowski-Sudarshan (GKS) operators. The coefficients γ jk are the elements of the positive semi-definite Kossakowski matrix γ ≥ 0. The generator of the quantum dynamical semigroup [1] driven by Eq. (1) is the Liouvillian L.
Interestingly, starting from the full microscopic model of system+environment and applying the Born-Markov approximations [1] is not the only way to obtain a GKLS master equation for the state of an open quantum system. For instance, the same dynamics may be recovered through a collision model, that is, through repeated interactions between the system and some external particles, which conceptually belong to the environment, that are called ancillas. If the duration of each collision (the timestep of the collision model) is infinitesimal, and if a new fresh independent ancilla is prepared at every timestep, then it has been shown that, tracing out the degrees of freedom of the ancillas, we recover a GKLS master equation for the state of the system [5,6].
The interest in quantum collision models has rapidly grown in the past few years. Indeed, this formalism can be successfully applied to different fields where open quantum systems play a key role, including quantum thermodynamics [7][8][9] and non-Markovian quantum evolutions [10][11][12]. A comprehensive review on quantum collision models has been recently published [13], as well as a special issue of the journal Entropy [14], a perspective article [15] and a tutorial [16].
The aim of this work is to explore the main collision models for the dynamics of a multipartite open quantum system, that is, an open system composed of multiple subsystems. The interest in this kind of open quantum systems is broad, owing, for instance, to their relevance for collective dissipative phenomena [17], correlated noise in quantum computers [18], decoherence-free subspaces [19], and quantum thermodynamics [20][21][22]. More specifically, in this paper, we focus on collision models whose description is based on elementary collisions between the subsystems and the ancillas 1 , and that can reproduce a (local or global [20][21][22][25][26][27]) GKLS master equation in the limit of infinitesimal timestep. We identify four main schemes of this kind, namely the multipartite collision model [28], the cascade collision model [29], the composite collision model [30], and the collision model with entangled ancillas [31]. After presenting the derivation of the GKLS master equation for a generic collision model in Section II, we will discuss the general idea and the possible limitations of the above-mentioned collision models for multipartite systems in Section III. We will also pictorially show how these models may be simulated on a quantum computer. As the properties of the collision model with entangled ancillas are relatively unexplored, we will devote Section IV to a new discussion about them, where we will also analyze the master equation it generates and its limitations. Some concluding remarks are drawn in Section V.

II. THE STANDARD DERIVATION OF A COLLISION MODEL FOR A SINGLE OPEN QUANTUM SYSTEM
In this section, we will provide the standard derivation of a collision model for an individual open quantum system 2 . The open system, labeled by "S", is considered as a single quantum system whose inner structure is not taken into account in the collision model. Our aim is to formalize the open Markovian evolution of S driven by n repeated collisions with different external ancillas labeled by j = 1, . . . , n, in addition to a free system dynamics that is independent of the ancillas.
We introduce the Hilbert space of the system H S , whereas the Hilbert space of jth ancilla is H Ej . Each ancilla j is initialized in the state ρ Ej , and we assume that the state of the system before the jth collision is ρ S ((j − 1)∆t), where we refer to ∆t as the timestep of the collision model. Then, after the collision, the system state becomes: The unitary operator describing the collision is given by: where H S is the free system Hamiltonian (acting on H S only) and H I,j is the interaction Hamiltonian, which for convenience are both dimensionless. The magnitudes of their energies are captured respectively by g S and g I . From now on, we will assume that all the ancillas are identical and that they are initialized in the same state ρ E . Moreover, we will employ the same interaction Hamiltonian H I for all collisions. Then, we can introduce the quantum map φ ∆t to represent the system evolution after a single collision between the system and a generic ancilla: The map acts on the density matrices of the system alone (represented by the dot in the above equation), E refers to the degrees of freedom of a single generic ancilla, while U C is still given by Eq. (3), dropping the dependency on j. Then, if the state of the system is initialized in ρ S (0), after n collisions we obtain: Our goal is to obtain a Markovian master equation from the dynamics driven by the quantum map in Eq. (4). Note that the fact that the ancillas are uncorrelated and collide only once against the system means that the system evolution during the jth collision does not depend on the previous history of ρ S (t), i.e., the Markovianity between the collisions is already enforced by Eq. (5). Then, what we are looking for is a "smooth" evolution of ρ S (j∆t) to mimic the GKLS dynamics in Eq. (1). Intuitively, we can obtain such a dynamics by assuming that the timestep is infinitesimal, i.e., ∆t → 0 + .
Assuming ∆t → 0 + , we employ the Baker-Campbell-Hausdorff formula to expand the collision operators in Eq. (4), and we obtain: We now make two crucial assumptions to simplify the above equation. The first one is Tr E [[H I , ρ S ⊗ ρ E ]] = 0 for all ρ S . This is always valid, for instance, if we decompose the interaction Hamiltonian into system operators A α and ancilla operators B α , as usually done in the derivation of a Markovian master equation from a microscopic model [1], i.e., and all the ancilla operators are such that Tr E [B α ρ E ] = 0. The second assumption reads g S g I (that is, the collision energy is much larger than the system energy) and g 2 I ∆t → γ, g 2 S ∆t → 0 for ∆t → 0 + , where γ is a finite constant with the units of the inverse of time. Under these assumptions, Eq. (6) becomes: For a single collision, ρ S (t + ∆t) = φ ∆t [ρ S (t)], so: where we have used Eqs. (7) and (8), and the assumptions on the magnitude of g S and g I . The above equation is in the non-diagonal GKLS form with GKS operators given by {A α }. Indeed, H † S = H S , and the correlation matrix 3 . Therefore, we can reproduce the action of a quantum dynamical semigroup driven by the Liouvillian L at any time t by performing n repeated collisions between system and ancillas, with n = t/∆t. That is to say: Finally, note that, assuming the same approximations, we obtain the same result and the same master equation if we decompose the collision operator as follows: or equivalently U C (∆t) → U S (∆t)U I (∆t). Indeed, the error due to the non-commutativity of H S and H I is of the order of g S g I ∆t 2 , and therefore it is absorbed by the corresponding remainder in Eq. (8). This means that we can first perform the collision and then, independently, the free system evolution (or vice versa). This result allows us to tune g I and g S by using different effective collision timesteps. 3 Indeed, let us introduce the matrix M αα = Tr E [B † α B α ρ E ]. Then, note that, given any complex vector v, we have αα v * α M αα v α = Tr E [C † Cρ E ] ≥ 0, because the state of the ancilla is positive semi-definite. We have introduced the operator C = α vαBα.

III. COLLISION MODELS FOR MULTIPARTITE OPEN QUANTUM SYSTEMS
In this section, we review the main collision models for multipartite open quantum systems that have been introduced to date. Our goal is to understand the crucial idea at the basis of each of them, their major limitations, and which situations they are suitable for. Moreover, to stress the connection between the formalism of collision models and ancilla-based quantum simulation algorithms, we will depict the elementary mechanism of each collision model through gates applied on the system and ancillary qubits of a quantum computer. We will not show here the mathematical details of the derivations of these collision models 4 , and we suggest checking the original papers for further details. Anyway, to get a grasp on these mathematical details, the reader may get familiar with the derivation of the collision model for a single open system we presented in Section II, which is at the basis of all the collision schemes we will discuss.
We will write here the expression of the GKLS master equation for a multipartite open quantum system we aim to simulate through a collision model. Let us suppose that the open system is composed of M (for simplicity identical) subsystems. Its Hilbert space can be written as For simplicity, we will focus on multipartite GKLS equations that can be structured such that their GKS operators F k in Eq. (1) act locally on a single subsystem only. Anyway, extensions to many-body GKS operators can be made for all the collision models we will analyze. Under the assumption of local GKS operators, we can write the multipartite Markovian master equation as: H S is the effective Hamiltonian driving the free system dynamics. The indexes m, m = 1, . . . , M run over the subsystems, while α, α run over the different local GKS operators of the open dynamics. If the Hilbert space of each subsystem has dimension d, then α, α = 1, . . . , d 2 − 1, according to the discussion after Eq. (1). F m,α is a GKS operator that is local on the subsystem m, while γ m,α,m ,α are the coefficients of the master equation 5 .
We point out that we will not discuss here the collision models for multipartite systems put forward in different papers by Ç akmak et al. to study entanglement generation [33], quantum synchronization [34,35] and quantum Darwinism [36], as it is not clear how to derive a GKLS master equation for these models (the authors do not consider the limit of infinitesimal timestep). However, these schemes are useful to study open multipartite dynamics via the formalism of quantum maps with a finite timestep, and we refer the readers to the original papers for further details. Moreover, a collision scheme leading to multipartite open dynamics driven by a GKLS master equation has been recently employed to study the charging of a quantum battery [37]. This collision model, which we will not discuss in this section, creates cross terms between the subsystems by engineering a many-body non-linear collision Hamiltonian that couples both the subsystems and the ancillas at the same time. In a similar spirit, a collective subsystems-ancilla interaction is engineered in recent collision models to study non-stationary phenomena of open spin systems [38] or to explore quantum thermodynamics with common reservoirs [39]. Finally, two recent papers [40,41] have shown that a quantum version of the cellular automata may implement the closed and open quantum dynamics of many-body systems. This formalism makes use of repeated gates that propagate the state of the system in time, and it is therefore analogous to quantum collision models.

A. Multipartite Collision Model (MCM)
Together with Gabriele De Chiara, we have recently introduced the Multipartite Collision Model (MCM), which is able to reproduce the most general Markovian dynamics of any open quantum system [28]. For this reason, the MCM may be considered the most general and complete scheme among the ones we are presenting here. The price to pay for such a generality lies in the fact that the MCM requires one ancillary qubit for each unordered pair of jump operators F m,α and F † m ,α in Eq. (12), i.e., for each quartet (m, α, m , α ). However, several shortcuts may be engineered depending on the symmetries of the multipartite system and of the master equation, and the number of necessary ancillas may be actually much lower. We refer the readers to the supplementary material of the original paper for further details on these shortcuts [28]. The action of a single ancilla of the MCM is depicted in Fig. 1, assuming, for simplicity, that the subsystems of the multipartite system are qubits. The MCM reproduces each single "brick" of the GKLS master equation, that is, each term coupling the jump operators F m,α and F † m ,α in Eq. (12) (including the ones with m = m and α = α ), through an interaction U p between the subsystems m, m and an ancilla labeled by p, with p = (m, α, m , α ), which can be decomposed into three elementary collisions via the following formula: Each elementary collision is given by: where σ + E,p is the operator that creates an excitation in the ancillary qubit, while λ m,α are dimensionless parameters that we can freely tune. We notice that the GKS operators are "encoded" in the interaction Hamiltonian of each elementary subsystem-ancilla collision.
The interaction operator U p (∆t) reproduces the corresponding term in the master equation by taking the standard limits ∆t → 0 + and g 2 I ∆t → γ, according to the discussion in Section II. The key idea behind this result lies in the structure of the decomposed ancilla-subsystems interaction, which is based on the second-order Suzuki-Trotter formula [42]: for any bounded operators A, B on the Hilbert space of the system. The reader can verify that inserting U p (∆t) in Eq. (4) (or analogously in Eq. 6), and then using Eq. (15), yields the correct result. Note that, clearly, switching A and B does not change the result, and the same holds for switching m and m in Eq. (13). We can build an interaction operator as in Eq. (13) for all possible ancillas p = (m, α, m , α ). Then, their composition in any order will produce the most general dissipator of Eq. (12), if the ancillas are initialized in the ground state. Moreover, the free system dynamics driven by H S can be added following the procedure discussed in Section II, and the most general GKLS equation for a multipartite system is simulated in the limit of infinitesimal timestep. That is, the quantum map for the MCM reproduces any quantum dynamical semigroup driven by a generic L, the latter given by Eq. (12), according to Eq. (10).
Ref. [28] shows that the MCM is efficiently simulable on a quantum computer under standard assumptions, and it also discusses how to extend this formalism to non-local GKS operators and how to compute the error due to a small but finite ∆t. Finally, the experimental implementation of the multipartite collision model on a near-term quantum computer has been recently presented [24].

B. Cascade collision model
To the best of our knowledge, the first collision model for multipartite open quantum systems was introduced by Giovannetti and Palma in 2012 [29,43]. This collision model reproduces a cascade interaction, that is, a single ancilla collides with the first subsystem of the open quantum system, then it goes away and collides with the second subsystem, then with the third one, and so and so forth. The interaction of the ancilla with all the subsystems (from m = 1 to m = M ) represents a single timestep ∆t of the cascade collision model. Additional timesteps may then be implemented by repeating the same collision pattern. This scheme is depicted in Fig. 2 for a single ancillary qubit (hence, for a single timestep) and four system qubits (M = 4).
Following the original papers, we can encode the jump operators of the master equation we aim to reproduce in the collision Hamiltonian between the ancilla and every single subsystem. Specifically, the operator describing the interaction between the jth ancilla and the mth subsystem in the cascade collision model is [29,43]: with where {B j,α } is a family of ancilla operators that satisfy the standard assumption (e.g., see the discussion related to Eq. (7)) Tr Ej [B j,α ρ Ej ] = 0, while ρ Ej is the initial state of the ancilla. This is valid, for instance, if the ancilla is a qubit initialized in the ground state and all the B j,α are σ + Ej , as for the multipartite collision model. The free system dynamics driven by H S can then be inserted following the procedure discussed in Section II [43]. It can be shown [29,43] that this collision model in the standard limits of infinitesimal timestep and g 2 I ∆t → γ leads to a master equation that can be written as follows: We observe a GKLS dissipator that acts locally on subsystem m and a global dissipator that has the following structure: By adding and subtracting the terms 1 2 ρ S F m,α F † m ,α and 1 2 F m ,α F † m,α ρ S , we can write D (glob) m,m in a manifestly GKLS form [44]: where we have introduced the Hermitian operator H (m,m ) LS that plays the role of an effective Lamb-shift Hamiltonian: In the original papers, the coefficientsγ m,α,m ,α are tuned by inserting a general quantum map, which acts locally on the state of the ancilla only, between two subsystem-ancilla collisions. The properties of this quantum map modify the autocorrelation functions Tr Ej [B † j,α B j,α ρ Ej ] on which the coefficients depend. This cannot be straightforwardly represented in Fig. 2, as reproducing a non-unitary dynamics on a quantum computer is not a trivial procedure and typically requires some additional ancillary qubits [28,[45][46][47][48]. Alternatively, one may freely tune these coefficients by implementing more collisions whose interaction Hamiltonian in Eq. (17) depends on controllable parameters, as for the MCM.
The form of the cascade master equation in Eq. (18) is not as general as Eq. (12) because of the condition m > m in the global dissipator. This structure is indeed capturing the "causality" of the cascade interaction, i.e., the fact that the ancilla is colliding with some qubits before than with some others. This can be seen by performing the partial trace over all the subsystems with m > 1 on Eq. (18), which would yield a Markovian local dynamics for the subsystem 1 only [29]. Curiously, this causal order is encoded in the Lamb-shift Hamiltonians defined by Eq. (22). Without these terms, the master equation would be equivalent to the one generated by the MCM. This means that, to obtain the master equation simulated by the MCM, one may also implement the cascade collision model and then simulate a free system dynamics driven by a counter-Hamiltonian that cancels out the Lamb-shift terms in Eq. (22). We clarify this point through an example in Appendix A.
To obtain the most general GKLS master equation for multipartite systems starting from Eq. (18), one may implement additional cascaded collisions in which the ancilla first interacts with the last qubit, and so on until the first qubit (i.e., we reverse the order of the collisions, and the requirement in Eq. (18) becomes m < m). To be really general, one should implement any possible permutation of the order of the subsystems in the global dissipator. Such a scheme, however, resembles the structure of the interaction of the multipartite collision model in Eq. (13), which is already optimal due to the second-order Suzuki-Trotter decomposition. Therefore, the MCM can be thought of as the optimal extension of the cascade collision model to the most general multipartite Markovian dynamics, where there is no "preferred time order" of the collisions between the ancillas and the subsystems. Intuitively, the two "spatial" directions of the MCM elementary collisions in Eq. (13) and Fig. 1 (first with subsystem m, then with subsystem m , and then again with subsystem m) are creating and then destroying the evolution generated by the Lamb-shift Hamiltonian in Eq. (22), thus removing any causal order. To have a better grasp on this issue, we refer the readers to the example in Appendix A.
The cascade open quantum dynamics has been proven useful to study, for instance, Landauer's principle in multipartite systems [49], the heat flux and quantum correlations produced by the cascade interaction [50], interferometry [44,51], giant emitters in quantum optics [52,53], and steady-state thermodynamics [54]. Interestingly, in 1991 Bergou and Hillery presented a study on a model of two masers pumped by the same beam of excited atoms that can be formalized exactly as a cascade collision model [55]. The master equation for the state of the two masers has the same structure as Eq. (18).

C. Composite collision model for local master equations
The composite collision model was introduced by Lorenzo, Ciccarello and Palma in 2017 [30]. Although this collision model is applied to multipartite open quantum systems, the dynamics it can generate corresponds to a local master equation [20][21][22][25][26][27], that is, it is not able to describe global master equations [21,22,26] that can capture collective dissipative phenomena [17]. First, the ancilla a (blue qubit) collides with the subsystem qubit 2 through the operator UI (∆t). Next, the two subsystems (red qubits) interact through the free system evolution US(∆t). The green and orange interactions represent two-qubit gates.
The idea at the basis of the composite collision model is quite simple but elegant 6 . We consider, for simplicity, a system of M = 2 subsystems. The subsystem m = 1 interacts only with the second subsystem, while it is not directly coupled to an external environment (i.e., it does not collide with any ancilla). In contrast, the subsystem m = 2 collides with a fresh ancilla at each timestep ∆t, and the magnitude of this interaction is given by g I . Moreover, a free system dynamics driven by a non-local Hamiltonian H S with magnitude g S that couples the two subsystems is implemented between two consecutive collisions (i.e., one free system evolution for each timestep). Therefore, the total evolution for a single timestep can be described by the operator in Eq. (11). A single timestep of the composite collision model is depicted in Fig. 3, where both the subsystems and the ancilla are qubits. Note that, even if the first subsystem is not directly interacting with the ancilla, its dynamics is anyway non-unitary because of the interaction with the second subsystem. For instance, the state of the first subsystem can thermalize.
In the framework of the composite collision model, we can recover a GKLS master equation in the standard limits ∆t → 0 + , g 2 I ∆t → γ, g S g I , as discussed in Section II. This master equation is local, that is, its dissipator can be decomposed as a sum of local terms whose jump operators act on a single subsystem only. By considering a generic number of subsystems and implementing all the possible local ancillary collisions in order to tune the coefficients at our will, the most general master equation generated by the composite collision model is: which is equivalent to Eq. (12) with m = m.
While the composite collision model is able to reproduce the most general local Markovian master equation, it was originally introduced to study non-Markovian dynamics by considering as the "open system" all the subsystems that are not directly colliding with the ancillas (e.g., subsystem m = 1 in Fig. 3), while the subsystems that are interacting with the environment particles (e.g., the second subsystem in Fig. 3) are considered as "pseudomodes". Tracing out the degrees of freedom of both the ancillas and the pseudomodes, we recover a non-Markovian dynamics for the open system only. This way to describe non-Markovianity has been proven to be particularly successful [56,57]. Furthermore, the formalism of the composite collision model has been employed to study local master equations, mainly for quantum thermodynamic analyses [9,[58][59][60], and quantum correlations in multipartite open dynamics [61]. Moreover, Ref. [62] compares the heat currents obtained either through the composite collision model with two system qubits or through a similar model with a collective term in the system Hamiltonian (these two models are the collisional analogue of the local and global master equations [21,22,26]), although the infinitesimal timestep limit and the corresponding GKLS master equation are not studied in this work. Finally, a similar collision model for local master equations (with possible extensions to global ones) can be found in Ref. [63]. In this work, the derivation of the master equation relies on the standard microscopic theory of open quantum systems [1] and not on the formalism described in Section II. First, the ancillas a and a (blue qubits) get entangled through the operator UE(∆t). Next, the ancilla a (a ) collides with the subsystem 1 (2) (red qubits) through the operator UI,1(∆t) (UI,2(∆t)). The green and pink interactions represent two-qubit gates.

D. Collision model with entangled ancillas
In this subsection, we will introduce the last type of collision model for multipartite open dynamics we consider in this work. Specifically, we will present a collision model that makes use of entangled ancillas at each timestep. Note that this is different from introducing quantum and/or classical correlations between ancillas prepared for the collisions at different timesteps, which is one of the ways to study non-Markovianity in collision models (see Ref. [13] and references therein). In contrast, here we are interested in collision models where pairs of entangled ancillas collide at the same time with different subsystems of a multipartite open quantum system. A single timestep of this collision model is depicted in Fig. 4, where for simplicity we consider a two-qubit open system colliding with two entangled ancillary qubits. To the best of our knowledge, this collision scheme has been introduced first by Daryanoosh et al. in 2018 [31], and then it has been employed to study the thermodynamics of thermal machines powered by correlated baths [64], the quantum Rayleigh problem [65] and quantum trajectories for correlated ancillas [66,67].
The properties of the collision model with entangled ancillas have not been studied as much as the ones of the other collision models we have introduced in this section. Specifically, Refs. [31,66] derive a master equation for a two-qubit environment whose ancillas are prepared in a linear composition of different Bell states and make some comments about the possible extensions to m qubits. It is our aim to characterize the collision model with entangled ancillas in the language of multipartite open quantum dynamics, so as to recover a GKLS master equation whose structure resembles Eq. (12). Therefore, we devote the following Section IV to a discussion on some properties of this collision model.

IV. PROPERTIES OF THE COLLISION MODEL WITH ENTANGLED ANCILLAS
In this section we will analyze some properties of the collision model with entangled ancillas that has been introduced in Section III D.

A. Model formalization
First, let us formalize the elementary timestep of this collision model that is shown in Fig. 4 for a simplified scenario. In general, we assume to have M subsystems and M ancillas that will collide locally with each of them. At time t = 0, the state of the open system is initialized in ρ S , while, to keep the analogy with the ancillary qubits of a quantum computer, we msy assume that the ancillas are all prepared in the ground state. Hence, we can write the total state of the ancillas as ρ E (0) = M j=1 |0 j 0|. Then, a single timestep of the collision model may be described as: where U tot (∆t) = U S (∆t)U I (∆t)U E (∆t), and: Note that we are keeping the same notation as in Section II, where U S describes a free system evolution (not shown in Fig. 4), while U I is the operator for the ancillas-subsystems interaction. Moreover, we have added the environment operator U E , driven by the Hamiltonian H E with magnitude g E , that entangles the ancillas. Our aim is to derive a GKLS master equation generated by the quantum map in Eq. (24) in the standard limit of an infinitesimal timestep. To do so, we assume to be in the regime of fast environment Hamiltonian 7 , that is: where µ is a finite dimensionless constant. Furthermore, we also assume the standard limits for g S and g I that were employed in Section II. Specifically, lim ∆t→0 + g 2 I ∆t = γ, where γ is a constant with the units of the inverse of time, and lim ∆t→0 + g 2 S ∆t = 0. Note that this implies g S g I g E , which clarifies the nomenclature "fast environment Hamiltonian". Therefore, using Eq. (26) and assuming that ∆t is very small, we can write U E (∆t) → U E (µ) = exp(−iµH E ), and rewrite Eq. (24) as: where we have introduced the entangled state of the ancillas: Note that, contrary to the decomposition of the operator in Eq. (11), we cannot commute U E with U I or U S , as µ is not an infinitesimal parameter. Therefore, U E (µ) must always be the first operator that is implemented during a single timestep. This is reasonable, as we need to prepare the ancillas in an entangled state before starting the multipartite collisions.
Finally, we point out that we may also choose a mixed entangled state of the ancillas by a linear combination of different states written as in Eq. (28). On a quantum computer, this may be realized either by a purification of the mixed state in an enlarged Hilbert space through additional ancillas, or by random sampling the preparation of different entangled pure states written as in Eq. (28). The most general state we can employ in the collision model can then be written as: where each ρ (k) E is given by Eq. (28) with a different environment Hamiltonian and a different µ k .

B. Derivation of the master equation and its properties
We can derive a GKLS master equation starting from Eq. (27) following the same procedure as in Section II, with an interaction Hamiltonian that can be decomposed as a sum of local ancilla-subsystem terms: where F m,α and B m,α are respectively local system and ancillary operators. Using Eqs. (8) and (9), in the limit of infinitesimal timestep we obtain:  7 We show in Appendix B the issues that arise if one does not assume to be in this regime. 8 Note that this assumption may not be trivial, contrary to the scenarios with separable ancillas (e.g., with σ + E j as ancilla operator, which has zero mean value on thermal ancillary states), and care must be taken to choose a suitable H I .
where the coefficients of this multipartite master equation are given by 9 : ]. (33) Note that the non-locality of Eq. (32) is due to the fact that the autocorrelation functions in Eq. (33) are evaluated on an entangled state of the environment. Eq. (32) resembles the GKLS master equation for multipartite systems we introduced in Eq. (12). Indeed, in principle we are able to encode any possible GKS operator F m,α in the interaction Hamiltonian in Eq. (30). However, we immediately find a limitation of Eq. (32) that is due to the structure of its coefficient, and that inevitably limits its generality. In particular, we observe that [B m,α , B m ,α ] = 0 if m = m , because they are operators acting on different ancillas. As an immediate consequence, we find: This is a serious limitation. For instance, suppose that the system is composed of two qubits (M = 2) and F m,1 = σ − m , F m,2 = σ + m . Then, Eq. (34) implies that, in the final GKLS master equation, the term σ − 1 ρ S σ + 2 will have the same coefficient as the term σ + 2 ρ S σ − 1 , and, as a consequence of the properties of the master equation, these coefficients have the same module as the ones of the terms σ + 1 ρ S σ − 2 and σ − 2 ρ S σ + 1 . This means that Eq. (32) is not able, for instance, to reproduce the action of a master equation for a global thermal bath at finite temperature (e.g., see Ref. [26]), as an equal decay rate for the emission and for the absorption in the cross terms is a signature of infinite global temperature.
Possible generalizations of Eq. (32), with the aim of tuning the coefficientsγ in the most general way, may be implemented by assigning two independent entangled ancillas to each possible pair of subsystems. Each ancilla would then collide locally with the corresponding subsystem. Each pair of collisions would generate a term of the GKLS master equation that couples two subsystems with tunable coefficients given by Eq. (33). Furthermore, more freedom may be granted by adding some tunable weights in the subsystem-ancilla interactions in Eq. (30). The number of necessary ancillas to achieve this generality is clearly exponential in the number of subsystems. Moreover, note that the constraint given by Eq. (34) is always valid for this type of collision model.
We have finally understood that the collision model with entangled ancillas gives rise to a limited variety of multipartite GKLS equations, and does not provide the generality that is granted by the multipartite collision model. However, quite intuitively, Eq. (32) is particularly convenient to simulate the Markovian dynamics that is microscopically generated by entangled separate baths [1,26]. Indeed, the structure of the coefficients of this type of master equation is quite similar to the one in Eq. (33). We will provide an example of such a master equation in the next subsection.

C. Example: two-mode squeezed ancillas
We provide here the example of a master equation generated by the collision model with entangled ancillas when the system is made of two qubits (M = 2) and the ancillas are two bosonic modes prepared in a two-mode squeezed thermal state. If S(ζ) = exp ζb † 1 b † 2 − ζ * b 1 b 2 , where b 1 and b 2 are the bosonic annihilation operator of respectively the first and second ancillary mode, and ζ is a complex variable, we have [68]: where ρ th (N j ) is a thermal state with average photon number equal to N j . Let us suppose that the subsystems-ancillas interaction Hamiltonian reads: Then, the coefficients of the master equation (Eq. (33)) are given byγ m,α,m , , where the labels α, α = †, · indicate whether or not the bath operators are daggered. By using the cyclic property of the trace, we obtain:γ According to the Baker-Campbell-Hausdorff formula, we have: where ζ = re iψ . Using the above relations, we observe that Tr Finally, we use the above results to derive the master equation for the collision model with entangled ancillas, given by Eq. (32) (for simplicity, we ignore the free system dynamics driven by H S ): The coefficients read: where for simplicity we have used a notation such that N 3 = N 1 . Note that the coefficient associated with the term σ − 1 ρ S σ − 2 must be equal to the one for σ − 2 ρ S σ − 1 , according to Eq. (34), and analogously for all the other cross terms. A microscopic model of an open quantum system whose master equation would resemble Eq. (39) may be two qubits interacting with two local, separate environments [26], where the initial state of the baths is such that the modes of each environment with the same frequency are entangled in two-mode squeezed thermal states. The same structure of the coefficients would then be recovered in the singular-coupling limit [1], corresponding to autocorrelation functions of the environments that decay instantaneously.
A master equation with the same structure as Eq. (39) can be simulated on a quantum computer by using ancillary qubits prepared in the entangled state b ee |ee + b gg |gg , where |e and |g are respectively the qubit excited and ground state, while b ee and b gg are normalization coefficients that we can freely choose [31]. To properly tune the values of the coefficients of the master equation, one may then add some weights in the interaction Hamiltonian Eq. (30), as done for the MCM in Eq. (14), and/or apply more collisions with new ancillas prepared in a different entangled state.

V. CONCLUSIONS
In this paper, we have analyzed four different collision models for multipartite open quantum systems that reproduce a GKLS master equation in the limit of an infinitesimal timestep. First, in Section II we have reviewed the derivation of such a master equation for the standard collision model for a single open system. Then, we have discussed the ideas and limitations of the collision models for multipartite open dynamics in Section III, also showing how to simulate them on a quantum computer for some very simple scenarios. We will briefly summarize these discussions here.
The multipartite collision model [28] (MCM) is the most general protocol we have discussed in this work. It can reproduce the Markovian master equation of any multipartite open quantum systems thanks to the structure of its collisions, which are based on second-order Suzuki-Trotter decomposition. To achieve this level of generality, however, it requires a large number of ancillas, which can be reduced in the presence of symmetries in the master equation.
The cascade collision model [29,43] can simulate master equations for cascade dynamics, where the environment ancillas collide with the subsystems in a time-ordered way (i.e., first with the first subsystems, then with the second one, and so on). However, it cannot efficiently reproduce any open dynamics that does not have this cascade structure.
The composite collision model [30] is a very simple and elegant collision model that simulates local master equations for multipartite systems, where each ancilla (say one for each subsystem) collides locally with a single subsystem only, while the different subsystems interact only through the free system dynamics driven by the Hamiltonian H S . However, it cannot simulate global master equations, that, for instance, display collective effects.
The collision model with entangled ancillas [31] can naturally simulate master equations for correlated environments (e.g., see the example we proposed in Section IV C). However, we have shown in Section IV that the coefficients of the master equations it generates are constrained by the condition expressed in Eq. (34). This limits the range of master equations it can simulate and, for instance, it is not suitable to reproduce the dynamics driven by a global bath at finite temperature.
Let us now analyze the cascade collision model. Due to the order of the elementary collisions in a single timestep, we immediately realize that the cross terms in the cascade Liouvillian will be (up to some constants) the same as i), ii), iii) and iv) for the MCM. However, v) and vi) are missing, because the ancilla first collides with subsystem 1 and then with subsystem 2, without bouncing back towards subsystem 1. Therefore, the symmetry between the two subsystems is broken. We thus understand that the second-order Suzuki-Trotter decomposition is the way to ensure that the subsystems in the expansion of the collision operators in the quantum map are treated on the same footing, i.e., all the possible combinations of GKS operators are present in the Liouvillian. This is also a requirement to obtain the most general GKLS equation for multipartite systems, expressed by Eq. (12). In contrast, the Liouvillian of the cascade collision model does not have two crucial cross terms: The absence of the terms ρ S F † 2 F 1 and F † 1 F 2 ρ S has important consequences on the dynamics of the two subsystems, which reflect the time order with which the ancilla collides against them. For instance, Eq. (A6) for two bosonic modes (specifically, with F j = a † j ) is exactly the master equation derived by Bergou and Hillery for two masers [55], which displays an unbalanced dynamics between the two cavities. For instance, if the initial state of the maser fields is the ground, the authors found that the photon number in the second cavity grows faster than in the first.
Following Ref. [44] and our discussion in Section III B, the cascade Liouvillian in Eq. (A6) can be brought into the symmetric structure of Eq. (A5), but then a Lamb-shift term appears in the master equation, which restores the time order of the collisions. This Lamb-shift term reads: . (A7)