Dynamical heat engines with non--Markovian reservoirs

We discuss whether, and under which conditions, it is possible to realize a heat engine simply by dynamically modulating the couplings between the quantum working medium and thermal reservoirs. For that purpose, we consider the paradigmatic model of a quantum harmonic oscillator, exposed to a minimal modulation, that is, a monochromatic driving of the coupling to only one of the thermal baths. We demonstrate, at any order in the system/bath coupling strength, that in this setup non--Markovianity of the bath is a necessary condition to obtain a heat engine. In addition, we identify suitable structured environments for the engine to approach the ideal Carnot efficiency. Our results open up new possibilities for the use of non--Markovian open quantum systems for the construction and optimization of quantum thermal machines.

In this work, we address this question for a minimal disturbance of the quantum system, that is, a monochromatic modulation of the coupling to one thermal bath. The same achievement of a heat engine in such a setup, without directly driving the system, is a * fabio.cavaliere@fisica.unige.it non trivial result. Indeed, modulation of the couplings is intuitively associated with dissipation, akin to friction induced by moving parts in strokes of a macroscopic heat engine. Non-Markovianity, associated to the spectral properties of the bath, is here investigated in the paradigmatic model of a quantum harmonic oscillator (QHO) [51,55,56], coupled to two bosonic thermal (hot and cold) baths. Such approach is quite versatile, since it is possible to study the QHO dynamics and thermodynamics, without resorting to any approximations, both in the quantum and in the classical regime, and for arbitrary spectral features of the environment.
We show that, as counterintuitive as it might seem, a dynamical heat engine can be obtained in the above configuration. To achieve such a result we demonstrate that non-Markovianity inherited from the reservoir that feels the driven contact is a necessary but not sufficient condition. Furthermore, we show that by taking advantage of a suitable structured environment, the engine can even approach the Carnot efficiency.
The paper is organized as follows: in Sec. II we introduce the model and outline the methods employed to evaluate the average power and heat currents. The connections between non-Markovianity and a working heat engine are discussed in Sec. III. Subsequently, exploiting the prototypical non-Markovian bath with a Lorentzian spectral density we discuss in Sec. IV the performances of the ensuing heat engine, both in the weak and in the strong coupling regime. An interpretation of the weak coupling limit in terms of quantum Otto cycles is also provided. Conclusions are finally drawn in Sec. V.
FIG. 1. Sketch of the setup: Jν represents the energy current flowing between the QHO and the ν-th contact. The two reservoirs are in equilibrium at a temperature Tν and described by a spectral density Jν (ω) with ν = 1, 2. The coupling of the bath ν = 1 is modulated by a monochromatic driving with frequency Ω, while the coupling to the bath ν = 2 is static.

A. Model
The working medium of the thermal machine is a QHO whose Hamiltonian reads (h = k B = 1) where m and ω 0 are its mass and characteristic frequency, respectively. The QHO is linearly coupled to two reservoirs, with the total Hamiltonian Each bath (ν = 1, 2) is modeled as an ensemble of harmonic oscillators in the usual Caldeira-Leggett [28-30, 57, 58] framework with Hamiltonians We assume that the system/baths couplings can be varied in time [23,47,59], described by the interaction contribution The interaction strengths are described by the parameter c k,ν , and the time-dependence of the couplings is in the dimensionless periodic functions g ν (t) satisfying g ν (t) = g ν (t + T ) with Fourier decomposition In this paper we consider the minimum modulation needed for the couplings in the search for a heat engine, that is a monochromatic drive at frequency Ω for the first contact while the second is kept constant see the sketch in Fig. 1. To model the bath properties we introduce their spectral densities [29] whose precise forms will be specified later. At the initial time t 0 → −∞ the baths are assumed in their thermal equilibrium at temperatures T ν , with the total density matrix, written in a factorized form ρ(t 0 ) = ρ QHO (t 0 ) ⊗ ρ 1 (t 0 ) ⊗ ρ 2 (t 0 ), where ρ QHO (t 0 ) is the initial system density matrix.
The out of equilibrium dynamic of the QHO obeys the generalized quantum Langevin equation [47,60,61] where overdots denote time derivatives and the damping kernels γ ν (t) are linked to the spectral function by with θ(t) the Heaviside step function. The fluctuating force operators ξ ν (t) to the right hand side of Eq. (8) explicitly depend on the initial values of the bath operators X k,ν (t 0 ) and P k,ν (t 0 ). Their expression is We recall that these operators have zero quantum aver-age ξ ν (t) ≡ Tr[ξ ν (t)ρ(t 0 )] = 0 and their time correlators are given by To solve the Langevin Equation (8) one needs the retarded Green's function G(t, t ), which obeys the integrodifferential equation At long times t, when the memory of the initial state for the QHO is lost, the time evolution of the position operator x(t) is directly expressed as a time integral of the retarded Green function as: As we will see shortly, the key relation Eq. (13) will allow us to evaluate all quantum correlation averages, associated to thermodynamic observables. Notice that at long times G(t, t ) acquires the following peculiar form: whereG µ (ω) are the so-called Floquet coefficients obeying the following set of algebraic equations [47]: Here we have introduced and withγ ν (ω) the Fourier transform of γ ν (t).

B. Average thermodynamic quantities
We study thermodynamic quantities in the long time limit, when a periodic steady state has been reached. Notice that the study of instantaneous quantities and their time dependence should be considered with care, especially at finite coupling strength [23][24][25]42]: this however is a problem outside the scope of our paper. We are in particular interested in quantities averaged over one period of the cycle, that are well defined both in the weak and in the strong coupling regime [23][24][25][26].
We start by considering the operatorial time evolution of the total Hamiltonian (2) in the Heisenberg representation: where int,ν (t) -see Eq. (4). Taking the quantum ensemble average, and performing the average over one period of the cycle, we obtain where we have introduced Tr ∂H Tr ∂H Here, P is the total power associated to the time evolution of the system/bath couplings, J ν is the current energy flow from the ν-th reservoir, and A represents the remaining contributions stemming from the QHO and the interaction term. It has been argued that, in general, the term A can be nonzero [25]. However, in our case we can show -see Appendix A -that A = 0. This important result implies that in the long time limit and after the cycling average the total power due to the coupling drives is totally balanced by the reservoir heat currents and fulfills the relation which can be interpreted as a manifestation of the first law of thermodynamics. To actually evaluate the average power P and heat currents J ν at periodic steady state, we resort to a non-equilibrium Green function formalism [47,[60][61][62][63]. Deferring all details to Appendix B, here we quote the final results for the power and the heat current via bath ν = 2 (recall that J 1 = −P − J 2 ) where the influence kernels have the form induced by the specific driving scheme considered in this work -see Eq. (6). We close noting that the Floquet coefficients possess only even µ components and have the following symmetry properties:G * µ (ω) =G −µ (−ω) andG µ (ω − µ 2 Ω) =G −µ (ω + µ 2 Ω) -see Appendix. C for the latter.

III. DYNAMICAL HEAT ENGINE VERSUS MARKOVIANITY
In general, a structured environment can induce memory effects and non-Markovian dynamics. To assess non-Markovianity, several estimators have been introduced recently [50,[64][65][66]. As suggested in Refs. [65,67], a proper criterion to quantify non-Markovianity in the asymptotic regime is the violation of divisibility property of the dynamical map. Moreover, a direct link between the non-divisibility notion of non-Markovianity and the form of the bath spectral density has been given. In particular, it is possible to show [65,67] that only a strictly Ohmic [68] spectral density J (ω) = mγω in the classic (high temperature) regime leads to a separable map, hence to a Markovian dynamics -see Appendix D for details.
We now inspect general and necessary conditions to reach a heat engine regime, that is P < 0. By studying Eq. (24), as discussed in Appendix B, one can realize that J 1 (ω), i.e. the bath spectral density linked to the driven contact, dictates the sign of the different contributions of the average power.
We begin considering a strictly Ohmic spectral function in the high temperature regime. We now show that this choice, which represents a Markovian dynamics at high temperature, cannot support a heat engine at any order in the system/bath interaction in the high temperature regime. As shown in Appendix B, the average power P can be decomposed as Inserting Eq. (26) into Eqs. (B9), (B12) and (B14) an explicit expression for the above three quantities is obtained. We thus arrive at and It is now easy to see that P (b,2) ≥ 0. It is also possible to show that P (b,1) ≥ 0. Indeed, by using Eq. (B7) we first note that The remaining part of P (b,1) can then be rewritten, after some algebra, as which is also manifestly positive. This implies that P (b) ≥ 0, and that the possibility to have a heat engine only depends on the sign of P (a) . Let us now consider the Markovian limit of high temperatures T 1 → ∞. Here coth ω 2T1 → 2T1 ω , then using also the odd parity property of Im[G 0 (ω)] we find P (a) = 0.
This finally proves that in the Markovian regime P > 0 to every order in the system/bath coupling strength, demonstrating that in this case it is not possible to obtain a working heat engine (P < 0). Thus, non-Markovianity is a necessary condition to achieve a working heat engine.
It is worth to stress that the above argument holds true independently from the shape of J 2 (ω), related to the static bath.
It is now natural to wonder if this is also a sufficient condition. However, this is not the case, as we now argue by providing a counterexample. To this end, we consider the case in which the system/bath coupling strength of the modulated in time reservoir (ν = 1) is weak. In this perturbative regime, simpler closed expressions are obtained, from which one can also get useful physical intuitions. Up to linear order in J 1 (ω) the average heat power can be written as (see Appendix E for details, in particular Eq. (E5)) Here, we introduced χ 0 (ω) = −[ω 2 − ω 2 0 + iωγ 2 (ω)] −1 the bare susceptivity and where we recall that ω ± = ω ± Ω and where n B (x) = (e x − 1) −1 is the Bose distribution function. Note that J 2 (ω) only enters into the expression of χ 0 (ω) through γ 2 (ω) [69]. Recalling that J 2 (ω) = mω Re[γ 2 (ω)], one can realize that Im[χ 0 (ω)] is positive for ω > 0. Therefore, the regions of a working heat engine are given by To look for a counterexample, for the modulated bath we consider a monotonically increasing spectral density of the form J 1 (ω) = mγ 1 ω| ω ω | s−1 , which describes a large class of spectral function: Ohmic behavior for s = 1, sub-Ohmic for 0 < s < 1, and super-Ohmic for s > 1 [29,68]. In this case, the last sum that appears in square brackets on the second line of Eq. (34) is always negative, then the most favorable requirement for a working heat engine is in the T 2 → 0 limit, where the second line vanishes. Then, taking advantage of the relationL 1 (−ω) = 2J 1 (ω)n B (ω/T 1 ), withL 1 (ω) Fourier transform of the fluctuating force correlator L 1 (t), the condition for a working heat engine reads L 1 (−ω − Ω) >L 1 (−ω + Ω). In passing, this confirm the necessary condition of non-Markovianity, since a Markovian bath has L 1 (t) ∝ δ(t), and it never satisfies the above constraint [70].
More importantly, this relation can be used to obtain the counterexample we are looking for: Indeed, after some lengthy calculations reported in Appendix F, it is possible to demonstrate that in the non-Markovian case of sub-Ohmic spectral density no heat engine can be achieved. Therefore, we conclude that non-Markovianity is a necessary but not sufficient condition to obtain a heat engine.

IV. NON-MARKOVIAN ENGINE WITH A LORENTZIAN BATH
Having established the importance of non-Markovianity for a dynamical heat engine, we now characterize its performance. To this end, hereafter we choose a strictly Ohmic [68] spectral density J 2 (ω) = mγ 2 ω for the ν = 2 static bath. For the modulated bath of interest, we focus on a paradigmatic example of structured non-Markovian environment, i.e. a Lorentzian spectral function [31,32,[71][72][73] with a peak centered at ω 1 , an amplitude governed by d 1 , and a width determined by γ 1 (parameter linked to the damping). Such environment can be physically realized with cavity architectures [74][75][76][77] (see the sketch in Fig.  2(a)). For instance, in cavity optomechanics [74] a mechanical oscillator (the QHO of frequency ω 0 ) is embedded in a optical cavity. Then, a laser detuning is imposed on the bare cavity, in order to adjust the resulting frequency resonance ω 1 close to ω 0 , like in sideband-resolved cooling experiments [78]. A possible way to implement a temporal modulation of the system/bath coupling g 1 (t) is to superimpose a modulation (by optical pulse or mechanical vibration) to one of the cavity mirrors of characteristic frequency Ω. For clarity we now introduce the dimensionless parameter κ ≡ d 1 /(ω 2 0 ω 2 1 ), which governs the coupling strength.
A. Performance in the weak coupling regime (κ 1) In the weak coupling regime (κ 1), one always finds the possibility for a heat engine. Indeed, under the assumption of a sufficiently sharp peak around ω 1 and looking at Eq. (34), one can argue that the dominant contributions are for either ω + ω 1 or ω − ω 1 . In the former case, when one can drop the off-resonance contribution due to J 1 (ω − ), to have P < 0 one needs while when ω − ≈ ω 1 with analogous reasonings one obtains the condition show the corresponding power and the efficiency η ≡ −P/J 1 for a representative temperature arrangement, T 1 /T 2 = 0.1. We note that, having chosen γ 2 ω 0 , one finds Im χ 0 (ω) peaked around the QHO frequency ω 0 . The behavior of the average power in panel (b) well agrees with the resonance condition ω 1 ω 0 −Ω, along which the maximum power occurs (see the dashed line in the plots) [79]. With regard to engine efficiency, using the results shown in Appendix F (see Eqs. (F2)-(F4)) along this resonance one has that P ≈ P (−1) and and from the expressions quoted above one immediately finds where η C = 1 − T 1 /T 2 is the efficiency of the Carnot machine. Indeed, from Eq. (38) one finds ω1 ω1+Ω > T1 T2 , showing that the efficiency is upper bound by the Carnot limit which can therefore be approached in a realistic parameter window. Note, however, that maximum power and maximum efficiency are reached at different values of ω 1 . In particular, inspecting Eq. (F2) one immediately sees that when ω1 ω1+Ω = T1 T2 the power vanishes.
Within the model of a cavity of Fig. 2(a), one can intuitively interpret the produced power as the unbalance between the energy flux flowing from the contact at temperature T 2 > T 1 and the energy flux that can be absorbed by the red-shifted cavity (ω 1 = ω 0 − Ω) at the colder temperature T 1 .
Similar results are found for the resonance condition ω + ω 1 for T 1 > T 2 . Following steps perfectly analogous to the ones outlined above, one finds the efficiency to be where now η C = 1 − T 2 /T 1 . Also in this case the Carnot limit is achieved when P → 0, when ω1−Ω ω1 = T2 T1 .

B. Effective model in terms of quantum otto engines
The above results can be interpreted in terms of an effective model in which the QHO is regarded as a thermodynamic substance performing a quantum Otto cycle [80,81]. To illustrate this fact observe that, in the weak coupling regime, the expressions for the power and heat currents are sums over two independent "channels" labeled by p = ±1 -see the last identifications in Eqs. (F2)-(F4). For a sharp Lorentzian spectral density, the resonance condition , which allows to focus only the p-th channel ignoring the negligible contribution of the channel −p, i.e. to treat the two channels separately.
Let us consider now the channel p = +1, i.e. the resonance ω 1 = ω 0 + Ω with T 1 > T 2 . We can build an effective model in terms of a Otto engine cycle C + represented schematically in Fig. 3 and composed as follows: • an isochoric transformation A→B along which the QHO is kept at frequency ω 0 + Ω and allowed to exchange heat with the hot bath at temperature T 1 over a characteristic time τ 1 ; • an isentropic expansion B→C, where the oscillator is decoupled from the baths and its frequency evolves adiabatically from ω 0 + Ω to ω 0 ; • an isochoric transformation C→D at frequency ω 0 exchanging heat with the cold bath at temperature T 2 over a characteristic time τ 2 ; • an isentropic compression D→A, where the frequency adiabatically turns back from ω 0 to ω 0 + Ω. .
The characteristic time spent along each isentropic branch is τ is > ω −1 0 , which we assume large enough so that average occupation number of the QHO (denoted here as N ) is conserved [80]: N B = N C , and N D = N A . Also, we assume that at the end of each isochor the QHO has thermalized to the corresponding bath: = N A . Thermalization along the isochors takes finite characteristic times τ 1,2 which, to lowest order, can be identified with the inverse rates [28] The heat exchanged with the contacts are given by [80] To complete the mapping we remind that we deal with a perturbative regime for the spectral density J 1 and thus τ 1 τ 2 , and that to obtain Eqs. (F2)-(F4) the regime γ 2 ω 0 has been considered -see Appendix E, in particular the argument leading to Eq. (F1) -so that τ 2 τ is . Thus the total time spent on the cycle is ≈ τ 1 which allows to estimate the average heat currents as Q ν /τ 1 ≡ J (+1) ν and the average power as W/τ 1 ≡ P (+1) -see Eq. (F2).
For the cycle to operate as a heat engine one needs W < 0 which implies as also found previously. The efficiency η + = −P (+1) J (+1) 1 of C + is given by is governed by the compression ratio of the QHO and is in accordance with the physics of a Otto cycle and with the results quoted in the previous section. If With similar arguments one can interpret the resonance ω 1 = ω 0 − Ω (channel p = −1), where the cycle C − is dominant. It is comprised of two isentrops operating between the frequencies ω 0 and ω 0 − Ω and two isochors where the QHO is kept, with fixed frequency ω 0 (or ω 0 − Ω), in contact with the hot bath at temperature T 2 (or the cold bath at temperature T 1 ). Identifying the time spent in contact with the bath at T 1 as τ ≈ 4mω0 J1(ω0−Ω) , which is also the longest time in the cycle, allows to evaluate the heat currents and the power with reasonings similar to those made for C + . The expressions again coincide with J For C − to operate as an engine, one needs in accordance to what discussed in the previous section.
In the heat engine regime, the efficiency η − = −P (−1) Also this result agrees with the ones reported in the previous section.
An important remark must be made: the effective model leading to C − breaks down when Ω → ω 0 : here the effective volume of the QHO and the cycle time diverge and correspondingly P (−1) → 0. When Ω > ω 0 , one can see that P (−1) > 0 and J (−1) ν < 0: Now the p = −1 channel acts as a "heater" absorbing work from the driving mechanism and discharging it into both baths [82].
The analysis conducted above can also be applied to spectral densities other than the Lorentzian one. However, when J 1 (ω) is not sharply peaked the contributions of the two channels p = ±1 cannot be clearly separated. In this case one can interpret the results as the action of two thermal machines running in parallel: heat currents J ν through the baths split/recombine into the two channels J (p) ν and the total power P is the net sum of the powers P (p) exchanged by each machine. When Ω < ω 0 the two machines perform the Otto cycles C ± discussed above, while for Ω > ω 0 only the channel p = +1 behaves as a Otto cycle, while p = −1 acts as a heater. From the discussion above it is clear that for given Ω and temperature ratio T 2 /T 1 only at most one channel can operate as a heat engine. Therefore, either T2 T1 > ω0 ω0−Ω or T2 T1 < ω0 ω0+Ω is a necessary but not sufficient condition for obtaining P < 0 and the precise balance between the power exchanged by the two channels must be studied [83].

C. Beyond the weak coupling regime
We conclude this section studying a regime beyond weak coupling. In Figs. 4(a,b) we show heat engine performance, obtained by numerically solving Eq. (15) and evaluating the expressions for average power in Eq. (24) and corresponding heat currents [84]. A clear broadening of the power resonance line can be observed (panel (a)). Indeed, the maximum power is no longer achieved along the "bare" resonance ω 1 = ω 0 − Ω (see the dashed line) but now appears detuned. Also the efficiency, reported in panel (b), displays an overall broadening.
To further inspect the strong-coupling regime we also show in Fig. 4(c,d) the average power and efficiency as a function of the coupling parameter κ and the driving frequency Ω near the maximum of Fig. 4(a). Operating the engine beyond weak coupling allows to achieve sensibly larger power outputs. Indeed, the maximum power occurs at moderate/strong coupling κ ≈ 0.2 and is over 10 times larger than the maximum power obtained at weak coupling κ = 10 −3 . However, the average power is a non-monotonic function of κ and for κ > ∼ 1 the heat engine is lost. Looking at Fig. 4(d), the maximum efficiency is achieved operating the engine in the weak coupling regime κ 1. Therefore, the parameter κ can be used to tune the trade-off between power and efficiency. As a final remark, we note that stronger coupling strengths induce a marked detuning of the resonance frequency, due to an energy renormalization which can be captured by solving self-consistently ω 2 − ω 2 0 − Re k 0 (ω) = 0 and ω 1 = ω − Ω. The solution is shown as a dashed line in Fig. 4c).

V. CONCLUSIONS AND OUTLOOK
We have shown that by properly modulating the system-bath coupling it is possible to obtain a working heat engine. Here, non-Markovianity is a useful resource for quantum thermodynamics, in that it allows for efficient dynamical heat engines, even approaching Carnot efficiency. Our results open up new possibilities for the exploitation of non-Markovianity. For instance, one could consider the combined effect of modulating couplings and driving the system, looking for a cooperative effect to enhance the performance of thermal machines. In such a quest, machine learning tools [85][86][87][88] could prove to be useful. Regarding possible implementations, our results on structured environment could be tested in the field of cavity optomechanics, which are emerging as an interesting platform for new quantum technologies [74,77,89,90]. We believe that these findings are not restricted to the investigated working medium of a QHO, and an interesting follow-up would be to investigate the role of non-Markovian contributions with different quantum system, for instance one or more qubits, that can be integrated in superconducting waveguide quantum electrodynamics architectures [75,76,[91][92][93].
Appendix A: Proof that A = 0 Here we prove that A convenient strategy is to rewrite the above quantity as Below, we will show that H (t) int (t) and H QHO (t) are periodic functions with period T , which eventually implies A = 0. The time dependence of the interaction term is given by Using the equations of motion for the position x(t) and momentum p(t) of the QHȮ withẍ(t) given in Eq. (8) and those for the position and momentum of the degrees of freedom of the bathṡ X k,ν (t) = P k,ν (t) m k,ν , P k,ν (t) = −m k,ν ω 2 k,ν X k,ν (t) + g ν (t)c k,ν x(t) we can rewrite H (t) int (t) in terms of the system position operators x(t) alone, as Similarly, for the QHO term in Eq. (1) we have The above expressions show that their quantum ensemble averages can be written as correlators of the QHO position operator only. Indeed we have Note that all these correlators can be represented in terms of the function as: We then focus on the evaluation of the quantum average of M (t, s). Plugging into Eq. (A1) the time evolution given in Eq. (13) we obtain where in the second equality we have inserted the bath correlator ξ ν (t)ξ ν (t ) = δ ν,ν L ν (t − t ) with L ν (t − t ) given in Eq. (11). The time integrals are performed using the Fourier representations in Eq. (14) for the Green functions G(t, t 1 ) and G(s, t 2 ), Eq. (5) for the driving g ν (t 1 ) and g ν (t 2 ), and for the bath correlator. We finally obtain Notice that now the times t and s only appear in the exponential factors: this means that after derivatives with respect to t and s as appropriate according to the definitions in Eqs. (A2)-(A4), the limit s → t implies a time dependent term always of the form e −iΩt(n1+n2+µ1+µ2) . This shows a clear periodicity with respect to the cycle time T . This result demonstrates that the required correlators are indeed periodic: M xx (t + T ) = M xx (t), M xẍ (t+T ) = M xẍ (t) and Mẋẋ(t+T ) = Mẋẋ(t). Since, as shown above, these correlators are the building blocks of H (t) int (t) and H QHO (t) we eventually arrive to the conclusion that the quantum averages of the interaction term and of the QHO part are periodic. This finally proves the key result that A = 0 and then that the cycling average of the total power is totally balanced only by the reservoir heat currents -see Eq. (23).

Appendix B: Average power and heat currents for a monochromatic drive
Here we derive a compact form for the power and heat currents by focusing on the case discussed in the main part, namely a constant coupling to the bath ν = 2 and a monochromatically modulated coupling to the bath ν = 1, see Eq. (6). It is useful to begin recalling the general expressions of the average total power P and heat currents J ν performed using Eqs. (20) and (21) as explained in Ref. [47]. We have: where n tot = n 1 + n 2 + n 3 + n 4 . We remind that the spectral densities J ν (ω) are odd functions of frequency. We now specialize to the driving considered in this work, with Fourier coefficients g n,1 = (δ n,1 + δ n,−1 )/2 and g n,2 = δ n,0 . Then the kernel in Eq. (17) reduces tõ with ω ± = ω ± Ω. This shows that only the kernels k 0,±2 (ω) are different from zero. We now plug the ex-pressions of g 0,ν into Eq. (B1) to write down the average power. Notice that only the ν = 1 term contributes to P , that can be conveniently decomposed into two contributions P = P (a) + P (b) . The first, stemming from the first line of Eq. (B1) reads while the second originates from the second line of the same equation and it is given by To obtain Eq. (B5), we have used that n 1 , n 2 = ±1, and the property of the Floquet coefficients Exploiting the symmetry property (see Appendix C for details)G one sees that the last term in the square brackets of Eq. (B5) has a null contribution upon integration giving then Using Eq. (B7) and renaming µ + n 2 + n 4 → µ we can then rewrite the P (b) term as We now insert the explicit form of g n,ν , and for notational convenience we write P (b) = P (b,1) +P (b,2) corresponding to the ν 1 = 1, 2 terms in the above expression. The former contribution reads Taking the real part of the sum in the last square bracket and using again Eq. (B7), one has Note that since J 1 (ω) coth ω 2T1 > 0 and since the last factor in the summand is positive, the sign of P (b,1) is governed by the behavior of J 1 (ω − µΩ) only. We now consider P (b,2) related to the ν 1 = 2 contribution in Eq. (B6), where n 1 , n 2 = ±1 and n 3 = n 4 = 0: Performing the sum over n 1 , n 2 , we arrive at Also here, since J 2 (ω) coth ω 2T2 > 0 and since the first factor in the summand is positive, the sign of P (b,2) is determined by the terms J 1 [ω + (µ ± 1)Ω]. Finally, we sum the three contributions in Eqs. (B9), (B12), (B14) to obtain the expression for the average total power P = P (a) + P (b,1) + P (b,2) as reported in Eq. (24).

Appendix C: A useful property of the Floquet coefficients
Here we prove that, for the dynamical couplings considered in this work, We begin by recalling the expression in Eq. (B4) for k ±2 (ω). In addition, note that for a generic bath on contact 1 withγ 1 (ω) = γ 1 φ(ω), Eq. (15) can be conveniently rewritten as where we have introduced We now write a formal series expansion ofG µ (ω) in powers of λ:G which is plugged into Eq. (C1). Matching order-by-order in λ a hierarchy of nested equations for the n-th contribu-tionG (n) µ (ω) is obtained. In particular one immediately sees thatG (0) 0 (ω) = D 0 (ω) and that From Eq. (C3) one can conclude that: (1)G  One can picture the set ofG (n) µ (ω) satisfying n ≥ 0 and |µ| ≤ 2n as a lattice of dots on a Pascal triangle, whose rows are labeled by n and whose columns are labeled by µ. This is represented in Fig. 5, where green (red) dots represent the non-zero (zero)G (n) µ (ω). Equations (C3) can be solved recursively. As an example, to first order one immediately finds G (1) ±2 (ω) = D ±2 (ω)J ±2,0 (ω)D 0 (ω). As another example, a nontrivial solution for the second order is Each term in the above equation can be interpreted as a path on the Pascal triangle, linking dots (µ, n) and (µ, n ) within the triangle according to the simple rule |µ − µ | = 2 and |n − n | = 1. To each dot a factor D µ (ω) is associated, to each link between dots a factor J µ,µ (ω) is associated. Explicitly, the two paths representing the terms in Eq. (C4) are (0, 2) → (2, 1) → (0, 0) and (0, 2) → (−2, 1) → (0, 0) respectively. Notice that the index n to the l.h.s. of gives the number n of links between the n + 1 dots.
Proceeding with the recursion one quickly realizes that the situation depicted above is general. Indeed, the term G terms whereG (n,j) µ (ω) is associated to one of all the N (n, µ) distinct paths P j (µ, n) (with 1 ≤ j ≤ N (n, µ)) that connect the dot (µ, n) with (0, 0) with n links that follow the rules |µ − µ | = 2 and |n − n | = 1 as stated above.
More formally, each path P j (µ, n) associated tõ G (n,j) µ (ω) can be represented by the sequence of dots n ≡ µ and µ (j) 0 = 0. One of such path for µ = 2 and n = 5 is shown as a blue line in Fig. 5. Then, to constructG (n,j) µ (ω) we associate a term D µν (ω) to each dot in the path and a term J µν ,µν+1 (ω) to each link between consecutive dots (with 0 ≤ ν < n) which leads tõ The set of all paths P(μ, n) = ∪ N (n,μ) j=1 P j contributing to Eq. (C5) lies within a rectangular region of the Pascal triangle (a representative example for the case µ = 2 and n = 5 is shown as the yellow region in Fig. 5). It is simple    to see that such rectangle has vertices To prove Eq. (B8) for the n-th order term we need to shift the argument ofG (n) µ (ω). To this end, it is useful to observe that D µ (ω + kΩ) = D µ+k (ω) and J µ,µ (ω + kΩ) = J µ+k,µ +k (ω), with k an integer. Geometrically, this means that shifting the argument of G (n) µ (ω) by kΩ is equivalent to shift all paths [94] that contribute to it (and hence the whole set P(μ, n)) by |k|Ω to the right or to the left according to Sgn(k).
According to what discussed above, let us denote with P + the region containing all the paths contributing tõ G (n) µ ω − µ 2 Ω , with vertices while the paths contributing toG (n) −2µ (ω + µΩ) belong to the region P − with vertices Observe that all the factors in Eq. (C6) actually depend only on the ordered set {µ (j) ν } but are invariant under any permutation of the second coordinate of each lattice point. This allows to re-order the vertices of P − in decreasing order of their second coordinate as This allows us to conclude that the two regions are actually identical. Since shifting the argument only amounts to a rigid translation of the paths and the actual reordering performed above corresponds to reading each term from right to left rather than left to right, it follows that to each path ofG (n) µ ω − µ 2 Ω identically corresponds one and only one term ofG 2 Ω . This allows to conclude thatG −µ ω + µ 2 Ω is valid ∀n. By virtue of the series expansion in Eq. (C2), the above property is valid also for the completeG µ (ω):

Appendix D: Non-Markovianity criterion
Memory effects and non-Markovian dynamics are the subject of many studies, and several notions of non-Markovianity have been introduced recently (see the review in Ref. 50 and references therein). Different estimators have been investigated to witness and to quantify non-Markovianity, and these not always are equivalent [50]. Moreover, many criteria, such as the ones related to the trace distance, are not well-suited in the study of asymptotic states, i.e. looking at properties in the long time limit (like in our case of interest). There, different estimators should be used, as discussed in Refs. [64][65][66][67]. In particular, it has been shown that non-Markovianity can be assessed through the violation of the divisibility condition: the evolution is Markovian if and only if it is described by a divisible, completely positive map. Importantly, this criterion is valid for both finite and asymptotic times. In Refs. [65][66][67] this criterion has been applied to the case of a QHO coupled to a thermal bath at temperature T , exploiting the exact solution derived in Ref. [56], and the divisibility estimator has been linked directly to the form of the bath spectral density J (ω). We now recall the definition of this estimator, details on the derivation can be found in Refs. [65,67]. The punctual non-Markovianity measure is given by where Γ(t), ∆(t), and Π(t) are the damping, direct and anomalous diffusion coefficients, respectively, and are completely determined by the form of the bath spectral density J (ω). Notice that we have indicated the damping coefficient with Γ(t), instead of γ(t) used in the original paper [67] to avoid confusion with the quantity introduced in the main text. Equation (D1) is valid at any time, and in particular also in the asymptotic t → ∞ regime. It is worth to note that this quantity is bounded 0 ≤ N p (t) ≤ 1. If N p (t) = 0, one has a Markovian dynamics, hence described by a separable dynamical map. On the other hand, if N p (t) > 0 the divisibility criterion is violated and the dynamics is non-Markovian. In general, the above three time-dependent coefficients are given by cumbersome integral expressions. However, at weak coupling the above coefficients can be evaluated in closed form, and in the asymptotic regime they read [65,67] and where the integral is taken as a principal value. We underline that, since ∆(∞) is positive definite, in the asymptotic regime N p (∞) is bounded between 0 (Markovian regime) and the maximal non-Markovianity 1/2. Hereafter, we analyze the behavior of N p (∞) for the case of baths with a spectral density ∝ ω s (s > 0) or with a structured Lorentzian spectral density, which are the examples discussed in the main text. We will show in particular that only the strictly Ohmic spectral density at high temperature (classical regime) leads to N p (∞) = 0 and thus to a Markovian dynamics.
We begin considering the spectral density where θ(x) is the Heaviside theta function, which models a sub-Ohmic (0 < s < 1), Ohmic (s = 1) or super-Ohmic (s > 1) bath with a hard cutoff ω c . Figure 6 that only for ω c ω 0 and T ω 0 one has N p → 0. Indeed, setting T ω 0 and evaluating the principal value in Eq. (D3) one finds Π(∞) ≡ 0, while in the same regime Γ(∞)/∆(∞) = ω0 2T → 0, which proves that N p ∝ ω T 2 → 0. Away from this regime deviations from the Markovianity occur, particularly in the T < ω 0 (quantum) regime. It is also worth investigating the suband super-Ohmic cases for large cutoff (ω c ω 0 ). The results are summarized in Fig. 6(b) where N p is shown as a function of s for three different temperatures. In all cases except the high-temperature Ohmic one, the dynamics is non-Markovian.
We now turn to the case of a Lorentzian spectral density with a hard cut-off discussed in Sec.IV in the large cut-off limit. Figure 7(a) shows the non-Markovianity measure as a function of the cutoff and the damping parameter γ 1 for typical values of ω 1 and T . In this case we found for N p a minimum value ≈ 0.36 for a broad Lorentzian peak (γ 1 ≈ ω 0 ) which increases to the maximum N p = 1/2 when the peak is very sharp (γ 1 ω 0 ), signaling a distinctly non-Markovian dynamics. Also, it is worth to notice that N p is essentially insensitive to the cutoff when ω c > ω 1 . Figure 7(b) confirms that the non-Markovian dynamics is stable against thermal effects.
when Ω → ω 0 . At least for Ω ≈ ω 0 and 0 < s ≤ 1 it is then clear that f s (Ω) ≤ g T1 (Ω) and thus no working engine can be obtained there. To prove that this is the case for any Ω we now inspect the general properties of f s (Ω) and g T1 (Ω), and their derivatives, to show that for 0 < s ≤ 1 Eq. (F6) cannot be verified. Observe that df s (Ω) dΩ = s 2ω 0 ω 2 0 − Ω 2 f s (Ω) for Ω = ω 0 , and where in the second passage we have used Eq. (F10). Thus we arrive at the following inequality Integrating Eq. (F14) from 0 to Ω one obtains f s (Ω) < g s T1 (Ω) or equivalently f s (Ω) g T1 (ω) < g s−1 T1 (Ω).