Finite-time quantum Stirling heat engine

We study the thermodynamic performance of the finite-time non-regenerative Stirling cycle used as a quantum heat engine. We consider specifically the case in which the working substance (WS) is a two-level system. The Stirling cycle is made of two isochoric transformations separated by a compression and an expansion stroke during which the working substance is in contact with a thermal reservoir. To describe these two strokes we derive a non-Markovian master equation which allows to study the dynamics of a driven open quantum system with arbitrary fast driving. We found that the finite-time dynamics and thermodynamics of the cycle depend non-trivially on the different time scales at play. In particular, driving the WS at a time scale comparable to the resonance time of the bath enhances the performance of the cycle and allows for an efficiency higher than the efficiency of the slow adiabatic cycle, but still below the Carnot bound. Interestingly, performance of the cycle is dependent on the compression and expansion speeds asymmetrically. This suggests new freedom in optimizing quantum heat engines. We further show that the maximum output power and the maximum efficiency can be achieved almost simultaneously, although the net extractable work declines by speeding up the drive.


Introduction
A flourishing research activity has developed recently around the understanding of the thermodynamic properties of quantum systems [1,2,3,4]. Special attention has been devoted to quantum heat engines and refrigerators triggered by both new theoretical questions and technological advancements in dynamical control of microscopic systems [5,6,7,8,9,10,11,12,13,14,15,16]. From the theoretical point of view a natural question is whether quantumness of the working substance (WS) can be exploited to achieve better performances over the classical systems. The role of quantum effects has been demonstrated for example in [17,18,19,20,21,22,23]. It has been also ascertained that the creation of coherence between energy levels leads to inner friction and reduction of the extractable work [24,25,26,27,28,29,30]. Thus, the sole use of a quantum WS does not in general guarantee superiority over the classical counterparts [31]. The study of quantum thermal machines has relied mostly on the Markovian (Lindblad) description of an open system dynamics, which guarantees non-negative entropy production rate and consistency with the second law of thermodynamics [5]. It has been shown that non-Markovian dynamics could lead to negative entropy production for the open system reduced state, however, the sum of the entropy change of the bath and the open system together is positive [32]. Besides these studies, non-Markovianity has been found to be influential in the performance of the quantum heat engines [33], and may also enhance the output power [34].
Quantum heat engines are composed of a series of strokes defined by completelypositive and trace preserving maps whose product forms the propagator over the full cycle [1]. Each stroke corresponds to coherent drive, dissipation to a heat bath, or simultaneous driving and dissipation. The maximum amount of extractable work produced by a cycle is obtained by ideal reversible processes which has the minimum entropy production. However, a true reversible process is infinitely slow and gives rise to zero output power. Therefore, to find a trade-off between the power and the efficiency one has to consider cycles running at finite times. Concerning the finite-time adiabatic strokes, i.e. coherent drive on the isolated WS, the shortcut-to-adiabaticity approach provides a way to mimic the adiabatic process [35,36] and has recently been demonstrated experimentally with superconducting circuits [37,38]. This technique has been employed, for example, to boost the performance of an Otto refrigerator [39] and an Otto heat engine [40]. The Otto cycle with a quantum WS has been extensively studied in the adiabatic as well as in the non-adiabatic case [28,12,19,13,16,21,40,34].
The quantum Carnot and Stirling cycles running in finite times have received less attention. The main reason is that a finite-time isothermal stroke is more tricky to study and optimize due to the simultaneous drive on the WS and its coupling to the heat bath. Usually, a slow drive whose effect falls within the validity of the adiabatic limit is assumed, allowing one to ignore the non-adiabatic effects [41,42,43]. This assumption is relaxed in the derivation of a time-dependent Markovian master equation to capture non-adiabatic effects but retaining the assumption that the time scales of the external drive and the ones of the coupling to the bath are still well-separated [44]. Using this master equation, it has been proposed to reverse-engineer the thermalization to find a corresponding driving protocol which provides shortcut to equilibration [45]. Alternatively, manipulating the coupling between the WS and the bath is also shown to speed-up isothermal strokes [46]. In particular, the Stirling cycle has been studied in the ideal adiabatic regime [47,48] and only very recently a finite-time scenario in an optomechanical implementation has been studied [49]. There the compression and expansion strokes have been treated in the Markovian regime and the adiabatic limit and, as the authors state, a deeper investigation to include non-Markovian and out-ofequilibrium effects is still lacking.
Here we fill this gap by studying the thermodynamics of the Stirling cycle, used as a quantum heat engine, in finite-times. A two-level system is considered as the WS and we investigate the role of different time scales involved in its dynamics. The description of the compression and expansion strokes in contact with the thermal baths requires solving the real-time dynamics of an open system beyond the adiabatic limit. In this work, we analyze the isothermal stroke in finite times without any restriction on the time scale of the drive, thus allowing for the dynamics to be non-Markovian. We note that in the non-adiabatic regime the quantum system, although in contact with a bath at fixed temperature, is brought out-of-equilibrium and its temperature is in general not defined. Hereafter by isothermal we shall refer to the fact that the WS is in contact with a bath in equilibrium at a well defined temperature. To study the dynamics of the WS during the isotherms we derive a non-Markovian master equation using the results presented in Refs [50,51]. We show that the master equation contains two time-dependent parts, the rotating (R) and the counter-rotating (CR) terms. Each part also includes a Lamb shift term with the important difference that the one coming from the rotating part commutes with the non-interacting Hamiltonian whereas the second one does not.
We observe that the efficiency of the cycle depends on the interplay between the driving time, the bath's correlation time and the resonance time of the hot and the cold baths. Interestingly, the efficiency exceeds that of an ideally slow cycle if we drive the qubit at a frequency comparable to the resonance frequency of the baths. The average output power also shows a similar behavior, allowing to get maximum power and maximum efficiency approximately at the same time scale for the drive. This is however not true for the net extractable work, which decreases as we speed up the drive. Our results also show that the performance of the cycle is non-trivially dependent on the individual speed in the compression and expansion strokes when the qubit is coupled to the hot and the cold bath respectively. As this dependence is in general asymmetric regarding the cold and the hot bath, it opens the possibility to optimize the performance of the cycle by choosing asymmetric compression and expansion speeds/protocols. The paper is organized as follows: in Section 2 we introduce the master equation, followed by the presentation of the Stirling engine in Section 3. The calculation of work and heat for the Stirling engine is done in Section 3. Section 4 deals with the evaluation of the thermodynamic performance of the engine. Finally, section 6 is devoted to concluding remarks.

The Master Equation
To study the dynamics of a driven WS in contact with a thermal bath we employ a non-Markovian master equation obtained by applying the approach developed in Refs. [50,51]. There, assuming weak coupling to the baths and the Born approximation, a general time-convolutionless non-Markovian master equation is derived using the Nakajima-Zwanzig method. Such a master equation is valid for any characteristic time scale of the drive, e.g. the period in a periodic drive or the ramping time in the case of a switching. The master equation retains both rotating and counter-rotating contributions, where the latter is specifically non-negligible at fast driving speeds. Here we introduce the operatorial form of the master equation and discuss its main features, leaving more details on the numerical implementation and how to recast the master equation in the adiabatic basis in appendix A.
Let us consider a quantum system subject to an external coherent driving field and weakly coupled to a thermal bath at an inverse temperature β. The total Hamiltonian whereĤ S (t) andĤ B (t) are the bare Hamitlonian of the open system and the bath respectively. We write the interacting Hamiltonian in the form withŜ(t) being a time-dependent operator acting on the open system, andB an operator acting on the bath. The non-Markovian master equation reads [50] Here Φ(t) is the correlation function of the bath and ρ B denotes the equilibrium state of the bath at an inverse temperature β. The correlation function is related to the bath's coupling spectrum G β (ω) via the Fourier transform G β (ω) = +∞ −∞ ds Φ(s)e iωs . By decomposing the operators in the master equation w.r.t. the instantaneous eigenvectors ofĤ S (t), denoted by {| i (t) }, we get (see appendix A for more details) whereĤ ef f (t) =Ĥ S (t) +Ĥ   [·] account for the exchange of energy with the bath and/or decoherence.
Note that the dissipators accounting for two different baths are additive by construction if one assumes that the baths are initially uncorrelated. Naturally, the specific expressions of the different terms appearing on the r.h.s. of Eq. (7) depend on the choice of the free Hamiltonian of the system and more importantly on the coupling HamiltonianĤ I (t).

Quantum Stirling heat engine
The Stirling cycle is composed of two isothermal strokes and two isochoric thermalizations. Classically it is common to supplement the cycle with two extra steps which involve the interaction of WS with the so-called regenerator. The latter is typically some substance with a very high heat capacity whose task is to absorb heat from the WS during the cooling isochoric stroke and transfer this heat back to the WS during the heating isochor to improve the overall efficiency and minimize the waste heat. In this work we do not consider a regenerative setup, i.e. the WS interacts directly with the heat baths instead of the regenerator. The diagrams for temperature and polarization versus level separation for the Stirling cycle as a heat engine are depicted respectively in the panels (a) and (b) in Fig. 1, where the polarization is given by n(t) = tr[Ĥ S (t)ρ(t)]/ω(t). Consistent with these diagrams and assuming a two-level system (TLS) as the WS, the cycle consists four strokes: 1-Isothermal compression, process a → b, with duration τ ab : the level separation of the TLS reduces from ω 2 to ω 1 while it is coupled to the hot bath at an inverse temperature β h .
2-Isochoric thermalization, process b → c, with duration τ bc : the TLS is disconnected from the hot bath and is brought to contact with the cold bath at an inverse temperature β c , with which it thermalizes while the external drive is off.
3-Isothermal expansion, process c → d, with duration τ cd : the level separation of the TLS increases from ω 1 back to ω 2 while it is still coupled to the cold bath.
4-Isochoric thermalization, process d → a, with duration τ da : the TLS is disconnected from the cold bath and is brought back to contact with the hot bath. The TLS thermalizes while driving is off.
Therefore, total duration of a full cycle is T = τ ab + τ bc + τ cd + τ da . We specifically consider a setup implementable with a superconducting circuit schematically shown in the panel (e) of Fig. 1. It is worth mentioning that a related design has been also put forward in [39] as a possible implementation of the Otto refrigerator. The free Hamiltonian of the TLS and the TLS-bath coupling Hamiltonian are respectively denoted byĤ S (t) andĤ The quantum Stirling heat cycle and its implementation using superconducting circuits. The β-ω and n-ω diagrams of the ideal Stirling cycle working between inverse temperatures β c and β h are respectively shown in panel (a) and panel is the instantaneous polarization of the two-level system. Panel (c) shows the coupling spectra of the two heat baths peaked at resonance frequency ω r and the range of frequency drive of the two-level system ([ω 1 , ω 2 ]). Panel (d) shows the piece-wise continuous coupling to the heat baths and the frequency drive of the two-level system as a function of time during a full cycle. Panel (e) shows a proposed circuit to implement the cycle using a superconducting two-level system (TLS) capacitively coupled to two RLC resonators, playing the role of the cold and hot baths. The energy levels of the two-level system is drived by tuning the external magnetic flux applied to the superconducting qubit.
Here ω 0 denotes a reference energy scale for the non-driven qubit. The operatorB α acts on the cold/hot bath, with α = c, h, and incorporates the coupling amplitudes between the WS and the corresponding bath. In order to realize the connection and disconnection from the baths required at steps b and d in the cycle, a tunable coupling element between the TLS and the resistor is required. Several types of tunable couplers have been proposed and studied, e.g. based on dressed states [52], additional qubits [53], additional single Josephson junctions with current bias [54], and using a SQUID junction whose effective inductance is modulated by a bias magnetic field [55,56]. As our main motivation here is not the realization of the setup, and for the sake of simplicity, we just assume an ideal connection/disconnection protocol described by a piece-wise continuous function λ α (t), as shown in the panel (d) of Fig. 1.
As depicted in the panel (d) of Fig. 1, we choose the driving protocol q(t) such that the instantaneous level separation ω(t) = 2ω 0 q(t) 2 + ∆ 2 of the TLS changes linearly Resonance time of the bath 2π/ω r τ C Correlation time of the bath Controlled by f β h Inverse temperature of the hot bath 2/ω 0 β c Inverse temperature of the cold bath 5/ω 0 (g c , g h ) Set of TLS-bath coupling amplitudes Resonance frequency of the baths Quality factor of the bath's resonators 2 or 3 τ D Unit of driving duration Duration of the isochoric strokes 6 × τ R (g 1 ) with time within the interval [ω 1 , ω 2 ] with a given constant speed. This requirement fixes unambiguously q(t) = ω(t) 2 /4 − ∆ 2 . A relevant coupling spectrum for the baths regarding the setup considered in this work is shown in the panel (c) of Fig. 1 and takes a specific expression given by [39] With i = c, h denoting again the cold and hot baths, the coupling strength to the bath is described by g i , the resonance frequency of the bath is denoted by is the quality factor of the resonators. We assume identical resonance frequency for the two baths denoted by ω r and set the values of coupling strengths g c and g h such that the corresponding spectra have the same amplitudes at ω r (see panel (c) of Fig. 1). All the relevant physical parameters and their values used in this work are reported in Tab. 1.
For the specific Hamiltonian given in Eq. (8), the instantaneous energy basis of the TLS reads | e (t) = cos θ t |e + sin θ t |g , with θ t = (1/2) cot −1 (q(t)/∆) and |e(g) as the eigenbasis ofσ z Pauli operator. By defining the transition operatorL(t) = | g (t) e (t)| between the instantaneous energy basis of the TLS, the master equation in Eq. (7) takes the explicit form Finite-time quantum Stirling heat engine where the exact expressions for the Lamb shift terms δ [·] is, however, too cumbersome to be included here. An interesting feature of the Lamb shifts is that while the Lamb shift contribution due to the rotating term is proportional toĤ S (t), the counter-rotating one does not commute withĤ S (t). Temporal behavior of the rates involved in the generator L t is shown in Fig. 2 for different compression and expansion speeds during the isothermal branches. In the rest of the paper, we consider a scale for the driving duration denoted by τ D . Recalling that the amplitude of the spectra of the cold and hot baths are set to be identical at the resonance frequency ω r , the value of τ D is fixed to the relaxation time of the TLS when the coupling amplitudes are (g c = 0.2, g h = 0.17), thus τ D := τ R (g 1 ). In addition, duration of the isochoric branches are always fixed at τ th = 6 × τ R (g 1 ).
By neglecting the counter-rotating terms in Eq. (12) we get a time-dependent master equation in the Linblad form, which we describe by L Consider this generator at a given fixed time t = τ denoted by L (R) τ , which means all the rates, Hamiltonian, and jump operators are set to their configuration at t = τ and remain unchanged for t > τ . We define the invariant state of this generator byρ (R) eq (τ ), such that L (R) τ [ρ (R) eq (τ )] = 0. It is straightforward to check that the invariant state is given bŷ with Γ(τ ) = γ (↑) (τ ) + γ (↓) (τ ). We stress that due to the explicit time dependency of the decay rates,ρ (R) eq (τ ) is not necessarily identical to a Gibbs state at the same temperature of the heat bath. As depicted in the two upper panels of Fig. 2, it is only for the asymptotic slow driving (adiabatic limit) that the decay rates γ (↑) (τ ) and γ (↓) (τ ) approach to their Markovian limits 2πG β (±ω(τ )) [51] and one consequently gets the equilibrium stateρ eq (β, τ ) = exp(−Ĥ S (τ )β)/tr[exp(−Ĥ S (τ )β)], where β is the inverse temperature of the bath with whom the TLS interacts. Asymptotic state of the full generator L t which includes the counter-rotating terms is, however, more complicated and does not depend solely on the rates γ (↑/↓) (t), especially for fast drives. Consider the full generator at a given fixed time t = τ denoted by L τ . We define its asymptotic state formally byρ * τ = lim t→∞ [exp(tL τ )ρ i ], whereρ i is some initial input state. We now proceed to utilize these tools to characterize the Stirling cycle used as a heat engine.
We set t a = T and exclude the first cycle 0 ≤ t < T to guarantee that the rates and state of the TLS all reset to their initial values at t = 2T , i.e.ρ(2T ) =ρ(T ) and L 2T = L T . Moreover, the duration of isochoric strokes are set sufficiently large (τ bc = τ da = 6 × τ R (g 1 )) such that the TLS can reach its asymptotic equilibrium states at the end of b → c and d → a branches. Accordingly, the two points a and c are always fixed in our analysis, as shown in the panel (a) of Fig. 3. Nonetheless, we consider arbitrary duration for the isothermal strokes. If the driving is sufficiently slow, the WS remains in an instantaneous equilibrium state with the bath during the whole process. This ideal case corresponds to the ab and cd trajectories in the panel (a) of Fig. 3. However, a faster drive kicks the WS out of the manifold of equilibrium states and, consequently, its trajectory deviates from the ideal isothermal(adiabatic) ones. The opposite regime is when the drive is so fast that the dynamics of the WS is essentially diabatic and its state remains unchanged during the process. Therefore, at the end of the diabatic process we end up at the pointb(d), instead of adiabatic points b(d).
Let us now define the actual target points of the WS at the end of the compression and expansion processes with some arbitrary speeds, respectively by b and d . The corresponding points of the asymptotic (equilibrium) statesρ * t b andρ * t d of the WS at the end of the processes are also denoted by b * and d * . Therefore, as it is shown in the panels (b), (c), and (d) of Fig. 3, by increasing the speed of driving one changes the trajectory of the WS within the two areas abb and cdd, moving from the ideal adiabatic trajectories towards the diabatic ones. Within this general picture, we now study in detail how different speeds affect the thermodynamic performance of the Stirling heat engine.

Work and heat
Studying the performance of the heat engine requires to calculate the work done as well as the energy exchanged with the baths during each stroke of a full cycle. In making the separation between work and heat it is important to include the Lamb shift in an effective Hamiltonian of the system [57], which readŝ where the summation is over the terms corresponding to the hot and the cold baths. Likewise, the full dissipator acting on the WS has two parts each corresponding to one of the baths: Having the Hamiltonian and the dissipator of the dynamics, we can calculate the average of work and heat transferred. For the average work done on the WS in the time interval [t 1 , t 2 ] one has which relates to the average output power P (t 1 , t 2 ) during this time interval via The fact that the Lamb shift enters the definition of the work has important consequences. Specifically, if the Lamb shifts vary in time (which indeed happens here due to presence of a memory kernel in the master equation) it is possible to have nonzero work even when the external drive is off. Furthermore, the average heat transferred into the WS in the time interval [t 1 , t 2 ] is given by By denoting W net as the average net extractable work during a full cycle, according to the first law of thermodynamics one has W net = − Q net , where Q net is the net average heat transferred. Consider the net positive heat transferred into the WS labeled by Q h , then the efficiency of the cycle is determined by Let us recall that in a regenerative classical Stirling heat engine the heat transferred during the isochoric branch d → a is not included in calculating the efficiency, since the regenerator is an internal component of the engine. However, here we let the WS interact directly with the hot bath during the stroke. Therefore, the net positive heat transferred into the WS has contributions both from the a → b and d → a branches. Moreover, in the classical regenerative Stirling heat engine, the regenerator helps to minimize the wasted heat and increase the efficiency closer to the Carnot bound: η C = 1 − β h /β c . As we do not resort to regeneration, we expect the efficiency to be well-below the Carnot bound η C = 0.6 (considering β c = 5 and β h = 2).
Before presenting our numerical results, we note that by considering the bare HamiltonianĤ S (t) (excluding the Lamb shifts) one can provide analytic expressions of the efficiency regarding four limiting cases. As depicted in the panel (a) of Fig. 3, these case are: (abcda) trajectory corresponding to the ideal adiabatic processes -(abcda) trajectory which follows the diabatic passage -(abcda) trajectory corresponding to the adiabatic compression and diabatic expansion processes -and finally (abcda) trajectory which follows the diabatic compression and adiabatic expansion processes. We call these cases respectively (ss), (f f ), (sf ), (f s), with s and f denoting slow and fast, respectively. The analytic analysis of efficiency for these cases is presented in Appendix B. Our aim is to examine the performance of the heat cycle for the situations between these four limiting cases and by considering the real-time evolution of the WS in finite times.

Thermodynamic performances
We first consider the situation in which the speed of compression and expansion processes are identical, i.e. τ ab = τ cd , which corresponds to the situations shown in the panels (b), (c), and (d) of Fig. 3. Efficiency of the cycle is plotted as a function of τ ab,cd using the solid curves in the upper panel of Fig. 4. We have plotted the efficiency calculated using both the effective HamiltonianĤ ef f (thick green curves) and the bare HamiltonianĤ S (thin green curves). The most striking observation is a peak in the efficiency at some values of τ ab(cd) which exceeds the efficiency of the adiabatic cycle shown by the dotted Regarding the efficiency, the solid curves correspond to the symmetric driving with τ ab = τ cd . The asymmetric cases are plotted with the large dashing red and small dashing blue, corresponding respectively to (1): when τ ab = τ D is fixed and τ cd changes and (2): when τ cd = τ D is fixed and τ ab changes. In the upper panel, the thick curves indicate the efficiency calculated w.r.t. the effective Hamiltonian and the thin lines correspond to the calculations considering the bare Hamiltonian. Note that, however, the power is calculated only w.r.t. the bare Hamiltonian. The efficiency of the asymptotic cases discussed in the appendix B are also marked, specifically the efficiency of the ideally slow cycle (ss) is plotted using the dotted black line. black line. To realize the relation between the observed efficiency enhancement and different physical time scales involved in the dynamics of the WS, to say relaxation time τ R , bath correlation time τ C , and bath resonance time scale τ B , we have plotted the efficiency for four different cases. These are considering two different values of the relaxation time, τ R (g 1 ) and τ R (g 2 ) = τ R (g 1 )/2, and two different values of the bath correlation time, τ C (f = 2) and τ C (f = 3) = 1.43 × τ C (f = 2). Moreover, we set the value of τ B fixed for all the four mentioned cases. Looking at Fig. 4, it is clear that the relevant parameter for the observed peak in the efficiency is the bath resonance time τ B , such that when τ ab(cd) are close to τ B we observe the enhancement in the efficiency. On the contrary, it is clear that when the time scale of the drive is close to the bath correlation time τ C the efficiency decreases. The output power of the cycle for the same settings is also plotted as a function of τ ab,cd in the lower panel of Fig. 4. We note that since we are interested in the work done by the external drive, the output power is only calculated with respect to the bare Hamiltonian. Interestingly, the average output power benefits from enhancement when τ ab,cd τ B as well. However, the peak in the power is happening at a slightly larger time scale than those for the efficiency. As expected, the output power decreases by increasing τ ab,cd . The same behavior also holds for very short time scales, when the extractable work diminishes at ultra-fast driving due to an increase in the irreversibility.
Studying the energy flow from or into the WS is essential to comprehend the observed boost in efficiency and power. Due to the limited space and without the loss of generality, we present the energetic results only for the case with f = 2 and (g c , g h ) = g 1 . The average heat transferred, average work and the net average work are plotted in Fig. 5 considering the effective HamiltonianĤ ef f in the left panel, and the bare HamiltonianĤ S in the right panel. Using the effective Hamiltonian to calculate the energy terms leads to some non-zero amount of average work for the isochoric strokes due to the time-dependent Lamb shifts. The corresponding terms are absent when we use the bare Hamiltonian as the drive is off during the isochoric strokes. Moreover, the average work in the expansion stroke is higher when we consider the effective Hamiltonian, which shows up also in the net extractable average work and consequently results into a higher efficiency in comparison to the case of using the bare Hamiltonian (see Fig. 4). This behavior is again due to the non-zero Lamb shift terms and the fact that by including them the effective frequency span of the WS is higher than the bare frequency span ∆ω = ω 2 − ω 1 , specifically for c → d process. The average net work approaches its adiabatic limit as we increase τ ab,cd and decreases by speeding up the drive. One can see that the average heat transferred during the compression and expansion processes goes to zero as we decrease τ ab,cd , because the WS does not have enough time to exchange energy with the baths. Moreover, the heat transferred during the isochoric strokes reaches its non-zero minimum by approaching the diabatic limit (pointsb andd in Fi.g 3).
Besides these asymptotic scenarios, we observe a dip in the heat transferred and the net average work at some values of τ ab,cd coinciding with the peak in the efficiency. The dip especially indicates some extent of suppression of heat transferred to the cold bath. With a given amount of heat absorbed from the hot bath, if the WS dissipates less to the cold bath it means that the work done is higher and thereby the efficiency as well. This may suggest that a faster a → b process in Fig. 3 is in general beneficial, as the state at the end point b gets closer to the equilibrium state at the point c and there would be less dissipated heat to the cold bath. However, the faster is a → b process, the less amount of heat is absorbed from the hot bath which restricts the amount of extractable work too. Note that a similar situation also happens for the c → d process considering the heat dissipated during the expansion and the heat absorbed during the thermalization d → a. Therefore, there must be some trade-off giving us the optimum efficiency in the intermediate situation.
To shed some light on the facts discussed above, we consider the distance between the states at some end points in Fig. 3. First, the distance between b (d ) and b * (d * ) allows us to figure how far we are from the instantaneous equilibrium states at the end of the compression and expansion processes. Second, the distance between the states at b (d ) and c(d) indicates how far the WS is from the thermal states at the end points of the isochoric strokes. To measure the distance between two states ρ 1 and ρ 2 we use the relative entropy between them defined by Looking at the left panel of Fig. 6, the distance between the end point stateρ b (d ) and the corresponding instantaneous asymptotic steady stateρ b * (d * ) decreases as we increase the driving duration, although there are some fluctuations in the process c → d . The interesting feature is that in the same time scale at which we observed the boost in the efficiency, the distance betweenρ b andρ b * is close to its maximum, whereas, two statesρ d andρ d are rather close. The distance to the instantaneous asymptotic steady stateρ b * (d * ) is an indicator of irreversibility of the process, i.e. a large distance indicates higher irreversibility and smaller amount of extractable work. Looking at the right panel of Fig. 6, we realize a dip at τ ab,cd τ B . In addition, we observe that distance for the a → b process is in general higher than c → d process.
An interesting feature about the results in Fig. 6 is the asymmetry in the behaviors of a → b and c → d processes. This fact rises the question whether considering different Here we set f = 2 and (g c , g h ) = g 1 Figure 7. The net average work and the net heat transferred into the WS as a function of τ ab,cd . The solid green curves correspond to the symmetric case with τ ab = τ cd , the large dashed red represents the asymmetric case (a), and the small dashed blue to the asymmetric case (b) discussed in the main text. Here we set f = 2 and (g c , g h ) = g 1 .
speeds for the two processes has some non-trivial effects on the performance of the heat engine. To make this point clear, we consider two different situations: (a) setting τ ab = τ D fixed while changing τ cd and (b) setting τ cd = τ D fixed while varying τ ab . The efficiency and the average output power of these cases are plotted in Fig. 4 using large dashing red and small dashing blue lines, respectively. Again the thick curves correspond to calculating the energies w.r.t. the effective Hamiltonian and the thin curves to the bare Hamiltonian. Interestingly, efficiency of the asymmetric cycles is always higher than the symmetric ones. However, superiority of the two asymmetric cases with respect to each other depends non-trivially on the time scale of the driving. Let us also examine the energetic of the asymmetric cycles in comparison to the symmetric ones depicted in Fig. 7. We note that the amount of net work is dependent on the total time of the expansion and compression process, τ tot = τ ab + τ cd , and in general decreases by decreasing τ tot due to irreversibility. However, within the two asymmetric cycles with the same value of τ tot , we notice slightly different values for the net average work, indicating again the importance of the finite-time effects in the performance of the heat engines.

Conclusions
In conclusion, we have studied the performances of a Stirling cycle when operated as a finite-time quantum heat engine. We first derived a non-Markovian master equation which allows us to study the dynamics of an open quantum system including counter rotating terms and avoids to make a clear distinction between the time scales of the system and of the environments. Thanks to this, we have been able to study the effect of the competing time scales, such as the typical time scale of the drive and the bath correlation/resonance time, on the performances of the heat engine. The main motivation of this work was to explore the performance of the heat engine operating in the non-adiabatic regime. Interestingly, we found that driving the WS at a time scale comparable to the resonance time of the bath, in addition to a boost in the output power, let us get an efficiency that is higher than the efficiency of the slow adiabatic cycle. One should note however that the net extractable work decreases by speeding up the cycle due to higher and higher degree of irreversibility. The other important finding in this work was the non-trivial dependency of the performance of the heat engine on the individual speed of the compression and expansion processes. Interestingly, one may achieve better performances by applying asymmetric compression and expansio speeds rather than a symmetric one. The latter opens new possibilities to optimize the performance of the quantum heat engines. An important aspect that is not covered in the current work is the reverse Stirling cycle working as a refrigerator. In [47], the authors have discussed that in general the revers cycle of a quantum Stirling heat engine might not be a refrigerator. As an outlook of our work, it is therefore interesting to explore how finite-time effects influence the operating range of the quantum Stirling as a refrigerator and also its performance. In addition, our results motivate for further studies aiming at optimization of quantum thermodynamic cycles with finite-time driving protocols. Here we briefly elaborate our numerical method to solve the master equation The unitary propagatorÛ (t, 0) = T exp −i t 0 dτĤ S (τ ) can be calculated numerically in a time interval [0, t max ] by solving the Schrdinger equation Then owing to the divisibility of the unitary propagator we get Inserting this solution inŜ(t, τ ) =Û (t, τ )Ŝ(τ )Û (t, τ ) † and decomposing the operators with respect to Pauli operator basis {σ 0 =Î,σ x ,σ y ,σ z } we get By introducing a similar decomposition for the operatorŜ(t) given byŜ(t) = Σ j λ j (t)σ j , the second term on the r.h.s. of Eq. (A.1) will be rewritten as with the time-dependent rates One can further arrange Eq. (A.9) into rotating (R) and counter-rotating (CR) parts with respect to the instantaneous energy basis. The rotating part takes the form while the counter-rotating part L (CR) t includes all the remaining terms. We focus now on the specific model considered in this paper given by the Hamiltonian in Eq. (8), and instantaneous energy basis labeled by | e (t) and | g (t) . HavingŜ(t) = λ(t)σ y , and the expression for the energy basis given in Eq. (10) and Eq. (11), we get (A.13) Considering the numerical solution in Eq. (A.4), a decomposition forŜ = n,m=e,gsnm (t, τ )Ê nm (t) in the energy basis is given bŷ S(t, τ ) = (A.14)  s 0 +s x ∆ +s z q(t) q(t) 2 + ∆ 2  Ê ee (t) +  s 0 −s x ∆ +s z q(t) q(t) 2 + ∆ 2  Ê gg (t) wheres i ≡s i (t, τ ). Note that one has R (↓) eg,ge (t) = R (↑) ge,eg (t) * and R (↓) ge,eg (t) = R (↑) eg,ge (t) * . Moreover, sinces 0 (t, τ ) ≡ 0 and all others i are real valued, we also have R (↓) ee,ge (t) = −R (↓) gg,eg (t) * and R (↑) ee,eg (t) = −R (↑) gg,ge (t) * . Accordingly, the rotating part L  (t) = Im[R (↓) ee,eg (t)]. However, the expression for the counter-rotating dissipator is so lengthy that does not fit here.

Appendix B. Analytic considerations for the asymptotic Stirling cycles
Apart from the dynamical approach of the paper, we provide some analytic analysis of the energy transferred for some asymptotic cases. Assume that at the end of the isochoric strokes b → c and d → a the TLS relaxes to the corresponding thermal stateŝ Then according to the first law of thermodynamics we get the net average work done on the qubit by W net = −( Q ab + Q bc + Q cd + Q da ).

The (fs) cycle: ideally fast compression and slow expansion
In this extreme, a → b process is done in a finite but very fast time scale, such that the process is diabatic. For a sufficiency fast process, TLS does not have time to exchange heat with the hot bath and Q ab = 0. Nonetheless, there is some non-zero average work done on the TLS that can be obtained by the change in its internal energy. Since the process is diabatic, state of the TLS remains at its initial configuration at time t a , therefore The remaining energy terms can be calculated similar to the (ss) case.

The (sf ) cycle: ideally slow compression and fast expansion
This is the opposite situation of the (f s) cycle, such that Q cd = 0 and W cd = tr (Ĥ S (t d ) −Ĥ S (t c ))ρ tc . (B.8) The (ff ) cycle: ideally fast compression and fast expansion Finally when both the processes are diabatic, one has Q ab = Q cd = 0 and the amounts of average work can be obtained as discuss in the two previous cases. Bibliography