Dynamically Emergent Quantum Thermodynamics: Non-Markovian Otto Cycle

Employing a recently developed approach to dynamically emergent quantum thermodynamics, we revisit the thermodynamic behavior of the quantum Otto cycle with a focus on memory effects and strong system-bath couplings. Our investigation is based on an exact treatment of non-Markovianity by means of an exact quantum master equation, modelling the dynamics through the Fano-Anderson model featuring a peaked environmental spectral density. By comparing the results to the standard Markovian case, we find that non-Markovian baths can induce work transfer to the system, and identify specific parameter regions which lead to enhanced work output and efficiency of the cycle. In particular, we demonstrate that these improvements arise when the cycle operates in a frequency interval which contains the peak of the spectral density. This can be understood from an analysis of the renormalized frequencies emerging through the system-baths couplings.


I. INTRODUCTION
Quantum thermodynamics is currently a prolific and exciting research field.Technical advancements, the precision of controlling quantum systems in the lab, and the chase for the quantum advantage, have led to the fast-paced development of quantum thermal machines, insofar as to already bring experimental realization [1][2][3][4][5][6].Yet, because of this accelerated growth, the field is also plagued with inconsistencies: quantum engines and other applications are being theoretized, but essential foundations of the field's theory are still under debate.
While the formulation and fundamentals of quantum thermodynamics seem to be well understood for ultraweak system-bath coupling and continuous information loss, the community is still not converging on basic definitions for thermodynamic quantities outside of this regime.There are many aspects that make the quest challenging, but we can perhaps identify two major contributors.First, the interaction energy between a quantum system and its environment becomes non-negligible at stronger coupling parameters.It therefore needs to be understood how much of this energy is accessible from the system and should be counted as part of the system internal energy, and whether (or how) it counts as work or heat.The second aspect is the complex information exchange between system and environment.It is now clear that information theory and thermodynamics share a deeper link and that irreversible behaviour is connected to information loss; thus, memory effects and information backflow from the bath to the system should also significantly influence the theory [7].
In the literature one may find numerous theoretical proposals addressing these questions [8][9][10][11][12][13][14][15][16].However, few of these disparate approaches made it so far as to enter the stage of theoretical applications.Perhaps the most popular strong coupling treatment of quantum thermodynamics is given by identifying the heat exchanged as minus the energy change of the bath, as proposed by Esposito, Lindenberg and Van den Broeck (ELB) [9].This corresponds to assigning the interaction energy fully to the system's internal energy, in the original formulation.Later treatments have instead kept the definition of internal energy as in the weak coupling formulation (i.e.excluding the interaction energy) while adopting the ELB definition for heat, imposing some work exchange ad hoc to maintain a first law of thermodynamics [16].This definition of heat exchange seems therefore too arbitrary a choice; moreover, it does not seem to reflect the impact of memory effects, which could arise from finite heat baths or structured spectral densities.The oscillations of entropy production according to this definition, i.e. temporary violations of the second law, seem in fact to have no relation with non-Markovian behavior [17].Further, this definition of heat exchange has been found to be incompatible with the weak coupling formulation in the limit of vanishing coupling parameters when the system is driven by an external time-dependent force [18].
Already a vast quantity of studies emerged in very recent years, that are concerned with investigating strong coupling and/or non-Markovian effects on a quantum Otto cycle [19][20][21][22][23][24][25][26][27][28][29].While these works explore the impact of different causes of non-Markovianity and find an array of interesting consequences, most of them either employ the weak-coupling treatment of quantum thermodynamics or the ELB approach within an exact treatment.Because of the reasons explained above, we believe some features might have been missed.We hope to help fill this gap by properly accounting for strong coupling and non-Markovian effects, starting by employing more suitable definitions of the main thermodynamic quantities.
We have recently proposed a different approach to the thermodynamic laws that is tailored to exactly resolve the consequences of coupling a system to an arbitrary reservoir [30].This arbitrariness can refer to the coupling strength, to the protocol for switching on and off the interaction, the finiteness of the bath or to its structural properties, including its initial state.This approach particularly was shown to link violations of the second law with information backflow, and thus represents in our opinion the best formulation in which to study the properties of a non-Markovian thermal machine.In Ref. [30] the thermodynamic quantities are given in terms of an effective system Hamiltonian K S emerging from the dynamics of the system when it is coupled with the environment.This emergent Hamiltonian contains details of the environment and total evolution, and describes how the system energy levels react to the interaction with the environment.Accordingly, the internal energy of the system can change either due to a change in the energy levels or eigenstates of the effective Hamiltonian (work), or due to a change of the populations (heat).It turns out that non-Markovian environments can in general induce a time-dependent effective Hamiltonian K S (t) even when the total system (system+environment) is an isolated system (following a unitary evolution according to a time-independent Hamiltonian) and even for weak coupling parameters.This reveals a major characteristic of non-Markovian baths: They can perform work on the system, which is interpreted as useful, coherent energy granted to (or taken away from) the system as a result of non-trivial information exchange with the environment.This can of course lead to crucial changes in the treatment of heat engines and is the main motivation behind this work.
The focus of this paper is the variation in work output and efficiency of a heat engine when the baths can be strongly coupled and exhibit memory effects.We do this by comparing the standard Otto cycle -a harmonic oscillator driven and coupled to Markovian baths [31] -with the corresponding cycle where the system-baths couplings are taken into account within an exact non-Markovian treatment.To this end, the system-baths dynamics are described by the Fano-Anderson model [32][33][34] featuring a structured spectral density.We note that in Refs.[35,36] an approach to quantum thermodynamics of the Fano-Anderson model has been developed which is closely related to our approach.
The paper is organized as follows.In Sec.II we collect the theoretical background needed for this study: We summarize in Sec.II A the approach for quantum thermodynamics developed in [30] and present some details of the microscopic model and its solution in Sec.II B. In Sec.II C, we show the theoretical results for the Otto cycle where the system is connected via arbitrary coupling parameters to two thermal baths of harmonic oscillators with peaked Lorentzian spectral density.Since the coupling with these baths leads in general to both heat and work exchange, we factor this in for the engine's work output and efficiency.In Sec.III we show numerical results in different parameter regimes and find that if the bath spectral density peaks at a frequency within the driving range of the system, then the cycle can be enhanced, showing higher efficiencies than its standard counterpart.Finally, we draw our conclusions in Sec.IV.

A. Nonequilibrium thermodynamics at strong couplings
We consider the nonequilibrium quantum thermodynamics of an open system S coupled to an environment E representing the heat baths.The quantum states of S and E are given by density matrices ρ S and ρ E , respectively, while states of the total composite system S + E are denoted by ρ SE .The microscopic Hamiltonian is taken to be of the form where H S and H E represent the system and the environmental Hamiltonian, respectively, and H I denotes the interaction Hamiltonian.Note that the system Hamiltonian and the interaction Hamiltonian may be time dependent, considering, e.g., an external driving of the system or a switching on and off of the system-environment interaction.
Our approach to nonequilibrium thermodynamics at arbitrary system-environment coupling developed in [30] is based on the following master equation for the density matrix ρ S (t) of the open system S: This is an exact time-local differential equation with an explicitly time-dependent generator L t which completely describes all memory effects of the open system dynamics [37,38].Since there is no time-convolution in this master equation, it is also known as time-convolutionless master equation.The generator of the master equation (2) consists of two parts.The first one is a Hamiltonian contribution given by the commutator with some effective Hamiltonian K S (t), while the second contribution is provided by the dissipator (3) with certain decay rates γ k (t) and corresponding Lindblad operators L k (t), which are in general timedependent.The above general structure of the time-local master equation can be derived from the requirements of the preservation of Hermiticity and trace of the density matrix [39,40].
The explicit form of the rate functions γ k (t) and of the Lindblad operators L k (t) for a given microscopic Hamiltonian can be found by means of a perturbation expansion in the system-environment coupling by applying the time-convolutionless projection operator technique [37,38,41].This expansion can also be viewed as an expansion of the generator in terms of so-called ordered cumulants [41][42][43].In this paper, we will investigate an analytically solvable system-reservoir model for which the exact master equation ( 2) is known [44][45][46][47] and can be constructed from the solution of the full Heisenberg equations of motion (see Sec. II B).
Given the master equation ( 2) it seems natural to interpret the Hermitian operator appearing in the commutator part as an effective system Hamiltonian.However, this interpretation has to be taken with care, because the splitting of the generator L t into Hamiltonian and dissipative part is not unique.In other words, one can modify the commutator part and, at the same time, the dissipator part in (2) in an infinite number of different ways without changing the generator itself and, hence, without changing the dynamics of the system.To fix a unique splitting of the generator into Hamiltonian contribution and dissipator we have proposed a principle of minimal dissipation [30] which is based on the approach developed in [48].This principle requires that the dissipator of the generator must be minimal with respect to a certain superoperator norm.Physically this means that the dissipative contribution of the master equation corresponds to that part of the total energy flow which cannot be interpreted as mechanical work due to a timedependent effective Hamiltonian of the system and, thus, must describe the flow of heat energy.As shown in [48] the dissipator is minimal if and only if the Lindblad generators L k (t) are traceless, a condition which leads to a unique expression for the effective system Hamiltonian K S (t).
In accordance with these considerations we define the internal energy of the system by the expectation value of the effective Hamiltonian, and formulate the first law on the change of internal energy as where we have defined work and heat contributions by Note that according to the master equation ( 2) the heat contribution Q S (t) can be expressed in terms of the dissipator of the master equation as

B. Microscopic model and effective Hamiltonian
Our goal is a study of the thermodynamics of the Otto cycle in which the full non-Markovian behavior of the system-baths coupling is taken into account, including all memory effects due to strong couplings and structured spectral densities.To keep the analysis as simple as possible we consider an analytically solvable model known as Fano-Anderson model [32][33][34], which can be described by the microscopic Hamiltonian where a † , a are the creation and annihilation operators of the open system mode with frequency ω 0 , and c † j , c j are the creation and annihilation operators of the environmental modes with frequencies ω j , satisfying the commutation relations [c j , c † k ] = δ jk .The interaction Hamiltonian H I describes the coupling between the system mode a and the environmental modes c j with corresponding coupling constants g j .As the total Hamiltonian is quadratic, the model can easily be solved analytically to obtain the reduced dynamics of the system mode a.To this end, we define the memory kernel which can also be written as where we have introduced the spectral density which will later be replaced by a continuous function.
Defining further the Green function G(t) as the solution of the integro-differential equation corresponding to the initial value G(0) = 1, we obtain the following exact solution of the Heisenberg equations of motion for the system mode, Here, a(0) = a represents the Schrödinger picture operator and we have defined which involves the mode operators c j (0) = c j of the environment in the Schrödinger picture.Employing Eq. ( 14) we can determine all relevant moments of the system mode a and also an exact time-local master equation for the open system's density matrix ρ S (t).For simplicity we assume a factorizing total initial state [49] ρ SE (0) = ρ S (0) ⊗ ρ E (0), (16) where the initial environmental state represents a Gibbs state with Z E the partition function and β = 1/k B T the inverse temperature.It follows from Eqs. ( 16) and ( 15) that the first and second moments of c(t) vanish, and, moreover, that the cross correlations between a(0) and c(t) are zero, With the help of ( 14) we then find the following expressions for the first and second moment of the open system mode, Here, we use the notation ⟨X⟩ t = Tr{Xρ S (t)} for any system operator X and define the noise integral I(t) through (23) From the expressions (20)- (22) for the moments of mode a one can derive an exact time-convolutionless master equation for the open system density matrix [44][45][46][47] which is of the form given in Eq. (2) (see Appendix A for details), where the effective Hamiltonian takes the form and the dissipator reads Thus we see that the effective Hamiltonian (24) represents a harmonic oscillator with a renormalized frequency which is given by The quantity γ(t) plays the role of a generalized decay rate defined by and, finally, the quantity N (t) is obtained from Note that we have written the effective Hamiltonian (24) and the dissipator (25) in a form which clearly reveals the Born-Markovian approximation of the exact master equation.In fact, within second order in the systemenvironment coupling one obtains the time-independent Markovian decay rate and the time-independent renormalized frequency while N (t) has to be replaced by the Planck distribution N 0 = 1/(e βω0 − 1).These replacements directly lead to the well-known Born-Markov approximation for the quantum master equation in Lindblad form: As is well known [41] a renormalized frequency ω M r even appears within the Born-Markov approximation, for which the frequency shift ∆ω(t → ∞) is given by the principal value integral in Eq. (30).However, when applying the master equation ( 31) e.g. in quantum optics the frequency shift is often very small for weak couplings and, hence, is commonly neglected completely.This is consistent in the sense that the stationary solution of the master equation ( 31) describes a Boltzmann distribution over the energy levels of the bare oscillator with frequency ω 0 .Therefore, in the following we define results of the Born-Markov approximation as results which are obtained by modelling the coupling to the baths by the master equation (31) with ω M r replaced by ω 0 .

C. Otto cycle
We study a quantum Otto cycle using a harmonic oscillator with time dependent frequency as working system.Our aim is to compare the results from the open system dynamics in the standard Born-Markov approximation with the exact approach in the non-Markovian regime.

Standard Otto cycle
The standard quantum Otto cycle [31] consists of four phases (see Fig. 1): I. Adiabatic compression: The frequency ω(t) of the harmonic oscillator is varied from the initial frequency ω 1 to the final frequency ω 2 > ω 1 .During this phase the system is decoupled from the baths such that the dynamics is unitary with Hamiltonian Consequently, there is no heat exchange, only work is performed on the system.For the Hamiltonian (32) the mean occupation number n = ⟨a † a⟩ is obviously a conserved quantity.At the initial time the oscillator is assumed to be in a thermal equilibrium state and, hence, we have where T 1 is the temperature of the cold bath and the index M serves to indicate that the quantity refers to the standard Born-Markov treatment.The work and heat transfer in phase I are thus given by Note that we are assuming here that a and a † are frequency independent.When starting from the expression 1 2 p 2 + 1 2 ω 2 (t)x 2 for the system Hamiltonian this assumption corresponds to fully adiabatic motion during this phase, which means that the adiabaticity parameter Q * is equal to 1 [50,51].We emphasize that this assumption is made here only to simplify the presentation since our main goal is a comparison between the standard Markovian and the non-Markovian treatment.
II. Isochoric heating: The oscillator is connected to a hot bath at temperature T 2 and evolves until thermalization is reached.The coupling to the bath is modelled through the Born-Markov master equation (31) which implies that the mean occupation number in thermal equilibrium becomes Note that the steady state does not depend on the initial conditions, but only on the temperature of the bath T 2 and on the frequency ω 2 .The frequency remains constant at ω 2 , such that all changes of internal energy corresponds to heat exchange Q H and that there is no work contribution.Hence, we have III. Adiabatic expansion: As in stroke I the oscillator is disconnected from the bath, and the frequency of the system is driven back from ω 2 down to ω 1 following a unitary dynamics given by the Hamiltonian (32).The occupation number remains constant at the value n M 2 .This leads to IV. Isochoric cooling: Finally, the oscillator is connected to a cold bath at temperature T 1 < T 2 , and evolves at constant frequency ω 1 until it thermalizes with mean occupation number n M 1 .Hence, there is only heat exchange with the cold bath and no work contribution which yields: We define the efficiency η of the cycle as the ratio of the total work output W out and the heat input Q H from the hot bath, when W out > 0, and as 0 otherwise.According to Eqs. ( 34)- (38) in the Markovian case the total work output takes the form while the heat input Q H from the hot bath is given in Eq. ( 36).Thus we obtain the well-known expression for the efficiency of the cycle within the Born-Markov approximation [31], Employing Eq. ( 40) we find that the condition on the frequencies and temperatures for the system to have net work output (W out > 0) and, thus, to operate as a heat engine, can be written as

Non-Markovian Otto cycle
In the following we will compare the exact approach, which is valid for arbitrary couplings and non-Markovian dynamics, with the standard treatment within the Born-Markovian approximation discussed in the previous section.The crucial differences to the standard approach are the time-dependent renormalization of the frequency of the system oscillator and the time dependence of the occupation number, as is visualized schematically in Fig. 2.
Here we assume the non-Markovian bath to be such that the system still thermalizes to a unique occupation number n 1,2 depending on the temperature T 1,2 of the bath, and to a unique renormalized frequency ω r 1,2 depending on the bare frequency ω 1,2 of the oscillator.The steady state also depends on the coupling and on the structure of the bath.The strokes of the cycle are then altered in the following way.

I. Adiabatic compression: This phase starts now at
frequency ω r 1 and occupation number n 1 which remains constant throughout the increase of the frequency up to ω 2 .Since the baths are decoupled we have again unitary evolution and only a work contribution which is given by Initially, the oscillator is assumed to be in a thermal equilibrium state at the temperature T 1 of the cold bath, i.e. we have II. Heating phase: The oscillator is coupled to the hot bath with temperature T 2 .The exact dynamics of the oscillator is given by the non-Markovian master equation (2) involving the effective Hamiltonian (24) and the dissipator (25), giving rise to a time dependent oscillator frequency, which results in a work contribution to the energy change.Hence, this is no longer an isochoric process, as there is generally both work and heat exchange between the system and the bath due to the non-Markovian behaviour.
The final frequency is ω r 2 = ω r 2 (t → ∞) and the final occupation number becomes Note that these quantities are generally different from ω 2 and n M 2 , respectively.Using Eqs. ( 6) and the effective Hamiltonian (24) we find the work contribution where the expressions for ω r 2 (t) and n 2 (t) = ⟨a † a⟩ T2 t are found in Eqs. ( 26) and ( 22).

III. Adiabatic expansion:
The oscillator is again disconnected from the baths, and the frequency is driven down to ω 1 , while the occupation number n 2 stays constant.The change of internal energy is fully due to work exchange which is given by with n 2 given by Eq. ( 45).
IV. Cooling phase: The oscillator is coupled to the cold heat bath with temperature T 1 which leads to relaxation to thermal equilibrium with frequency ω r 1 = ω r 1 (t → ∞) and occupation number n 1 given by Eq. (44).As in phase II we have both heat and work exchange, the latter given by The efficiency is again defined by (39).However, for the non-Markovian treatment the total work output consists of contributions from all four phases of the cycle: The expression for the heat input from the hot bath is given by but we can also make use of the first law of thermodynamics to determine Q H by means of In our numerical simulations we consider two baths at temperatures T 1 and T 2 with identical structured spectral densities described by a Lorentzian where λ represents the spectral width and ω c = ω 0 − ∆ describes the position of the peak with ω 0 the bare frequency of the oscillator and ∆ the detuning.We note that we have modified the Lorentzian shape of the spectral density in the regime of very low frequency modes in order for the noise integral (23) to be finite; see appendix B for details.Depending on the phase of the cycle we have ω 0 = ω 1,2 and, considering a fixed ω c , the detuning is ∆ = ∆ 2 = ω 2 − ω c for the heating phase, and ∆ = ∆ 1 = ω 1 − ω c for the cooling phase.Using Eq. ( 29) we see that the parameter γ 0 , which is used as a measure for the strength of the system-bath coupling, corresponds to the Markovian decay rate at resonance (∆ = 0).With the spectral density (52) Eq. ( 13) can easily be solved by Laplace transformation if one extends the range of the integral (11) to the entire real axis, leading to the solution where µ 1,2 are the roots of the quadratic equation For this particular case, the coupling to the bath induces a time-dependent renormalization of the oscillator frequency which, according to Eq. ( 26), is given by Using Eq. ( 53) one can show that in this model the system always relaxes to a unique steady state representing a thermal equilibrium state.Moreover, if we take the long time limit of (55) we find that ω r (t → ∞) > ω 0 for ∆ > 0 and ω r (t → ∞) < ω 0 for ∆ < 0. This behavior can be seen explicitly for the Born approximation of the renormalized frequency which is given by showing that the sign of the frequency shift is given by the sign of ∆.
Recently, the impact of the switching on and off of the coupling to the baths has been studied in detail investigating the system-bath interaction energy [52,53].We note that in our approach this effect is taken into account through the time dependence of the renormalized frequency of the system oscillator.

III. SIMULATION RESULTS
Let us now discuss some numerical simulations of important quantities characterizing the performance of the Otto cycle.All quantities and parameters presented in the following are either dimensionless or have the dimension of energy (or frequency, since we have set ℏ = 1).Therefore, all dimensionful quantities will be represented in terms of an arbitrary energy or frequency unit.
In Fig. 3 we show the behaviour of a) the total work output W out , b) the net heat input from the hot bath Q H and c) the efficiency η when varying ω c (the position of the peak of the spectral densities of both the hot and the cold bath).The structure of the spectral density does not play a role in the Born-Markov approximation, therefore this parameter will only affect the exact dynamics.The pink dashed lines show the Born-Markov case for comparison, and we see that they remain constant for all ω c .The dots show the results from the exact dynamics for different values of γ 0 , for the weak coupling case (blue) and the strong coupling case (yellow).As expected, we observe that in the weak coupling regime the results are generally closer to the Born-Markov approximation than in the strong coupling regime.
We also look at the behaviour of work output, heat input and efficiency with varying ω c for different values of the spectral width λ in Fig. 4. As can be seen from the second-order Born approximation in Equation (56), the value of λ also affects the coupling strength.This fact is reflected in the graph: For a narrow spectral density (blue lines) the results come closer to those from the Born-Markov approximation than for a wider one (yellow lines).
We note that in both cases shown in Figs. 3 and 4 when changing ω c the work output deviates more strongly from the Born-Markov approximation than the heat input.In our numerical simulations we find this to be a general behaviour when varying parameters of the spectral density.Therefore, the behaviour of work with changing characteristics of the spectral density has a greater influence on the behaviour of the efficiency than that of heat.
Figures 3 and 4 clearly show a characteristic behaviour: Work output and efficiency are greater than in the Born-Markov case in the region ω 1 < ω c < ω 2 , and smaller outside this region.The first intuition for this fact comes from Fig. 5, where we represent, for different values of ω c , the Otto cycle in the (ω, n)-plane, where ω is the frequency of the oscillator and n the mean occupation number.The colored solid curves represent the exact cycle, while the dashed black lines show the standard Born-Markov approximation.The area enclosed by the curves corresponds to the total work output.We can already see in this picture that the work output is higher for the case ω 1 < ω c < ω 2 than in the Born-Markov approximation (Fig. 5b).This case corresponds to ∆ 2 > 0 and ∆ 1 < 0 and, consequently, to ω r 1 (t → ∞) < ω 1 and ω r 2 (t → ∞) > ω 2 , which results in the cycle taking place within a larger range of frequencies.For the cases ω 1 > ω c (a) and ω 2 < ω c (c) we see that the work output is smaller than in the standard approach, and that here the frequency renormalization reduces the work output, instead of enhancing it.
This behaviour of the work output translates into efficiency, which is also larger than the Born-Markov efficiency within the region ω 1 < ω c < ω 2 .A simple approximation can be derived if we assume that within second order ω r (t) relaxes much more rapidly than n(t).We see from (56) that ω r (t) relaxes on the time scale λ −1 , while n(t) relaxes on the time scale γ −1 M , where γ M is the Markovian relaxation rate given by Eq. ( 29), which yields in the present case γ M = γ 0 λ 2 /(λ 2 + ∆ 2 ).Thus, the assumption λ −1 ≪ γ −1 M leads to the condition Under this condition we can use ( 46) and ( 48) to get the approximate expressions such that the efficiency of the cycle may be approximated by the renormalized effciency η r : Note that η r is obtained from the Markovian efficiency (41) by replacing the bare frequencies with the renormalized frequencies.The renormalized efficiency (59) is represented by the solid lines in Figs.3c and 4c.For the blue curves condition (57) is satisfied outside the resonances at ω c = ω 1,2 , and we see that (59) indeed provides a good approximation for the cycle efficiency.Note that in the vicinity of the resonances (∆ ≈ 0) condition (57) is violated which explains why the approximation (59) is less accurate close to resonance.To further understand the curves shown in Figs. 3 and  4 we give an interpretation based on the height and the slope of the spectral density at the initial frequencies of the harmonic oscillator when coupled to each bath.In Figs. 6 and 7 we represented the height (top) and the slope (bottom) of the spectral density at ω = ω 1 (solid lines) and ω = ω 2 (dashed lines), i.e., J(ω 1,2 ) and d dω J(ω 1,2 ), with varying ω c on the horizontal axis.
First, we will analyze the different behaviour of work output, heat input and efficiency for different values of the coupling strength γ 0 .In Fig. 6 the lines show γ 0 = 1.0  (blue) and γ 0 = 10.0 (yellow), and all other parameters are the same as in Fig. 3.When we look at the curves in Fig. 3, we find, first of all, that the closer to resonance with the corresponding baths (ω c → ω 1,2 ), the more the results deviate from the standard Born-Markov approximation (except at exact resonance, see below), while far away from resonance the approximation predicts the results accurately.We can relate this to the fact that the height and the slope of the spectral density at ω 1,2 close to resonance take significant values, while far away from resonance they have decayed to negligible quantities, and the bath looks flat in the vicinity of the initial frequency of the central harmonic oscillator.Additionally, we see that, for weak coupling (blue lines), W out presents a dip in the region ω 1 < ω c < ω 2 , while for strong coupling (yellow lines) work output and efficiency have no inflection point in this region.If we look at the spectral densities for these two values of γ 0 , we see that in the weak coupling regime the spectral densities and their slopes have decayed to negligible values for ω c centrally located between ω 1 and ω 2 , while for strong coupling the spectral densities still take high values in this region, with significant slopes, and therefore strong effects in both the heating and the cooling phases combine.This is the reason why work output (and consequently efficiency) is maximally enhanced in this region for strong coupling.We also see that at resonance with either of the baths, the quantities seem to come close again to the Born-Markov approximation.If we look at the slope of the spectral densities in Fig. 6, we see that it vanishes as it approaches resonance (ω c = ω 1,2 ), resembling a flat bath.
In Fig. 7 the height (a) and slopes (b) of the spectral density at ω 1,2 can be found for different values of λ, corresponding to the results shown in Fig. 4. The faster decay of the spectral density and its slope for small λ (blue) can account for the higher resemblance of the results to those in the Born-Markov approximation far from resonance.On the contrary, for large λ (yellow) the cycle is most altered for intermediate ω c between ω 1 and ω 2 , when the height and slope of both J(ω 1 ) and J(ω 2 ) are still significant, meaning that both the heating and the cooling phases feel the effects of a strongly coupled, structured bath.
We now turn to look at the different modes of operation of the cycle in different parameter ranges.In the standard Born-Markov approximation this distinction is clear: The cycle operates as a heat engine (W out > 0 and Q H > 0) for the configurations of frequencies and temperatures that satisfy (42), and as a refrigerator (W out < 0 and Q H < 0) otherwise (due to the fact that work output and heat input change sign at the same points).In the exact approach we can only determine the behaviour of the cycle numerically, and it will additionally depend on the specific parameters of the spectral density.We find that the exact approach can give rise to three different types of behaviour: As a heat engine or refrigerator as before, but also as a heater (W out < 0 and Q H > 0), where the cycle just uses work to amplify the flow of heat from the hot to the cold reservoir.
In Fig. 8 we depict the regions in which the cycle operates as a heat engine (green), as a refrigerator (blue) or as a heater (red).The dashed black lines represent the limit between the heat engine and the refrigerator regions in the standard Born-Markov approximation.On the horizontal axis we change k B ∆T , where we have chosen a reference temperature k B T 0 and vary the temperatures of the cold and the hot bath symmetrically such that k B T 1,2 = k B T 0 ∓ k B ∆T /2.On the vertical axis we change ∆ω, where we have chosen a ω c as the reference frequency and vary ω 1,2 = ω c ∓ ∆ω/2.It can be observed that in the exact approach the region in which the cycle operates as a heat engine is limited compared  to the Born-Markov case, with work output being positive for a smaller range of frequencies and temperatures.For small γ 0 and λ this difference is almost negligible, but it increases as we increase γ 0 and λ.Additionally, a new type of behaviour appears and expands as γ 0 and λ increase: The cycle operates as a heater, where a work input is used to dissipate heat from the hot reservoir into the cold reservoir.In spite of this, for large values of γ 0 and λ (Fig. 8d) the refrigerator region broadens.
Finally, we analyze the values of γ 0 and λ for which we can obtain the highest efficiencies.From Figs. 3  and 4 it may seem that, for a spectral density satisfy-ing ω 1 < ω c < ω 2 , the higher the parameters γ 0 or λ are, the higher is the efficiency.However, in Fig. 9 we show the efficiency as a function of both γ 0 and λ for a spectral density centered between ω 1 and ω 2 , and we see that the behaviour of η is not monotonic with increasing γ 0 or λ.We observe that, fixing one parameter, there exists an optimal value of the other parameter such that the efficiency is maximal.We show the results for two different temperatures of the hot bath k B T 2 , and it can be seen that the values of γ 0 and λ for which the efficiency is maximal change with temperature.Therefore, stronger coupling can enhance the efficiency of the cycle, but only up to a limit, at which point it starts being a hindrance.Larger spectral width λ can also start improving the efficiency, as the spectral density gets wider and starts having relevant effects on both the heating and cooling phases.However, if λ becomes too large, the slope of the spectral density decreases, resembling a flat, Markovian bath and resulting in a lower efficiency for the cycle.

IV. CONCLUSIONS
We have investigated the thermodynamics of the quantum Otto cycle in the presence of non-Markovian effects and strong system-bath coupling.Our analysis was based on a novel approach to quantum thermodynamics, which allowed us to exactly resolve the consequences of coupling a system to an arbitrary reservoir with memory effects.By comparing the standard Markovian Otto cycle with the exact non-Markovian treatment, we have identified feasible conditions under which the cycle exhibits enhanced work output and efficiency.These enhancements occur when the spectral density of the bath peaks at a frequency which is located within the driving range of the system, leading to an effective broadening of the range of frequencies in which the cycle operates.Conversely, outside these regions, the cycle performance either is lower or remains close to the standard Born-Markov approximation.When taking baths with a spectral density that enhances the efficiency, we encounter a small trade-off in the parameter region of frequencies and temperatures for which the cycle is producing work output, but this effect becomes noticeable only at very large coupling parameters.We have also shown that the behaviour of work output and efficiency strongly depends on the strength of the system-bath coupling and the width of the spectral density, showcasing how much the nuances of the bath spectral density can influence the performance of the cycle.
It is important to acknowledge that our analysis focused solely on the quantum Otto cycle.Further research is required to investigate the implications of memory effects in other quantum thermal machine models and explore their influence on diverse thermodynamic processes.Still, the present study highlights the role of memory effects and structured environments as a promising avenue for the enhancement of current quantum thermal machines -one notable implication of our findings being the potential for exploiting strong system-bath coupling to achieve more efficient quantum heat engines.
While in this work we have not concerned ourselves with the duration of the cycle -allowing fully adiabatic driving and infinitely long thermalization timesthis could in principle represent another positive aspect of considering reservoirs outside of the weak-coupling regimes.A promising path would be to investigate how the coupling strength can influence thermalization rates, possibly leading to shorter relaxation times and, in turn, higher efficiencies at maximum power.Moreover, while we have demonstrated enhancements in performance for specific scenarios, it still warrants exploration whether surpassing the Carnot efficiency is feasible with non-Markovian baths.

Figure 1 .
Figure 1.Standard quantum Otto cycle of a harmonic oscillator.Phases I and III are adiabatic compression and expansion processes, respectively, while phases II and IV consist of isochoric heating and cooling by coupling to heat baths.

Figure 2 .
Figure 2. Non-Markovian quantum Otto cycle of a harmonic oscillator.In this case, phases II and IV entail a timedependent renormalization of the frequency due to the interaction with the bath.

Figure 3 .
Figure 3. a) Total work output, b) net heat input in phase II and c) efficiency as a function of the central frequency in the spectral density, for different values of γ0.The dots are calculated through the exact approach, the pink dashed lines indicate the results from the Born-Markov approximation and the vertical black dotted lines indicate where ωc = ω1,2.The solid line in the efficiency graph shows the renormalized efficiency (59).The parameters used are ω1 = 14.0, ω2 = 16.0,k B T1 = 15.0,k B T2 = 20.0 and λ = 0.2.

Figure 4 .
Figure 4. a) Total work output, b) net heat input in phase II and c) efficiency as a function of the central frequency in the spectral density, for different values of λ.The points are calculated through the exact approach, the pink dashed lines indicate the results from the Born-Markov approximation and the vertical black dotted lines indicate where ωc = ω1,2.The solid line in the efficiency graph shows the renormalized efficiency (59).The parameters used are ω1 = 14.0, ω2 = 16.0,k B T1 = 15.0,k B T2 = 20.0 and γ0 = 1.0.

Figure 5 .
Figure 5.The Otto cycle represented in the (ω, n)-plane for different values of ωc.The dashed black lines show the Markovian Otto cycle, while the colored solid lines stand for the exact results.The area enclosed by the curves gives the total work output.The parameters used are ω1 = 14.0, ω2 = 16.0,k B T1 = 15.0,k B T2 = 20.0,γ0 = 5.0 and λ = 0.2.We see that the area is bigger than in the Born-Markov approximation for the case ω1 < ωc < ω2 (b), and smaller otherwise (a and c).

Figure 6 .
Figure 6.a) Value of the spectral density and b) slope of the spectral density at the initial frequencies of the central harmonic oscillator: J(ω1) (solid line) and J(ω2) (dashed line) as a function of ωc, for different values of γ0.The frequencies, temperatures and spectral width λ are the same as in Fig. 3.

Figure 7 .
Figure 7. a) Value of the spectral density and b) slope of the spectral density at the initial frequencies of the central harmonic oscillator: J(ω1) (solid line) and J(ω2) (dashed line) as a function of ωc, for different values of λ.The frequencies, temperatures and coupling strength γ0 are the same as in Fig. 4.

Figure 8 .
Figure 8. Representation of the parameter regions for different behaviours of the cycle: In the green region the cycle behaves as a heat engine (Wout > 0 and QH > 0), in the blue region as a refrigerator (Wout < 0 and QH < 0) and in the red region as a heater (Wout < 0 and QH > 0).The black dashed lines represent the limit between the heat engine and the refrigerator regions in the standard Born-Markov approximation (the cycle cannot function as a heater in this case).We choose ωc = 15.0 as a reference frequency and kBT0 = 17.5 as a reference temperature, and we vary ω1,2 and kBT1,2 symmetrically (such that ω1,2 = ωc ∓ ∆ω/2 and kBT1,2 = kBT0 ∓ kB∆T /2).

Figure 9 .
Figure 9. Efficiency of the Otto cycle as a function of γ0 on the horizontal axis and λ on the vertical axis, for different values of the temperature of the hot reservoir: a) k B T2 = 20.0 and b) k B T2 = 22.0.The other parameters are: ω1 = 14.0, ω2 = 16.0,ωc = 15.0 and k B T1 = 15.0.