The multilevel four-stroke swap engine and its environment

A multilevel four-stroke engine where the thermalization stokes are generated by unitary collisions with bath particles is analyzed. Our model is solvable even when the engine operates far from thermal equilibrium and in the strong system-bath coupling. Necessary operation conditions for the heat machine to perform as an engine or a refrigerator are derived. We relate the work and efficiency of the device to local and non-local statistical properties of the baths (purity, mutual coincidence etc.). In particular, we relate the Clausius inequality to the symmetrized relative entropy of the baths (Jefferys divergence). Other Clausius-like inequalities are derived as well. Finally, in the ultra-hot regime we optimize the work of the multilevel engine and obtain simpler forms for the work and efficiency.


Introduction
Present day technology is on the verge of enabling different realizations of quantum heat machines where the engine core ("working substance") comprises of a single particle in a discrete level system (in particular a qubit). In analogy to their classical counterparts quantum heat machines can be used to cool or to produce work. Furthermore, similarly to classical thermodynamics, performance analysis of heat machines without going into the details of a specific realization is a key theme in quantum thermodynamics (QT). There are several approaches to QT . The more recent one comes from quantum information resource theory. In this approach thermal states are considered free while non-equilibrium states are considered as a resource ( [1] and references therein). An energy conserving unitary is used to couple the system and the bath. A different viewpoint to QT is dynamical, and uses the framework of quantum generators of open systems [2]. In this approach, models of quantum heat machine are constructed and analyzed. To a certain extent the present article combine the two approaches by analyzing the action of a heat engine constructed from discrete elementary unitary interactions with the reservoirs. A third approach called typicality aims to find the condition when complex interactions lead to thermal behavior for a typical set of states (see [3] and references therein).
The goal of the model studied here is to lead to a more comprehensive understanding of multilevel quantum heat engine operation and to study the byproducts of their operation. In addition, the open system approach is typically restricted to a weak coupling to the baths. Even this weak coupling should be carefully handled in order to avoid conflicts with the second law [2]. Our model is both solvable and consistent with the second law in all coupling strengths.
Various quantum engine models have been suggested -some are continuous (where the baths are always connected to the engine) and some are reciprocating (where in different strokes different processes takes place). The equivalence between a laser and the Carnot engine [4] has inspired the study of quantum heat engines. For example studying maximum power operation [5,6,7,8], the role of coherence and entanglement [9,10,11,12,13,14,15], quantum refrigerators and the third law of thermodynamics [16,17,18,19,20,21,22] and the connection to quantum information [23,24].
Different quantum working mediums have been considered, two-level systems, Nlevel systems or harmonic oscillators [25,26,27,28]. In most previous studies the system bath dynamics was modeled by a reduced description where the dynamics of the system is described explicitly by a Master equation with a generator cast into the Lindblad form [29,30]. In the reduced formalism the bath is only considered implicitly, therefore the effect of the engine on the baths is almost always ignored. In this paper we take the opposite point of view and study the heat machine operation using the properties of the baths. One of the main difficulties in the Lindblad approach is that the reduced dynamics generators are not uniquely defined. Various choices of the Lindblad generator lead to the same reduced dynamics of the working medium but to different entropy fluxes in the baths. As it turns out the Lindblad operators have to be carefully chosen in order to be consistent with the second law of thermodynamics [2]. The problem becomes more severe when the time dependent external fields drive the system [31,32].
The analysis carried out in this paper overcomes this problem by considering a very simple yet physically plausible bath model that is based on collisions [33,34]. In this model the bath consists of non-interacting particles which are initially in the same thermal state. The engine interacts with the bath particles via a two-particle collision described by a unitary operation. We choose to model the collision by the unitary swap operation but other unitary operations such as CNOT and random unitaries have been considered as well [35,36] even though they were not applied to the study of heat machines. As will be discussed later on, the swap operation has two interesting physical limits. The first is the full swap that describes the exchange of particles and the other is the weak partial swap that mimics quasi-static evolution.
Typically the system that interacts with the bath is taken to be a single qubit but qudits and a chain of coupled qubits have been considered as well [37,38]. The collision model can also represent a general non-Markovian dynamics [39,40]. If, however, the particle that interacted with the system is extremely unlikely to interact again with the system or the bath particles, then the dynamics is Markovian. See [41,42] and reference therein for more studies on collision-generated Markovian dynamics. In this article we consider this type of Markovian collisions in order to study the most basic features of an engine driven by collision baths.
The importance of studying multilevel heat machine is two-fold. First, experimentally speaking, in some systems it is not possible or not practical to interact with only two levels. Secondly, two-level systems may contain some non-generic features and it important to understand the more general framework of quantum thermodynamics. For example, in a two-level system temperature is always well defined if the density matrix is diagonal in the energy basis. Clearly this is not true if there are more than two levels. While in a two-level engine it is fairly simple to write down the map the parameter space where the system operates as an engine and as refrigerator, this situation is considerable more complicated when there multiple levels. The different operation regimes are determined by the temperature and the gap energy associated with the cold and hot strokes. It is not straightforward what energy scale replaces the two-level energy gap in the multilevel case (the variance? the maximal gap?). Furthermore it is not clear at this point if necessary and sufficient conditions for engine (or refrigerator) operation can be formulated without explicitly solving the full dynamics of the system.
The present analysis can be related to "finite time thermodynamics" [43,44]although there is no explicit reference to time. Implicitly the cycle can be related to the collision timscale and in addition there is no assumption of thermal equilibration of the device with the baths upon contact.
The article is organized as follows: Section 2 contains the description of our engine and bath model. In Sec. 3 we analyze the steady state operation of the engine and discuss the difference between quantum swap and classical random swap. Next, Sec. 4 explores various aspects of the engine evolution: Clausius and, generalized Clausius inequality and the baths purity reduction. In Sec. 5 we set upper bounds on the work and efficiency in term of statistical quantities like purity, mutual coincidence and Wooters distance. Finally, Sec. 6 presents two possible regimes of operation: the ultrahot bath regime and the quasi-static regime.

The bath
Our bath model is based on the "quantum homogenizer" introduced in Ref. [33,34,45]. This bath contains a very large amount of single species non-interacting multilevel particles. At the beginning all the particles are in the same thermal state. However, each time a bath particle interact with the engine the population of the particle and of the engine changes according to the simple partial swap rule that will be described shortly. Every time the system (engine) interacts with the bath it collides with a new thermal particle in the bath. If the system is coupled only to one bath then repeated collision would lead to thermalization of the engine particle at as exponential rate. In this work, however, we will assume that there is only one collision or none at all at each thermal stroke of the engine. Therefore, the engine will typically not be in a Gibbs state. The extension to multiple collisions is straightforward.
The bath is assumed to be large enough so that the probability to re-collide with a bath particle that already interacted with the system is negligible. Consequently, the resulting dynamics is Markovian. We assume that the particles do no interact with each other in the bath but only with the system particle at the bath-system interface. That is, in our model there is no thermalization inside the bath. In principle, for large baths, this has no impact on the engine's operation. If the bath is large enough it does not matter if the scattered particle get thermalized again or not as the engine will (almost) always interact with a new thermal particle. Nevertheless, it is interesting to explore the operation of small baths with and without intra-bath thermalization.
There is a different class of quantum engines where the baths are always connected to the engine but different levels are connected to different baths [16]. This implies that if only one bath is connected the relative probabilities of the specific coupled levels will be thermal but the engine as a whole will not be in a Gibbs state. Alternatively stated, these bath models have a continuum of steady states. In our model the Gibbs state is the only single bath steady state.

The engine's cycle
Our engine comprises of a multilevel system that is driven to a four stroke Otto cycle [46]. Figure 1a illustrate the engines operation for a two-level system. In stroke A the gap spacing is increased by applying an external field. The engine is decoupled from the baths at this stage. The state evolution is such that the system's density is diagonal in the energy basis at the beginning and at the end of the stroke and the populations do not change. Adiabatic change of the Hamiltonian will achieve such a transformation but other faster options exist (see Sec. 2.5). Nonetheless, this part of the cycle is termed "adiabatic".
In stroke B the gap is kept fixed in time, but the engine is allowed to interact with the hot bath. During this time interval there may be one collision, multiple collisions or none at all. The average collision rate is a parameter that characterizes the bath and it's coupling to the engine. It encapsulates within the particle density in the bath, the engine cross-section for collision etc. The collision entangles the engine and the particles in the bath. However since the interacting bath particles will not interact with the the engine again (the top spheres with temperature denoted by primes), we consider the reduced dynamics of the engine by taking a partial trace over the scattered particle of the bath. This stage changes the entropy of the engine. In stroke C the external field is changed again so that the gap returns to its initial value. Finally, in stroke D, the system is coupled to the cold bath.
In the multilevel case the evolution can be considerable more intricate. Firstly, we do not assume that all the levels are increase or decreased by the same ratio (Fig. 1c). Secondly, we allow the levels to cross (Fig. 1d).

The collision model
We assume that at the initial state the system and the bath particles are in a product of thermal Gibbs state: The energy of the ground state and excited state of the two level system that comprises four-stroke, two-level, partial swap collision driven engine. (1b) Dynamic of the ground state population. In strokes A and B the population is fixed and the energy gap changes and in stroke C and D the gap is fixed and the population changes via collisions with the baths particles. In general the collision stage is not an isotherm. The work is done during stroke A and C and heat exchange takes place during strokes B and D. (1c) In a multilevel the level dynamics can be more intricate.
The levels do not have to be compressed by the same factor. Furthermore they can even cross each other (1d) as long as it does not conflict with energy population invariance in stages B and C.
where 'b' stands for 'c' (cold bath) or 'h' (hot bath). The free energy F b takes care of normalization. Let U be some two-particle unitary operation that conserves the total energy of the two particles. The reduced density matrices after the interaction are: When the collision Hamiltonian that generates U is invariant under the transformation b ↔ s one can show that: This condition, together with the total energy conservation, implies that the energy levels of the bath and the system have to be equal.

Density matrix swap and energy population swap
In this work we focus on collision that induce a single parameter convex transformation for reduce density matrices or for the energy population. The density swap rule is: where x is the swap parameter. For x = 0 nothing happens and for x = 1 the density matrix of the bath and the particle completely swap. It is easy to verify that this transformation satisfies (6) and that the only steady state of this transformation is Re-colliding the system particles with new bath particles in a thermal state will lead to exponential decay of ρ s towards ρ b . Yet, in this work, we do not assume that equilibrium is reached. In fact we will mostly consider a single collision at most in a given stroke (B or D).
Although the density swap rule has a very simple structure our framework is valid for a more general type of swap operation. Let ρ b and ρ s be diagonal matrices in the energy basis, where p b and p s denote the energy level populations (the diagonal terms in the density matrix) of the bath and system particle. We shall use bold face characters to denote quantities that have a level index such as probability or energies of a bath E b . The energy population swap rule is: where as before the prime relates to the traced-out outcome after the collision. We do not claim that the swap operation in its present form is general enough to describe many typical thermal baths used in experiments. Yet, such bath could in principle be built. One advantage of this bath model is that it leads to a solvable dynamics that can be used as a reference for other more complex model. Another advantage, as will be shown later on, is that it helps to formulate new questions and new points of view concerning the operation of heat machines. Nevertheless, in a two-level system an energy population swap follows immediately from (6) since only one parameter is needed to describe the population transfer in a two-level system. For example, the Einstein rate equation for an interaction of two levels with thermal radiation describes such an energy population swap interaction.
In appendix I, we show examples of two (partial) swap Hamiltonians. One generates a density matrix swap and the other an energy population swap.

The adiabatic evolution step
Stroke A and C are termed "adiabatic" as we require that the evolution will be diagonal in the energy basis at the beginning and the end of the stroke. The coherences will remain zero and the energy populations will remain as they are. Such a process does not have to be slow.
For example, consider the Hamiltonian that commutes with itself at all times such as: A density matrix that is diagonal in the energy basis |i will be invariant under this type of Hamiltonian regardless of how fast the energy levels are changing. In a two-level spin system this Hamiltonian will be H adiab = B z (t)σ z , where σ z is the z Pauli matrix. The population and coherences evolution can be obtained by more general time dependent Hamiltonians that do not commute at different times [H(t 1 ), H(t 2 = t 1 )] = 0. Nevertheless the possibility of an adiabatic transformation is guaranteed by the coherent control theorem [47,48]. Alternatively it is possible to use a method known as "quantum driving" or "shortcut to adiabaticity" [49,50,51,52] to generate an evolution that preserves the energy-basis diagonal form of the density matrix at the end of the process.
In principle, we allow the energy levels to cross ‡. However in such a case it important to remember that the energy index is a level index and not some order index that indicates how the levels are ordered. For example, in Fig 1b-

The average populations in steady state operation
The dynamics of the system (engine) reduced density ρ s and its statistical features are explored. We assume that during stroke B a single particle interact with the engine with probability R. The description can easily be generalized to include the possibility of more than one collision during stroke B or D. The probability R appears naturally in a collisions model as there is no certainty that a particle from the a bath particle will be available to interact with the engine. For simplicity we assume that the bath parameters are identical for both baths: they have same R and the same swap parameter x. A generalization to different x and R for different baths is straightforward using the same methods.
Even though for some applications transient behavior may be of interest, here we focus on the steady state operation. The average engine population at stage C, ρ C s , is related to that of stage A, ρ A s , via: The first term describes a no-collision event and the other describes a collision with a swap parameter x. Equation (14) simplifies to: Here x and R are inseparable since both of them have the same effect on the average population. In the same manner we can write equation for ρ A s : By combing (15) and (16) we get: Note that a combination of Gibbs states is not a thermal state if there are more than two levels. Hence, as expected in a finite time thermodynamic framework, the multilevel swap engine is never in a Gibbs state for xR < 1. The expectation value of the population change is: Inspection of (17) and (18) shows that even if the initial density matrix of the engine has some coherences they have no impact on the energy diagonal steady state. Thus for all energy observables computed locally (i.e., with the reduced density matrices of the baths or the engine) it is sufficient to consider the population vector p i = diag(ρ i ) where i = s , c , h . Hence, in this notation: In general, the population change under a swap operation in steady state is always proportional to p h − p c . This result will have a large impact later on. We point out that if the baths are only coupled to some of the system levels as in [4], then the result is different. When xR = 1 in a single collision, or xR < 1 with a vast number of collisions a complete thermalization of the engine's population takes place. Therefore, many of the results in the limit xR = 1, apply to a much more general setup than the swap model described above. They apply wherever the engine is connected long enough to effectively reach a Gibbs state (and not some other state).
For work and efficiency investigations, only the populations in the energy basis are important. Hence, it is not necessary that the whole density matrix is swapped as in (7)-(9), rather it is sufficient that an energy population swap takes place (see Eqs. (10)-(12)).

Quantum vs. classical swap and entropy generation
In Eq. (20) x and R are lumped together in a product form xR. Yet, they are not physically equivalent. R controls how deterministic the system is, and x determines the interaction strength and the "quantumness" (coherences and entanglement before the partial trace). Note that the quantum behavior is not monotonically increasing with x, since the system is classical when x = 1.
In the quantum case, where R = 1 x < 1, x can generate coherences and entanglement in the joint density matrix of two interacting particles. The entropy increase in the sum of entropies of the reduced density matrices is the result of ignoring entanglement and classical correlations.
In contrast, in the classical case where x = 1, and R < 1 , the particles are either fully swapped or left as they are. Therefore entanglement cannot be produced from the initial product state. Here, the sum of the entropies of the individual particles also increases since the information if a collision took place or not is discarded.
In both cases the total entropy increase of the reduced entropies is encapsulated in the mutual information of the colliding particles. In the quantum case, the mutual information contains contribution from entanglement. The separation to quantum and classical correlations can be studied using the quantum discord tool [53]. This however is outside the scope of this work. Furthermore in our model the observables we studied depend only on xR so they can equally be obtained from a classical or a quantum realization. However, we do no expect this to hold for observables that are not functions of the steady state population.

First law
In the present model the expectation value of the energy is: On average in a complete cycle U initial = U f inal we have: where the first term correspond to stroke A, the second to B, and so forth. Regrouping we identify: where we used the "positive heat in, positive work out" sign convention. Note that these quantities are invariant to any constant shift of the level in one or two of the bath:  Note that the energy at the end of the cycle is equal to the energy at the beginning of the cycle only on average. In a specific cycle it is not zero. Since the engine's Hamiltonian is time dependent its energy need not be conserved at each instant. Using (20) and (26) we obtain: As mentioned before, the limit xR = 1, (28) holds for any interaction that leads to a complete thermalization (or very close to it) and not just for a swap interaction. The modes of operation of a two-level swap machine are shown in Fig. 2. In general it is not straightforward to write analytically the condition for engine or refrigerator operation using the energy levels and the temperature due to the exponential dependence of the energy difference on the same parameters. Yet, in some temperature regimes this is considerably simpler as discussed in section 7.

Second law and the Clausius number
In cyclic processes in classical thermodynamics the second law can be expressed in terms of the Clausius inequality¸δ Q T ≥ 0. In our case we calculate: where, as before, the bracket denotes the average value in steady state. Strictly speaking, this is not the classical Clausius inequality. The swap collision process is not an isotherm at all. For a multilevel system a temperature cannot be assigned to the bath or engine particles since they are not in a Gibbs state after (or during) the collision. In this section we find that R ≥ 0 holds in steady state for any multilevel swap engine and that it has an information theoretic interpretation. Furthermore, it is shown that R ≥ 0 belongs to a family of more general inequalities. It is easy to verify that R ≥ 0 entails within the second law. For T c = T h we get W ≤ 0. That is, in steady state, no work can be extracted from a single bath. Alternatively, when calculation the efficiency, , R ≥ 0 ensures that the efficiency is smaller than the Carnot efficiency (see Sec. 6.2).
Using the expressions for heat we want to show that for swap heat machines: In a specific cycle this does not have to be true. Our aim is to show that it holds on average when the system is in steady state. To obtain the average we use (20) and get: which leads to the following result for swap engines: where J is the Jefferys divergence and D KL (p|q) = i p i ln p i q i is the Kullback-Leibler divergence or the relative entropy. A similar result is known for isothermal processes [54]. Since D KL (p|q) ≥ 0 for any two probability vectors (also known as Klein's inequality), R ≥ 0 follows naturally as information theory inequality.

Generalized Clausius inequality
In this section we show that the Clausius number inequality is a special case of a family of inequalities that hold for swap heat machines: Clearly, if this equality holds, it holds for any odd analytic and monotonically increasing function of Tc as well. The proof is straightforward. While R ≡ R 1 can be understood in thermodynamics terms of heat, temperature and entropy, for m > 1 higher energy powers are involved and there is no straightforward thermodynamic interpretation.
Denoting ε b,i = E b,i /T b and D i = ε h,i − ε c,i we get that inequality (35) is equivalent to: Next we write R 2m−1 = 1 2 R 2m−1 + 1 2 R 2m−1 , exchange the indexes names in the second term and obtain: The term D 2m−1 i −D 2m−1 j has the same sign as D i −D j and the same sign as (e −D j −e −D i ), hence the product of the two differences that appears in (37) is always positive. The rest of the multipliers are positive and symmetric under i ↔ j and therefore R 2m−1 ≥ 0.

Clausius dominated level
Immediate consequence of this generalized Clausius number inequality follows from considering m → ∞. In this case, the term with largest Tc becomes enormously larger than all the other terms. Therefore, this single term in R ∞ must be positive. If this term is positive, it is positive also in R 1 and therefore the Clausius dominated level i max defined by: satisfy: where we assume that has a single maximum and that dp A→C imax = 0. Equation (39) allows us to immediately determine the direction of heat flow into the baths for this level. Of course, this can be done explicitly by evaluating numerically the probabilities. This equality is not trivial since the Clausius dominated level is not necessarily the largest element in the Clausius number sum (30).

Local and non-local quantities
We call a quantity 'local scalar' if the scalar quantity can be written in terms of the single-bath parameters E b and T b (i.e. f (E b , T b )). The mean energy p b · E b of the bath particle, its purity p b · p b , its entropy, its energy variance and its free energy, are all examples local scalars. A local function is a function of local scalars.
Many other quantities of prime importance cannot be written in this form: the mutual coincidence p c · p h , the fidelity of the two baths i √ p c,i √ p h,i , the Jefferys and Kullback-Leibler divergence, the Wootters distance. Hence these scalars are 'non-local scalars'. More importantly heat, work and efficiency are non-local scalars. In our article thermodynamics amounts to finding relations between non-local scalar (e.g. work and efficiency) and local scalars (temperature entropy etc.).

Purity reduction
Let us define the purity of the bath after the collision by the purity of the reduced density formed by tracing out the engine and the other bath. Although the purity does not have all the appealing properties of the von Neumann entropy, it provides a simple and convenient measure of impurity. Furthermore, in contrast to entropy, it is not sensitive to whether we keep track of the particles position or not (i.e. there is no mixing entropy issue). In addition the purity will naturally emerge in the derivation of bound on maximal work and efficiency.
The purity change in one cycle in one bath is: where the absolute value of a vectors refers to the standard Using dp c = −dp h for steady state operation we finally obtain: In a refrigerator the cold bath the purity increases while the hot bath becomes more mixed. Yet, the change in the hot bath is larger so the purity of the whole system decreases. For xR < 1 the whole system (engine+baths) becomes increasingly mixed in all modes of operations.
In the spirit of section we want to express in the purity reduction using local functions. This can be achieved by applying the inverse triangular inequality: Interestingly, in contrast to the Clausius inequality, we did not assume a thermal distribution. In fact, in steady state this device will decrease the total purity for any bath population that is diagonal in the energy basis. While the Clausius is not well defined (there is no notion of temperature for non thermal baths) the purity decrease still holds. Finally, we note that using the Jensen inequality for the ln function it is straightforward to show that the purity is related to the entropy through: where S = − i p i log p i is the standard entropy function. This inequality can be interpreted in the following way: the Chebyshev sum inequality yields: P ≥ 1/N where the equality holds for uniform distributions. e S is the Shanon's effective number of degrees of freedom. The right hand side of (45) is also equal 1/N for uniform distribution. Hence (45) expresses the relation between two measures that count the degrees of freedom in the system.

Mutual coincidence necessary conditions for engines and refrigerators
The engine regime is defined by the condition W = xR Using the mean bath energy E b = p b · E b and the free energy defined we get: After rearranging and using the Jensen inequality to get ln p c · p h ≥ i p c,i ln p h,i , ln p c · p h ≥ i p h,i ln p c,i we get that the work is smaller than: where we used: If the left hand side is larger than zero so will W . Thus, demanding that the left hand side of (47) is positive we get a necessary (but not sufficient) condition for the heat machine to perform as an engine: where P ch is the mutual coincidence p c · p h . While the purity describes the probability that the same result will appear when measuring the energy of two particles in the same bath, P ch describes the probability that particles from different baths will be in the same level. P ch is often used in cryptography and code ciphering. Note, that it does not imply a correlation but simply the likelihood of identical events in both baths. Note that (48) actually imposes a restriction also on the individual purities since: For the refrigerator the work is not a good indicator since it is possible to apply work to the system without cooling the cold bath. Cooling occurs when the cold bath gives away heat to the system: Repeating the same procedure as before we get a necessary refrigerator condition: In contrast to other measures like fidelity, P ch has a very simple statistical interpretation which also make it very easy to measure in practice. i p c,i p h,i is the probability that two particles chosen from different baths will be in the same state (like getting the same result with two different unbalanced dices). To evaluate it, there is no need to keep track of the exact result and then to estimate the probabilities p c,i and p h,i through their frequencies. One only needs to keep record if the results are the same or not. Thus, P ch corresponds to a binomial random variable. Like the purity P ch also satisfies P ch ≥ 1/N (if the levels don't cross) but this is true for any two monotonic distributions and does not give any indication of the heat machine's functionality.
6. Upper bounds on work and efficiency

Bounds on the maximal work production
The first work upper bound can be obtained from (47) where the inequality ln p c · p h ≤ 1 2 ln P c P h is used to get to bring it to the local form: A different bound can be obtained by writing the work in a different form: The last two terms are always negative (including the minus sign). Therefore the first term must be positive for the machine to perform as an engine. S h ≥ S c is expected for engines as engine transfer entropy from the hot bath to the cold bath but (53) quantifies how much large the entropy difference must be: The work expression (53) contain contain the non-local quantity D KL . To obtain a non-local upper bound we use the inequality: and get: Different approaches and different inequalities may lead to work upper bound that may perform better in certain regimes. In addition, further reasonable restrictions on the system may lead to better results. In appendix III we derive another work bound under some assumptions about the energy levels structure.

Upper bounds on the efficiency
In the two level case the efficiency is simply given by 1 − and the population change cancels out. In the multilevel engine the expression is more complicated: Using the Clausius number defined in Sec. : This result is still exact. To get a simpler bound we use again the inequality KL(p c |p h ) ≥ 1 2 |p c − p h | 2 and Cauchy-Schwarz in the numerator to obtain: where we used E h to denote the centered energy vector: The replacement E h → E h is highly important. It keeps the bound invariant to energy shifts and at the same time makes the bound tighter.
In order to obtain a local-quantities bound we use the same procedure we use P ch ≤ √ P c √ P h and obtain a weaker yet local upper bound : This inequality is often very close to the Carnot efficiency. Yet, it still shows that the baths purity imposes some limitations on the efficiency. A different upper bound can be written down in term of the Wootters statistical distance [55] between the hot and cold probability distribution: From [54] it follows that: The efficiency bound obtained from using (64) is not always smaller or larger than (59) and it uses the non-local quantity N i=1 √ p c,i √ p h,i (Fidelity). Yet, it introduces another relation between the difference in the statistics of the baths and the efficiency.

The Ultra hot baths regime
In atomic physics often the cold bath temperature is more difficult to produce. In contrast it is fairly simple to produce a heat source so that, roughly speaking, E h /T h , E c /T c 1 (a more rigorous condition will be given later) . In this section we study the multilevel heat engine operation when the hot bath is so hot that we consider only the first order correction in 1/T b to the T b → ∞ limit where the p b,i → 1/N . The ultra-hot regime is defined by the small parameters: This time the centered energy form E b emerges naturally from the free energy normalization factor in the ultra-hot limit. W > 0 yields the necessary and sufficient engine condition: Applying the Cauchy-Schwartz we get a necessary condition for ultra-hot swap engines: To optimize the work output, let us first assume that the norms of the energies |E c | and |E h | are fixed. Under these two constraints the last two terms in (66) are fixed. To maximize the first term, and consequently the work, E c must be to be parallel to E h : where C is the compression ratio and (70) follows from (68). For "uniform compression" (69) the device works as a refrigerator when C < 1, and as a heater when C ≥ T h /T c . In the regime (70) it performs as an engine. The uniform compression can be studied for any temperature but in the ultra hot regime it is found that the uniform compression maximizes the work when the energy norms are fixed. Now we remove the restriction that |E h | is fixed and instead use (69) in (66). Imposing ∂ C W = 0 we get: where |E h | in the subscript signifies the |E h | = const optimization constraint. The efficiency is: Where η c is the Carnot efficiency. We note that this efficiency is always tighter (smaller) than the Novikov, Curzon and Ahlborn (NCA) efficiency bound η CA = 1 − T c /T h [56,44]. The NCA derivation assumption assumes the working substance is in some well-defined temperature. As Eqs. (17) and (18)  and η ultra hot max,|Ec| = T h −Tc T h +Tc = ηc 2−ηc . Interestingly, this efficiency is always larger than the NCA result. The same expression for the efficiency was obtained in [57] for a different type of engine in the ultra hot regime.
The expression for maximal work for a given |E h | is: The expression looks asymmetric in T c , T h and E c , E h since we already applied the optimal ratio between E c and E h . Bringing back E c to the equation will restore symmetry but will hide the choice already made for E c by the optimization procedure. At first sight, it seems that the work becomes smaller when the number of levels increases. Yet, E 2 h may also depend on N . For example if the levels are degenerate so that there are N/2 replicas of two levels, then we get that the work does not depend on the number of replicas, as expected.
It is interesting to look on the Clausius number in this limit: The constant term can be dropped out since it contributes zero to the total sum. After the constant term is removed, all the terms in the Clausius sum are positive. This is not true for below the ultra hot baths regime. In this regime, the Clausius inequality holds since the probability difference and the Clausius factor Tc have the exact same form.

Almost "Quasi-static evolution" at finite time
Consider the case where the swap parameter is small enough so that after a collision the change in a bath particle population, dp b,i , is small with respect to the original bath population so that δp b,i p b,i . At the end of appendix I it is shown that the energy diagonal form of the density matrices is conserved in an energy population swap (this trivially holds in a density swap interaction). Consequently, it is possible to use just p b and δp b to calculate the von Neumann entropy change. To first order in δp b,i the change in the bath particle's entropy is given by: Note that this equation holds only for the bath particles and for tiny deviations from the Gibbs state. The engine is not in a thermal state at all, and cannot be assigned a temperature. However, in a complete cycle, on average, the engine returns to its initial state so the engine does not contribute to the average total entropy production of the total system. Therefore the total entropy increase in both baths and the system satisfies: by virtue of the Clausius inequality for swap heat machines proven earlier.
Assuming that on average there are n cold collisions and n hot collision, the work can be expressed in terms of the entropy changes: Although (75) seems like a plausible property for a bath, it is not mandatory and in our model it emerges only in the weak coupling limit (small population change per collision).
Even though we assumed the change in a single collision is small it does not mean that the heat exchange in strokes B or C must be small in our model. If we allow multiple collisions, it is possible to have large changes while still satisfying In our entropy considerations we did not take into account mixing entropy so it is not claimed that (75) is the change in the entropy of the system but just the part of entropy change that is responsible for the heat exchange.

Conclusion
We have presented an analysis of a multilevel heat machine that is driven by a partial swap interaction sequentially with a hot and cold baths. Our approach emphasizes the back reaction of the engine on the baths. By studying the bath's thermodynamics observables we were able to relate the work, efficiency and entropy production to a number of measurable bath properties like purity and mutual coincidence. The equivalent of Clausius inequality in this system was generalized to higher order energy moments. We identified a quasi-static regime in finite time evolution and an ultra-hot regime where stronger statements can be made.
While our findings are valid for any multilevel engine and collision that can bring the system to a Gibbs state, it is interesting to consider other bath models like the ones used in continuous engines where different levels interact with different baths or where the thermalization rate is different for different levels.

Appendix 1 -The two-level swap Hamiltonian
The partial swap unitary is given by: where σ i are the Pauli matrices and φ is the swap angle for φ = π/2 a complete swap takes place. Let ρ 1,2 be diagonal matrices. The reduced density matrix of particle 1after applying U is: One can verify that the density swap rule is satisfied: if ρ 1,2 are diagonal then ρ 1,2 are also diagonal.
Since the Hamiltonian is written in terms of the Pauli matrices it is not straightforward to generalize it to a multilevel swap operation.
For energy population swap, the generalization is simpler when the two-particle states are considered. A complete swap operation yields: U cs |i, j = U cs |j, i . The states |i, i are invariant under this operation. The complete swap unitary in the twoparticle state can be written as: 11  22  33  12  21  13  31 . . .
The column to the left shows how the states are ordered in the two-particle density matrix. In this form the unitary has a clear block diagonal structure. In space of each non identical state {|ij , |ji } a simple spin flip takes place: Replacing the spin flip by a more general rotation in this Hilbert subspace, it is easy to deduce the partial swap time-independent Hamiltonian: The φ ii terms just contribute a phase to the invariant states. A complete swap takes place when φ i,j =i = ±π/2. In general, there is no particular reason why the rotation rate, φ ij should be the same for all pairs of states. However, when it is the same, φ i =j = φ can show that the energy population swap probability swap rule follows: p 1 = diag(tr 2 (U ρ 1 ⊗ ρ 2 U † )) = cos 2 φp 1 + sin 2 φp 2 This swap model is specially designed for the energy basis. Consequently, for partial swap, in general, it does not satisfy the "density swap" rule but the energy basis probability rule. Yet, a proper choice of φ ii leads to the density swap Hamiltonian that appears in the exponent of (78). We conclude by noting the important fact that the energy population Hamiltonian preserves the energy diagonal form after the partial trace. That is, if the input state is ρ 1 ⊗ ρ 2 where ρ 1 , ρ 2 are diagonal in the energy basis then ρ 1(2) = tr 2(1) (U ρ 1 ⊗ ρ 2 U † ) will be diagonal in the energy basis. For simplicity we illustrate it for the φ i =j = φ case. The density matrix after the collision and before the trace is: where f ij are some complex coefficients. The third term vanishes when taking partial trace on either particle. Due to the block structure (80) this holds even when φ i =j depends on i and j.

Appendix II -Markovian swap formulation
As shown in Sec. 3 in steady state all coherences of the density matrix vanish thus only the diagonal elements are important and the Markov chain formalism can be used to describe the evolution of the diagonal elements (probabilities). In this appendix it will be more convenient to use Dirac's "Bra-Ket" notation for the probabilities vectors in order to distinguish between right vectors and left vectors. However the normalization is still the regular probability normalization. The probabilities before and after the collision with the hot bath satisfy: wherex is an effective swap parameter that may contain "classical" contribution from the collision probability R and quantum contribution from a quantum partial swap (x < 1 or φ < π/2). I N ×N is the identity operator. One can verify that the above equations lead to the population swap rule p s = (1 −x)p s +xp h . A Markov chain with a single steady state vector is always associated with a left eigenvector of the form 1, 1, 1..|. Thus, |p T h is a right eigenvector of K h and it has a unity eigenvalue: Note that: P c + P h − 2 N ≤ √ 2 N −1 N . If all the levels are compressed (not necessarily by the same factor) so that: In this case, we can use (a − b) 2 ≤ |a − b| |a + b| = |a 2 − b 2 | ∀ ab > 0 to get: Finally we can write down a bound based only on individual bath properties: We cannot use here, as it will typically violate (95). Furthermore, although (97) seems to imply that W → 0 if |E h | = |E c |, the norms cannot be equal and satisfy C i ≥ 1 (with the exception of the uninteresting case where all the compression factors C i are equal to one in which the work is trivially equal to zero). The advantage of bound (97) is that it uses only local scalars.