Quantum homogenization in non-Markovian collisional model

Collisional models are a category of microscopic framework designed to study open quantum systems. The framework involves a system sequentially interacting with a bath comprised of identically prepared units. In this regard, quantum homogenization is a process where the system state approaches the identically prepared state of bath unit in the asymptotic limit. Here, we study the homogenization process for a single qubit in the non-Markovian collisional model framework generated via additional bath-bath interaction. With partial swap operation as both system-bath and bath-bath unitary, we numerically demonstrate that homogenization is achieved irrespective of the initial states of the system or bath units. This is reminiscent of the Markovian scenario, where partial swap is the unique operation for a universal quantum homogenizer. On the other hand, we observe that the rate of homogenization is slower than its Markovian counter part. Interestingly, a different choice of bath-bath unitary speeds up the homogenization process but loses the universality, being dependent on the initial states of the bath units.


Introduction
Over the last decade or so, exceptional advancement in quantum technologies [1] has inspired an extensive interest in the emerging field of Quantum Thermodynamics [2,3].Dealing with thermodynamic processes in out of equilibrium scenario leads to dynamical considerations [4] and a proper description of open quantum system becomes immensely crucial.Process of thermalization is one of the most fundamental problems of non-equilibrium thermodynamics.To put it simply, when a system is kept in contact with a bath, eventually it reaches thermal equilibrium.Microscopic description of this phenomenon has a rich history starting with Boltzmann's transport equation [5] in the context of classical statistical mechanics.In quantum domain one needs to describe the evolution of system density matrix which is interacting with bath degrees of freedom.But in practice, first principle derivation of an exact master equation for system density matrix is extremely difficult due to the lack of precise knowledge about the environment and the interaction.To circumvent this problem usually the Markovian approximation is made to pose the master equation in the so called Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form [6] neglecting all the memory effects.One benefit of this approach is that the system dynamics describes a legitimate physical evolution being completely positive and trace preserving (CPT).But, this leaves out a vast variety of realistic phenomena involving non-Markovian [7,8,9] effects.One approach to tackle these problems is based on the time-nonlocal Nakajima-Zwanzig equation [10,11,12,13], which is notoriously hard to solve.Even the analysis of complete positivity condition of the dynamics is a highly non-trivial task [14].Alternatively, collisional models (CMs) [15,16] offer a simplified approach in a totally controllable manner to model an open quantum system dynamics, where complete positivity is ensured by design.Historically, first introduced in a paper by J. Rau [17] in 60's, CMs gained a renewed interest in early 2000's [18,19,20,21,22].The most basic building block of CMs consists of two ingredients: firstly, the bath (environment) is made up of identically prepared sub-units (ancillas) and secondly, the system interacts sequentially with one ancilla at a time via a unitary.With initially uncorrelated bath and in the absence of ancilla-ancilla interaction, CMs automatically lead to a GKSL type master equation for the system density matrix in the continuous time limit [20,21] without performing any approximation.Such simplicity makes the CMs advantageous to study non-Markovian dynamics by some modifications in the basic model outlined above.By introducing ancilla-ancilla interaction [23,24,25,26,27,28,29], using initially correlated bath [30,31,32,33,34], through a composite collisional model [35] or through multiple collisions of the system with each ancilla [36,37] can give rise to non-Markovian dynamics.With rapid developments in CMs, a number of applications in (but not limited to) quantum thermodynamics [38,39,40,41,42,43,44,45,46,47,48,49,50,51,52] have been explored.
In the setting of Markovian CMs, above mentioned thermalization problem has been addressed [18,19] through a more general process called homogenization.It is a process through which the system density matrix is transformed to the same state of the identically prepared ancillas in the limit of infinitely many collisions.In particular, authors in Ref. [18] showed that for a qubit system interacting with qubit ancillas, partial swap (PSWAP) is the unique system-ancilla unitary, for which homogenization is achieved regardless of the initial states of the system and the ancillas.Recently, the problem of homogenization has been studied [53] in a specific non-Markovian scenario, when the ancillas are locally identical but globally correlated.Interestingly the authors found that some initial correlations do not allow the system to homogenize.This prompts the question what happens if introduce the non-Markovianity through ancilla-ancilla interaction.Does the homogenization occur?In this paper, we study this particular scenario answering this question affirmatively.Specifically, we take a qubit system and infinitely many identical qubit ancillas, which constitute the bath.When system-ancilla and ancilla-ancilla unitaries are both PSWAP, we numerically give evidence that the system homogenizes with the bath irrespective of the initial states we choose for the system or the ancillas albeit slowly, compared to the case when there is no ancilla-ancilla interaction.The decrease in the rate of homogenization can be understood from the fact that due to the additional ancilla-ancilla interaction, system does not interact with a fresh ancilla each time, which takes more collisions for the homogenization to happen.Interestingly, for a specific choice of the ancilla-ancilla interaction unitary, which is not PSWAP, homogenization process gets almost as fast as the Markovian counterpart.But as PSWAP is the unique universal homogenizer [18], in this case homogenization occurs only for the ancilla states which are diagonal with respect to the computational basis.If system and ancilla Hamiltonian are explicitly mentioned this can be regarded as thermalization [54] by assigning a temperature to the ancillas.
The paper is organized as follows, In Sec. 2, we recapitulate the basics of collisional model.In Sec. 3, we employ a new technique to show that homogenization happens universally for PSWAP operation in the Markovian case.Subsequently, in Sec. 4 and 5, we discuss the homogenization problem for non-Markovian scenario taking ancillaancilla interaction PSWAP and modified PSWAP respectively.Finally, in Sec.6, we conclude.

Framework of collisional model
In this section we briefly discuss the basics of collisional model that we will follow in the subsequent sections.First, we start with the Markovian scenario.Consider a quantum system S in contact with a bath B, which is assumed to be a collection of smaller and identical sub-units {B n }.There is no restriction on the dimension of the Hilbert space for system or the ancillas.With initial system state ρ S 0 , and ancilla state η, the joint state is taken to be the product state, By design, the dynamics takes place through successive pairwise collisions between the system and the ancillas.In Fig. 1a, we provide a schematic diagram to visualize this.First, the system interacts with the ancilla B 1 , then the system is disconnected from the bath.Next the system interacts with a fresh ancilla B 2 and so on.Each ancilla interacts with the system only once.The system-ancilla interaction is modelled by a unitary.Assuming all the collisions have same duration ∆t, the unitary for the n-th collision is given by, where, H S is the system Hamiltonian, H Bn is the Hamiltonian for n-th ancilla, and V n is the interaction Hamiltonian between the system and the n-th ancilla.After n collisions the joint state of system and bath is given by, The system state is obtained by tracing out the bath degrees of freedom, where, ξ[ρ] = Tr Bn {U n (ρ ⊗ η)U † n } denotes a completely positive map (known as dynamical map), ensured by the fact that there is no initial correlation between the system and the bath.It follows immediately form Eq. ( 4) that ρ S n = ξ n [ρ S 0 ].To see the Markovian property, we introduce the dynamical map, It is easy to see that Λ n obeys the well known semi group property, Λ n = Λ n−m Λ m , for any integer m, between 0 to n.This justifies the name Markovian collisional model.Next, we discuss the non-Markovian collisional model introducing the ancilla-ancilla interaction.Similarly as before, we take product system and bath state without any initial correlation between the ancillas.We denote the unitary between the system and n-th ancilla by U n as before, and the unitary between (n − 1)-th and n-th ancillas by U B n−1 Bn .In Fig. 1b, we present a schematic diagram of the model.The dynamics is described as follows.The system first interacts with the ancilla B 1 , then the system is detached from the bath and ancillas B 1 and B 2 now interact with each other.The ancilla B 1 is now detached from B 2 .Next the system interacts with the ancilla B 2 and the process goes on.After the collision between the system and the n-th ancilla, the joint state is given by, where, U ′ n = U n U B n−1 Bn and the system state is given by, Here, we have not specified the strategy for tracing out the bath degrees of freedom.
Originally it means to trace out the bath degrees of freedom all at once at the very end of the process which deals with the total state of the system and n ancillas.But computationally it is very hard when n is large.But as long as only the dynamics of the system is concerned we have a nice way around.Note that, to work out the dynamics of the system only its reduced state is required before the interaction with an ancilla B i .It can be understood from the following argument [55].For a density matrix ρ 12 ∈ H 1 ⊗ H 2 , where H 1 and H 2 are two Hilbert spaces, the following is true with U a unitary operation.
Although the above equation is valid for any CPTP map, for our scenario unitary is enough.This means that, before the system interacts with let's say B n , all the previous correlations with n − 1 ancillas can be destroyed without affecting the dynamics.Consequently we do not have to deal with the joint state σ n .Now the question is, before the system interacts with B n , when to destroy the correlations.Thus giving rise to three different ways [25,45,55] to do this retaining or destroying the correlations in different stages.Below we list them and give a short description.
Scheme 1 : Tracing out is performed after the system interacts with the ancilla B n−1 , before B n−1 interacts with the ancilla B n followed by system's interaction with B n .After the collision of the system and the ancilla B n−1 , the reduced state of the system and the ancilla B n−1 are respectively given as the following, Here, ρ S n−2 is the system state after the collision of the system and the ancilla B n−2 .Similarly, ρ ′ B n−1 is the reduced state of the ancilla B n−1 after its collision with the ancilla B n−2 .After the collision of the ancillas B n−1 and B n , the reduced state of the ancilla B n is given by, Finally, state of the system after it interacts with the ancilla B n is given as, Scheme 2 : In this scheme, bath degrees of freedom of the ancilla B n−1 are traced out after it interacts with the ancilla B n .Joint state of the system and the ancilla B n after its collision with B n−1 is given by, where, State of the system after it interacts with the ancilla B n is, Scheme 3 : Here, tracing out of the ancilla B n−1 is done after the system completes its interaction with the ancilla B n .The correlation is kept until then.Denoting Bn , the system state after its collision with the ancilla B n is given as, Clearly scheme 2 and scheme 3 will give rise to the same dynamics.Because, as mentioned earlier, the reduced state of the system is the only thing we need before its interaction with B n .
A simple example can be given now to show that semi-group property does not hold in the presence of ancilla-ancilla interaction.We take the ancilla-ancilla unitary to be the swap operation: U B n−1 Bn ≡ S n−1,n (lower indexes indicate that swap operation takes place between B n−1 and B n ancillas), defined as, for any two states |ϕ⟩ and |ψ⟩.One can show that, the system state after its collision with n-th ancilla is [15], This can be contrasted with the scenario with no ancilla-ancilla interaction, where the map Λ n can be realized as operating ξ for n times.Clearly, the semi-group property does not hold here implying the presence of memory effect.This concludes our discussion about collisional model.

Homogenization in Markovian collisional model
Quantum homogenization in the context of Markovian collisonal model was first introduced in the Ref. [18] for qubits.In the limit of infinite collisions, system state homogenizes with the initially identically prepared ancilla states.In the language of previously introduced notations this means, Homogenization process is called universal if this holds for any states ρ S 0 and η.Clearly, this process is a generalization of the thermalization process [19], where the ancilla states η are thermal with a fixed temperature.The unitary interaction U between the system and the ancilla will allow homogenization iff the following conditions hold, Taking a unitary interaction that satisfies the above conditions, next step is to show that the system actually reaches the initial ancilla state asymptotically such that further interactions do not change the states of the system or the ancillas and the homogenization is achieved.In Ref. [18], it was shown that for qubit scenario, PSWAP is the unique operation for which homogenization is achieved irrespective of the initial states of the system or the ancillas.In the following, using a new technique, we show that with PSWAP one can achieve homogenization universally for qubits.We do not provide here the proof for the uniqueness of the PSWAP.In the subsequent sections, we use this technique to discuss the non-Markovian scenario.
We take an arbitrary initial system state in the Bloch vector notation, Initial states of the ancilla are taken to be identical as, where ⃗ l j = ⃗ l, for all j.Though all the initial ancilla states are same, we use the subscript for easy reading.We take the system-ancilla unitary to be PSWAP, which obviously satisfies the conditions of Eq. ( 19).We denote the unitary between the system and the n-th ancilla as following, where, the SWAP operator S 4×4 is given by, Such PSWAP operation can be generated by two-qubit Hamiltonian S 4×4 .Such a twoqubit Hamiltonian can, in principle, be realized in laboratory using a Nitrogen vacancy centre (which is a spin-1/2 system) in the environment of a bath of individual spin-1/2 impurities where the centre, in fact, interacts with the individual spin-1/2 impurities in succession -reminiscent of the collisional model [56].In fact, in the context of a spin-bath [57] corresponding to a spin-star model (where a single central spin interacts with a bath of several spins), it may be possible to realize the corresponding global unitary evolution as a concatenation of unitaries corresponding to the interaction of the central spin and the 1st spin of the bath, followed by a 1st bath spin and 2nd bath spin interaction, followed by the central spin and 2nd bath spin interaction, and so on.After the first collision, states of the system and the ancilla B 1 are respectively given by, where, Superscript in the above equations denotes the number of collisions a state encounters.
After the first collision with the ancilla B 1 , the system interacts with a fresh ancilla B 2 , and the process goes on.Proceeding similarly as before, we get the following relations after the system interacts with the n-th ancilla B n , where, Now to show that Homogenization happened, we have to show that in the limit, n → ∞, From the expression of ⃗ k (n) in Eq. ( 29), it is straightforward to see that, where, K = (cos The above inequalities hold as | ⃗ l| 2 ≤ 1 and 0 < α < π.Now, it follows (from Eq. ( 32) -using the ratio test) that, ⃗ k (n) n→∞ − −− → ⃗ l.Putting this in Eq. ( 30), one immediately obtains ⃗ l (1) This shows that with PSWAP operation, homogenization occurs for all initial states of the system or the ancillas.Next we will use this technique and some other tools to study the homogenization in the presence of ancilla-ancilla interaction.

Homogenization in non-markovian collisional model with Partial SWAP
We now introduce the ancilla-ancilla interaction to model a non-Markovian dynamics as described in Sec. 2. In this section, we consider both system-ancilla as well as ancillaancilla interaction to be PSWAP, guaranteeing the conditions of Eq. ( 19) are satisfied for both of these unitaries for any state ρ.This means if homogenization occurs that will happen for all initial system or ancilla states.As ancilla-ancilla interaction is present in this scenario, we also have to consider whether the ancilla states reach its initial states asymptotically along with the system state.Because only then the conditions of Eq. ( 19) will be satisfied implying that no further interaction will change the state of the system or the ancillas.Similar to the previous section, we take the system and ancilla states in Bloch vector representation given by the same expressions, as in Eq. ( 20) and Eq. ( 21).Additionally, system-ancilla interaction is given by Eq. ( 22) and the ancilla-ancilla interaction between the ancillas B n−1 and B n is given by, From now onward, we will be adopting the scheme 1 for carrying out all the analytical calculations in the subsequent discussions.However, only numerical analysis have been performed for scheme 2, leaving the analytical part for future investigations.After the collision between the system and the first ancilla B 1 , system state is given by, ρ S 1 = 1 2 (1 + ⃗ k (1) .⃗σ) (see Eq. ( 24)) and the ancilla state is given by, η 1 .⃗σ) (see Eq. ( 25)), where ⃗ k (1) , and ⃗ l (1) 1 are given by Eq. ( 26) and Eq. ( 27) respectively.The next step is new compared to the previous section.First ancilla B 1 now interacts with the second ancilla B 2 via the unitary of the Eq.(33).After that, the system interacts with the ancilla B 2 and the process goes on.Now, the state of B 1 , after the B 1 − B 2 collision is given by, Similarly, in this case, the state of B 2 is given by, Here,

⃗ l
(2) 1 = cos 2 δ ⃗ l (1) 1 Now, the system and the ancilla B 2 will interact with each other and the states of the system and B 2 after this collision are given by, where, One can carry on the similar calculations and get the following recursive relations, For the first three relations, Eq. ( 42), Eq. ( 43), and Eq. ( 44), n ≥ 2, while, ⃗ k (1) and ⃗ l (1) 1 are given by the Eq. ( 26) and Eq. ( 27) respectively.Now, to show that homogenization occurs we have to show that in the limit n → ∞, the above four recurrence relations converge to ⃗ l.The proof is not straightforward like the Markovian scenario.We are not in position to prove homogenization from the above set of recurrence relations following the ratio test without any assumption.For example, let us assume that, ⃗ l With this assumption we can show that other three vectors also converge to ⃗ l in the limit n → ∞.This is shown by the ratio test -as done for the Markovian scenario.In fact, from Eq. ( 42) one can show that (for large n), where, (given in Eq. ( 45)) converge to ⃗ l asymptotically then it is easy to show that ⃗ l (0) n n→∞ − −− → ⃗ l.To establish the justification behind these assumptions (namely, for large n, ⃗ l one way is to give numerical evidence.But in principle it is never complete with a finite number of states.That is why we now follow a different method comprising of two steps to establish the phenomena of homogenization.This method needs numerical support only for six ( in fact, it reduces to five since, for ⃗ l = ẑ, the aforesaid convergences are immediate) states.In the following we elaborate upon the steps of the method.

First step
In this step we show (with numerical help) that for a particular initial state of the ancillas the convergence of the set of recurrence relations occurs for all initial system states.Recall the dynamical map where, n , as given in Eq. ( 6).Here, again we will be applying the scheme 1 for tracing out the bath degrees of freedom.First we choose the following specific initial ancilla state, This means, from previous notation (Eq.( 21)), ⃗ l = ẑ.We show that for this particular ancilla state, any initial system state approaches to it in the asymptotic limit.For this we take the numerical support for the initial system states whose Bloch vectors are along x, y, or z directions and then generalize it to arbitrary system states.This means we start with the following, The validity of the above equation is supported by the numerical evidence as shown in the plots of Fig. 2a, 2b and 2c, where we plot the fidelity of the states ρ S n and η (given in Eq. ( 48)) with n, for ρ S 0 = 1 2 (1 ± σ i ), (i = x, y, z).Fidelity between two quantum states ρ and σ is given as F (ρ, σ) = Tr{ρ 1/2 σρ 1/2 } [58].The plot clearly shows that the initial system state reaches the initial ancilla state in the large n limit.In the figures we have taken different values α for better demonstration, as for same values of α, the curves almost overlap and therefore indistinguishable.This has been followed in later plots also.One may notice that behavior of the plots is similar for x, and y direction but different for z direction.When initial system state Bloch vectors are along either ±x or ±y direction, all of them make same angles with the Bloch vector of the initial state of the ancilla qubit.So it gives us same behavior of the system qubit for the x and y directions.But when the initial system state is along −z direction the behavior of the fidelity is different due to the difference in angle with the x and y directions.Furthermore, +z direction initial state of the system is already homogenized so the fidelity is at the value 1 from the very beginning.Though we have shown the plots for some particular values of α, and δ, we have numerically checked for a large number of randomly generated values of α, and δ, and in each case homogenization is achieved.This is illustrated in Fig. 3, where we plot the final fidelity of the initial system state ρ S 0 = 1/2(1 + σ x ) after n = 2000 collisions, for a 10 5 set of randomly generated values of α and δ.
This plot shows that in the large n limit, the initial system state converges to the initial bath state for all values of α and δ.Similar conclusions are obtained for other five (essentially four as the state 1/2(1 + σ z ) is already homogenized) states.We have not shown these latter numerical simulations for brevity.
Nevertheless, one can ask the question, whether for some parameter values homogenization is not achieved owing to different behavior of fidelity.To justify that it is not the case, we have provided three plots in Fig. 4. In Fig. 4(a) we plot the fidelity with number of collisions for five different values of α spread across the range of α (0 to π) with a fixed value of δ.We notice that if we increase the value of α, the rate of homogenization increases upto α ∼ π/2 and then decreases.This can be understood from the fact that for α around π/2, the probability of swap (cf.Eq. ( 22)) in the system-ancilla interaction operator is maximum causing fast homogenization.Basically, the homogenization rate roughly follows the periodic trend of a sine curve with varying α.So the nature of the fidelity as a function of α is expected to be regular indicating no exceptional behavior.Similar but opposite nature is observed for the parameter δ, as shown in Fig. 4(b), where we plot the fidelity with no. of collisions for four different values of δ with fixed α.One can notice the opposite behavior that upto δ ∼ π/2, the rate of homogenization decreases and after that it increases.This can be understood from the fact that from δ ∼ 0 upto δ ∼ π/2 ancilla-ancilla interaction (cf.Eq. ( 33)) increases and as a result the rate of homogenization decreases.This fact is also explained in the introduction.Hence, the rate of homogenization is essentially a competition between these two parameters.As shown in Fig. 4(c), we have plotted two extreme cases with small α, large δ, and large α, small δ.As expected we notice extreme slow and extreme fast homegenization respectively.If we take values of α and δ in between, the rate of homegenization falls between the last two curves also shown in the same plot.We have also checked this fact for different system and ancilla states and similar pattern is observed.These numerical results provide justification for Eq.(49).Now to generalize this result for an arbitrary initial state, first observe, In the last step, we have used the linearity of the map Φ n and the Eq. ( 49).Equipped with this, for any initial system state, ρ S 0 = 1 2 (1 + ⃗ k (0) .⃗σ), we can write, In the second equality we have used Eq. ( 50), and in the next line we have again used the linearity of the map Φ n along with Eq. ( 49).Now, from our previous notation, Using this and the fact that sin α ̸ = 0, from Eq. ( 42), in the asymptotic limit, we get, As ( ⃗ l = ẑ.Moreover, it is easy to note that, using lim n→∞ ⃗ l (0) n = ẑ, and lim n→∞ ⃗ k (n) = ẑ, Eq. ( 44) and Eq. ( 45) give lim n→∞ ⃗ l (1) n = ẑ, and lim n→∞ ⃗ l (2) n = ẑ.Evidently, the last step remains is to generalize this for any initial ancilla states and then the proof will be complete.Thus, we have proved here analytically -with the support of the numerical results provided in figures Fig. 2a, 2b, 2c, and 3 -that homogenization happens for every initial system state whenever the Bloch vector of the initial ancilla state is considered to be ẑ.

Second and final step
As mentioned before, the final step is to generalize the conclusion of the first step to an arbitrary initial ancilla state.In the previous step we took ⃗ l = ẑ and to generate any Bloch vectors ⃗ l from this with | ⃗ l| ≤ 1, we need two operations: rotation and scaling.First we apply a rotation matrix R on the initial ancilla vector ⃗ l = ẑ.In the set of recurrence relations, Eq. ( 42) to Eq. ( 45), we replace ⃗ l with R ⃗ l.This gives us, with, and, Operating R −1 from left on both sides of the above set of recurrence relations, and defining, n , and, ⃗ l n , we have the following set of relations, with, and, The above set of recurrence relations are exactly same as the one from Eq. ( 42) to (45).Now, we have already shown that with ⃗ l = ẑ, and any ⃗ k ′(0) , homogenization occurs, meaning, Next, we show that similar conclusion can be drawn by performing the scaling operation on the initial ancilla state Bloch vector ⃗ l = ẑ.Proceeding along the same lines as before, we replace ⃗ l by λ ⃗ l (with |λ| ≤ 1), in the set of recurrence relations given in Eq. ( 42) to Eq. ( 45).We define ⃗ k Then we get the following set of relations, Again, we have already shown that, for ⃗ l = ẑ, and any ⃗ k ′′(0) , ⃗ k ′′(n) n→∞ − −− → ⃗ l.From this, following the similar reasoning used previously, it is straightforward to show that, ⃗ l The above two results, one for rotation and another for scaling operation finally prove that for any initial system state and any initial ancilla state we achieve homogenization.
At this point we give some plots to demonstrate how fast or slow the homogenization takes place in comparison to the Markovian scenario.In Fig. 5, we take the initial system state and ancilla state both to be a pure state, such that there is coherence with respect to the computational basis |0⟩ and |1⟩.Here, we have taken the eigenstates of σ z (|0⟩,|1⟩) to be the computational basis, where, σ z |0⟩ = |0⟩, and σ z |1⟩ = − |1⟩.With explicit reference to system Hamiltonian, the computational basis can be taken to be the eigenstates of the Hamiltonian.We notice that scheme 1 and scheme 2 for non-Markovian dynamics give almost same homogenization profile, although the homogenization rate for non-Markovian scenario is much lower than the Markovian scenario.We have checked (numerically) that the rate is always slower than the Markovian counterpart irrespective of the values of α, and δ.In the next section we show that for a slightly modified ancilla-ancilla interaction we can achieve the homogenization rate almost as fast as the Markovian scenario.

Non-Markovian Homogenization with different kind of interaction
Keeping the system-bath interaction same as before (given in Eq. ( 22)), we introduce a different ancilla-ancilla interaction between the ancillas B n−1 and B n , which we call Partial S θ,ϕ (P-S θ,ϕ ), given as below, where, Note that for θ = π/2 and ϕ = 0, S θ,ϕ = S 4×4 -the SWAP operator.
Like in the previous sections we take the system and ancilla states in Bloch vector representation given by the expressions, Eq. ( 20) and Eq. ( 21).After the collision between the system and the first ancilla B 1 , reduced states of the system and the ancilla B 1 are given by ρ S 1 = 1 2 (1 + ⃗ k (1) .⃗σ) (see Eq. ( 24)) and η (1) 1 .⃗σ) (see Eq. ( 25)) respectively, where ⃗ k (1) , and ⃗ l 1 are given by Eq. ( 26) and Eq. ( 27) respectively.Next, the ancilla B 1 interacts with a fresh ancilla B 2 through the proposed interaction in Eq. (73).After this collision the reduced state of the ancilla B 1 is given by, Similarly, the state of B 2 is given by, The detailed expressions of ⃗ l (2) 1 and ⃗ l (0) 2 are quite long and are given in the Appendix A. Similarly, continuing the process iteratively as before, we end up with following set of recurrence relations, where, Also we have, with, ⃗ l (1) 1 Explicit expressions of the components of the vectors ⃗ l (2) n and ⃗ l (0) n are given in the Appendix A. Now, to show that the homogenization occurs for this process, we have to demonstrate that in the asymptotic limit, each vector on the left side of the above recurrence relations must converge to ⃗ l.To prove this statement, once again ratio test like Markovian scenario is not an ideal choice as we have to assume one or three vectors' convergence.Instead, like in the last two sections, we proceed through two steps.First step is to show for all initial system states and a particular ancilla state homogenization occurs and then in the second step, we generalize it to arbitrary ancilla states.
Before proceeding, we note an important observation.For this new type of ancillaancilla interaction P-S θ,ϕ of Eq. ( 73), the conditions for homogenization Eq. ( 19) are not satisfied for all states ρ (here the two states should be regarded as the states of two ancillas) unlike the PSWAP operation.If we take the states |0⟩ and |1⟩ (eigenstates of σ z ) to be the computational basis as before, only those states, which are diagonal with respect to this basis will satisfy the conditions of Eq. ( 19).Therefore, it is necessary that we have to take the initial ancilla state to be diagonal in this computational basis for homogenization to occur.So, in our subsequent proofs, we always take the initial ancilla states to be diagonal in the eigenbasis of σ z .Now to show the homogenization, we proceed in two steps as mentioned.

First step
This step is exactly similar as the step followed for homogenization with PSWAP.We choose the initial state of each ancilla to be η = 1 2 (1 + σ z ), meaning, in our notation ⃗ l = ẑ.If we take the initial system state whose Bloch vectors are along x, y or z direction then the system state approaches the initial ancilla state in the asymptotic limit.To put it mathematically, this means, Where Φ n is the dynamical map given in Eq. (7).Fig. 6a, Fig. 6b, and Fig. 6c, plotted for two different pairs of values of θ and ϕ, show the validity of the expression in Eq. (85).We have also checked numerically for a large number of randomly generated values of α, δ, θ, and ϕ, the system state approaches to the initial ancilla state for large n.In Fig. 7, we plotted the final fidelity for the initial system state ρ S 0 = 1/2(1 + σ x ) after n = 2000 collisions for 10 5 sets of randomly generated values of α, δ, θ, and ϕ.The plot illustrates that for all values of α, δ, θ, and ϕ, homogenization is achieved.Same conclusion is obtained for other five states also.
Nature of fidelity with different θ and ϕ is explained later.Next, similarly to the case of PSWAP scenario, in Fig. 8, we provide two plots for fidelity with no. of collisions.In Fig. 8(a), we vary α with fixed δ, and in Fig. 8(b), we vary δ with fixed α.Similar characteristics as PSWAP is also observed here, which provides a sound justification about no exceptional behavior of fidelity for any parameter values.Having justified this for six initial system states numerically, the next step is to generalize it for an arbitrary initial state of the system.For this, one proceeds along exactly the similar lines as in previous section and finally obtains, where, ρ S 0 is an initial system state.Therefore from our notation, lim n→∞ ⃗ k (n) = ⃗ l = ẑ.From which we can easily obtain lim n→∞ ⃗ l (0) n = ẑ, lim n→∞ ⃗ l (1) n = ẑ, and lim n→∞ ⃗ l (2) n = ẑ, just from the recurrence relations in a similar way as in the previous section.Next step is to generalize it for any initial ancilla state -diagonal in the computational basis.

Second and final step
To generalize the result of the very previous subsection to any initial ancilla states which are diagonal in the computational basis, we only need to consider the scaling operation on the Bloch vector as rotation operation will inevitably introduce off-diagonal elements.Similarly as before, we replace ⃗ l by λ ⃗ l (with |λ| ≤ 1) in the set of recurrence relations given in Eq. (81-84) and Eq.(A.7-A.12) of the Appendix A and then define, n , and ⃗ l With these notations, the recurrence relations now become, with, The detailed expressions for the components of ⃗ l * (2) n and ⃗ l * (0) n are long and provided in the Appendix B. Now, we have numerically demonstrated that for any initial system state Bloch vector ⃗ k * (0) , and initial ancilla state Bloch vector ⃗ l = ẑ, ⃗ k * (n) n→∞ − −− → ⃗ l.From this, following the similar reasoning used previously, it is straightforward to show that, (91) So, we have shown that with the new ancilla-ancilla interaction P-S θ,ϕ , homogenization occurs provided we take the initial ancilla state to be diagonal in the computational basis.
As explained earlier, with explicit reference to the system and ancilla Hamiltonians, this scenario can be regarded as thermalization, a special case of homogenization.Note that thermalization under some non-Markovian open quantum system dynamics, turns out to be a resource for achieving better performance of some quantum heat engines (see, for example, ref. [59,60]).Now, to see how fast the homogenization occurs comapred to the Markovian or the non-Markovian scenario with PSWAP, we provide plots of Fig. 9a and 9b for illustration.We notice that the rate of homogenization depends on the values of θ, ϕ, and the intial ancilla states η.We keep the value of α and δ unchanged for a fair comparison with the PSWAP case.In Fig. 9a, for the initial ancilla states η = 3 5 |0⟩ ⟨0| + 2 5 |1⟩ ⟨1|, two situations arise for two different pairs of θ and ϕ.For θ = 0.40, and ϕ = 0.15, the rate of homogenization for P-S θ,ϕ is higher than PSWAP and almost same as Markovian scenario.Whereas for θ = 1.60, and ϕ = 1.15, the homogenization rate of P-S θ,ϕ is still higher than the PSWAP but visibly lower than the Markovian scenario.We have checked that for this initial ancilla state, no values of θ and ϕ can give rise to other scenarios like rate of homogenization for P-S θ,ϕ is higher than Markovian or lower than PSWAP.This changes when we take a different ancilla state.In Fig. 9b, we take the initial ancilla state to be η = |0⟩ ⟨0| and the same two pairs of θ and ϕ.Now, for θ = 0.40, and ϕ = 0.15, the homogenization rate for P-S θ,ϕ is visibly greater than Markovian scenario.On the other hand for θ = 1.60, and ϕ = 1.15, the rate of homogenization for P-S θ,ϕ is marginally higher than the PSWAP in the asymptotic limit.Again, no values of θ and ϕ can give rise to the situations of plot 9a.To summarize, homogenization rate of P-S θ,ϕ is never less than the PSWAP, but depending upon the values of θ, ϕ, and η, it can be higher, almost same or lower than the Markovian counterpart.In conclusion, with this second type of ancilla-ancilla interaction, we can achieve faster homogenization rate than that of the PSWAP.Of course, to compare between Markovian and non-Markovian scenario we have to choose same α, and δ for both the cases.

Conclusion
In this paper, we study homogenization for non-Markovian collisional model.Non-Markovianity in the dynamics has been incorporated by introducing an additional ancilla-ancilla interaction in the standard collisional model.We start with Markovian scenario and show that homogenization is achieved by employing a new technique which relies upon the ratio test for convergence.For Markovian case, PSWAP is the unique system-ancilla unitary for universal homogenization, meaning, it occurs for any initial system and ancilla states.Subsequently we consider the non-Markovian scenario where both system-ancilla and ancilla-ancilla interactions are taken to be PSWAP.We numerically illustrate that homogenization occures in this situation too.Additionally, due to the use of PSWAP operation, homogenization is also universal here.But, interestingly, the rate at which homogenization is achieved is much slower than that in the Markovian scenario.To improve the rate we consider a new ancillaancilla interaction P-S θ,ϕ which is a modified version of PSWAP.First, demonstrating homogenization numerically, we notice that the rate of homogenization is always higher (in the asymptotic limit) than the PSWAP.Depending upon the values of θ, ϕ and the initial ancilla states η, the rate can be higher, almost same or lower than the Markovian counterpart.But the price we pay for this is that it loses the universality.Homogenization occurs only for those initial ancilla states which are diagonal in the computation basis.With explicit mention to the system and ancilla Hamiltonians, in a particular scenario, we may see this as the phenomenon of thermalization.