Optimal cycles for low-dissipation heat engines

We consider the optimization of a finite-time Carnot engine characterized by small dissipations. We bound the power with a simple inequality and show that the optimal strategy is to perform small cycles around a given working point, which can be thus chosen optimally. Remarkably, this optimal point is independent of the figure of merit combining power and efficiency that is being maximized. Furthermore, for a general class of dynamics the power output becomes proportional to the heat capacity of the working substance. Since the heat capacity can scale supra-extensively with the number of constituents of the engine, this enables us to design optimal many-body Carnot engines reaching maximum efficiency at finite power per constituent in the thermodynamic limit.

The Carnot engine has a pivotal role in thermodynamics, both from a fundamental and applied perspective, being the reference point for other engines in terms of efficiency [1,2]. It is thus of paramount importance to understand its limits and strategies for its best utilization. In this article, we consider the optimization of a finite-time Carnot cycle within the so called low-dissipation (LD) regime [3][4][5][6][7][8][9][10][11][12][13][14], where the driving of the control parameters is slow but finite. Previous studies of Carnot engines in the LD regime have considered bounds on the reachable efficiencies [3], tradeoffs between efficiency and power [7][8][9]15], the coefficient of performance of refrigerators [12,13], the impact of the spectral density of the thermal baths [14], and other thermodynamic figures of merit [10,11]. Despite this remarkable progress, the following crucial question has remained unaddressed: given a certain level of control on the working substance (e.g. some parameters of the Hamiltonian, or some macroscopic variables such as volume or pressure), what is the optimal cyclic modulation of the control parameters to maximise the power output (or, more generally, any figure of merit involving power and efficiency [7][8][9]15]) of a finite-time Carnot engine? Such an optimal cycle has been designed for a single-qubit engine in [16,17], but a general understanding is lacking.
Using recent insights on a geometrical approach to thermodynamics [18][19][20][21][22][23][24][25][26][27][28][29], we show that, given any reasonable figure of merit involving power and efficiency, the optimal control strategy is to perform infinitesimal cycles around a fixed point. Furthermore, when the thermalization of the relevant quantities can be described by a single time-scale τ eq (see details below), the optimal power output becomes proportional to C/τ eq , where C is the heat capacity of the working substance (WS). Hence, the optimization of the heat engine cycle becomes intimately related to the maximization of C of the WS (interestingly, maximizing C is also crucial in thermometry [30][31][32][33]). * paolo.abiuso@icfo.eu We then use these insights to design many-body heat engines that can operate at Carnot efficiency with finite power per constituent of the WS through a supraextensive scaling of C/τ eq (e.g. in a phase transition), in the spirit of [34,35] (see also [35][36][37][38]). Despite differences w.r.t. previous proposals, which were based on Otto engines [34,35], we find the same asymptotic scalings for performance, while proving by construction their optimality in the slow-driving regime. Other recent proposals towards the possibility of reaching Carnot efficiency at finite power have been developed in [34][35][36][39][40][41][42][43][44][45] (see also [15,46,47] for no-go results [48]).
Finite-time Carnot cycle. We consider a quantum working substance (WS) with a driven Hamiltonian H(t), which interacts alternatively with a cold (B c ) and a hot (B h ) heat bath at temperature T c and T h , respectively (the results presented in this article can be naturally extended to classical systems). The Carnot cycle consists of four steps: (i) While being coupled to B c , H(t) is modified continuously from H(0) = H (X) to H(τ − c ) = H (Y ) for a time τ c . (ii) With the system isolated from the reservoirs, an adiabatic process is performed taking H (Y ) → H (Y ) T h /T c . During this process H(t) satisfies [H(t), H(t )] = 0 ∀t, t and commutes with the initial state (a Gibbs state w.r.t H(τ − c )), so it is possible to perform it arbitrarily quickly without affecting the state (hence corresponding to Hamiltonian quenches). (iii) While being coupled to B h , H(t) is modified back from H(τ + c ) = H (Y ) T h /T c to H(τ c + τ h ) = H (X) T h /T c in a time τ h . (iv) Finally a second quench is performed to restore H (X) T h /T c → H (X) , closing the cycle.
It is convenient to introduce the adimensional Hamiltonian G(t) := βH(t), where β = 1/k B T is the inverse temperature of the bath that the WS is coupled to (we shall set the Boltzmann constant k B equal to k B = 1). Note that the Carnot cycle becomes smooth with respect to G(t), and in what follows, we take the driving protocol to be time-reversal symmetric, more precisely that G(sτ c ) = G(τ c + τ h (1 − s)) with s ∈ [0, 1]. This property is always satisfied by optimal heat engines in the LD regime as long as the arXiv:1907.02939v3 [quant-ph] 5 Mar 2020 two baths have the same spectral density [14], which is the subject of interest of this work. By expressing G(t) as where λ j (s) are the control parameters and X j the conjugated forces, the cycle control can be characterised by τ c , τ h and its shape λ(s) (notice that τ c = τ h in general; the time-reversal symmetry is intended in the adimensional unit s). We will not write explicitly the dependence on s, which will be clear from the context, but will indicate with a dot the time derivative w.r.t. s, i.e.Ġ ≡ ∂ ∂s G ≡ τ x d dt G, x = (h, c). We now assume the slow driving or low dissipation approximation, which is crucial in this work. That is, d dt G is finite but small, so that we can expand the relevant quantities and keep only leading terms (formally, the small parameter of the expansion is τ eq /τ x , where τ eq is the relaxation time of the dynamics). In this regime, the state of the WS is always close to thermal equilibrium, and the heat exchanged during the isotherms (steps (i) and (iii)) can be divided as where ∆S x is the reversible contribution obtained in the quasistatic limit τ x → ∞, which is given by ∆S ≡ ∆S h = −∆S c = S(ω(τ c )) − S(ω(0)) with ω(t) = e −G(t) /Tr(e −G(t) ) and where S(ρ) is the Von Neumann entropy. The irreversible term W diss x can be described in this regime by the so called thermodynamic length [20][21][22][23][24][25]28] W diss where m ij is given, when the driven observables X j relax to equilibrium with the same time-scale τ eq , by [20][21][22][23][24][25]28]: where Z = Tr(e −G ) is the partition function. Given the time-reversal symmetry of the driving protocol, it follows that Σ h = Σ c (symmetric low-dissipation regime SLD). Importantly, while the results presented in the main text are based upon the standard thermodynamic metric (4), generalizations (including τ eq depending on the trajectory [24,25], the possibility of having several relaxation time-scales, general Lindbladian dynamics [24,28], and protocols in which Σ h = Σ c ) are developed in the Supp. Mat. [49].
Optimisation of the cycle. We now optimise the power (and efficiency) of the Carnot engine over τ c , τ h and its shape λ j (s), which are all the possible degrees of freedom. This enables us to obtain a fundamental upper bound on the power in the slow-driving regime, as well as the corresponding optimal control.
The work extracted during a cycle is given by W = Q h + Q c , and the total time is τ = τ c + τ h . The power hence reads P = (Q h + Q c )/τ , and the efficiency η = (Q h +Q c )/Q h . By substituting the expressions (2)- (3) and appropriately setting τ c and τ h , one can maximize the power of the engine (∂P/∂τ j = 0) obtaining [3,50] and the corresponding efficiency at maximum power (EMP) is given by the Curzon-Ahlborn EMP, η CA = 1 − T c /T h [51]. In the most general case one might seek, in order to not sacrifice completely the efficiency optimization over the power, to maximize a hybrid figure of merit [7][8][9]15]. The maximum efficiency for any given power output of the engine has been derived in [9] (see also [7]). Analogously, we can express the best power for a given efficiency, fixed to be a fraction of the maximum one: where η C = 1 − T c /T h is the Carnot efficiency. In the SLD regime this leads to a maximum power (cf. Supp.Mat. [49]), obtained by setting Essentially, by tuning γ in τ c and τ h , one can interpolate between a maximally powerful engine with power (5) at efficiency η CA , and a Carnot engine with maximal efficiency and null power.
At this point, we note a crucial observation: after the optimisation of P over τ c and τ h , the remaining figure of merit is always (∆S) 2 /Σ, independently of the value of γ. In fact, this is a property that can be argued to be general given any figure of merit combining power and efficiency [49]. We now show how to maximize it by an opportune use of a Cauchy-Schwarz inequality. First, using the formula for the derivative of an exponential [52]: ∂e −G /∂λ j = − 1 0 e −G X j e −(1−s)G , we can express again Σ in (3) in the more compact form where cov ω (A, B) is the Kubo-Mori-Boguliobov inner product: cov ω (A, B) = 1 0 ds Tr(Aω 1−s (B − Tr(ωB))ω s ). Next, we note that the Von Neumann entropy S = −Tr[ω ln ω] satisfies:Ṡ = −Tr(ω ln ω) = −cov ω (ln ω,Ġ) = cov ω (G,Ġ) where we used again the formulaω = − 1 0 dx ω 1−x (Ġ − Tr(ωĠ))ω x , as well as Tr(ω) = 0 and cov ω (A, 1) = 0. Hence, we can write the total change in entropy ∆S as: Crucially, both ∆S and Σ can be expressed through an infinite-dimensional scalar product given by: , that depends on the path {λ j (s)}. Using the Cauchy-Schwarz (C-S) inequality | A, B | 2 ≤ A, A B, B , the ratio (∆S) 2 /Σ can then be bounded as: where C is the heat capacity of the WS, i.e. the variance of the adimensional Hamiltonian G for its thermal state. The C-S inequality is saturated forĠ ∝ G, which means that creating quantum coherence cannot favour the power output in the slow driving regime, in agreement with Refs. [53,54]. More importantly, it shows that optimal thermodynamic protocols take the simple form G(t) = λ(t)G(0). We can further maximize (63) by noting that 1 0 ds C ≤ max s C. To saturate this inequality in practice, one needs to consider cycles where the modulation of G is small around an optimal point G(0) where C/τ eq is maximised (as long as is small enough so that G(sτ c ) with s ∈ [0, 1] does not change substantially along the cycle, the precise form of G(sτ c ) is not important; see the SM for examples of explicit cycles). In this case, in the limit → 0 the maximal power (7) for a given efficiency γη C is given by where C is the heat capacity at G(0). We stress that (13) has been obtained after maximising P (for a fixed η) over all degrees of freedom: τ c , τ h and the protocol {λ j (s)}. This result shows a fundamental link between maximal power of a finite-time Carnot cycle and the heat capacity of the WS. The simplicity of (13) can be contrasted to exact optimisations of heat engines [45,[55][56][57][58], where the full solution easily becomes too complex or not even computable with the size of the WS; and with other geometric optimisations which require solving geodesic equations to design optimal paths in the parameter space [24,26,28,29,[59][60][61]. In our approach, from the point of view of optimization all that is left to do in (13) is to maximise C over the control parameters {λ j } to identify the optimal working point G(0) in (12). In Fig. 6 explicit results are reported for the value of maximum C, for different paradigmatic levels of control on the same system of N qubits.
Crucially, this approach can be generalized to any metric m ij in (3) describing dissipation: in the SM [49] we show that the optimal control problem is always reduced to infinitesimal cycles, and the optimal working point can be found by a scalar maximisation problem. We work out as well examples of standard microscopical dynamics in open quantum systems: a 2) For an Ising chain Given complete control over the spectrum with 2 N levels, Cmax N 2 /4. Details can be found in [49].
qubit, a qutrit, or an harmonic oscillator as a WS in contact with bosonic thermal baths with different spectral densities.
It is important to point out that (13) should be understood as a theoretical ultimate upper bound on power, obtained by taking → 0 in (12). In practice, any experimental or realistic protocol will have finite . In this case, as long as C is sufficiently smooth along the thermodynamic cycle (12), the power output P (given by the integrated C in (63)) will not change considerably, so that realistic cycles will provide a similar P than the theoretically maximal one. In practice, keeping finite is also important to ensure the consistency of the slow-driving approximation τ eq /τ 1, given that for the optimized protocols that lead to the power (7), one has τ ∝ , more precisely Note that this can always be guaranteed for high efficiencies γ → 1. From the same equation we note incidentally that engines whose maximum efficiency is constrained to be low (T c /T h 1), i.e. arguably those engines that mostly need optimization in the high efficiency regime, show better convergence to the absolute bound.
Reaching Carnot efficiency at finite power. As an example of application of the previous results, we now use the designed optimal finite-time Carnot cycles (Eqs. (12) where G(0) will depend on each model of interest) to explore the possibility of reaching Carnot efficiency η C at finite power in the macroscopic limit. We follow the approach put forward in Refs. [34,35]: considering a N -particle WS, we aim at approaching Carnot efficiency in the macroscopic limit N → ∞ without giving up power per constituent.
To reach Carnot efficiency, we need γ = 1 in (6), and hence we take 1−γ = N −ξ , where ξ > 0 can be chosen at will. On the other hand, the maximal power P (max) γ in (13) depends only on C and τ eq ; we then assume C = c 0 N 1+a and τ eq = τ 0 N b , where the meaning of the different constants will later be described for each model of interest. Expanding the relevant quantities for N 1, we obtain at leading order in N : where σ 2 w = w 2 − ( w ) 2 is the variance of the work distribution, which measures the work fluctuations per cycle of the engine (see Supp.Mat. [49] for details on the calculation). Let us now discuss two separate cases, inspired by [34] and [35], respectively.
(a) Control on the engine and the enginebath interaction. We first assume full control over the engine Hamiltonian with 2 N levels: that is, all levels can be modified at will by the experimentalist. While this is extremelly challenging in practice, it is useful to obtain fundamental upper bounds on the maximal power. The optimal Hamiltonian maximising C then consists in a ground level and a 2 N − 1 degenerate level (see [30,62]) which, as shown in Fig.6, leads to C ∝ N 2 , i.e., a = 1. This supralinear scaling is obtained in an increasingly smaller region of the parameter space, which requires taking ∝ 1/N in (12) (see details in Supp. Mat. [49]), and this constraints from Eq. (14) also 1−γ to scale accordingly, i.e. ξ = 1. Furthermore, it is possible to reach in realistic collisional scenarios τ eq ∝ √ N (i.e. b = 1/2), or constant τ (b = 0) if one is allowed to fine-tune the interaction between the WS and the baths [34] [63], in agreement with Ref. [34]. In the Supp.Mat. [49], we solve exactly this proposal for a feasible driving protocol close to the optimal one.
(b) Engine working on a phase transition point. A promising platform to obtain supralinear scaling of power with realistic control is by choosing the engine to work in a phase transition point of the many-body WS. For a finite WS operating close to the critical point, finite size scaling theory tells us that C develops a peak of height C ∝ N 1+α/(νd) and width δ ∝ N −1/(dν) , while τ eq ∝ N z/d (here α, ν and z correspond to the specific heat, correlation length and dynamical critical exponents, while d is the spatial dimension of the engine [64,65]). In order to exploit the critical scaling of the WS, we need to perform the cycle (12) where C becomes peaked, and hence ∝ δ ∝ N −1/(dν) , which implies ξ = 1/(dν) from Eq. (14). Then, from (15) with a = α/(νd) and b = z/d, supralinear scaling of P This condition is the same found for the Otto cycle proposed in [35]. Examples of physical systems where (16) is satisfied are also provided in [35], particularly in the presence of critical speed-ups of thermalisation (z < 0 [66][67][68]). Besides efficiency and power, another crucial aspect of a heat engine is its reliability, i.e. the fluctuations in the output power. In fact, it has been recently pointed out in [36] that the Otto-cycle of [35] suffers from macroscopic fluctuations in the thermodynamic limit. For the Carnot-cycle considered here, from (15) the relative work fluctuations read (15), one has f w ∼ O(1) in the macroscopic limit. A similar situation takes place for the critical heat engine (b) as one simultaneously has ∝ N −1/(dν) and a = α/(dν), and hence f w ∝ N (−2+dν+α)/(dν) . Using the relation dν = 2 − α [69], we hence obtain f w ∼ O(1). Therefore, for both proposals f w ∼ O(1) in the thermodynamic limit, hence hindering their reliability, which is the same result found in [36] for the Otto cycle. Despite f w ∼ O(1), these fluctuations can be suitably avoided when the number M of cycles is large (the argument below applies to the Otto and Carnot cycle). Given M cycles, we have that f w ∝ 1/ √ M as the average work W ∝ M whereas the work fluctuations σ w ∝ √ M (think of a biased random walk). Therefore, the ratio between the fluctuations per unit time and the power goes to zero as M grows even when the fluctuations per single cycle are large. Since we have that M ∝ τ tot /τ eq , for a total time τ tot of observation, fluctuations can be suitably avoided e.g. for critical speed-ups where τ eq ∝ N b with b < 0.
In actual implementations, the technical requirements to realise such optimised Carnot cycles are: global control of the WS, H(t) = λ(t)H(0), and enough precision to engineer small cycles in the region where C/τ eq has supralinear scaling with N . Since the width of this region shrinks with N , in a realistic implementation the supralinear scaling will be eventually lost as the control precision is limited [70]. We remark that even when the experimental control may be limited, our considerations provide upper bounds on the maximal power of finite-time Carnot engines.
Conclusions. We have characterized the optimal cycle of a finite-time Carnot engine in the lowdissipation regime. The dissipation has been characterized by the thermodynamic metric (4), which is justified when the thermalization of the working substance (WS) is well described by a single timescale τ eq . In this case, the optimal cycle turns out to be remarkably simple: it consists of modulations in the form λ(t)H(0), where H(0) is the Hamiltonian of the WS. The power output is then proportional to the heat capacity C of the WS, linking the optimal performance to the nature of the WS: as an application we showed how the critical scaling of C can enable the design of optimal engines with extensive power reaching Carnot efficiency. These results have been generalised to general metrics in the Supp. Mat. [49], which we have used to derive the optimal cycle and corresponding power output of different WS (qubit, 3level system, or harmonic oscillator) interacting with a bosonic thermal bath. Putting everything together, our reults provide a general framework to efficiently optimise the control of slowly driven Carnot engines.
In particular our results hint at the answer for two open problems: 1) small cycles are optimal for engine performance in all regimes [17,45,57], and 2) the performance of the proposals [34,35] cannot be improved.

SUPPLEMENTAL MATERIAL
This Supplemental Material (SM) contains a generalisation of the results presented in the main text to open quantum systems, as well as technical details. We start in Sec. I with a review the quantum Carnot cycle in the quasistatic limit. In Sec. II we introduce finite-time corrections for open quantum systems dynamics and show to characterise them geometrically. In Sec. III we generalise the main result (the optimal cycle and corresponding maximal power) for general open quantum system dynamics by giving the solution as the maximisation of a scalar function. In Sec. IV, we present the maximisation of the heat capacity discussed in the main text. In Sec. IV we discuss the technical details of the asymptotic expansions of power and efficiency presented in the main text. Finally, in Sec. VI we solve analytically an illustrative solvable case. To help the reader find information in this Supplemental Material we introduced Table I.
Section I Description of a quasistatic, reversible Carnot cycle for a general quantum system. Section II Introduction of finite-time Carnot cycle and thermodynamic length. IIA Thermodynamic metric for exponential relaxations of the observables. IIB Thermodynamic metric from microscopical models. IIB1 Example: Qubit with bosonic baths. IIB2 Example: Harmonic oscillator with bosonic baths. IIB3 Example: 3-level system with detailed balance evolution. Section III Generalizations of cycle optimization to general metrics and protocols. IIIA Generalization for asymmetric dissipation and symmetric protocols.

IIIA1
Optimization for systems with point-dependent thermalization timescale.

IIIA2
Optimization for general metrics. IIIA3 Solution of the optimization for the models presented in Section IIB1-2-3.

IIIB
Bounds on completely asymmetric protocols. Section IV Maximization of the heat capacity for systems with different degree of control. Section V Time-tuning optimization of a low-dissipation Carnot engine. VA Critical scalings of power, efficiency and fluctuations of optimal protocols. VB Comparison with Otto cycle of Ref. [36] Section VI Explicit analytical control protocol for a D − 1 degenerate model.

I. QUASISTATIC QUANTUM CARNOT CYCLE
For completeness, in this first section we review the quantum Carnot cycle in the quasistatic limit (i.e. without including finite-time corrections).
The internal energy of a system with Hamiltonian H in the state ρ is defined as Considering the variation dU , it is possible to identify [79][80][81][82] the work and heat contributions To simulate sensible restraints on the system, external control is assumed on dynamical parameters of the local Hamiltonian of the system H t = H( l(t)). When in contact with a reservoir at temperature T = 1/β (we use units in which the Boltzmann's constant is k B = 1), such a system relaxes to the the Gibbs state is the partition function of the system in the canonical ensemble). A Carnot Cycle [14,17,57,82] is identified with a 4 steps process, that is two isothermal strokes alternated with two isoentropic (adiabatic) strokes (cf. Fig. 2). Consider a system with a controlled Hamiltonian H t which can be coupled independently to two reservoirs with temperature T h > T c . In the ideal quasistatic limit the operations are performed slowly enough to allow the system to be in thermal equilibrium ρ(t) ≡ ω β ( l(t)) at every instant. The 4 steps are: 1) while being coupled to the cold reservoir, the Hamiltonian is modified continuously from H (X) to H (Y ) such that Tr[ω β H t ] is negative, in order for heat to be released to the cold source.
2) with the system isolated from the reservoirs, a quench is performed taking 3) while being coupled to the hot reservoir, the Hamiltonian is modified continuously from Note the factors T h Tc and Tc T h are chosen in order for the state to be continuous during the quenches. In fact the thermal state uniquely depends on βH T (cf. (20)); i.e. for example during the quench 2) the relation . We shall thus define then the adimensional Hamiltonian at temperature 1/β (note that the temperature takes only two values, that depend on the respective isotherms) G t := βH t τ (21) so that the thermal state is ω t = e −G t /Tr[e −G t ] on both the cold and hot isotherm. The time reparametrization is conceived in order to isolate the shape of the control G t , with 0 ≤ t ≤ 1. Note that G (and hence ω) is continuous also on the quenches; that is, we can consider the cold isotherm consisting in a transformation (ω (X) , G (X) ) → (ω (Y ) , G (Y ) ) with duration τ c and the hot isotherm the opposite (ω (Y ) , G (Y ) ) → (ω (X) , G (X) ), with duration τ h . In a Carnot cycle heat is exchanged only during the 1), 3) steps (absorbed from the hot source, released to the cold one), hence we can compute the efficiency using the observation that over a cycle ∆Q + ∆W = 0 and the heat is absorbed from the hot bath Q (0) h (we use the superscript (0) to indicate quantities in the quasistatic regime) where we have used in the last step that in the quasistatic limit the above mentioned considerations together with Eq. (19) give

II. FINITE-TIME CARNOT CYCLE AND THERMODYNAMIC LENGTH
We now consider finite-time corrections on the above quasistatic Carnot cycle through a quantum open system approach, where the dissipation in linear response can be described in a geometric form by following [28]. We consider an isothermal process where the Hamiltonian H(t) of the WS is driven for t ∈ [0, τ ] in contact with a thermal bath at temperature T = 1/β. In order to characterise the process beyond the quasistatic limit, we need to assume some structure on the thermalization processes of the working substance (WS) induced by the reservoirs. In a rather generic scenario we consider that the relaxation of the WS can be described by a master equation with the thermal state as a unique fixed point: where we note that L ≡ L(t) is also time-dependent. Following [14], the solution of this master equation can be found perturbatively around 1/τ . Using the renormalised time s ∈ [0, 1] with the convection H s ≡ H(sτ ) and ρ s ≡ ρ(sτ ), the solution is given by (indicating with a dot˙the time derivative w.r.t. s): with ω s = e −βHs /Tr(e −βHs ), and where L −1 s is the inverse of L s within the traceless subspace of density matrices, the so called Drazin inverse (see e.g. [14,17,28,92]). Plugging this expression into Q = 1 0 ds Tr(ρH) and using integration by parts we obtain: with the first order correction to the quasistatic limit given by: In terms of the adimensional Hamiltonian G s ≡ βH s , and using the formula for the derivative of an exponential, we can write the variation of entropy as while Σ can be reexpressed in the convenient quadratic form: Now, expanding G s as we can conveniently write ∆S and Σ as with The matrix m ij is symmetric, positive-definite due to the second law dΣ ≥ 0, and it depends smoothly on the base point ω; hence it defines a metric. Finally, we note that the linear expansion (32) can also be obtained in exact treatments of the system-bath Hamiltonian dynamics through linear-response theory [25,[83][84][85]. As a rule, the expansion (32) can be argued to be general for any system with dissipations that are linear (at the lowest order) in the speed of the driving. Suppose indeed that the dissipation along an infinitesimal segment of the trajectory depends only on the local point and the local driving. As a consequence it must be in the form dΣ = dλ i f i (λ,λ,λ, ...), but the 1/τ scaling implies that the first derivative terms enter linearly in the product while higher orders are suppressed, which implies dΣ = dλ i m ij (λ)λ j , which is equivalent to (32).
A. Thermodynamic metric for single or multiple time scales We first show how to obtain the Kubo-Mori-Boguliobov metric used in the main text (Eqs. (4) and (8)). This can be easily done by taking an exponential relaxation of ρ to equilibrium with a single time-scale, as described by the Lindbladian: which has the Drazin inverse L −1 s [.] = τ eq (ω s Tr(.) − I). In this case, using that Tr(J ω [Ġ]) = 0 one finds that Σ in (30) is given by: that is Eq. It is important to keep in mind that these considerations are only relevant close to equilibrium, where the metric (4) in the main text becomes a good thermodynamic description. Furthermore, although we have assumed that the whole state ρ converges to equilibrium as in (34), strictly speaking it is only necessary that the driven observables (the X j in (31)) converge to equilibrium with a single time-scale, i.e., Ẋ i ρ = τ −1 eq ( X i ω − X i ρ ), with X ρ = Tr(Xρ). This is especially relevant in complex systems (e.g. many-body systems), where the full dissipative dynamics can be extremely complex but the equilibration of some macroscopic observables can be well described by an exponential relaxation with a suitable time-scale. In this sense, it is also worth pointing out that if each generalised observable X i decays with a different time-scale τ i then the metric Eq. (4) of the main text can be easily generalised as [28]: where we have absorbed the dependence on τ i in the metric and where τ i can in principle depend on the point of the trajectory. As a final remark, we note that while we have derived the metrics (35) and (37) with an heuristic approach based on exponential relaxation near equilibrium, one can derive them from a microscopic derivation based upon linear-response [24,25], in which case τ i will in general depend on the point of the trajectory, i.e., G. This case will be treated in Sec. III A 2.

B. Thermodynamic metric on standard microscopical models
In this section, we apply our general considerations to derive thermodynamic metrics for systems described by quantum master equations. Before, let us discuss here some generic properties. Following [14,86], we have in standard scenarios where a quantum system is coupled to a bosonic bath: where γ 0 ν α is the spectral density of the bath (α defines the Ohmicity), and is the Bose-Einstein distribution (one can consider ferimonic baths by replacing N (βν) by the Fermi distribution). The general form (38) of L depends on H s , β and the spectral density J(ν) = γ 0 ν α . Consider now a Carnot cycle with two baths at temperature β c and β h , and the action of each bath being described by a L(H S , β, J(ν)). Assuming that the spectral density of both baths is the same, then both Lindbladians are related by the transformation β → λ −1 β and H S → λH S , with β ≡ β c and λ ≡ β c /β h . In terms of the Lindbladian, note that [14]: which shows how the generators of the dynamics are related between the cold and hot isotherm, i.e.
After noticing that the dissipation is related to L through Eq. (30), we can write a simple proportionality relation for the metric describing the dissipation on the hot and cold isotherm, which implies, for time-reversal symmetric protocols described in the main text, While from this relation we see that the considerations of the main text (for symmetric protocols Σ c = Σ h ) in principle only apply for flat spectral densities α = 0, we will show in Sec. III A that the same figure of merit ∆S 2 /Σ c (or equivalently ∆S 2 /Σ h ) needs to be maximised to obtain maximal power. An important comment is now in order. Whereas the dissipator (38) is usually derived for time-independent Hamiltonians [86], here we are interested in slowly driven Hamiltonians. Nevertheless, the same form for the dissipator is justified as long as the bath dynamics are fast compared to the driving rate of the system Hamiltonian, which leads to the well-known adiabatic master equation [87][88][89]].

Two-level system with a bosonic bath
The well-known optical master equation [86] can describe a two-level system with Hamiltonian G ≡ βH = w 2 σ z relaxing in a bosonic thermal bath at temperature 1/β, and from (38) and working in the interaction picture it takes the formρ where the rate Γ ∝ w α depends on the ohmicity of the bath (for a standard atomic-optical field interaction α = 3 [86]) and σ ± = (σ x ± iσ y )/2. It is easy to translate this equation on the single Bloch-vector components This equations being in the form (36), imply that the thermodynamic metric is indeed in the form (37). Moreover the covariances of different cartesian components decouple as cov(σ i , σ j ) ∝ δ ij (with i, j, = x, y, z) hence implying Consider now the external control on the qubit, i.e. (n(t) is a unit vector) Assuming that the driving is much slower than the internal dynamics of the bath (i.e. the adiabatic master equation [87][88][89]), the same master equation (45) and metric (47) instantaneously holds in the rotated basis where σ z (t) =n(t) · σ and with ω → ω(t) becoming time dependent. The same metric is expressed in polar coordinates in [28].

Harmonic oscillator with a bosonic bath
As in the previous example, the optical master equation can be written similarly for an harmonic oscillator with gap control βH(t) = w(t)a † a (again working on the interaction picture and assuming an adiabatic master equation)ρ Considering a modulation of the level splitting w, we have for the average occupation number Thus the metric for the single control parameter w(t) is in this case simply Note that a frequency-controlled harmonic oscillator can be realized in a different set up, i.e.
(to see how this control is not equivalent to (49) it is sufficient to notice [Ḣ, H] = 0 in general). The metric arising with this control is being worked out in [90].

Three-level system with a bosonic bath
For a final example we consider here a three-level system satisfying a detailed balance master equation. We assume control on each of the energy levels, while keeping the basis fixed (this is motivated by the observation that creating coherence does not increase power in the linear-response regime [78]). Without loss of generality we can assume the ground state energy to be zero (βE 0 ≡ 0) and the two excited states with energies β −1 E 1 (t) ≤ β −1 E 2 (t), hence the thermal state being characterized by a probability vector The evolution of the probability vector is described by the Markovian master equation [86] p The rate matrix Γ ij can be obtained from (38), yielding where Γ is a rate, N (w) = 1/(e w −1) for bosonic baths, and α defines the spectral density of the bath. Note that these dynamics cannot be described by the simple exponential decay (36), so that we need to use the more general approach (33) to find the metric. Indeed, we use Eq.
whereΓ −1 ij is the Drazin inverse of Γ ij (see e.g. [14,17,28,92] for details on the Drazin inverse) while ∂ j ω k is the variation of ω k due toĖ j , i.e. ∂ j ω k = −ω k δ jk + ω j ω k . The analytic form of m ij can be computed from this expression even for larger systems with more than 3 levels using symbolic computation software.

III. GENERALISATION OF THE MAIN RESULT: OPTIMAL CYCLES FOR GENERAL METRICS AND ASYMMETRIC DISSIPATIONS
We now have the necessary tools to generalise the optimisation to more generic metrics and thermodynamic processes. In the most general case, we can consider the Carnot cycle where the cold isotherm is characterised by the control parameters {λ ij are simply proportional to each other (this class includes e.g. all the standard microscopic scenarios with same spectral density of the baths, as showed in II B), the time-reversal property is not an assumption but simply follows from the optimization process [14]. This implies, as we show in the following III A, that the relevant figure of merit to be maximized is still ∆S 2 /Σ. Therefore, we show how to accomplish this task for different structures of the metric m h,c ij , up to the most general case where the dependence on the temperature, spectral density, etc. is encoded in L −1 in the expression (33).
Finally in Section III B we will consider generalizations to cases where m

A. Simple asymmetric dissipations
We consider here cases in which the dissipations along the hot isotherm are proportional to the dissipations along the cold isotherm via a constant factor that does not depend on the specific control protocol: where σ is a number. Importantly this class of cases includes the scenario where the two baths have the same (non-flat) spectral density, as explained in Sec. II B and in [14], where we obtained that m (c) ij ) is the metric associated to the dissipation of the cold (hot) bath), which implies Σ c = Starting from (57), the power to maximize is given by under the efficiency constraint As the only time unit is Σ/∆S ≡ τ * and the temperature ratio r = Tc T h the problem can be rephrased as where x j = τ j /τ * . Here the maximization with its constraint is expressed in full adimensional terms, meaning that after solving it the resulting power will be in the form T h (∆S) 2 Σ f (γ, r, σ) for some scalar function f . The general form of the function f can be obtained analytically but it is in general non-trivial and very lengthy (so we do not show it here), but in the limit of high efficiency it simplifies to At this point the optimisation of the power boils down to the maximisation of (∆S) 2 /Σ, as in the main text,. In the following we show how to maximise ∆S 2 /Σ in general scenarios, and prove that asymptotically infinitesimal cycles are optimal. Before, let us note that given any reasonable figure of merit between efficiency and power, the relevant figure of merit will always be (∆S) 2 /Σ. Given as objective any figure of merit f (η, P ) expressed in terms of the efficiency and the power, after optimization on τ h , τ c the resulting valuef := f (η,P ) (we indicate withq uantities after time optimization) will be given, byη =η(T h , T C , σ) andP =P (T h , T c , σ)(∆S) 2 /Σ for some adimensional functionη and homogeneous functionP of the temperatures, by dimensional analysis. Moreover, any reasonable f will be monotonously increasing in both its parameters singularly (i.e. for a given efficiency we wish to enhance the power and vice-versa), therefore after speed optimization it will still be possible to improve the engine performance by increasing the ratio (∆S) 2 /Σ. The role of the above mentioned ratio as a characteristic scale defining the performance of both engines and refrigerators was already noticed in [6] without further analysis.

Dependence of the time scale in G
As a simple first extension of our results, we consider the case of the standard thermodynamic metric with a time-scale that can depend on the point of the parameter space G = βH. In other words, we consider the heuristic model of exponential thermalisation (34) with a G-dependent τ eq (see also [25] for a microscopic derivation). Following the reasoning of the main text, we can define an infinite-dimensional scalar product given by with ω s = e −Gs /Tr(e −Gs ), and where we note that the scalar product is defined for a fixed trajectory G s with s ∈ [0, 1]. We then have, so we conclude that the results of the main text can be naturally extended to time-scales that depend on βH, meaning that the performance is upper-bounded by small cycles with proportional modulation of the HamiltonianḢ ∝ H performed on the point where the ratio between heat capacity and relaxation time of the system is maximum.

General metrics
Let us consider the maximisation of ∆S 2 /Σ for more general protocols where the thermodynamic metric is given by the general expression (33) which depends on the particular Lindbladian describing the thermalisation dynamics. Following the considerations of the main text, to maximise (∆S) 2 /Σ, we can expand them as In contrast to the main text, where m ij is given by Eq. (4) in the main text, here we a general metric of the form (33). The Cauchy-Schwarz inequality can be applied by considering the following: Inequality. Consider two vectors a, b and a quadratic invertible form g > 0 defined on their vector space. Then the standard C-S inequality applied to g 1/2 a and g −1/2 b tells If we now consider this inequality applied to two vectors and the metric where m(t) is a positive time-dependent quadratic form in R k , we have and thus it is possible to write the C-S inequality as This inequality allows us to bound the figure of merit for the maximum power as where we recall that s i = Tr[(Tr[Gω]ω − Gω)X i ], and the X i are given by the expansion (64). Hence, the initial optimisation over all possible trajectories {λ j (s)} boils down to maximising a single real function given the control parameters. In other words, our result shows that in order to design optimal finite-time Carnot cycles one does not need to optimise over all trajectories {λ j (s)} with s ∈ [0, 1], but it is enough to search a suitable point in the thermodynamic space {λ j }. To saturate it in practice, one needs to maximise sm −1 s over the control parameters λ, and consider infinitesimal variations around this optimal point: where we have used the more compact vector notation: s = {s 1 , s 2 , ...}, etc. The direction of the modulations is defined, from the Cauchy-Schwarz saturation condition, by the vector m −1 s, while the normalization of the modulation is in principle to be taken infinitesimal. The abstract expression for µ corresponds, in the case of the analysed in the main text, to λ * itself, i.e. the modulations of the control parameters are proportional to the control themselves in that case.

Examples
a. Qubit with optical master equation. We first consider here a slowly driven qubit with full hamiltonian control (here G is the adimensional Hamiltonian in temperature energy units) in contact with a bosonic bath, so that the evolution is described by the adiabatic master equation (in the interaction picture) Here, σ ± ≡ σ ± (t) = σ x (t) ± σ y (t) in the basis where σ z (t) =n(t) · σ, and ω ≡ ω(t), as in Section II B 1. Given the metric (47) and considering that the variation of entropy is non-zero only along the instantaneous z direction n(t) it is evident how, to maximize ∆S/Σ 2 , changing the direction ofn is useless (as noted in [28], generation of coherences is detrimental, and the eigenvectors of m ij (47) are larger along the tangential directions x and y), hence the optimal trajectories can be recognized to be in the form of simple gap modulation w(t), reducing the metric to Therefore we fall into the category of systems analysed in Section III A 1, for which the optimal control consists in modulations on the point where the ratio between the heat capacity Var ω (G) and the relaxation timescale is maximum. In this case τ eq ≡ Γ −1 (2N (w) + 1) −1 and C Q = w 2 e −w /(1 + e −w ) − (we −w /(1 + e −w )) 2 . In Fig. 3 we show the maximization of C Q Γ(2N (w) + 1), corresponding to the maximum value of (∆S) 2 /Σ in this case, choosing Γ constant (i.e. flat spectral density, α = 0). Allowing α = 0 would not change significantly the difficulty of the maximization problem, here α = 0 has been chosen to allow a comparison with the case studied in the main text, where the thermalization timescale does not depend on the value of the Hamiltonian. To sum up, by making use of equation (61), the full maximization of power for a qubit in contact with bosonic thermal baths, leads to where the maximization result depends on α.  b. Harmonic oscillator The above analysis for the qubit can be replicated thoroughly for the controlled harmonic oscillator described in II B 2, and the final result where heat capacity of the harmonic oscillator is equal to C HO = (2 sinh(w/2)) −2 . In Fig. 4 an example is shown for α = 3. c. 3-level system with master equation Finally we consider a 3-level system as described in II B 3, controlled by the energy levels of the two excited states E 1 and E 2 , and whose evolution is given by a standard detailed balanced master equation (54)- (55). By following the general approach described in III A 2, we bound ∆S 2 /Σ by ij s i m −1 ij s j , see (71). We symbolically compute the metric (56) and the vector {s i } (see details in Sec. II B 3) in order to express ij s i m −1 ij s j as a function of {E 1 , E 2 }, and then maximise it over {E 1 , E 2 } numerically. An example is shown in Fig. 5, where we plot ij s i m −1 ij s j for a linear spectral density α = 1. Note that the maximum is obtained for E 1 = E 2 . We also find that in this case the vector µ (73) is proportional to 1 1 on the bisector, meaning the modulation on the optimal working point are symmetric on E 1 and E 2 .
The maximum power results in such a case Furthermore, in Fig. 4 (right), we also compare ij s i m −1 ij s j with the heat capacity C divided by the characteristic timescale τ eq = Γ −1 w −1 (but remember the thermalization process is not simply exponential in this case), and it can be seen how this ratio provides a rather good approximation of the maximal power. That is, while C/τ eq is not exactly the correct figure of merit in this case, it can be used as a relevant figure of merit to determine where the optimal power is. This example also illustrates that our approach can be used to deal with more complex dissipative dynamics than the standard exponential decays discussed in the main text.

B. General asymmetric protocols
Finally, we consider the most general setting (time-reversal asymmetric protocol, different baths, generic dependence on the temperature of L): Σ c = Σ h . Then, the maximum power expression after isotherms-duration tuning (with no efficiency constraints) is given by [9] which can be thus written as from Eq.s (30)-(33) where now two different metrics m (j) appear in the denominator. In this case we lose the simple structure (scalar product) 2 /(quadratic form), hence it is not possible to apply directly the C-S inequality; nevertheless it is possible to upper-bound Eq. (81) using simple inequalities e.g.
These can be now optimized using C-S as we did in the previous section, and depending on the relative size (in the interval 1 3 ≤ a b ≤ 3 the former is tighter, otherwise the latter) they will give a bound (not tight in general). Note that while inequalities (82)(83) are useful to give upperbounds to the maximum power theoretically obtainable, in practical terms it is also possible to give lower bounds, which are useful to certify that it is possible to reach at least a given value of the power, and maximizing it by the same methods we used in the main work. Namely the bounds can be written As a final observation, we note how all of the above can be applied also to the power in the high efficiency regime, which is (cf. Eq. (61))

IV. MAXIMIZING THE HEAT CAPACITY VS. DEGREE OF CONTROL
The optimal cycles described in the main text (induced by the standard thermodynamic metric (4)) make the maximal power proportional to the heat capacity (13) of the chosen working point of the infinitesimal Carnot cycle. The task of maximizing the power is then translated in finding the point in the control parameter space that maximizes the heat capacity, i.e. the variance of the adimensional Hamiltonian G for a thermal state ω = e −G /Tr[e − G]. We show in the following 3 paradigmatic cases that are summed up in Fig.6 showing how different degrees of control offer different possible performance on the same system, which is made of N qubits.
(a) Full control over the spectrum. We first assume total control on the Hamiltonian of the N qubits, which is made up of D = 2 N energy levels. That is, any desired (possibly long-range) interaction can be engineered, so that the D-level spectrum can be controlled at will. While this might be extremely challenging to realise in practice, it is useful to consider this situation to obtain a fundamental upper bound on the maximal power. Indeed, the maximization of heat capacity C of a general D-dimensional system at thermal equilibrium has been carried out in [30,62]. The optimal Hamiltonian consists of a ground level and a D − 1 degenerate level, with an optimal gap x in adimensional units (i.e. rescaled by the temperature) defined by e x = (D − 1)(x + 2)/(x − 2) and the corresponding C is C max = x 2 e x (D − 1)/(D − 1 + e x ) 2 . This expression gives in the asymptotic regime (D → ∞) x ln D, hence C max (ln D) 2 /4; which in terms of the particle number N ∝ ln D means C max N 2 /4 for N 1, i.e., a quadratic scaling. On the robustness of result (a). Suppose the control of the optimal energy gap x has some imprecision (or simply needs to be modulated to perform the cycle). We show here that as far as the modulation/imprecision doesn't scale with the dimension, the heat capacity behaves smoothly. Suppose indeed that x = (1+ ) ln(D −1). Then the adimensional variance of the flat Hamiltonian with d ≡ D − 1 degenerate excited states and gap x is which, to keep the scaling as (ln d) 2 ∼ N 2 needs the rest to stay finite, i.e. to scale as 1/ ln d = 1/N . (b) N independent qubits. As another extreme case, corresponding to almost no control on the spectrum, we can consider N independent qubits, that is an Hamiltonian H (N ) = N i=1 λ i (t)σ z i . Note that given the qubits do not interact the thermal heat capacity will result additive, and the optimal gap will be the same for all the qubits λ i = λ j . This means it is enough to solve the previous case for D = 2 (i.e. N = 1), to find C max ≈ 0.439N (with an optimal gap λ * i ≈ 2.40), in agreement with [17]. (c) Ising chain. Finally, we consider an Ising chain ( , and assume control over λ 1,2 (t). Numerical results in Fig. 6 show how the interactions allow for substantially improving on (b), although asymptotically we obtain a linear scaling in N , C max 0.59N for N 1. This linear scaling is in fact expected for systems away of a phase transition point (e.g., a rigorous proof that C is extensive for translationally invariant gapped systems in the low-temperature regime is provided in Appendix A of [32]). The optimal controls are found to be λ * 2 ≈ 5, λ * 1 ≈ 8.

V. ASYMPTOTIC EXPANSIONS AND CRITICAL SCALING
Inspired by Ref. [9] we can find the best power for a given fixed efficiency, i.e. η = γη C = γ 1 − Tc T h , which in the symmetric low-dissipation regime means to fix (we identify ∆S ≡ ∆ h S) which relates τ c and τ h , after which it is possible to maximize the power by enforcing ∂ τx P = 0 (here τ x can be either τ c (τ h ) or τ h (τ c ) via Eq. (88)), to obtain in adimensional units), for a system of N qubits with different degrees of control. We obtain asymptotically: Cmax 0.44N for N independent qubits with a single control parameter, Cmax 0.59N for an Ising chain with two control parameters, and Cmax N 2 /4 for full control on the spectrum. by choosing The correspondent work is A. Critical scaling Now we consider the optimal engine described in the main text whose cycle consists of G(sτ c ) = (1+ g(s))G(0) with s ∈ (0, 1), 1 but finite, and G(τ c + sτ h ) = G((1 − s)τ c ). G is chosen appropriately to maximize the heat capacity given the allowed control. In the limit 1 the shape g(s) becomes irrelevant as far as it is a smooth function with g(0) = 0 and g(1) = 1, which implies ġ = 1. This leads to We also consider that we scale up the engine with N while approaching γ → 1, by setting (using the notation for critical exponents of a second order phase transition as in [35,36]) Then, expanding the relevant quantities for N 1, we obtain at leading order in N : Now we look at the fluctuations of work. The cycle consists of four processes: Two (quasi)-isothermal processes and two quenches. We note that as N → ∞, we have that W → (T h − T c )∆S, so that the isothermal processes become exact at leading order in N which implies that they become fluctuationless [29]; note that this is expected as τ c → ∞ with N → ∞, which makes the low-dissipation assumption also more and more exact as we increase N . Secondly, regarding the two quenches: H → HT h /T c and the inverse one H ← HT h /T c ; it is easy to see that the work fluctuations are given by and thus Hence we have that As explained in the main text, if we want to exploit the critical scaling in the phase transition of the system must satisfy ≤ N −1/(2−α) . Then, using dν = 2 − α [69], we obtain which is the same result found in [36] for the Otto engine proposal. More details on the necessary critical exponents, and on the fluctuations of the engine, are provided in the main text.

B. Comparison with Otto cycle
Let us check in more detail how our results compare to Ref. [35], where an Otto cycle is considered. Following [35], let us define the internal energy and the heat capacity The work output of a Otto cycle working between Hamiltonians λ h H ← λ c H for a fixed efficiency η = γη C , ∆η = η C (1 − γ) is given by Expanding for low ∆η which corresponds to λ c H c ≈ λ h H h , using (102) and keeping only leading terms in ∆η we obtain where we defined It is important to note that the linear expansion (104) of (103) is only justified close to the phase transition point (105), whose width scales as δ ∝ N −1/(dν) . This sets an extra requirement on ∆η: The power is simply W/τ where τ is the time of the cycle, which involves two thermalization processes [35]. Let us take τ = 2κτ eq , where κ measures how exact the thermalisation process is (the error being exponentially small with κ if one assumes a standard exponential relaxation). Then the power reads Let us now consider the Carnot cycle. After maximisation over (∆S) 2 /Σ we have shown in the main text that Expanding around γ ≈ 1 , we have which is correct for a fixed γ close to 1.

VI. EXPLICIT PROTOCOL FOR A 2 N − 1 DEGENERATE SYSTEM
For the sake of completeness, we show here how the results of the paper apply explicitly for a feasible driving protocol close to the optimal one, in the case of a D-level system with a ground state and D − 1 degenerate excited states. As we explained in Sec. IV case (a), this Hamiltonian is motivated as it maximises the heat capacity given a D-level system. For comparison purposes, we will sometimes use N satisfying D = 2 N , so that the Hamiltonian can be thought of N suitably interacting particles. Note that the case N = 1 corresponds to the driving of a qubit [17,57].
Consider a D-level system in a thermal state at temperature 1/β with an engineered Hamiltonian with a ground state and a d ≡ D − 1 degenerate excited states with gap E/β, i.e.
such that the thermal state is where we defined the driving ground state population Consider the dynamics given by the equation d dt ρ(t) = Γ(ω(t) − ρ(t)), Γ ≡ 1 τ eq (113) (which induces the metric (4) of the main text in the context of continuous time dynamics, as explained in Section II A or Ref. [28]). We can then consider a solution in the form characterized in terms of the ground state probability, which by the slow-driving approximation [14] can be solved to be that is, expressing quantities in terms of the adimensional time unit ρ(t) = ρ(sτ ) (here τ is the duration of the driving), and using the notationȦ = ∂ s A = τ d dt A, The heat exchange can be computed at all the orders βQ (j) = The point of maximum heat capacity that optimizes the power output of this system, as shown in Sec. IV case (a), is asymptotically E ∼ ln d which implies q ∼ 1 2 . We consider thus a modulation near this point, i.e. a driving protocol in the form q(t) = 1 2 1 + ε cos(πt/τ ) , with a small abuse of notation q(s) = 1 2 1 + ε cos(πs) , or equivalently E(s) = ln( dq(s) 1−q(s) ). Notice here that we chose to use the letter ε instead of used in the main text and in Eq. (87), because the latter multiplies the Hamiltonian, which scales linearly in N on the optimal point, while here the modulation of E(s) is of order ε with no pre-factor. Hence for scaling comparisons it is important to remember ∼ ε/N .
The low-dissipation we considered in the text consists in neglecting all the terms of order O( 1 τ 2 Γ 2 ) on (which in this specific case would consist only in a renormalization of the total result), hence identifying βQ = ∆S − Σ/τ : In the high efficiency regime we can then write from the asymptotic expansions of Eq. (90) and (91) In Figure 7 we show the evolution of the ground state probability p(s) (which characterizes the whole state solution) according to different approximations, as well as the actual performance of the engine compared to the maximal one (which is obtained in the limit ε → 0), for feasible choices of the parameters. The slow-driving approximation gives results that almost coincide with the full analytical series for this particular protocol, and the performance is close to the optimal one.