Extraction of ergotropy: free energy bound and application to open cycle engines

,


Introduction
The second law of thermodynamics governs what state transformations are allowed if work and heat are exchanged with the environment. These transformations are described in terms of the thermodynamic potentials, which, regardless of the microscopic details, provides fundamental limitations for energy transformations in thermodynamic processes. In particular, these limitations set the upper limits on power and efficiency. One of the most important thermodynamic potentials is the free energy, whose non-equilibrium generalisation is given as the so-called thermal operations, where an explicit work storage is introduced [4][5][6]. Remarkably, the inequality (2) can also be formulated as the (Jarzynski) equality if the free energy is treated as the fluctuating quantity [7][8][9][10].
On the other hand, if the system is isolated from a thermal bath, or the bath itself is finite, the inequality (2) cannot be saturated for processes leading to maximal change in free energy (reached in the Gibbs state), unless the initial state is Gibbs. Considering an isolated system driven by a cyclic force, Allahverdyan et al. [11] have replaced the free energy with a new thermodynamic quantity: ergotropy, defined as where the minimum is taken over all unitary operatorsÛ . If the stateρ is diagonal one can restrict to minimization over the set of permutations [12]. The minimum is achieved for the stateρ P referred to as the passive state of the system, and E(ρ P ) is called the passive energy of the stateρ. As it was shown in the seminal paper [11], ergotropy can be directly related to the free energy (cf. Section 2 below). The concept of the ergotropy appears in many areas of quantum thermodynamics. Firstly, it characterizes the optimal energy extracted from the so-called quantum batteries (see, e.g., [13][14][15][16][17][18]). Recently, it was related to measures of entanglement [19][20][21], coherences [22] and relative entropy [23]. It is worth noticing that ergotropy is in general not additive. This leads to the concept of complete passivity: the state is completely passive if the state composed of arbitrary number of its copies is passive. Remarkably, this links the idea of ergotropy with the equilibrium state, since any completely passive state has necessarily the Gibbs form [24]. Moreover, the complete passivity was recently related to the inability of extracting ergotropy in a presence of a catalyst [25].
Ergotropy is also connected with an explicit model of the translationally-invariant energy-storage, called the ideal weight. The ideal weight model was introduced in [5] to formulate a thermodynamic framework which satisfies the Second Law. In particular, the authors proved that work W performed via lifting the weight (increasing its average energy) is always smaller than the non-equilibrium free energy of the quantum stateρ (remaining in contact with a thermal reservoir): whereγ β is the Gibbs state at inverse temperature of the bath β and relative entropy S(ρ γ β ) is given by S(ρ γ β ) = Tr ρ(logρ − logγ) .
This result was subsequently generalized to an arbitrary correlated initial state of the weight and the system, and the extracted work was expressed as [26]: where ∆R is a change in ergotropy of the so-called effective control-marginal state of the systemσ i,f (with indices i and f standing for initial and final states, respectively), containing information about coherences and correlations influencing the work extraction process [26]. This reveals the connection between ergotropy and translational symmetry of the ideal weight model. Eq. (6) was further used in a formulation of the "tight Second Law", involving the ultimate bound given in Eq. (4), but lowered by corrections expressing the effect of initial coherences and finite-size of the heat bath [27].
Since ergotropy can be understood as work done by the system, it is worth to explore its relevance in the operation of quantum thermal machines. Quantum thermal machines were characterized from multiple perspectives: starting from quantum optics [28][29][30][31], thermalisation hypothesis in many-body systems [32,33], dynamical description in quantum open systems [34], and, most recently, theory of information [35] and resource theories [4,[36][37][38][39]. From the theoretical perspective, the family of heat engines in quantum regime can be categorized into two fundamental classes: autonomous and minimal coupling engines. Autonomous engines are the ones where the working body is simultaneously coupled to both heat baths, as well as to a system storing energy [40,41]. On the other hand, minimal coupling engines act in discrete strokes [26]. A simple construction of minimal coupling heat engines would enable us to emphasise and investigate the role of ergotropy extraction. On the contrary, for autonomous engines, the difficulty in optimizing work production and efficiency lies in the necessity for taking into account all possible simultaneous interactions between working body and heat baths, thus restricting the analysis only to particular physical models [42][43][44][45][46][47].
In this paper, we address the problem of ergotropy extraction via an energy conserving process. While recent interest in storing energy in quantum systems (quantum batteries) focuses mostly on cases when ergotropy is extracted via the time-dependent driving [48,49], here we concentrate solely on a fundamental, thermal contact scenario. In other words, we ask how much ergotropy can be extracted via a heat flow through the system coupled to the thermal bath. Since the ergotropy quantifies the amount of inversion of the population, one should expect that this is only possible via the non-Markovian evolution. Consequently, we use thermal operation framework [37,50,51] to optimize over such evolutions, emerging from energy conserving interactions. Our goal is to establish a connection between the ergotropy extraction and the Second Law of thermodynamics, as well as to characterize quantitatively the process of ergotropy extraction for finite and infinite dimensional systems. Alternatively, the problem of battery charging was recently formulated in terms of repeated interactions with a subsystem of environment [52], where a dissipative process leading to an active steady state of the battery was identified. Secondly, we focus on the application of ergotropy extraction process in a design of thermal machines. We introduce a simple model of quantum heat engines under a name of open-cycle engines. Open-cycle engines form a subclass of the aforementioned minimal coupling engines. The simple operational idea of the open-cycle engine is to use the gradient of temperatures to extract ergotropy from the hot bath using the state taken from the cold reservoir. The name open-cycle is motivated by the demand to discard or thermalize the used resource after the energy is stored in the work reservoir. The complete analytical analysis of these machines in terms of efficiency and power production per cycle is performed for two-and three-dimensional working body systems.
Here are the most important results that we obtain: • We formulate the Second Law of thermodynamics for a task of ergotropy extraction via any complete positive trace preserving (CPTP) map which preserves the Gibbs state of the system. The formulation is through an upper bound on the value of extracted ergotropy, expressed solely in terms of the non-equilibrium free energy; • We construct the process which saturates the bound for the quantum harmonic oscillator in the zero-temperature limit.
• On the opposite side of the spectrum of the dimensionality of the system, we identify optimal operations for ergotropy extraction for qubits and qutrits. We obtain this desrciption within the toy model of open-cycle engines, in which ergotropy extraction powers work production.
• We prove that optimal work production and efficiency of open-cycle engines is achieved for extremal thermal processes, disregarding the dimension of the working body.
This paper is organized as follows. We start by introducing ergotropy of an isolated system and a system coupled with a bath in Section 2. We derive upper bounds on ergotropy in terms free energy and relative entropy, and discuss the cases when these bounds can be saturated. In Section 3 we define the amount of ergotropy that can be extracted via any CPTP map that preserves the Gibbs state of the system and derive the upper bound in terms of quantum relative entropy (non-equilibrium free energy difference). We prove that the upper bound can be saturated for a harmonic oscillator system initialized in the ground state. We also provide the numerical simulation for systems with finite number of levels which shows how the value of ergotropy extraction tends to the general bounds. We proceed to describe the open cycle heat engines in Section 4, where we characterize a set of finite number of states over which optimal work production and efficiency can be achieved. Finally, we provide a complete analysis of the work production and efficiency for an open cycle heat engine with a qubit and qutrit working body. We conclude with the discussion in Section 5.

Ergotropy and free energy
We start from discussing the established bounds on ergotropy of an isolated system [13] and system coupled with a bath [27] and their saturation.
From the definition of free energy and ergotropy given in Eq.
(1) and Eq. (3), and using the fact that entropy is invariant under unitary evolution, we immediately have the following relation whereρ P is the passive state of the system as given in Eq. (3), and relative entropy is defined in Eq. (5). Note that the above holds for β taking arbitrary values. Among states with fixed entropy, energy of a system is minimized for the Gibbs state. This can be used in Eq. (7) to minimize the free energy component: where the inverse temperature β * is now uniquely specified by the relation S(ρ) = S(γ β * ). We see that the bound cannot be achieved in a single shot scenario if spectrum ofρ is different than that of γ β * , as in this case it is impossible to achieve γ β * fromρ P using an unitary transformation. However, Alicki and Fannes showed in [13] that this bound can be achieved by processing asymptotically many number of copies of the stateρ, i.e., We proceed by reporting a bound on the ergotropy of a system coupled with a bath at inverse temperature β [27]. Consider the system in the stateρ coupled with a thermal bath in Gibbs statê τ β at inverse temperature β. The composite state of the system and bath is expressed as a product ρ ⊗τ β . Based on Eq. (7), one can provide the following upper bound on the ergotropy for a system in contact with a bath: whereγ β is the Gibbs state of the system in the inverse temperature of the bath β. This is similar to the bound given in Eq. (8), with the crucial difference that the inverse temperature β is now fixed and inherited from the bath. In [5,27], the saturation of the bound is established for a quasi-reversible process consisting of infinitely many steps conserving the average energy, under the assumption of presence of an infinitely large heat bath.

Free energy bound for ergotropy extraction
In the previous section, we have discussed bounds on ergotropy of an isolated system and a system coupled with a bath in equilibrium. These bounds and their saturation relate the ergotropy of the system with the non-equilibrium free energy (see Eq. (8) and Eq. (10)). This suggest that the ergotropy can be understood as a generalization of the free energy for finite-dimensional quantum systems [27].
In this section, we move to the question of how much ergotropy can be extracted through system's coupling with a heat bath. We define ergotropy extraction for a state ρ evolving via a CPTP map Φ as the change of ergotropy in the process, i.e.: We are interested in inducing ergotropy on a system through an interaction with a single heat bath, and we assume the total energy of the system and bath to be conserved. As a consequence, we consider CPTP maps Φ of the form: such that [U, H S + H B ] = 0, where H S and H B are the Hamiltonian of the system and the bath, respectively. As before,τ β is the Gibbs state of the bath at inverse temperature β. These maps are known as thermal operations [37,50,51]. Let us stress that the defining property (12) is not restrictive, and serves only as a formulation of the first law of thermodynamics in the spirit of collision theory. Namely, one should consider the above unitaries U as generated by an arbitrary interaction Hamiltonian, with the interaction negligible at first and last stages of the process, when the systems remain spatially separated. Thermal operations constitute a subclass of a bigger set of transformations called Gibbs-preserving maps, which preserve the Gibbs state of the system: Φ(γ β ) =γ β . The sets of Gibbs preserving maps and thermal operations are equal when restricted to statesρ which are diagonal in the energy eigenbasis i.e., [ρ, H S ] = 0. Before addressing the question of ergotropy extraction via thermal operations, we propose a bound on the extraction by a general CPTP map.

Proposition 1.
Consider an arbitrary stateρ that is evolving via a CPTP map Φ. Then, the ergotropy of the final state Φ(ρ) can be decomposed into: whereγ β is a Gibbs state of the system at some inverse temperature β.
The formula (13) follows from the application of the identities given in Eq. (7) : In Eq. (13) we identify three different contributions to the ergotropy of the final state Φ(ρ). The first one 1 β S(ρ γ β ) is the free energy difference of the initial stateρ and the Gibbs stateγ β . This expression has appeared in the derivations of the bounds on ergotropy (see Eq. (8) and (10)), and provides the ultimate bound for the work extraction from a single heat bath (see Eq. (4)). The second contribution − 1 β S(Φ(ρ) P γ β ) is always non-positive (due to the non-negativity of the relative entropy) and characterizes a thermodynamic distance of the final passive state to the equilibrium. In particular, this term is zero if and only if the final state Φ(ρ) has the same set of eigenvalues as the Gibbs state, so that Φ(ρ) P =γ β . Finally, the last term: can be interpreted as the entropy production for a CPTP map Φ. Indeed, for the maps Φ being thermal operations as given by (12), the change of the system energy can be identified with heat, i.e., ∆E = Q, and then F (Φ(ρ)) − F (ρ) = Q − 1 β ∆S. Considering the map Φ to be Gibbs-preserving leads to the crucial corollary Corollary 2. If the CPTP map Φ is Gibbs-preserving, i.e., Φ(γ β ) =γ β , then ergotropy of the final state Φ(ρ) is upper-bounded as Proof. The proof of the corollary immediately follows from two properties of the relative entropy and Gibbs-preserving map: the non-negativity of the relative entropy and monotonicity of the relative entropy under Gibbs preserving maps: Inserting these Eq. (17) and Eq. (18) in Eq. (13), we obtain the bound.
The inequality (16) shows that the final ergotropy cannot be higher than the initial thermodynamic resource given by the non-equilibrium free energy. This remains in agreement with the bound (4), derived for a model which defined work through changes of the avarage energy of the quantum weight. As mentioned earlier, in case of no quantum correlations between the ideal weight and the system, its ergotropy can be associated with maximal change of average energy of the weight. Then the bound (16) is just another statement of the Second Law.
This result can be generalised to encapsulate ergotropy extraction for non-passive initial states via Gibbs-preserving maps.
Proof. We begin by writing the change in ergotropy using Eq. (7) as where the first inequality follows from the monotonicity of free energy under the Gibbs-preserving map and the second inequality follows from non-negativity of relative entropy.
Note that the bound given in Eq. (19) does not depend on a particular CPTP Gibbs-preserving map Φ, but only on the initial stateρ. This bound is consistent with the previous bound in Eq. (16), since S(ρ P γ β ) ≤ S(ρ γ β ). If the initial stateρ is already passive R(ρ) = 0, then inequality in Eq. (19) boils down to Eq. (16). Therefore, we can interpret Eq. (19) as the second law for ergotropy extraction that sets the upper bound on the amount of ergotropy that can be extracted via a Gibbs preserving transformation.
In the next section, we design an example showing saturation of the bound for Φ given by a thermal operation (12).

The maximal-energy thermal process and saturation of the bound on ergotropy extraction
In order to identify a scenario in which the general bound (19) is saturated, we restrict ourselves to a subset of Gibbs preserving maps, generated by thermal operations. We recall the notion of thermal processes [37], which completely characterise the action of a thermal operation on the diagonal of a quantum state.

Lemma 4.
Consider a given stateρ that is diagonal in the energy basis, i.e., where {| } d i=1 are energy eigenvectors, and a thermal operation Φ which transforms the stateρ toσ. This is equivalent to existence of a stochastic matrix A Φ such that where p, q and γ β are probability vectors generated from diagonal entries ofρ,σ, andγ β respectively, i.e., The stochastic matrix A Φ associated with thermal operation Φ is called thermal process.
The bound on ergotropy extraction given by (19) is saturated for a thermal operation Φ max acting on a ground state of a harmonic oscillator system,Ĥ = n nω|n n|. Its associated thermal process A Φmax is given by where with the choice Figure 1: Ergotropy extraction for the d-dimensional qudit with frequency ω prepared in the ground state, for the thermal process at temperature β that maximizes the final energy. Horizontal dashed lines correspond to the ergotropy extraction bound (19). It can be seen that the final ergotropy for the maximal-energy process saturates with growing dimension d, and it approaches the optimal ergotropy if δ → 0.
taking integer values for 1 βω log Z ∈ N, where Z = Tr e −βĤ . One can immediately see that A Φmax is stochastic, and that it preserves the probability vector The action of A Φmax on the probability vector (1 0 . . . 0) T determines the transformation of the ground state: It can be seen that the stateĜ is the Gibbs stateγ β with populations shifted by ωL. In order to calculate the ergotropy of the final state, we find the corresponding passive stateĜ P , i.e., the minimal energy state obtained via the permutation of occupation probabilities. For this infinite-dimensional system, the passive state is found in the limit where all "zeros", which occupy the first L energy levels, are "sent to infinity" (as in the Hilbert's Hotel paradox). In this limit, we haveĜ P =γ β , i.e. the passive state ofĜ is the Gibbs state. The proposed construction of the passive state does not apply in finite dimension, since the Gibbs stateγ β in this case does not have any "zeros" within its spectrum. Putting G P =γ β turns the LHS of the inequality (19) to while the RHS of the inequality is equal to: We see that the harmonic oscillator in the zero-temperature limit (i.e., in the ground state) and with frequency ω satisfying 1 βω log Z ∈ N saturates the bound of the ergotropy extraction (19). The thermal process which lead to the saturation maximizes the final energy of the state. This is visible from the fact that all the action of the process acting on the ground state |0 0| is described by its first column. It can be easily verified that for the n row of the first column of Ω, the value of the corresponding matrix element is upper bounded by e −βω(L+n−2) (otherwise the map is not Gibbs-preserving). However, from (27) we read the value of this element to be equal to it: e −βω(n−1) Z = e −βω(L+n−2) . Therefore, the map leads to maximal possible occupation of high energy states, at the expense of the lowest L − 1 energy levels.
Below, we extend to a finite dimension the protocol maximazing energy of the final state, and calculate the amount of extracted ergotropy. Our numerical calculations are based on the fact that the state with maximal energy which can be obtained via a thermal operation, can be associated with a so-called thermomajorozation curve which is tightly-thermomajorized by the curve associated with the initial state, and which has the so-called descending β-order. We refer the reader to Lemma 16 of Supplementary Material for the detailed description of the protocol.
Let us define the detuning parameter: In Fig. 1 we present the results of the numerical simulation for the final ergotropy of the d-dimensional system with different frequencies ω and non-zero value of the detuning parameter δ, with the ergotropy extraction protocol maximizing the energy of the final state. In the regime δ → 0, extractable ergotropy approaches the bound already for moderate dimension d, while the non-zero detuning parameter prohibits the strict saturation even in the limit d → ∞. Note that in the limit β → 0, thermal processes become permutations, and the protocol maximizing energy of the final state is the one maximizing its ergotropy. Therefore, the numerical results presented in Fig. 1 indicate that no thermal operation is able to strictly saturate the bound (19) for non-zero detuning δ and high temperatures of the bath, with harmonic oscillator prepared in the ground state.
In the next section, we move to the case of small dimensional systems, in which the bound for ergotropy extraction can be far from being saturated. We analyze these cases in the context of utilizing small dimensional system as working bodies in thermal machines which extract ergotropy from the environment. We will show that for finite temperatures of the hot bath, optimal protocols for ergotropy extraction dependents on the spectrum of the Hamiltonian and bath temperatures, even if we relax cyclicity constraints by demanding that the initial state is in thermal equilibrium with the cold bath.

The open-cycle heat engines
In this section, we consider the model of the open-cycle heat engine and describe the role of ergotropy extraction from the bath as work that can be stored in a battery modelled by ideal weight [5]. We shall analyze the optimal performance of such an engine and characterize the states for which optimal work production and efficiency can be achieved. We consider open cycle engine within the family of minimal coupling engines [27], where the working body alternately couples to the heat baths and the battery, resembling traditional stroke engines.

Description of the open cycle heat engine
An open cycle heat engine is a minimal coupling thermal machine introduced in [26], consisting of two baths in equilibrium associated with two different inverse temperatures β H and β C , and a battery. We denote by β H the inverse temperature of the hot heat bath, while β C stands for the inverse temperature of the cold bath, and β H < β C . The working body is initially in a thermal state with respect to the cold bath. In subsequent stages, it alternately couples to the hot bath and the battery through discrete two-body energy conserving operations referred to as strokes. Initially, the working medium is thermalized with respect to cold bath. In stroke 1, working body extract ergotropy from hot bath, that induces non-passivity in the working body. In stroke 2, ergotropy of the working body is stored in the battery as work. This interaction makes the state of the working body passive. Finally, in stroke 3, working body can be discarded or can be thermalized with the cold bath.
The design of the open-cycle engine is schematically presented in Fig. 2. As in the minimal coupling quantum heat engines, open-cycle engine in its first stroke utilizes the interaction between the hot bath and the working body to extract ergotropy. In the second stroke, extracted ergotropy is stored as work on the battery. In the final stroke, the working body gets thermalized with the cold bath (or discarded). We proceed with the detailed description of the engine.

Working body thermalized with respect to cold bath.
At the beginning of the procedure, we initialize both of the baths in equilibrium i.e., whereĤ H andĤ C are Hamiltonians of hot and cold reservoir, respectively. We assume that heat capacities of the baths are infinite, thus the temperatures baths remain unchanged through the operation of the engine. The working body is a d-dimensional system initially in equilibrium with cold bath. The Hamiltonian associated with the working body is given bŷ while state can be written asρ The initial state of the engine can be written aŝ where ρ B is a state of the battery. The free HamiltonianĤ 0 of the combined system is given bŷ whereĤ S ,Ĥ H ,Ĥ B are Hamiltonian of the system, hot bath, and battery respectively. If the opencycle engine runs iteratively, such that after storing the ergotropy of the working body in the battery it has been thermalized again with cold bath, then the engine at the beginning of each of the runs is affected only by the change of the battery state ρ B . The role of the working medium is to steer the energy flow from hot bath to the battery, as described below.

Interaction of working body with the hot heat bath and battery via discrete two body energy conserving strokes.
Following the definition of stroke operations in [26], we allow only for the interactions among any two bodies at a given moment. This implies that the unitary evolution of an engine can be decomposed into a product of two unitaries: whereÛ H describes the interaction between the hot bath and the working body,Û B describes the interaction between the system and the battery. In order to account for energy in the engine, we demand that the unitaries conserve the total energy. Therefore, we impose the constraint The purpose of a unitaryÛ H is to extract ergotropy from the hot bath to the working body (ergotropy extraction step) in the first stroke, whileÛ B aims at charging the battery by converting the ergotropy that has been extracted from the hot bath into work in the second stroke. This justifies the ordering of the unitaries given in Eq. (38).
The interaction between the hot bath and the working body generatesÛ H , which induces a thermal operation E H given in Eq. (12) which leads to extraction of ergotropy from the hot bath: As the initial stateρ S of the working body is in equilibrium with cold bath, we have R(ρ S ) = 0. Therefore, for a work storage modeled by the ideal weight (see discussion in Supplementary Material A.3), work stored at it in the second stroke is equal to ergotropy extracted in the first stroke [5,26]: Since the working body interacts with the hot bath, the heat flow from the hot bath to the working body increases its average energy. The transferred amount of heat Q is defined as the change of average energy of the system due to interaction with the hot bath, i.e., In agreement with Clausius's statement of second law of thermodynamics, in the absence of external work source, heat flows from hot to cold body. The definition (42) abides this, which ensures that the transferred amount of heat is always non-negative: Lemma 5 (Directionality of heat flow from hot bath to a working body that is thermalized with cold bath). Consider a stateρ S that is initially thermalized with cold bath at temperature β C , transforming via any thermal operation E H with respect to the hot bath at temperature β H , then the amount of the exchanged heat, given in Eq. (42), is always non-negative i.e., Proof sketch: From the linearity of the function Q(E H (ρ S )) we conclude that minimal value of the exchanged amount of heat is obtained at one of the extremal points of the set of states T (ρ S ) -the set of states which can be achieved by applying a thermal operations onρ S . One of those extremal points is the stateρ S itself. Thus we have to show that minimum of the average energy of extremal points of T (ρ S ) is never smaller than average energy ofρ S . To this end we shall first note that there is only one extremal element of this set with β-order (for definition of β-order see section A.1 of appendix) given by (1, 2, . . . , d). Next we show, that for any state with different β-order, there exists a thermal operation A which decreases its average energy. This implies that for extremal states of T (ρ S ) which are different thanρ S , their average energy is strictly higher than average energy of some other states from T (ρ S ) (i.e. those resulting from action of A). Therefore the minimum must be achieved onρ S , which gives the minimal value of Q equal to zero. The complete proof can be found in section A.2 of the Supplementary Material.
Finally, we define the efficiency of open cycle heat engine as the ratio work stored in the battery with the amount of heat exchanged from the hot bath i.e., In the next section we shall discuss in detail the optimal performance of the open cycle heat engine. Precisely, how one can deposit the maximum amount of work in the battery and perform it the most efficiently? This question is equivalent to extraction of maximum ergotropy from the hot heat bath as one can see from Eq. (41).

Optimal performance for the open cycle heat engine
Below, we optimize the work and efficiency defined in Eq. (41) and Eq. (44) over the set of states that can be achieved from the state of working bodyρ S via a thermal operation. As a set of thermal processes forms a polytope, it immediately follows from lemma 4 that the set of states achievable via thermal operations from any given diagonal state also forms a polytope. We denote by T (ρ S ) the set of states which can be achieved fromρ S via a thermal operation. This set is called thermal polytope of the stateρ S . Therefore, we are interested in maximizing work and efficiency over all the states in The theorem below characterises the states in the T (ρ S ) which lead to optimal work and efficiency. Proof. Optimal work is obtained by maximizing the ergotropy function over the polytope T (ρ S ) (see Eq. (45)). From the definition of ergotropy given in Eq.
(3), we can prove the convexity for ergotropy as follows, whereρ andσ are any arbitary states and α ∈ [0, 1]. We exploit the fact that if we maximize a convex function over a polytope, the function takes maximal values at the extremal points. This implies that the maximum value of work production is achieved at an extremal point of T (ρ S ). For maximizing efficiency, we notice that an arbitrary stateσ ∈ T (ρ S ) can be written asσ = where the first inequality exploits convexity of ergotropy and linearity of the function Q(·), and the second inequality uses lemma 5 (Q(σ j ) is non-negative for all extremal states, and therefore Therefore, maximization of work and efficiency of the open-cycle heat engine is tantamount to optimisation over all extremal points of T (ρ S ). Note that these two may not be achieved simultaneously. Now we shall discuss the open cycle heat engine with qubit and qutrit working body and analyze their optimal performance, in this way addresing the ergotropy extraction process for low dimensional systems. We use the results of [53] to characterize the extremal points of the set T (ρ S ).

Open cycle heat engine with qubit working body
Consider a qubit working bodyρ S with Hamiltonian ω|1 1| (48) in thermal state at inverse temperature β C given bŷ From [37], it is straightforward to characterize T (ρ S ) as convex hull with two extremal points as follows whereρ S is given in Eq. (49). From Theorem 6, we optimise the performance of open cycle heat engine over extremal points of T (ρ S ) as given in Eq. (50). Therefore, we evaluate work production and efficiency of the engine at which is the only extremal point that gives non-trivial work production and efficiency. We calculate the optimal work output defined in Eq. (41) as Therefore, the engine produces non-zero work iff 2e −β H ω − 1 > e −β C ω . To calculate the efficiency, we proceed by calculating the amount of heat exchanged given in Eq. (42) as Therefore, the optimal efficiency is given by Let us now compare with the more general case of minimal coupling quantum heat engines [26], for which we have Thus W mc ≥ W and η mc ≥ η, with the equality holding in the limit β C → ∞. This relation is in agreement with the fact that open cycle engines are minimal coupling engines where working body goes to its initial state (after storing work) by thermalization with the cold bath. From (52) and (55) it is also visible that the range of temperatures in which the machine can work (non-zero work is produced) is larger for minimal coupling heat engines. Now we move to the qutrit case of open-cycle heat engines, which, to the best of our knowledge, for minimal coupling heat engines remains uncharacterized.

Open cycle heat engine with a qutrit working body
As a working body we select a qutrit working body with Hamiltonian that is initialized in thermal state given bŷ Like before, by employing Theorem 6 we optimise work output and efficiency over all extremal points of T (ρ S ). To do so, we employ a result from [53] which characterizes the set T (ρ S ) as a convex hull of its extremal points. We are going to state the lemma by using the following notation where ω 0 = 0. We rewriteρ S given in Eq. (58) aŝ where Z C = 1 + q C 10 + q C 20 . The application of [53] (see Supplementary Material, Lemma 17 for details) results in the following: Lemma 7 (Extremal points of the qutrit thermal polytope). The set of achievable states via thermal operation in presence of bath at inverse temperature β H , from the initial stateρ S given in Eq. (61), is given by the convex hull of the extremal points: where β 0 is specified by the relation and extremal points are defined aŝ Extremal states Amount of ergotropy extracted where Z C = 1 + q C 10 + q C 20 .
The maximization over these extremal points turns out to be difficult as it dependents on the values of β H , ω 1 and ω 2 . We start from evaluating work produced by the engine as the amount of ergotropy extracted from the hot bath, according to (41). We enlist the values in Table 1. Values of heat exchanged with the hot bath for each of the extremal states are given in Table 2, and efficiency simply follows from taking the ratio of the amount of ergotropy extracted with corresponding amount of heat exchanged.
As optimal values of ergotropy and efficiency dependent on the temperature regime and system Hamiltonian, we analyze the work production and efficiency for a qutrit Hamiltonian given by We identify optimal protocols for work production. These can be characterized by different states of the working body (64)-(69)) obtained by application of a thermal operation to the initial state (see Fig. 3). It should be stressed that, in contrast to the infinite dimensional case investigated before, for low dimensional systems the optimal protocols for ergotropy extraction may not maximize energy of the working body. In Fig. 4 we show the values of extractable ergotropy and efficiency of the engine. We see that the region in which the ergotropy is extracted shrinks with increasing ω. In the limit ω → ∞ and β C → ∞ we recover the condition β H < log 2, bounding the operational regime of a qubit minimal step engine [26]. Indeed, for high values of ω, extremal thermal processes acting on a qutrit Gibbs state cannot lead to higher population inversion than the ones acting on a qubit system (see Supplementary Material Sec. A.4 eq. (92-93)).
Maximization of energy (provided by a process leading to the stateρ 6 S ) is the optimal strategy in the limit β H → 0, but otherwise the landscape of optimal protocols is very complex. In particular, close to critical values of β H satisfying 1 = e −β H + e −β H ω , efficiency and extracted work may strongly depend   Figure 3: Protocols for optimal ergotropy extraction from a qutrit Gibbs state with inverse temperature βC , subjected to interaction with thermal bath with inverse temperature βH . Hamiltonian has spectrum (0, 1, ω).

Extremal states
Amount of Heat Exchanged on β H . This should be attributed to different structures of the set of extremal thermal processes above and below this limit [53]. We expect that diversity of optimal strategies increases with the dimension of the working body.

Open cycle heat engine with qudit working body
As it was seen in the previous section, the structure of non-passive states becomes more complicated while going from the two-level to three-level systems. Thus, in order to characterize the higher ddimensional qudits, we compute optimal work production and efficiency via the numerical simulation. The method of computing the optimal quantities is based on the thermomajorization diagrams (which we discuss in Supplementary Material A.1): For a particular initial stateρ S , we numerically search for all extremal points in the set T (ρ S ).
The results are presented in Fig. 5. Due to the fact that number of extremal points in a thermal polytope grows with factorial of the dimension, we were not able to perform these calculations for moderate d > 8. Nevertheless, optimizing over all possible protocols for ergotropy extraction suggests that work production and efficiency increase with the growing dimension of the working body, while the bound given in (19) cannot be saturated for small dimensions. These conclusions are similar to those valid for a restricted class of protocols aiming at optimizing energy (Fig. 1), and match the intuition that low dimensionality of the working body enforces separation (in terms of relative entropy) between the Gibbs state and the state emerging from ergotropy extraction, prohibiting saturation of the bound. Therefore, using low dimensional working body may require careful optimization.

Outlook
We have proposed a bound on ergotropy extraction, showed that it can be saturated for infinite dimensional systems, and, for finite dimensional systems, we have investigated the gap between the bound and the value of extractable ergotropy. We use the idea of ergotropy extraction to design a stroke engine in the quantum regime, called open cycle engine. We prove that optimal work production and efficiency for these engines are achieved for extremal states of the thermal polytope of the initial state. We have analytically characterized the optimal protocols for work production and efficiency for such engines in two and three dimensions and analyze higher dimensional cases numerically.
In the construction of open cycle engines we have assumed that the state of the working body is initialised in the Gibbs state, hence performance of the engine was not affected by coherences in the local Hamiltonian basis. However, understanding the role of coherence in ergotropy extraction in general remains an open problem. We may formalise it as follow: consider a working body in an initial stateρ (with coherences). Is it possible to achieve a stateσ such that R(σ)−R(ρ) > R(D(σ))−R(D(ρ)) where D(·) denotes dephasing ? If it was so, it would constitute a quantum advantage for ergotropy extraction, and may lead to constructions of heat engines with higher work yield per cycle than their classical counterpart.
Secondly, it remains an open problem to investigate the relation between extractable ergotropy (optimized over all energy conserving operations) and the dimension of the working body. We have numerically attested this in Fig. 1 for a qudit system initialized in the ground state being transformed via a thermal operation that exchanges maximum energy from bath.
On the general note, the bound (9) attests for the ergotropy of a composed system in the asymptotic limit of copies, singling out the Gibbs state as a completely passive state. Analogously, one could ask about extractable ergotropy for composed systems coupled with a heat bath, investigating the effect in which the emerging entanglement within the system influences the extractable ergotropy. Finally, to complete the characterization of the closed cycle engines [26], protocols for higher dimensional working body systems would need to be developed by optimizing ergotropy extraction, with cyclicity constraints taken into account.
Finally, a problem of finding physical systems saturating bounds derived with the help of thermal operations boils down in our case to the question if strong dependence of optimal values of efficiency and work production on system Hamiltonian and bath temperatures is present within a fixed physical implementation of an open-cycle heat engine. We leave as a subject of future research to check if this dependency results from identifying different energy preserving interactions as optimal in different temperature regimes.

A Supplementary information A.1 Mathematical preliminaries of Thermal operation
In this section, we shall discuss the required mathematical properties of thermal operations and Gibbs preserving operations which have been used to analyse the ergotropy extraction protocols and the performance of open-cycle heat engines.

Definition 8. Thermal operation:
A complete positive trace-preserving (CPTP) map Φ is a thermal operation on the system in the stateρ if it transforms the product state of system and bath at equilibrium via an energy conserving unitary, followed by tracing out the bath: where U is a joint unitary commuting with the total Hamiltonian of the system and bath [U,Ĥ S +Ĥ B ] = 0, andτ β is a thermal state of the bath at some fixed inverse temperature β.
Conservation of energy by the unitary can be understood as the first law of thermodynamics. The operation preserves the Gibbs state: Φ(γ β ) =γ β . Next, we shall state the theorem which establishes the equivalence between thermal operation and Gibbs-preserving operation when they acts on a state that is diagonal in energy eigenbasis.

Theorem 9 (Equivalence between thermal operation and Gibbs-preserving operation for a state that is diagonal in energy eigenbasis ([37]). Consider a stateρ is diagonal in eigen basis of
are eigenvector ofĤ S with corresponding eigenvalues i . Then the set of achievable states from a given stateρ via thermal operation in presence of a bath at inverse temperature β, and Gibbs-preserving map that preservesγ β are identical.
Since in this work we want to analyze the extraction of ergotropy and optimal performance of open cycle engines, we focus on transformation of states which are diagonal in the energy eigenbasis. Setting up this framework, the central question to ask is what are the conditions on two statesρ andσ that are diagonal in the eigenbasis of HamiltonianĤ S , such that they are related via a thermal operation Φ: Φ(ρ) =σ.
By Lemma 4 from the main text of the paper, the transformation can be described by a stochastic matrix A Φ such that probability vectors p and q generated from the diagonal entries To describe the necessary and sufficient conditions for the existence of a stochastic matrix satisfying A Φ p = q and A Φ γ β = γ β , we proceed by introducing the concept of β-ordering of a population vector and thermo-majorization curves. These concepts are illustrated in Fig. 6.
k=0 will be called elbows of the curve β(p). The curve is concave due to the non-increasing order of elements in d. Let us note that there might be more than one permutation leading to a concave curve β(p). of ρ is defined as a vector π(1 . . . d).

Definition 11. β-order: β-order
β-order shows the order of segments assuring convexity of β(p). See Fig. 6 for examples of curves with different β-orders for a d = 3 case.
Thermomajorization curves are used to characterize possible transitions between states under Thermal Operations [37] : Note that the set of all stochastic matrices which preserves the probability vector γ β forms a polytope (see [54] for its complete characterization), thus the set of all achievable states also forms a polytope. We call this polytope a thermal polytope ofρ, and denoted by T (ρ). Its extremal points are defined with the help of the notion of tight thermomajorization:   Figure 6: Thermomajorization diagram for a three dimensional system with HamiltonianĤ = 1| 1 1| + 2| 2 2|, 2 > 1. Elbows of curves are indicated by circles. Curves β(s) and β(q) are thermomajorized by curve β(p). Curve β(r) does not thermomajorize β(p), nor β(p) thermomajorizes β(r). Curve β(s) is tightly thermomajorized by  curve β(p). p and r have β-order (3, 2, 1), while s, and q have β-order (1, 2, 3).
Finally, we will show that a thermal process which leads to descending β-order of the final state, i. e. (d, d − 1, . . . , 1), maximizes the final energy, providedσ belongs to the set of extremal states of T (ρ):

Lemma 16. For every initial stateρ with probability vector p in dimension d, among all statesσ, with probability vector q, which can be obtained from it under thermal operations, the state with the highest energy corresponds to thermomajorization curve β(q) which is tightly thermomajorized by β(p)
and has β-order (d, d − 1, . . . , 1).
Proof. As energy of the state depends solely on its diagonal, it can be deduced from the thermomajorization curve of the state (even in a case when some non-diagonal elements are initially present). A set of states achievable by thermal operations from a given initial state is a convex set. As energy is linear with respect to mixture of states, it is maximized for an extremal state in this set. From Lemma 16 and Lemma 13 we see that thermomajorozation curves of these extremal states are tightly thermomajorized by the curve of the initial state. Therefore, what remains to be shown is that the β-order of this curve should be (d, d − 1, . . . , 1). Assume that the β-order of the final state is different, i.e. that in (α 1 , . . . , α k , α k+1 , . . . , α d ) there exist α k and α k+1 such that α k < α k+1 . It is straightforward to see that swapping α k with α k+1 and transitioning instead to a final state with βorder (α 1 , . . . , α k+1 , α k , . . . , α d ) leads to a state with non-lower energy, as the modified β-order implies non-decreased occupation on level α k+1 , and non-increased occupation on the level α k , compared to the previous final state, while level α k+1 is associated with higher energy than α k , according to the assumption α k+1 > α k . We can use the above argument to transit in steps from an arbitrary β-order to (d, d − 1, . . . , 1), while having the energy of the final state unchanged or increasing in each step.
We use the above lemma to construct states of maximal energy in the thermal cone of the ground state of finite dimensional system, and analyze their ergotropy in Fig. 1.

A.2 Proof of Lemma 5
We begin by re-writing the definition of the amount of heat that is exchanged with hot reservoir given in Eq. (42) whereρ S is the state of the working body that has been thermalized with cold reservoir at inverse temperature β C as given in Eq. (35). We see that function Q(·) is linear. Therefore, when we minimise the function Q(·) over all the states that can be achieved fromρ S via thermal operation i.e., where the minimum is achieved at a extremal point labelled byσ * S of the set T (ρ S ). Therefore, it is enough to prove Q(σ * S ) ≥ 0 which we shall do by showingσ * S =ρ S . From lemma 4, we see that applying thermal operation on the stateρ S in presence of a hot bath at inverse temperature β H , is equivalent to applying a stochastic matrix on the probability vector such that the stochastic matrix preserves Note that the probability vector p S has β-order (12 . . . d). From lemma 13, we see that If s * S has β-order (12 . . . d) then using lemma 15, immediately gives s * S = p S . Now consider the situation when s * S has β-order (α 1 . . . α d ) = (12 . . . d). Then, there will be at least two conseceutive entries α k and α l in the β-order (α 1 . . . α d ) such that α k > α l . Let's call α k = i and α l = j where i , j ∈ {1, . . . , d}. Since i > j, therefore i > j since conventionally we take 1 ≤ 2 ≤ . . . d . Now we shall construct the following stochastic matrix where I Rest is the identity matrix on the subspace orthogonal to the subspace spanned by {| i , | j }.
It is easy to see that A(γ β H ) = γ β H where γ β H is defined in Eq. (78). Next, from the definition of β-order given in A.1, we can write since i = α k appears before j = α l in the β-order (α 1 . . . α d ). Note that we can always choose α k and α l in such a way that Eq. (81) is a strict inequality. If it is not true, this will lead us to conclude that s * S = γ β H which is not the case, because σ * S has to be the extremal point of T (ρ S ). From Eq. (81) we have, Now, applying the stochastic matrix A on s * S transforms the probability on the energy level | j and | i as and leaves the probability on the other level unchanged. As i > j , employing (82), we immediately have a state with diagonal entries given by As * S has strictly less average energy thanσ * S , that can be via a thermal operationσ * S . To put it simply , using thermal process A that acts non-trivially only on the energy levels | i and | j , we have decreased the population on the higher energy level | i , and increased the population in the lower energy level | j . Therefore, the state with diagonal entries given by As * S , has strictly less energy thanσ * S which is contradicting with the assumption thatσ * is the sate with lower energy in the set T (ρ S ). Therefore, we conclude thatσ * S has to be equals toρ S .

A.3 Explicit work storage battery modelled by translational invariant weight
In this section, we introduce a battery model that can be use to store the ergotropy present in a system as work. This battery is modelled by an ideal weight [5,38] with translational invariant symmetry. It is described by the HamiltonianĤ while the unitariesÛ B are restricted to these satisfying translational invariant symmetry given by whereΓ is a shift operator which displaces the energy spectrum of the weight, i.e.Γ † Ĥ BΓ =Ĥ B + , and is an arbitrary real constant. The above assumption implies that the resulting transformation on the battery does not depend on its initial state, a welcomed property for its iterative use. It was shown that work defined as the change of average energy of the battery, resulting from the unitaries with the above constraint, satisfies the Second Law of Thermodynamics [5]. Moreover, for the case of the battery state having no coherences in the eigenbasis ofĤ B , it was shown in in [1,10,27] that the average energy gain corresponds to ergotropy of the working body.

A.4 Detailed analysis of open cycle heat engine with qutrit working body.
In this section we provide the detailed calculation of the amount of extractable ergotropy from the hot bath for the initial state that is given bŷ where Z C = 1 + q C 10 + q C 20 where and q C ij is introduced in Eq. (59) as q C ij = e −β C (ωi−ωj ) . We write the probability vector generated from the diagonal ofρ S as We can immediately see the β-order of the of the probability vector at inverse temperature β H is p S is (123) since Now we shall employ the lemma from [53], that characterizes the extremal thermal processes that act on any probability vector p with β-order (123) and give the spectrum of the extremal states of T (ρ S ). Table 1 of [53]] For any probability vector p with β-order (123), thermal processes at the inverse temperature β H that gives the spectrum of non-trivial extremal states in T (ρ S ) are given by

Lemma 17. [Based on
where β 0 is defined by the following relation and Moreover, for the probability vector with β-order (123), the final β-order for the thermal process given in Eq. (90) is given by Thermal processes β-order A 1 (213) A 2 (132) A 5 (312) A 9 (231) or (321) A 12 (231) A 13 (321) Using this lemma, we shall analyze the optimal performance of open cycle heat engine with initial state having spectrum p S i.e., calculate the optimal ergotropy that can be extracted from the hot bath via thermal process A i and stored it as work. We denote it as R Ai (p S ). Thus we can calculate the efficiency of such heat engine.
Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 1 : Note that thermal process A 1 transforms the initial probability vector given in Eq. (88) in the following way where the last inequality follows from the fact 0 ≤ q H 20 ≤ 1. Therefore, we have the following order of population that is compatible with Eq. (96) is given by Note that if the order of the population given by Eq. (97), then no ergotropy extraction is possible from the final state because it is passive. Finally, if the order of the population of final probability vector in Eq. (95) is given by Eq. (98) then the amount ergotropy that can be extracted is given by Thus, a condition for positive ergotropy extraction is given by Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 2 : Thermal process A 2 transforms the initial probability vector given in Eq. (88) in the following way From the table 17, we have the β-order of the final probability vector in Eq. (101) is (132). That means where the last inequalities follows from 0 ≤ q H 10 ≤ q H 20 ≤ 1. Thus, order of population that is compatible with Eq. (102) is given by If the order of the population of the final probability vector in Eq. (101) is given by Eq. (103), then non-zero ergotropy extraction is not possible due to passivity of final probability vector. Thus, amount of ergotropy that can be extracted from the final probability vector in Eq. (101) using thermal operation A 2 is given by Therefore, a condition for positive ergotropy extraction can be expressed as Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 5 : Thermal process A 5 transforms the initial probability vector given in Eq. (88) as follows From the where we write the last inequality using the fact 0 ≤ q H 10 ≤ 1. Therefore, the order of populations that are compatible with Eq. (107) is given by If the order of the population for the final probability vector in Eq. (106) is given by Eq. (108), then no ergotropy extraction is possible since final probability vector is already passive. A non-trivial ergotropy extraction is possible if the ordering is given by Eq. (109) and Eq. (110). Therefore, we calculate the amount of extractable ergotropy as Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 9 : Thermal process A 9 transforms the initial probability vector given in Eq. Since q H 10 ≥ q H 20 is true for any plausible values of β H , therefore the possible orders of population is given by If the ordering of the populations in the final probability vectors are given by Eq. (113), then no ergotropy extraction is possible since initial state is already passive. Thus non-trivial ergotropy extraction is possible if the ordering of the population in the final probability vector in Eq. (112) is given by Eq. (114) and Eq. (115). Thus, the amount of extractable ergotropy is given by Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 12 : Thermal process A 12 transforms the initial probability vector given in Eq.
where the last inequality follows from the fact q H 10 ≥ q H 20 . Hence the ordering of the populations of the final probability vector Eq. (117) that are compatible with Eq. (118) is given by If the ordering of population of the final probability vector in Eq. (117) is given by Eq. (119), then it is not possible to extract ergotropy from final probability vector due to passivity. Therefore, non-trivial ergotropy extraction is possible only if the ordering of populations in final probability vector is given by Eq. (120) or Eq. (121). Thus, the amount of ergotropy that can be extracted using thermal operation A 12 is Calculation of the amount of transferred heat from Hot bath and extractable ergotropy from p S using thermal process A 13 : Thermal process A 13 transforms the initial probability vector given in Eq.
(88) as follows From the table 17, we see that β-ordering of the final probability vector is given by (321). It is straightforward to see that all ordering of the populations in the final probability vector given in Eq.
(123) is compatible with this β-order. Thus ergotropy can be calculated as R A13 (p S ) = max 0,