Imperfect Thermalizations Allow for Optimal Thermodynamic Processes

Optimal (reversible) processes in thermodynamics can be modelled as step-by-step processes, where the system is successively thermalized with respect to different Hamiltonians by an external thermal bath. However, in practice interactions between system and thermal bath will take finite time, and precise control of their interaction is usually out of reach. Motivated by this observation, we consider finite-time and uncontrolled operations between system and bath, which result in thermalizations that are only partial in each step. Our main result is that optimal processes can still be achieved for any non-trivial partial thermalizations, at the price of increasing the number of operations. We focus on work extraction protocols and show our results in two different frameworks: A collision model and a model where the Hamiltonian of the working system is controlled over time and the system can be brought into contact with a heat bath. Our results show that optimal processes are robust to noise and imperfections in small quantum systems, and can be achieved by a large set of interactions between system and bath.

Optimal (reversible) processes in thermodynamics can be modelled as step-by-step processes, where the system is successively thermalized with respect to different Hamiltonians by an external thermal bath. However, in practice interactions between system and thermal bath will take finite time, and precise control of their interaction is usually out of reach. Motivated by this observation, we consider finite-time and uncontrolled operations between system and bath, which result in thermalizations that are only partial in each step. Our main result is that optimal processes can still be achieved for any non-trivial partial thermalizations, at the price of increasing the number of operations. We focus on work extraction protocols and show our results in two different frameworks: A collision model and a model where the Hamiltonian of the working system is controlled over time and the system can be brought into contact with a heat bath. Our results show that optimal processes are robust to noise and imperfections in small quantum systems, and can be achieved by a large set of interactions between system and bath.

I. INTRODUCTION
Recent years have experienced a renewed interest in understanding thermodynamics in the quantum regime [1,2]. A common assumption in the field is that a system becomes a Gibbs state when put in contact with a thermal bath. While this being perfectly reasonable, it requires some implicit assumptions: For example, that a sufficiently long time is available to ensure full thermalization of the system [3]; or, within approaches based on repeated interactions between system-bath, that control is available on the system-bath interaction [4]. The goal of this article is to study how optimal thermodynamic processes are modified when one does not have access to such "perfect thermalizations". Our results thereby contribute to a recent research line which aims at understanding the level of control required in thermodynamic protocols [5][6][7][8][9][10][11]. To ensure that our findings are not particular to a specific model for thermodynamic transformations, we consider two different frameworks in this article.
We start by considering a framework where thermalization takes place due to repeated interactions between the system (S) and elements of a thermal bath (B), in the spirit of collision models [4,12]. We study work extraction from a qubit system S when having access to B, which consists of a collection of qubit systems in a thermal state. Following [13][14][15], the optimal protocol for work extraction can be constructed as a concatenation of unitary swap operations between S and qubits of B. In practice, there are many possible sources of errors for this protocol; for example, the experimentalist might not be able to control precisely the time of the interaction between S and B, resulting in a partial thermalization, where α ∈ [0, 1], ρ S is the state of S, τ B = e −βH B / Tr(e −βH B ) a thermal state, and H B is the internal Hamiltonian of B. For such partial thermalizations, we show that maximal work extraction can still be achieved for any α = 1, at the price of requiring more operations for the same amount of extracted work. Next, we prove that this result does not only hold for qubit systems, but in fact can be generalized to arbitrary systems. Furthermore, we also study work fluctuations, showing that for a large number of steps the final work distribution using such uncontrolled interactions converges to the same form as for the optimal protocol [14,16]. Optimal thermodynamic protocols can therefore be achieved through protocols that only use partial thermalizations in the form (1), which arise due to uncontrolled interactions between SB. In a second part of the article, we study the same question in a different and widely used framework, namely when B is a macroscopic bath and work is extracted by changing the Hamiltonian of S over time. In this case, partial thermalizations naturally appear due to finitetime contacts with B. However, we significantly generalize from the partial thermalizations presented above, by describing the action of B on S as essentially any nontrivial quantum channel with the Gibbs state as a unique fixed point. We then show that isothermal processes arise in the limit of arbitrarily slow transformations, independently of the thermal bath in consideration. This confirms the intuition that any slow process where the Hamiltonian of S is modified while being in contact with B tends to an isothermal process in the limit of arbitrarily slow transformations. This result also allows us to link discrete [13-15, 17, 18] and continuous [19,20] optimal thermodynamic processes.
The paper is structured as follows. In Section II A, after introducing the framework, we describe the ideal protocol for work extraction and later introduce a protocol which also works for imperfect operations. In Section II B, we show that the noisy protocol, like the ideal protocol, allows for maximal work extraction in the case where system and bath consist of qubits, and calculate the corresponding work fluctuations. These results are later generalized in Section II B 3 for generic qudit systems. In section III, we discuss the second framework and conclude in Section IV.

II. COLLISION MODEL
Collision models have been extensively used to describe open quantum systems [4], in particular thermalization processes [12,[21][22][23][24] and non-equilibrium scenarios in different thermodynamic contexts [25][26][27][28][29][30]. The main assumption is that when the system interacts with a bath, it collides and thereby interacts with only one particle of the bath at a time. This means that we can model a thermalization process as successive interactions of the system with different parts of the bath. Recently isothermal processes have also been described in the framework of collision models as a concatenation of interactions between the system S and elements of the bath B [13][14][15]. A difference between the latter and the standard description of thermalization via collision models [12] is that in the latter the system S interacts with non-identical Gibbs states, in order to construct an isothermal process. We follow this approach in the first part of the paper. We do so in an inclusive way,where all elements entering in the process are described explicitly as quantum systems, in the spirit of the resource theory of thermodynamics [5,[31][32][33].

A. Framework
The model we consider for describing optimal processes consists of a system S, which interacts with a thermal bath B at fixed temperature T , and a work battery W, in which the extracted work can be stored. Most of our results concern the scenario in which S is a twolevel system, a qubit -later we will consider generic qudit systems. The Hamiltonian of S is hence taken to be H S = S |1 1| S . Similarly, the bath consists of a collection of N + 1 qubits, with internal Hamiltonians H (k) For W, we consider a weight that stores the work as potential energy and take the Hamiltonian H W = mgx, wherex is the continuous position operator,x = dx x |x x|, and for simplicity we choose mg = 1 J/m. The total Hamiltonian on SBW (which is a shorthand for S ⊗ B ⊗ W ) hence reads The initial state is taken to be a product state, We take for S a diagonal state, the case of coherent states is discussed in [16,34], where it is shown that the coherent content of the free energy can be extracted by using an external source of coherence, a process where B does not intervene. Hence, without loss of generality, here we focus on diagonal states. The states of the bath τ (k) B are Gibbs states at temperature Since we are dealing with qubits, we can also write where we note the relation In order to construct an optimal process, the excitation probabilities {q k } N k=0 are chosen to start at q 0 = p 0 and increase/decrease strictly monotic to q N = e −β S /(1 + e −β S ). Hence, τ (0) B = ρ S and in the case of maximal work extraction τ Finally, the initial state for W is taken to be ρ (0) W = |0 0| W . As we now argue, this choice is in fact irrelevant for the extracted work.
Thermodynamic operations consist of unitary operations U on SBW. The average extracted work is then defined as the energy gained by W, where ρ W = Tr SB [ρ] is the initial state, and ρ W = Tr SB [U ρU † ] is the state after the transformation. Importantly, we impose two conditions on U : 1. The total energy must be preserved by U , more concretely [5], This condition ensures that all energy fluxes are accounted for, as there cannot be energy flowing out of SBW, and also justifies definition (8).
2. U must commute with the displacement operator on W . If ρ W is diagonal, this implies that the extracted work W is independent of the initial state ρ W of W. Since we will only consider diagonal states for ρ W , this condition ensures that the state of the weight cannot be used as a source of free energy, so that the extracted work comes solely from SB.
Thermal Bath

System
Work Battery Figure 1. The setting described in this paper consists of three subsystems: The system S initially in a non-thermal state, a work battery W in which the work is stored and a thermal bath B that provides thermal energy to be converted to work. They interact sequentially by means of U (k) SBW such that in each step another part of the thermal bath is involved. The total extracted work can then be computed as the sum of the extracted work in each step, W (k) .

Maximal Work Extraction
The maximum extractable work of ρ (0) S is given by the change of (non-equilibrium) free energy where F (ρ, H) = Tr(ρH) − T S(ρ), with S(ρ) the von Neumann entropy, S(ρ) = − Tr ρ ln ρ [35,36]. Following Refs. [13,14], we now briefly discuss an explicit protocol leading to (10). The total process is described by the following unitary operation on SBW, where each U (k) SBW acts on S, W and the k th bath qubit, see Figure 1. Using the short hand notation, |a, b, c := |a S ⊗ |b B ⊗ |c W and ω k : Let us also write this expression in matrix form using the basis The action of U (k) SBW is illustrated in Figure 2. Note that SBW swaps the energy degenerate states |01x SBW and |10x + ω k SBW for all x, where only the k th bath qubit is involved. The states |00x SBW and |11x SBW are mapped to themselves.
• U (k) SBW acts upon energy degenerate states, and hence condition 1., as given by (9), is satisfied. It is straightforward to see that condition 2. is satisfied as well.
• At the global level, the final state after applying U (k) SBW becomes correlated. At the level of the reduced states, it leads to a swap between the states ρ S ↔ τ (k) B , • As k increases, the state of S progressively changes from ρ (0) S to τ S , which mimics an isothermal process.
Furthermore, the total extracted work W = Tr H W (U ρ (0) U † − ρ (0) ) from (11) satisfies, hence becoming maximal in the limit N → ∞ [13,14]. As we shall see later, this expression is only valid under the assumption that ρ S is initially not in a pure state.

Imperfect Operations
Let us now describe the set of imperfect operations we consider. For clarity, we start by presenting a simple model to describe lack of control on the operations U (k) SBW . For that, note that the unitary (12) can be implemented by turning on the interaction gσ x in the subspace |0, 1, x , |1, 0, x + ω x for a time 2π/ g. If the timing is not precise, the corresponding unitaryŨ (k) SBW would take the form, where α k ∈ [0, 1] quantifies the error in the k th step: for α k = 0 we apply the desired full rotation and for α k = 1 we do nothing on the state. After (15), the reduced states of S and B, given by ρ Hence, here the lack of control essentially leads to a partial thermalization process.
In fact, there are more general processes that lead to (16). In Appendix A, we consider arbitrary thermal operations [5] on the relevant degenerate subspaces |0, 1, x , |1, 0, x + ω x . By this, we aim to model the situation where the relevant subspace where the experimentalist is acting upon is not thermally isolated, e.g., there is decoherence in the unitary operations being implemented. We then show that all operations of this form lead to the same model of noise (16), with the α k 's being complex functions of the applied specific thermal operation. Therefore, we conclude that (16) is generally valid when we restrict ourselves to the operations described above which map states from the energy degenerate subspaces |0, 1, x , |1, 0, x + ω x to states on that same subspace while acting arbitrarily on an auxiliary system. This is the model we will consider for the rest of our work. As a remark, note that (16) naturally allows for generalizations when considering systems beyond qubits.
Before moving further, it is important to keep in mind that the set of imperfect operations we described is of course not the most general noisy, or thermal, operation. The motivation behind our choice is that, as we will now show, thermodynamic processes can still be constructed despite the lack of control. Conversely, if we lift the restrictions made above, it is not difficult to convince oneself that other types of noise, e.g., thermal noise acting on S or W only, are detrimental for thermodynamic purposes and cannot be compensated for. This is discussed in Appendix B.

B. Results
This section contains the main results for the collision model. In Section II B 1, we present several extensions of (14): (i) we provide a bound on the N -dependence of the error term as O(1/N ), (ii) we show that when ρ S is initially a pure state, the scaling of the error has instead the form O(ln N/N ); and, most importantly, (iii) we show that (14) can be obtained by any thermal operation acting non-trivially in the degenerate subspace |0, 1, x , |1, 0, x + ω k . In other words, very generic operations lead to (14). In Section II B 2, we characterize the fluctuations of these optimal processes. Finally, in Section II B 3, we generalize our results to generic qudit systems by considering partial thermalizations.
We now proceed to derive our main results, taking (16) as a starting point. We recall that the α k 's characterize the noise level of the operation in the k th step. One may think of the α k as being equal to a fixed value α. For purpose of generalization, we treat the α k 's as independent random variables with average value α. More precisely, let the average of a function f over {α j } N j=1 be where P i (α i ) are some unknown probability distributions that satisfy dα i P i (α i )α i = α for all i. Since the α j 's are independent, in particular It is then straightforward to obtain from (16) that where we used that τ where the corresponding results should be understood as the average over uncontrolled parameters {α j }. Now, starting from the recursive definition of ρ (k) S in (20), we can determine its explicit form as Furthermore, since all τ S , and hence the same relation holds for the excitation probabilities p k , with ρ (k)

Average Work Extraction
We first show that it is possible to obtain W → ∆F for any α < 1 in the limit of N → ∞. For clarity of the presentation, we find it convenient to present here the case where ρ S = |0 0| S and S = 0, whereas the general case for arbitrary initial states of S is treated in Appendix C. Note that the case ρ S = |0 0| S and S = 0 corresponds to well-known information to energy conversion, where ∆F = k B T ln 2 work can be extracted from one bit of information. The opposite direction, resetting one bit of information while investing work ∆F = k B T ln 2, the socalled erasure, works analogously since we are looking at reversible processes.
To compute the average energy gain of the weight W, we first notice that the work extracted in step k is which follows from energy conservation. Thus, the total average work extracted after N steps is the sum Let us now use our choice of initial conditions S = 0 and p 0 = 0, i.e., initially ρ S = |0 0| S . The distribution of bath qubits is taken to be so that τ = τ S . Note that this specific distribution is chosen here for simplicity and to give an explicit upper bound; in Appendix C it is proven that any set of excitation probabilities {q k } 1...N with q k = f k N and f monotonically increasing/decreasing and of class C 2 yields (14). Thus, there is some freedom in the choice of Hamiltonians, as long as they are close enough. Using (22), we get an explicit form for the excitation probabilities of ρ where we can already see that they only differ by O (1/N ) compared to the excitation probabilities q k we would get using error-free unitaries. That is, for large N our state ρ S is at every step close to the corresponding thermal state, and hence it follows approximately the isothermal trajectory of the state transformation. We now show that for large N this also happens at an optimal work cost, i.e., that the work extracted becomes arbitrary close to the one in the error-free case. Using (26), the average extracted work simplifies to where W max := N k=1 E k ∆q k is the maximum average work that can be extracted with N steps in the error-free protocol, i.e., with α = 0. As we show in (C7)-(C10) of Appendix C, the loss in the optimal protocol due to the finite number of steps, goes as γ = O (ln N/N ), so for large N the W max converges to the difference in free energy. The loss due to noise ε gives the difference between the work extracted in the optimal protocol (α = 0) and the one with faulty unitaries (α = 0), and it can be upper bounded by This provides the desired result: In the limit N → ∞, ε → 0 and hence W → ∆F for any α < 1. Note also that the value of α sets the speed of the decay of ε. Roughly speaking, one needs approximately N α/(1 − α) operations to compensate for the presence of non-zero α.
In Figure 3 we plot the numerical value of loss ε and the upper bound (29). We observe that the deviation is small and the bound almost tight. We also plot ε together with the total loss ε + γ for finitely many steps as a function of α.
We now consider the more general case where ρ S = (1 − p 0 )|0 0| + p 0 |1 1| and general S . In Appendix C we show that W = ∆F − O(E max /N ) for any α < 1 and with E max = k B T ln((1 − p 1 )/p 1 ), i.e., for p 1 = O (1/N ) the error scales as O(ln N/N ) and for p 1 = O (1) as O(1/N ). As the protocols can be reversed, this scaling difference implies that pure states are harder to reach than mixed states, as expected. Importantly, this result is derived without assuming a particular form for q k , holding for arbitrary bath qubits which have a Gibbs distribution that satisfies q 0 = p 0 , strictly increases/decreases monotonically to q N = e −β S /(1 + e −β S ) and is drawn from a function f (k/N ) which is of class C 2 , i.e., twice differentiable with a continuous second derivative. Putting everything together, we thus conclude that a big freedom exists in both the evolution of SB and the distribution of B when constructing an isothermal process, which is used here to extract maximal work. The loss ε due to noise as a function of the total number of steps N for α = 1/2, p k = k/2N , kBT ln 2 = 1. The green line corresponds to the exact loss ε computed numerically with N steps, whereas the orange line corresponds to the upper bound (29). Inset figure: Errors as a function of α for N = 1000, p k = k/2N , kBT ln 2 = 1. The green line corresponds to ε, while the purple one is the total loss γ + ε (see (28)).

Fluctuations
The fluctuations of work depend on the probability distribution That is, the work fluctuations are mapped to energy fluctuations of the weight, which initially starts in a welldefined energy 1 .
In Appendix D, a lengthy calculation shows that for ρ S = |0 0| S , S = 0, α k ≤ α < 1 ∀k and distribution (25), the variance of P (W ) decreases as O (ln(N )/N ), and hence disappears in the limit N → ∞. This result implies that an amount k B T ln 2 of work can be extracted with certainty from one bit of information, even when dealing with partial thermalizations as in (20). The work probability distribution of this process is illustrated in Figure 4.
Given this result, the case of ρ S = (1 − p 0 )|0 0| S + p 0 |1 1| S and S = 0 can be easily understood qualitatively. First of all, recall that for general ρ S the optimal protocol is such that H is the probability distribution starting the protocol with |0 0| S and |1 1| S , respectively. From our previous result, it follows that P 0 (w) tends to a peaked distribution, which is now centered at (recall that we consider the protocol for maximal work extraction from ρ S ). Similarly, P 1 (w) tends also to a peaked distribution centered at Putting everything together, and noting that w 0 = k B T ln((1 − p 0 )/(1 − p (eq) S )) and w 1 = k B T ln(p 0 /p (eq) S ), we can write the work probability distribution in the limit N → ∞ as where δ(w − x) is the Delta function. Crucially, this distribution is independent of α, and hence we recover previous results where perfect thermalization operations (α = 0) were assumed [14,18].

Generalizations and Discussion
Although our results have been derived for qubit systems, they can be naturally extended to qudits, which is carried out in Appendix E. In order to get an intuition, consider the following distribution of bath elements, which can be seen as an extension of (25) when dealing with qudits. We then obtain from (21) Therefore, we see that the state ρ (k) S is thermal up to an error O(1/N ). Hence, when computing the extracted work, one expects W → ∆F with an error O(1/N ). We confirm this intuition in Appendix E, and in fact we show that W → ∆F for any α < 1 and for any smooth choice of Hamiltonians H = τ S . We note that this provides a big freedom in the choice of H (k) B , as any continuous function that interpolates between the initial and final points is expected to satisfy this condition. Summarizing, as long as the Hamiltonians of the bath are smoothly distributed, and the model of partial thermalizations (20) is satisfied, optimal processes can be achieved in the limit of infinitely many thermalization steps. It is worth pointing out that while the proof in Appendix E shows that the error scales as O(1/N ) (for non-pure initial states of S), it does not provide an explicit bound on the scaling, unlike previous results such as (29). That is, the price to pay for the generality of the result is the lack of an easily computable bound for the loss due to noise.
Another point worth stressing is that, even when dealing with qubits, there exist uncontrolled interactions that are thermodynamically detrimental and cannot be compensated for. For example, if B and W exchange energy, so that W tends to a thermal state, this clearly will render the extraction of useful work impossible. Similarly, letting S simply thermalize with B without raising the energy in W can only decrease the extractable work. In order to build an optimal thermodynamic process, it is needed that B continuously thermalizes S (and not W) while W is extracting work from it. A more detailed description of interactions that can not be tolerated can be found in Appendix B.
Let us also briefly comment on our assumption that the α k in (16) are random variables with identical expectation value α, which is independent of N . This assumption is in fact not crucial for the main results -although it is very convenient for the computations-, the reason being that our result holds for any α < 1. Note that we can always decrease the individual values α k without compromising the work extraction protocol. As a consequence, even if the expectation value of the α k would change with k, one can guarantee convergence to an optimal thermodynamic process as long as α k ≤ α < 1, ∀k. Furthermore, even if α = 1 − C/N for some constant C > 0, one still obtains the optimal expected work in the limit N → ∞ (also see section III). Similarly, we strongly expect that our results also hold for typical realizations of the α k , provided that each α k is chosen independently with expectation values bounded away from 1. Indeed, we expect that our results still hold even if a diverging number of α k is equal to 1 provided that these α k are a vanishing fraction of all the α k in the limit N → ∞.
Finally, it is important to realize that, although we have focused our attention on optimal work extraction, our results also allow for optimal state transformations in the other direction, i.e., they can be used for converting a thermal state to a given final state at optimal work cost. More generally, they apply to the interconversion of any two states at minimal work cost. Furthermore, we note that although our results have been derived within the resource theory of thermodynamics, it is not difficult to obtain essentially the same results in other frameworks, such as [15,17,18,38,39]. This is illustrated in the next section.

III. TIME-DEPENDENT HAMILTONIANS AND THERMALIZATION MAPS
The collision model presented above has the advantage of being very explicit, allowing us to derive precise results both about average work extraction and fluctuations. However, it has the drawback of requiring a relatively fine-tuned bath. We therefore now show similar results in a different framework, in which work is extracted by changing the Hamiltonian of S over time and contacts with the bath are described by a thermalizing channel, which can essentially be any channel that brings the system closer to the thermal state. The purpose of this section is two-fold (i) this model is experimentally more relevant and realistic, as no control over B is assumed (in the spirit of open quantum systems), and (ii) in this model we are able to generalize our results to arbitrary interactions with the bath that bring the system closer to the thermal state, with partial thermalizations as in the previous section being a specific instance.

A. Framework
In this framework, a protocol consists of a choice of trajectory of Hamiltonians H(t) and a choice of N times {t i } N i=1 at which the system is brought in contact with B (the continuum limit, where (t i − t i+1 ) → 0 so that S is essentially always in contact with B, will be discussed below). For simplicity, in the following we only consider cyclic protocols, for which the initial Hamiltonian coincides with the final Hamiltonian and denote by t 0 the time at which the protocol starts. We further assume that while the system is in contact with the bath, the Hamiltonian is not changed and the system thus (partially) thermalizes to the Hamiltonian H (i) = H(t i ). The contact to the bath at time t i is effectively described by a quantum channel G (i) and results in a state where ρ (i) is the state of S right before the thermal contact. Inbetween the contacts to the bath, the system evolves according to the time-dependent Hamiltonian H(t) from time t i to time t i+1 , represented by the unitary U ti+1,ti . The states ρ (i) and σ (i−1) are therefore related by the time-dependent evolution under H(t) as This time-evolution is associated with a work-gain with σ (0) = ρ (0) being the initial state. Given an initial state ρ (0) and trajectory of Hamiltonians H(t), we denote the total work in an N -step protocol by Let us finally discuss the assumption that we make on the thermalizing maps G (i) . Clearly, the thermalizing maps should have the thermal state τ (i) of the Hamiltonian H (i) as a fixed-point. It then follows in general from the monotonicity of the trace-distance under quantum channels that In the following we want to model that the thermalizing maps have some non-trivial effect on any non-thermal initial state. We will therefore first assume that there exists a fixed number κ < 1 such that This expresses formally the physical requirement that the interaction with the bath brings the system closer to the thermal state by at least some, even arbitrarily small, amount. In the following we call κ the thermalizationparameter of the protocol. For example, in the case of the partial thermalizations considered in the previous sections, we simply have κ = α i . Later we also discuss what happens when κ → 1 as the number of steps of the protocol is increased.

B. Results
Having introduced the framework, we can now easily state our result. Let us mention before, though, that the maximum amount of work that can be extracted using a cyclic protocol of the above form is bounded as [6,18] Furthermore, it has been known that for perfectly thermalizing maps (κ = 0), this bound can be achieved to arbitrary accuracy [17,18,39]. The optimal protocol consists of a first rapid change of the Hamiltonian to a HamiltonianH (0) such that (to arbitrary accuracy) followed by an isothermal process back to the initial Hamiltonian, meaning a process where t i − t i−1 = 1/N and N → ∞, so that the system stays in the thermal state of the current Hamiltonian throughout the protocol. The first quench yields a work-gain The optimality of the above protocol follows from the fact that an optimal isothermal process between two Hamiltonians H(0) and H(1) with thermal states τ (0) and τ (1) has work-gain Therefore, the total work-gain in the optimal protocol is given by (using H(0) =H (0) , H(1) = H (0) ):

Main result
Our following result shows that the protocol discussed above remains optimal for any thermalizing parameter κ < 1, since isothermal processes remain optimal in the limit N → ∞: In the above formulation, consider a protocol with a smooth curve of Hamiltonians H(t), thermalization parameter κ < 1 and initial state ρ (0) . Choose the time-steps as t i = i/N . Then the corresponding workgain is given by independent of the initial state ρ (0) . The independence of the initial state follows from the fact that the state converges exponentially to a thermal state within the first few steps of the protocol. The proof of (45) is given in Appendix F. There, we also discuss what happens if the thermalization parameter may depend on the number of steps N as for some characteristic parameters τ, γ > 0. We find that for γ ≤ 1, Thus, even if κ approaches unity with γ < 1 the error vanishes. In the case γ = 1 we obtain a constant error of the order of 1/τ . If γ > 1, on the other hand, we cannot bound the error terms in the limit N → ∞. This is reasonable since the case γ > 1 can be seen to correspond to scenarios where the total time of all the thermalization steps in the protocol goes to zero. Finally, as in the collision model, we expect that our results also generalize to cases where each of the κ (i) is chosen independently at random, since all that we require is that most of them are bounded away from 1. Finally, we connect these results to the continuum limit where S is continuously in contact with B. In our framework, this corresponds to the limit N → ∞ (i.e. t i − t i+1 → 0) and simultaneously κ → 1. A common scenario in open systems dynamics is an exponential convergence to a thermal state in time. We can model such a scenario by setting κ = e −t th /τF , where t th denotes the time that the system spends in contact with the bath in each step and τ F is some characteristic time-scale of thermalization. In the continuum limit, we have that t th → 0 and simultaneously N → ∞, which we obtain by focusing on the case t th = T /N , where T is the total time of the protocol. Using (47) and assuming that the model of exponential convergence remains valid in the limit t th → 0, we then obtain In this way, we rigorously obtain a convergence to ∆F iso with an error that is inversely proportional to the total time of the protocol, a behaviour that is often assumed on phenomenological grounds [40], or obtained by means of Markovian master equations [19,20]. In this way, our results make a connection between isothermal processes in discrete processes, where perfect thermalization is taken for granted [6,17,18,39,41,42], and in continuous ones [19,20].

Fluctuations
So far, we have only discussed average work-extraction. This has a good reason: In the current framework it is in general not possible to associate a classical random variable that measures the fluctuating work to the process due to the coherences that can be created by the unitaries U ti,ti−1 and the channels G (i) [43]. Let us emphasize, however, that in the limit N → ∞ in which the extracted work becomes optimal, these coherences become arbitrarily small. It then follows that we can associate a distribution of fluctuating work to the limiting process, which coincides with the limiting process when perfect thermalization and quenches of the Hamiltonian are assumed. Due to a result byÅberg [18], it follows that this limiting distribution is universal, i.e., independent of the trajectory H(t) and coincides with the one obtained from the collision model.

IV. CONCLUSION
In this work we have considered thermodynamic processes in which contacts between S and B correspond to partial thermalizations, which capture large classes of possible error models. The degree of thermalization is quantified by a parameter α (for the collision model) or κ (for the second framework), such that for α, κ = 0 we recover the standard case of full thermalization and for α, κ = 1 there is no interaction between S and B (and W in the collision model). Our main result is that optimal processes can be constructed for any α, κ < 1 with protocols consisting of sufficiently many steps. Furthermore, in the collision model and in the case of qubit systems, we have provided bounds for errors in terms of the number of steps, and characterized the fluctuations of work. In a second part of the paper, we have considered a different framework where work is extracted by changing the energy levels of S and interactions with B are described through quantum channels that provide an effctive model for thermalization. In this case, we have shown that any interaction with the bath that brings S closer to a thermal state can generate an isothermal process in the limit of slow transformations of S.

V. ACKNOWLEDGMENTS
We thank Guillem Perarnau for interesting discussions. We also thank Matteo Lostaglio and Chris Perry for useful comments on the manuscript. We acknowledge contributions from the Swiss National Science Foundation via the NCCR QSIT as well as project No. 200020 165843. M.P.-L. acknowledges support from the Alexander von Humboldt Foundation. All authors are grateful for support from the EU COST Action MP1209 on Thermodynamics in the Quantum Regime.  [44] C. Jarzynski, "Nonequilibrium equality for free energy differences," Phys. Rev. Lett. 78, 2690-2693 (1997).

APPENDICES Appendix A: Models of Noise
In the k th step of the protocol, the joint state of S, W and the k th bath qubit can be written as, where r k := q k (1 − p k−1 ), s k := (1 − q k )p k−1 , and note that r k > s k . p where τ A = e −βH A / Tr(e −βH A ) is an arbitrary thermal state. From the second law of thermodynamics, we know that the free energy of φ (deg) X can only decrease which here implies that the entropy of each φ (deg) X must increase -recall that φ deg X has a trivial Hamiltonian as it lives in an energy degenerate subspace. Hence we can write without loss of generality where z k (x) and α k (x) are in general functions of the specific τ A and V k , and α k (x) ∈ [0, 1]. On the other hand, the total state after the thermal operation reads For future convenience, let us also define Note that, for a specific choice of V , where z k (x) = α k (x)(1 − α k (x))(r k + s k ), we would get exactly the same result as after the application of (15). Another "extreme" case would be for z k (x) = 0, where the resulting state on SBW would be diagonal, corresponding to the resulting state of the map That is: with a given probability (1 − α k ), we perform the desired operation, and with an error rate α k we do nothing.
At the level of S, i.e., after tracing out also the bath qubit B and W, the coherence term z k (x) becomes irrelevant, and analogous to (16) the final state of S is given by Hence, for any unitary V k this more general model gives the same description of the resulting system state as the specific model we considered before (15), where V k only determines the value of α k . In any non-trivial unitary operation V k , i.e., for any α < 1, S and B get mixed, meaning that it brings the state of S closer to the thermal state of B. Analogously, when tracing out the system S and work battery W, we get the same evolution for the k th bath qubit, Here, ρ In particular, we note that in our previous considerations S interacts successively with each bath qubit and always on all three systems S, B and W together. This sequence is kept as we are interested here in quasistatic isothermal processes.

Appendix B: Discussion of noise types that the protocol can not tolerate
Here we can briefly discuss other imperfect operations which will lead to non-ideality of the process. The error models discussed previously are not the most general energy conserving unitaries. Indeed, since we assumed that W has a continuous Hamiltonian, there are many other degenerate subspaces. Therefore, there are in principle other energy-preserving operations, namely • First of all, note that we assume that in each step the correct bath qubit is chosen, i.e., the system interacts successively with each of the N bath qubits, such that we a have quasistatic isothermal process. Changing considerably the order would lead to energy-preserving dissipative processes.
• There can be energy-preserving interactions between S and W only, which would on average excite the system and therefore decrease the extracted work. The same would happen for random interactions between B and W, where in addition the correlations between S and W would be destroyed leading to uncontrolled fluctuations.
• When considering all three systems S, B and W together, there is not only the degeneracy between |0 S |1 B |x W ↔ |1 S |0 B |x + E B − S W , but there would also be a degeneracy between |0 S |0 B |x W ↔ |1 S |1 B |x − E B − S W due to the continuous Hamiltonian H W . Since starting in |0 S |0 B |x W is much more likely than starting in |1 S |1 B |x − E B − S , we would also loose work on average during this transition.
We can see that errors on all other kind of degeneracies prevent us from extracting maximal work. The only degeneracy on which errors can be compensated for is essentially the one on which the ideal protocol is interacting. While these considerations provide a qualitative analysis on the kind of errors that can/cannot be tolerated, we leave as an interesting task for the future to characterise a scenario where both kind of errors are present simultaneously.

Appendix C: Proof for Qubit Systems
In Section II B 1 we proved our result for simplicity for a specific distribution of excitation probabilities of the bath qubit. Here, we will show the general proof.
Let us choose a function f , s.t. ∀k ∈ N the excitation probabilities q k of the bath qubits are given by q k := f ( k N ). W.l.o.g. we assume that ρ S is initially purer than a thermal state, s.t. we can use its information to extract work. This function should be defined on the interval [0,1], starting at f (0) = q 0 and then be monotonically increasing until f (1) = q N . In addition, it should satisfy ∀k: We can prove that the requirements are satisfied for all functions that are of class C 2 on [0, 1], i.e., functions whose first and second derivatives both exist and are continuous on that interval. Due to the Weierstrass' boundedness theorem, we know that if the first derivative is continuous in a closed interval, then it is also bounded in that interval. This is equivalent to f being Lipschitz continuous with constant L on [0, 1], which already satisfies the first requirement (C1). For functions of class C 2 we can use the same reasoning one more time to argue that f (x) is bounded by a constant L'. Then we get Comment: If f is continuous only on (0, 1) and L if of the order of N γ with γ < 1, the proof still works the same, but results in a slower convergence. Here, we consider Lipschitz continuity on [0, 1] and therefore L = O(1).
This is sufficient to satisfy (C2), where we defined whose upper bound only depends on j, but not on k.
It can be shown that all such probability distributions {q k } k=1...N yield the optimal average work after N steps in the limit N → ∞ for any α > 0. Using (24) with those general excitation probabilities q k and the corresponding energy where E max := k B T ln 1−pmin pmin is the maximal energy gap of a bath qubit that the system is interacting with (in the case of work extraction it is E 1 ), which is O (ln N ) if ρ S is initially pure and O (1) otherwise. The first term can be identified as a Riemann sum, for which holds Here, f (q k ) = k B T ln 1−q k q k . The error γ for a monotone function with finite N is upper bounded by Thus, we can rewrite (C7):

Appendix D: Fluctuations
In order to calculate the fluctuations, we need to list each possible path with its corresponding probability and energy. We use vectors p [1] k that list the corresponding energies. With each step, the number of possible paths doubles, such that after k steps, we will have 2 k−1 -dimensional vectors.
As we have shown before, any of our considered noise models can be seen as "apply a swap" with probability (1 − α) and "apply identity" with probability α. Using (20), we can find a recursive way to describe the vectors: where q k corresponds to the excitation probability of the k th bath qubit and E k = ln 1−q k q k to the corresponding energy gap. We can calculate the average work after the m th step, as well as the square of it, since this is needed for calculating the variance, Next, we calculate the sum of all elements of p k , which corresponds to the excitation probability p k of the system qubit after the k th step: Using this recursive form, we can determine the explicit form: which gives us the same result as (22). In addition, we need the average value of the squared energy after m steps, where in the following W 2 [1] m denotes the entry-wise square product of each vector element Now we can calculate the variance after the total N steps as Writing the last term explicitly, which determines the average work after m steps when ending in the excited state |1 S , we get To simplify these expressions and show that the variance indeed goes to zero for N → ∞, we need to choose a probability distribution q k for the excitation probabilities of the bath qubits. For this proof we will use q k = k 2N and the corresponding energy levels are then given by E k = k B T ln 1−q k q k = k B T ln 2N −k k , but looking at simulations, it should also work for all the other distributions allowed in this paper (i.e., those mentioned in Appendix C). First, we will simplify some expressions that will help us later on: and estimate the sum: Using (D10) and (D12) we can already see that the last term in the first line of (D7) goes as O ln 2 N N and thus disappears for N → ∞. Now we rewrite the other terms, where the change of indices does not matter since E N = p 0 = 0. Writing the last two terms in the brackets explicitly, we get Next, we separately calculate the sum, using − ln(1 + which we can insert into (D17), Inserting this into the expression for the variance (D16) and rewriting the first line, we get In the last step we need to calculate these two sums, Thus, we get for the variance, where we can see that for infinitely many steps N the fluctuations vanish. This is in accordance with the Jarzynski equality e β∆F = e −βW [44], as for N → ∞ we are looking at a quasistatic process with ∆F = − W = F (τ Appendix E: Generalization to Arbitrary Systems Now we extend our previous considerations to generic qudit systems. We take as a starting point the model of noise (20). The derivation up to the second line of (24) still follows naturally. To move on, we notice the identity, where τ β = e −βHB / Tr(e −βHB ). Let us also define τ (k) β = e −βH (k) B / Tr(e −βH (k) B ). Then we can rewrite the second line of (24) as, Let us now define, We want to show that, = e −βH (k) B /Z k . Those we choose in such a way, that F (λ) := F (τ k β , H S ) is a monotonically increasing/decreasing function of class C 2 , i.e., is twice continuously differentiable in λ = k N on λ ∈ [0, 1]. Let us remark that that the free energy F (ρ, H S ) is not differentiable if ρ is a pure state because of the entropic term. Therefore, this derivation can not be applied when the initial state is pure. However, in practice this limitation is irrelevant, as any pure state can be approximated arbitrarily well by a mixed one. Coming back to our choice of H (k) B , note that where L is the Lipschitz constant of F (λ). In addition, we get where L is the Lipschitz constant of F (λ). This implies with Now we can compute W * : (1 − α) 2 α j j(j + 1) + N −1

k=0
(1 − α)α k k(k + 1) Note the similarity between this expression and previous expressions obtained for arbitrary choices of H (k) 's in the qubit case. Now we move to the more involved computation of E. Recall that we choose the Hamiltonians H (k) B so that F (λ) := F (τ k β , H S ) is a monotonic function that is twice continuously differentiable in λ = k N on λ ∈ [0, 1]. Similarly, we can define the family of Gibbs states τ (λ), with τ (j/N ) = τ (j) β . Now, note that S (τ (λ + x) τ (λ)) = x 2 2 which follows because S(ρ τ ) has a minimum at ρ = τ . Let us define, and the maximum of this function, with λ ∈ [0, 1]. On the other hand, for large N , we can write, Importantly, f (j/N ) in (E13) is a function independent of k and is bounded by f max , independently of the value of N . Now we have all the terms necessary to perform the computation, Hence we conclude that, which indeed goes to zero as N goes to infinity.
Note however, that the unitaries U i still appear implicitly, since This will be taken care of later. We now approximate each state σ (i) by the instantaneous thermal state τ (i) . We then get The first term can be interpreted as a Riemann sum and replaced by an integral in the limit N → ∞. This yields (see for example [39]) Tr ∆H (i) τ (i−1) = − where H (t) denotes the derivative of H(t) at t. The error in this approximation is again of order 1/N . We thus obtain It remains to show that the second term is of order 1/N . To do this, we bound the term σ (i−1) − τ (i−1) 1 . We first use the fact that the thermalizing maps have a thermalization parameter κ to write Since [τ (i) ,Ũ i ] = 0, the last term can be estimated as Adding and subtracting zero and using the triangle inequality again, we then arrive at Below, we show that the second term in the parenthesis is of order 1/N . We thus obtain for some constant K . Iterating the recursion, we hence obtain Since ∆H (i) Thus, if κ is independent of N , we obtain Note that in this case, the final work-value in the limit N → ∞ is independent of κ. Let us also discuss what happens if we let κ depend on the number of steps N . To do that, let us assume that for some characteristic parameters τ > 0 and γ > 0. For large N , this means that In particular, we have that Furthermore κ N → 0 for γ < 1, κ N → 1 for γ > 1 and κ N → exp(−τ ) for γ = 1 as N → ∞. As a consequence, we find that if γ ≤ 1: For γ < 1 the error vanishes in the asymptotic limit and for γ = 1 the error is of order 1/τ in this limit. If γ > 1, we find that so that the bound becomes useless for sufficiently large N . What is left to prove is that To see this, we use that (derived in E13) S(τ (t + δt) τ (t)) = O(δt 2 ), where S(· ·) is the quantum relative entropy. The result then follows using Pinsker's inequality ρ − σ 1 ≤ 2S(ρ σ) with δt = 1/N .