Thermalization of Finite Many-Body Systems by a Collision Model

We construct a collision model description of the thermalization of a finite many-body system by using careful derivation of the corresponding Lindblad-type master equation in the weak coupling regime. Using the example of a two-level target system, we show that collision model thermalization is crucially dependent on the various relevant system and bath timescales and on ensuring that the environment is composed of ancillae which are resonant with the system transition frequencies. Using this, we extend our analysis to show that our collision model can lead to thermalization for certain classes of many-body systems. We establish that for classically correlated systems our approach is effective, while we also highlight its shortcomings, in particular with regards to reaching entangled thermal states.


Introduction
Computer simulations of finite many-body systems have been challenging and expanding predictions of statistical mechanics since their first application to test equilibration of an anharmonic crystal modeled by a chain of masses with fixed-ends [1]. While standard methods to investigate equilibration and thermalization of quantum systems are based upon master equations [2], so called quantum collision models are introduced as versatile computational tools for simulating and studying open quantum systems [3,4]. The simplest collision model consists of a two-level system undergoing repeated collisions with environment, or ancilla, two-level systems. It is equivalent to a discrete time Markovian master equation in Lindblad form for the dynamics of the system, for short collision times [5]. Here, we address the question of how to generalize the collision models to finite quantum many-body systems for illuminating their thermalization dynamics.
Intuitively, it is reasonable to obtain a Markovian dynamics of the system using collisions if the colliding ancillae do not interact with any other degrees of freedom since such short time interactions should not allow any significant memory effects. However, the often implicit assumption of stronger interaction than the system Hamiltonian and the neglecting of the bath Hamiltonian are not always valid. Furthermore, using the typical formalism, e.g., [5,6] where energy preserving exchange interactions are considered, results in a dynamics which drives the system to the same state as ancillae, meaning that the result is independent from the system Hamiltonian and homogenization, rather than thermalization, is achieved [7,8]. This problem persists and is compounded for the generalization of collision models for many-body systems [8,9]. Interestingly, [10] derives a Lindblad type master equation for collisions with arbitrary interaction strengths and collision times and establishes that the thermal state of a system at the environment temperature with respect to the HamiltonianĤ 0 is an equilibrium state if [Û,Ĥ 0 +Ĥ bath ] whereÛ is the unitary evolution operator under the total Hamiltonian. However, settingĤ 0 =Ĥ system and finding the necessary interaction type and strength to satisfy this commutation property remains as a challenging open problem so far. At variance with this and other works that study collision models starting from a "global" unitary picture [9,10], in this work we propose a master equation derivation inspired by the well-known derivation for a time-independent system-bath interaction in the weak coupling regime [2].
Despite its drawbacks in describing Markovian open system dynamics, quantum collision models are still a good candidate for understanding the quantum thermodynamical phenomena from a microscopic perspective [11]. For example, the microscopic Markovian master equation derivation in [2] does not account for the information loss of the system about its initial state, while it is evident using the collision model that the lost information is kept by the entanglement between the system and ancillae [12]. Another study analyzes the entropy generation and distribution in a collision model and proves the asymptotic factorization of the total density matrix of system and environment into the density matrices of the system and the environment for a two level system in the strong coupling regime [13]. More complex collision models involving ancilla-ancilla collisions allow for the derivation of completely positive non-Markovian dynamics [14][15][16]. The controllable degree of non-Markovianity and its effect on the dynamics of quantum coherence has been examined [17]. Further attempts to study non-Markovian dynamics are made by using initially entangled ancillae [18], introducing time overlap of two consecutive collisions [6,19,20] and using a two-spin system in which only one of the pair interacts with the environment resulting in a Markovian dynamics for the composite system, while tracing out the spin interacting with the bath gives a non-Markovian dynamics for the remaining spin [21]. The versatility of collision models has resulted in other interesting research directions, such as the introduction of collisions with non-thermalized ancillae to study non-equilibrium effects in quantum thermodynamics [8,11,[22][23][24][25] and the generation of multi-qubit entanglement via a shuttle qubit colliding with disjoint qubit registers [26].
This work aims to examine the conditions required for thermalization in a Markovian collision model using two level ancillae. To this end, we will first carefully examine the microscopic derivation of a Lindblad master equation for a two level system in the weak coupling regime from [2] and introduce a time dependent interaction Hamiltonian in Section 2, where we also assess each assumption made for the derivation and examine their validity. Section 3 examines our collision model for many-body systems for both non-entangled and entangled energy eigenstates with an example for each of these cases illustrating how our proposed collisional route to many-body thermalization works. Finally, we conclude in Section 4.

Derivation and Validity of Lindblad Master Equation
We begin by following the microscopic derivation of the Lindblad master equation given in [2], however allowing for a time-dependent interaction Hamiltonian instead of using the second order approximation of the unitary evolution operator for the system and the ancilla with respect to the collision time [21]. The dynamics of the system and the bath is governed by the Liouville-von Neumann equation Integrating this equation with respect to time and plugging in the expression for ρ(t) in the commutator twice with the assumption of Tr B ([Ĥ I (t), ρ(0)]) = 0 we arrive at Applying the Born approximation by neglecting system-bath entanglement and the effect of the system on the bath allows us to write an equation for the dynamics of the system by tracing over the bath degrees of freedom At this point, the dynamics of the system is still, in general, non-Markovian and we have not made any explicit assumptions about the nature of the interaction. However, the finite time of a given collision may serve to justify the constancy of the bath density matrix along with the weak interaction assumption. Putting aside the validity of Born approximation, we need to explicitly assume that the density matrix of the system does not change significantly during the interaction with a single ancilla, which is justifiable for short collision times, in order to replace the past states of the system with its present state and to obtain the Redfield equation The standard master equation derivation in [2] for a time-independent interaction Hamiltonian continues with the assumption that the integrand above vanishes quickly enough to extend the integral to infinity with negligible difference on the system dynamics. In our case of short time collisions starting after t = 0, this extension is not an assumption to be checked as the integrand is explicitly zeroed out for s > t by the time-dependent strength of the interaction Hamiltonian. For simplicity, we assume that each ancilla interacts with the system once and these collisions start with a period of τ p and a duration of τ c . After explicitly defining our collision model, we can investigate the effects of the finite time interactions on the dynamics. As in the derivation in [2], we will introduce the interaction Hamiltonian in the Schrödinger pictureĤ where the Hermitian operatorsÂ α andB α act on the system and the bath respectively. After decomposing the operatorsÂ α into operatorsÂ α (ω) based on the energy transitions with frequency ω generated on the eigenstates of the system Hamiltonian and plugging the interaction picture interaction Hamiltonian in Equation (4), we obtain d dt where Γ αβ (ω) is the one-sided Fourier transform of the reservoir correlation functions where operators are defined in the interaction picture. For the evaluation of bath correlation spectra, we must specify our open system setup which consists of the same basic ingredients as [6,19,20,25]. We first consider a two-level system with time-independent HamiltonianĤ The reservoir consists of, an in principle infinite number of, two-level systems prepared at an inverse temperature β b =1/(k B T) for a bath Hamiltonian where the index n indicates that the operator acts on the n-th spin of the reservoir. The time-dependent interaction Hamiltonian in the Schrödinger picture is given bŷ where the operator without index acts on the system. For simplicity, we assume that the interaction strength is exactly zero before and after the interaction, remains constant during the collision, and has the same magnitude for all collisions. It should be noted that the interaction in Equation (10) is different from the often considered partial swap case which is known to lead to homogenization [4] rather than thermalization [7]. Knowing the collision period and duration, we can now define the time-dependent interaction strengths as where θ denotes the Heaviside step function. Before explicitly calculating the bath correlation spectra, we can make some simplifications. As each ancilla has one interaction component in the form of Equation (5), the indices α and β in fact denote the index of the corresponding ancilla. Also, knowing that all ancillae are prepared in a thermal state, it is easy to prove that cross correlations vanish and we can arrange Equation (6) After some manipulation, we find the explicit form of reservoir correlation spectra where ρ n ee and ρ n gg are excited and ground populations of n-th ancilla. It is clear that the bath correlation spectra are time-dependent and they are zeroed out by the step functions before or after the collision. We must evaluate this expression for the cases ω = ±2h b and ω = ±2h b separately, If ω = ±2h b , one of the complex exponentials in the integrand simplifies and gives a linearly growing term The final step of the derivation of a Lindblad type master equation is the decomposition of bath correlation spectra into its real and imaginary parts. The imaginary part results in an additional Hamiltonian term, the Lamb shift acting on the system. However, as this is not relevant to the equilibration of the system, we will neglect it in what follows. As we explicitly show in Figure 1 it is also reasonable to neglect situations where ancillae spins are not on resonance with the system, i.e., we only consider h b = h s . In this case, the bath correlation spectra consists of a real and linearly growing term and a rotating term with real and complex parts. The linearly growing term generates a dynamics similar to a Lindblad master equation with time-independent interactions, while the real part of the rotating term can be neglected assuming that the relaxation of the system is much slower than the dynamics of the closed system. The master equation in Lindblad form can be obtained after applying these assumptions to Equation (12) together with the secular approximation resulting in (16) where the Γ function contains the information about all of the collisions (17) where N denotes the number of ancilla spins. The function δ (ω) is defined as one for ω = 0 and zero elsewhere, not to be confused with Dirac delta function. Note that this equation neglects the case where the ancilla is not in resonance with the system and it is used throughout Section 3. However, the off-resonance effects in Figure 1 need to be interpreted using the bath correlation spectrum described in Equation (14). Simulation results for the thermalization of a single two-level system using our collision model. We show the fidelity of the system with the target thermal state as a function of the bath ancilla frequency and number of collisions. We clearly see that thermalization occurs when the system interacts with bath frequencies that are on resonance. We have fixed T = 10 mK, g = 1 MHz, h s = 1 GHz, t = 200 ns, and ρ s (0) = |1 1|.
The transition from Equation (12) to Equation (16) takes the secular approximation for granted, however it can be justified by some assumptions relating three different time scales of the open system dynamics: The natural evolution times of the system and ancillae and the duration of the collision, all of which play a critical role in constraining the validity of the derived master equation. We assume that the interaction vanishes before any significant change on the density matrix of the ancilla can happen. We also assume that the variation of the system state during one collision is small, which further constrains the maximum collision time. On the other hand, we also want to eliminate the rotating terms of the bath correlation spectra by averaging them over multiple periods of the system dynamics with a slow relaxation of the system which leads to a lower bound of the collision duration.
After justifying the derivation of Equation (16), it is straightforward to find the Kubo-Martin-Schwinger (KMS) condition for n-th collision exploiting the fact that the ancillae are prepared in a thermal state, resulting in vanishing cross bath correlations.
The interpretation of this equation is obvious: The thermal state of the system at the inverse bath temperature β b is the unique steady state of the Markovian dynamics generated by collisions with ancillae prepared in thermal state [2]. This result was also predicted in complementary works on collision models [5,25] derived using different parameter regimes.
In Figure 1 we simulate our collision model sweeping through a range of frequencies for the bath ancillae and show the final state fidelity between the system and its target thermal state. The simulation consists of the unitary evolution of the system and ancillae during the collision time with the sum of system, bath, and interaction Hamiltonians described above and the ancillae are traced out after each collision without interacting again with the system or other with ancillae. We clearly see that when the ancillae are close to resonance the collision model leads to thermalization of the system. Conversely, when the ancillae are far detuned from h s we find the system dynamics are almost frozen. This result can be predicted theoretically by calculating the real part of bath correlation spectrum without assuming resonance. Equation (14) has two terms which are inversely proportional to the difference between the transition frequency ω and ±2h b . Assuming a small detuning from either 2h b or −2h b , the other term becomes negligibly small. After dropping the small term, evaluating the real part for the other part gives ignoring the Heaviside step functions and taking the beginning of each collision as t = 0. Its limit for δ → 0 recovers the case of resonance. The off-resonance dynamics depend heavily on the product δτ c . As the average of sine function over a period is zero, we can conclude that the effect of the dissipative term should be negligible if the product δτ c = 2kπ where k is an integer and the dynamics is slow enough. On the other hand, in the case where the product is an odd multiple of π, the average of sine function is not zeroed out and we observe thermalization as seen in Figure 1. Furthermore, it is straightforward to prove that the fastest thermalization is achieved in the case of resonance using the identity sin(δt) The results in Figure 1 confirm the range of validity of our master equation and are in keeping with other results in the literature [25]. Furthermore, the clear importance of on-resonance ancillae indicates that, under suitable constraints, only particular bath frequencies are relevant for ensuring the system thermalizes. Thus, we can exploit this feature to explore the requirements for achieving thermalization for many-body systems.

Classically Correlated Systems
Let us consider the 1D Ising chain described by the Hamiltonian As stressed in the previous section, to achieve thermalization we require the driving frequency of the system and the ancillae to be the same. In the case of interacting many-body systems, it should be clear that there will be a range of frequencies, each of which will be related to the various transition frequencies of the many-body system. Thus, to examine the requirements to reach thermalization we use the expression of the interaction Hamiltonian in the form where sum over ω denotes the decomposition of each spin-ancilla collision operator into the different energy transitions it generates. We can make a temporary simplification to make the illustration of many-body system thermalization easier by replacing the ancillae with a set of harmonic oscillators forming a continuous spectrum prepared at an inverse temperature β b = 1/(k B T). In this case, we can find the energy transitions generated by each term of the interaction Hamiltonian with a partition of the Hilbert space of the whole system based on each nearest neighbor configuration with respect to a reference spin denoted as i. We can write all terms of the system Hamiltonian involving i-th spin aŝ where i = 1, N as the first and last spins of the Ising chain do not have a left and right neighbor, respectively. The Hamiltonian at the end points i = 1, N can be found by omitting the term corresponding to the lacking neighbors i = 0, N + 1 in the above equation.
We can now define the transition frequencies generated by flipping the i-th spin in terms of the state of neighbor spins Decomposing the operatorσ x asσ we obtain two dissipators for each term of the interaction Hamiltonian. The frequencies in Equation (24) correspond to the transitions generated byσ −i while their negatives correspond toσ +i . Expressing the frequencies as a function of nearest neighbor configuration for each spin results in the master equation Here, {s i } is a short hand notation for the respective arguments of the frequencies in Equation (24), corresponding to the set of basis vectors of the Hilbert space of the neighbor spins of i-th spin. The notationσ s i ±i implies that this operator can be decomposed aŝ and D(ρ,ô) is defined by Although Equation (26) is derived for a continuous set of harmonic oscillators, each term appearing in the double sum is similar to the master equation for a two-level system with the driving frequency depending on the nearest neighbor configuration. Therefore, the implementation of a similar master equation with collisions generating one spin flip operations with ancillae driven at the frequencies of single spin transitions, as illustrated in Figure 2, is possible if the secular approximation is valid such that the ancillae cannot generate any transitions other than those corresponding to its driving frequency. The results of Section 2 on the KMS conditions for the bath correlation spectra can be generalized for the master equation of 1D Ising model and this ensures that if all ancillae are prepared at an inverse temperature β b , the thermal state of the system at the same temperature is a steady state of the master equation [2]. However, the uniqueness of the stationary solution requires additional constraints. A sufficient condition for the uniqueness can be stated as follows [27,28]: Condition 1. Let L be the Lindblad superoperator describing the time derivative of the density matrix and σ ±i (ω(s i )) operators the generators of L. The dynamical semigroup generated by L is relaxing in the sense that it drives the density matrix to a unique final state as time tends to infinity regardless of the initial state if the linear span of the generators is an adjoint set and the bicommutant of the generators is the set of all bounded operators acting on the Hilbert space of the system. In order to check the applicability of Condition 1 to the thermal bath with local system-bath interactions, we start by checking the adjoint property of the linear span of generators. As established in [27], this follows sinceσ +i (ω(s i )) =σ † −i (ω(s i )), meaning that the adjoint of each generator is also a generator. The second property is easy to prove using the fact thatσ ± operators only commute with themselves and the identity operator and the only operator commuting with allσ ±i (ω(s i )) for all i and s i is the identity operator.
To simulate thermalization for a two-site Ising model, Equation (21) with N = 2, we require collisions corresponding to the one-spin flip transition frequencies as illustrated in Figure 2a. As each of the spin has a single neighbor, there are two nearest neighbor configurations, resulting in a total of four energy transitions for the whole system. For larger systems, each spin in the bulk of the chain has four different energy transitions and requires more ancillae to successfully thermalize, as shown in Figure 2b.
We implement our collision model for the two-site Ising chain, considering when the collisions with the various ancillae happen "sequentially", i.e., the whole system collides with one of the ancillae corresponding to one of the energy transitions at a time and the colliding ancilla is subsequently traced out before the next collision occurs. We also consider "simultaneously" occurring collisions where the whole system interacts with all of the four ancillae corresponding to different energy transitions at once, after which they are traced out. The minimum energy states are up-down and down-up states and these states cannot be prepared by a local master equation as the collisions are identical, verifying the effect of the system Hamiltonian on open system dynamics resulting in a global master equation. In Figure 3 we show that our collision model gives rise to thermalization for interacting systems. Furthermore, as the cross bath correlations vanish for a thermal bath, we expect that a time overlap between the collisions (such as that which occurs in the simultaneous collision case) does not change the form of the equation, and our numerical results confirm that both approaches generate an almost identical evolution. Figure 3. Simulation of a 2-spin Ising model with parameters h 1 = h 2 = 500 MHz and J = 1 GHz and corresponding transition frequencies ω 1 = ω 3 = 1.5 GHz and ω 2 = ω 4 = 0.5 GHz. All ancilla-system spin coupling strengths are set as 1 MHz and the collision times are fixed as 400 ns. Fidelity after each step consisting of one collision for each one spin transition frequency with respect to the thermal state of the system at the temperature of ancilla spins T b = 10 mK. The initial state is the thermal state of the system at infinite temperature.

Entangled Systems
The Ising model considered in the previous analysis has eigenvectors which are product states without any entanglement among the spin sites. In this section we elaborate on the validity of our collision model for realizing thermalization in more generic many-body systems, particularly those that exhibit entanglement. Addressing such an issue in full generality is a formidable task. Indeed, unlike in the case of non-entangled eigenstates where the generation of single-spin transitions for each interacting neighbor configuration was sufficient, even determining the minimum necessary number of collisions for the uniqueness of the equilibrium state is difficult for entangled states. As such we will restrict to a specific example in this section.
We begin our discussion by reminding that the matrix representation of any Hamiltonian has an eigenvalue decomposition in the form where D is a diagonal matrix with the values of eigenenergies on the diagonal, U is a unitary matrix such that its columns are the eigenstates of the Hamiltonian. Following our master equation derivation, each term of the interaction Hamiltonian is decomposed into different energy transitions, giving rise to operators in the formÂ where ψ k denotes the k-th eigenstate of the Hamiltonian, which is also denoted by the k-th column of the matrix U. This simple form of the energy transition operators can also be expressed in the basis consisting of the Kronecker product of the bases of the subsystems aŝ where N is the dimension of the Hilbert space of the system and states i and j are selected from the basis constructed by the Kronecker product of the subsystems, therefore these states are not entangled.
Knowing that the k-th column of the matrix U is equal to ψ k , we can write where U ab denotes the element of U at the a-th row and b-th column.
The existence of coefficients a ij,kl indicates that there is a one-to-one linear map from the vectors in the basis of eigenstates to the vectors in the Kronecker product basis. Furthermore, we can vectorize the indices i and j into one index u and the indices k and l into another index v. By these reductions, we can express our linear map in the form of a matrix M such that where the vectors A and x are the vectorized representations of an operator in the basis of eigenstates and in the Kronecker product basis, respectively, with the sum running over the repeated index v. As the matrix M represents a one-to-one linear map, its inverse exists and any vector A u can be expressed in the form Using this expression, we can conclude that a single subsystem transition generated by the interaction Hamiltonian, which can be expressed with a vector x v , with one non-zero element can generate multiple energy transitions by the multiplication by the inverse of the matrix M used for conversion into eigenstate basis. In this case, the thermalization conditions depend on the structure of the matrix M, however the thermalization of any many-body system is in principle possible with a sufficient number of energy transitions generated by the collisions with two-level ancillae driven at the corresponding transition frequencies, the appropriate choice of interaction Hamiltonian, and the validity of our assumptions for the master equation. As a concrete example consider a two-spin anisotropic XY-model with Dzyaloshinskii-Moriya (DM) interaction in z-direction with the eigenstates and eigenenergies [29] |ψ 1,2 = |↓↓ ± |↑↑ √ 2 , |ψ 3,4 = |↓↑ ∓ i |↑↓ √ 2 ; Using the definition of operatorsÂ kl , we can express the one spin flip operators as We can then describe the spin ladder operators acting on the first site aŝ It is clear that theÂ 41 andÂ 32 terms of the ladder operators and their Hermitians generate state transitions with non-zero energy difference. Other state transitions are between the states with the same energy which cannot be generated via with collisions with ancilla spins which do not have internal energy as we have assumed h s , h b >> g. If the zero energy transitions were allowed, we could make transitions from any state of the system to another state using intermediate transitions, impling the uniqueness of the thermal state as the equilibrium point of the dynamics [30]. In our case, this condition is not satisfied, and this leads to the equilibrium state of the system exhibiting an initial state dependence. Our numerical simulations in Figure 4 show that if the system is initially prepared in some thermal state, but not in equilibrium with the bath, a Gibbsian thermal state at the environment temperature is achieved. However, it is not guaranteed for generic non-equilibrium initial states, such as |1 1|. We understand this as follows: the choice of initial state as a thermal state at some temperature guarantees that the population of the states having the same energy is equal, thus implying that the zero frequency transition terms will not contribute to the dynamics of the system even if they are generated by the collisions. This means we can assume that the zero frequency terms exist and consequently the equilibrium state is the thermal state at the environment temperature. Another possible issue regarding thermalization of entangled many-body systems by our collision model is the additional terms of the master equation due to the non-vanishing bath cross correlations arising due to the decomposition of each term of the interaction Hamiltonian acting on a single subsystem into multiple energy transition terms, which implies that the bath operator of those energy transition terms are the same. For this reason, the positive definiteness of the bath correlation matrix for every frequency needs to be asserted for the uniqueness of the equilibrium state [27].
In summary, our example of two spin anisotropic XY model shows that our collision model can generate multiple energy transitions without the explicit calculation of the M matrix. Although thermalization is not guaranteed, this analysis nevertheless provides insight about how an entangled many-body system with non-degenerate energy levels can be thermalized as long as the secular approximation used in the master equation derivation remains valid and the bath correlation matrix is positive definite.

Conclusions
In this work we have presented a collision model using two level ancillae that leads to thermalization in the weak coupling regime, even for certain finite many-body systems. By carefully assessing the relevant timescales present, we showed that when the ancillae are tuned inline with the transition frequencies of the system, thermalization can be achieved. This is at variance with other schemes commonly examined in the literature where system and environment interact via a partial swap [6,20]. Our master equation derivation for 1D Ising model can be straightforwardly generalized to N-dimensional spin lattices by redefining the sums over the Hilbert space of neighbor spins. In the case of Ising spin lattices with more than one dimension, the system can be tuned to be an integrable or non-integrable system depending whether the external magnetic fields are turned off or on respectively [31] and our collision model achieves thermalization in both of the cases. If the eigenstates of the system Hamiltonian are entangled, our collision model gives valuable insight on the dependence of equilibrium state on the initial condition; in particular reveals the conditions to engineer Gibbsian thermal state at the environment temperature. Remarkably, for entangled eigenstates, the decomposition of single-spin transition operators into multiple energy transition operators may remove the necessity of bath interaction with each spin in the system.
Beyond the clear interest in understanding the phenomenology of thermalization using a collision model and its possible extensions to non-Markovian and non-equilibrium dynamics, our collision model also can be viewed as a versatile and implementable artificial environment acting as a temperature knob, as similarly considered in [30,32]. Contrary to the artificial temperature knob proposal in [30], our proposal satisfies the KMS condition for thermalization instead of an optimized approximation depending on tunable system parameters and it is promising to be scalable for large many body systems. The proposal in [32] relies on a similar idea to our proposal; its authors propose to sweep all possible energy transitions of the system with a slowly varying bath Hamiltonian strength, which can be considered as a different way of obtaining the effect of ancillae colliding to a subsystem with different bath Hamiltonian strength. Obviously, making use of only relevant transition frequencies leads to much faster thermalization and it is possible to get rid of some timescale constraints of [32] as the ancillae are supposed to be prepared in a thermal state for a time independent bath Hamiltonian before the collision in our proposal.
Our proposal can also lead to the cooling of the target system if it is possible to keep ancilla spins colder than the environment temperature. Indeed we mention two possible methods of spin cooling for the preparation of a cold environment that our scheme is well suited to. The first one is the use of frequent measurements on a two-level system interacting with a non-Markovian environment which brings the mean energy of interaction Hamiltonian to zero in order to reduce the total energy of the two-level system and its environment [33]. The application of this idea may suffer from the challenges posed by the necessary minimum frequency of the measurements. Another idea is to use quantum coherent or entangled two-level systems [34][35][36] to engineer the temperature of a two-level target system, which can then be used as an ancilla for the many-body system to be thermalized.
Our results can have practical significance for suggesting design principles of quantum thermalizing machines for finite many-body systems. Such devices would be compact as they can consist of few ancillae as artificial environment; they would be fast as they can engineer the target thermal state with high fidelity after a small number of collisions describing a unitary route to thermalization. These properties can be valuable for quantum thermal annealing [30] and quantum simulation applications [37], for example using superconducting circuits.

Author Contributions:
Ö.E.M. and O.A. are responsible for the conceptualization, the formulation of methodology, and the theoretical analysis of the ideas introduced in the paper. S.C., Ö.E.M., and O.A. equally contributed to the discussions throughout the work, the numerical simulations, and the writing of the paper.