Non-Markovian dynamics of a quantum heat engine: out-of-equilibrium operation and thermal coupling control

Real quantum heat engines lack the separation of time and length scales that is characteristic for classical engines. They must be understood as open quantum systems in non-equilibrium with time-controlled coupling to thermal reservoirs as integral part. Here, we present a systematic approach to describe a broad class of engines and protocols beyond conventional weak coupling treatments starting from a microscopic modeling. For the four stroke Otto engine the full dynamical range down to low temperatures is explored and the crucial role of the work associated with the coupling/de-coupling to/from reservoirs as an integral part in the energy balance is revealed. Quantum correlations turn out to be instrumental to enhance the efficiency which opens new ways for optimal control techniques.


Introduction
Macroscopic thermodynamics was developed for very practical reasons, namely, to understand and describe the fundamental limits of converting heat into useful work. In ideal heat engines, components are always in perfect thermal contact or perfectly insulated, resulting in reversible operation. The work medium of real macroscopic engines is typically between these limits, but internally equilibrated, providing finite power at a reduced efficiency. Any reduction of engine size to microscopic dimensions calls even this assumption into doubt.
At atomic scales and low temperatures, quantum mechanics takes over, and concepts of classical thermodynamics may need to be modified [1][2][3][4][5]. This is not only of pure theoretical interest but has immediate consequences in the context of recent progress in fabricating and controlling thermal quantum devices [6][7][8][9]. While the first heat engines implemented with trapped ions [10,11] or solid state circuits [12] still operated in the classical regime, more recent experiments entered the quantum domain [13][14][15][16]. In the extreme limit, the work medium may even consist of only a single quantum object [17].
Theoretically, one is thus faced with the fundamental challenge that a separation of time and length scales on which conventional descriptions of thermal engines is based, may no longer apply. This has crucial consequences: first, the engine's operation must be understood as a specific mode of the cyclic dynamics of an open quantum system with the coupling/de-coupling processes to/from thermal reservoirs being integral parts of the time evolution; second, thermal coupling strength and thermal times  k T B at low temperatures T may match characteristic scales of the work medium. The latter requires a non-perturbative treatment beyond standard weak-coupling approaches [17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] to include medium-reservoir quantum correlations and non-Markovian effects [32][33][34][35][36][37][38]. The former implies the introduction of two distinct sources of work. While typically only exchanged heat Q and work for compression/expansion of the medium (driving work) W d is addressed, a complete dynamical operation of any quantum thermal machine requires finite-time coupling/de-coupling processes to/from the respective thermal reservoirs that inevitably lead to an additional time-dependence of the system-reservoir interaction and, hence, to an additional source of work, see figure 1. This coupling work W I may turn into an essential ingredient in the energy balance as will be discussed in section 5 of this manuscript.
, for a quantum device this separation of scales may fail andW W d I | | | |. This is particularly true for scenarios to approach the quantum speed limit in cyclic operation [39][40][41]. Can a quantum engine under these conditions be operated at all?
To address this question, in this paper we provide a non-perturbative treatment and apply it to a finite-time generalization of the Otto cycle (figure 1). It is based on an exact mapping of the Feynman-Vernon path integral formulation [42,43] onto a SLN equation [44] which has been successfully applied before [45][46][47][48][49][50]. Here, we extend it to accommodate time-controlled thermal contact between medium and reservoirs and thus to arrive at a systematic treatment of quantum heat (QH) engines at low temperatures, stronger coupling and driving. Work media with either a single harmonic or anharmonic degree of freedom are discussed to make contact with current experiments. We demonstrate the decisive role of the coupling work W I which inevitably must enter the energy balance. In general, its impact turns out to be detrimental for the efficiency, however, its dependence on quantum correlations opens ways to improve the situation, e.g. using optimal control techniques [46].

Modeling
A quantum thermodynamic device with cyclic operation involving external work and two thermal reservoirs is described by the generic Hamiltonian . Not only the working medium is subject to external control, but also the couplings-this is required in a full dynamical description of the compound according to specific engine protocols. We consider a particle in a one-dimensional potential, ( ) ( ) ( ), also motivated by recent ion-trap experiments [10,11]. Reservoirs are characterized not only by their temperatures T c <T h , but also by spectral densities w J c h ( ) [43,51], related to the dynamical response functions c t c h ( ) of the reservoirs through Assuming the free fluctuations of each reservoir to be Gaussian, these can be modeled in a standard way [43] as a large collection of independent effective bosonic modes with bilinear coupling terms of the form Figure 1. Top: energy-frequency diagram of the work medium in a quantum Otto heat engine with frequency ω(t) varying around ω 0 . A cycle includes two isochore (A→B, C→D) and two isentropic strokes (B→C, D→A). Bottom: thermal contact to hot (cold) reservoirs is controlled by λ h (t) [λ c (t)] and expansion (compression) is due to ω(t). The cycle is specified by three characteristic time with coupling constants c k c h , and dimensionless coupling functions l t c h ( ) which are varied between zero and one to close and open thermal contact to the respective reservoirs. In terms of creation and annihilation operators of the bosonic excitations, the harmonic reservoir Hamiltonians for the cold and hot bath, respectively, read The terms ò ò are conventionally chosen such that no static force arises from the medium-reservoir coupling for constant q and λ [51]. Equivalently, this guarantees translationally invariant medium-reservoir interactions and the dynamical stability of the model for V(q)=0. In a quasicontinuum limit the reservoirs become infinite in size; thermal initial conditions are therefore sufficient to ascertain their roles as heat baths.
The dynamics of this setting will be explored over sufficiently long times such that a regime of periodic operation is reached, without limitations on the ranges of temperature, driving frequency, and system-reservoir coupling strength. The nature of the quantum states encountered, either as QH engine or refrigerator (QR), is not known a priori.
In order to tackle this formidable task, we start from the Feynman-Vernon path integral formulation [42,43,51]. It provides a formally exact expression for the reduced density operator r r = t t Tr are memory kernels of a non-local action functional, representing the effective impact of the reservoir dynamics on the distinguished system as a retarded self-interaction. They can be recovered from the spectral density through This formulation can be exactly mapped onto a stochastic Liouville-von Neumann equation (SLN) [44], an approach which remains consistent in the regimes of strong coupling, fast driving, and low temperatures [45][46][47][48][49][50], where master equations become speculative or inaccurate.
Without detailed microscopic knowledge of the reservoirs, the assumption of an ohmic reservoir w w µ J [ ( ) ] is commonly made; in a classical limit of equation (4), this results in white thermal noise and memoryless friction. For a quantum reservoir, however, ohmic dissipation results in colored noise with a long- for finite temperatures and an algebraic decay L(t)∝1/t 2 at zero temperature. Thus, non-Markovian dynamics inevitably appears for the ohmic case at lower temperatures ω cut ÿβ?1 with a typical reservoir high frequency cut-off ω cut significantly larger than any other frequency scale of the problem.
Here we substantially extend an SLN-type method for ohmic dissipation [52], i.e. spectral densities of the form w gw w w = + J m 1 2 cut 2 2 ( ) ( ) , to the highly non-trivial time-dependent control of the system-reservoir couplings. This includes a concise mapping of key reservoir observables and system-reservoir correlations to their respective stochastic representation. The resulting dynamics is given by which contains terms known from the conventional SLN equation and further parts that result from a careful handling of the explicit time-dependence of the system-reservoir coupling control and the finite-memory quantum noise x c h g l r r Averaging over samples of the operator-valued process ρ ξ (t) yields the physical reduced density The independent noise sources ξ α (t) are related to the reservoir correlation functions through . The time local equations (5) and (6) thus provide a non-perturbative, non-Markovian simulation platform for quantum engines with working media consisting of single or few continuous or discrete (spin) degrees of freedom; different protocols can be applied with unambiguous identification of per-cycle energy transfers to work or heat reservoirs. Next, we will apply it to a four stroke Otto cycle.

Engine cycle
For this purpose, steering of both the time-dependent potential V q t , ( )and time-dependent couplings l t c h ( ) in an alternate mode is implemented, see figure 1. For simplicity, ohmic reservoirs with equal damping rate γ are assumed. A single oscillator degree of freedom represents the working medium as a particle moving in with a parametric-type of driving ω(t) and anharmonicity parameter κ0. We consider ω(t) varying around a center frequency ω 0 between w w  D > w D , 0 0 2 ( ) , within the time τ d during the isentropic strokes of expansion (B→C) and compression (D→A); it is kept constant along the hot and cold isochores (A→B and C→D) (see figure 1). The isochore strokes are divided into an initial phase raising the coupling parameter l c h from zero to one with duration τ I , a relaxation phase of duration τ R , and a final phase with l  0 c h , also of duration τ I .
The cycle period is thus figure 1. The total simulation time covers a sufficiently large number of cycles to approach a periodic steady state (PSS) with r r ). Conventionally, one neglects what happens during τ I ; one assumes that modulating the thermal interaction has no effect on the energy balance (see also [32]). In the quantum regime, such effects may, on the contrary, play a crucial role as will be revealed in the sequel.

Periodic steady state
In figures 2(a)-(c) results are shown for a purely harmonic system, for which analytical results have been derived in limiting cases [10,33,53]. We use it as a starting point to refer to the situation in ion trap experiments [11] and to identify in (d) the role of anharmonicities. After an interval of transient dynamics (a), the elements of the covariance matrix settle into a time-periodic pattern with damped oscillations near frequencies w   time-reversibility which implies that a description in terms of stationary distributions with effective temperatures is not possible. Indeed, the PSS substantially deviates from a mere sequence of equilibrium states. Further insight is gained by taking the oscillator at its mean frequency w 0 as a reference and employ the corresponding Fock state basis to monitor populations r = á ñ p t n t n n ( ) | ( )| and coherences Population from the ground and the first excited state is transferred to (from) higher lying ones during contact with the hot (cold) bath. In parallel, off-diagonal elements ρ nm (t) are maintained. These are dominated by ρ 02 (t) contributions according to the parametric-type of driving during the isentropic strokes. While this Fock state picture has to be taken with some care for dissipative systems, it clearly indicates the presence of coherences associated with qp-correlations in the medium [26,54]. The impact of anharmonicities for stiffer potentials in (7) is depicted in (d). In comparison to the harmonic case, dynamical features display smoother traces with enhanced (reduced) variations in á ñ p 2 (á ñ q 2 ). Non-equidistant energy level spacings may in turn influence the efficiency (see below).

Work and heat
The key thermodynamic quantities of a QH are work and heat per cycle. Note that even though we operate the model with a medium far from equilibrium, these quantities have a sound and unique definition in the context of fully Hamiltonian dynamics involving reservoirs of infinite size. An assignment of separate contributions of each stroke to heat and work is not needed in this context. Moreover, any such assignment in a system with finite coupling would raise difficult conceptual questions due to system-reservoir correlations [36,55].
In the context of full system-reservoir dynamics, heat per steady-state cycle is uniquely defined as the energy change of the reservoir Within the SLN the integrand can be transformed into computable expressions involving system-reservoir correlation functions. We can thus determine the transferred heat without the use of generating functionals [56,57]. Similarly, work is obtained as injected power, i.e.
where separate driving and coupling work contributions W d and W I affect power output and efficiency of the cyclic operation.

Energy balance and first law
We first consider the energy balance of our engine in the context of the full system-reservoir model. Any changes in the energy of the global system as defined in equation (1) are due to work terms it is easily verified that these quantities obey the first law of thermodynamics It is instructive to note that equation (12) can be used to give an alternative expression for the heat terms per period The first term of (15) is an energy flow from the system due to the coupling; the second is the work performed through changes of l c h . Viewing this equation as a continuity equation, we conclude that the work described by the second term is completely dissipated (which has been previously shown for a Markovian classical heat engine model [58]).

Work and heat in the probabilistic SLN context
The work terms W d and W I as well as the expression (15) are needed. In simpler SLN approaches [44], the noise variable ξ can serve as a direct substitute for X [48]. Here we need to treat the terms which have been contracted to delta terms in the more complicated SLN equation (5) separately. A careful consideration of short-time dynamics on timescales of order w 1 cut before taking the limit of large ω cut leads to the results r x l m l g l g = á ñ+ á ñ-á + ñ -á ñ where α=c, h. Table 1 summarizes expressions for work and heat using either the full Hamiltonian dynamics or the SLN framework. All quantities appearing in the SLN column can be extracted from simulation data. The expression for Q c/h is based on equation (15), thus avoiding expressions involving momenta or velocities of the reservoir. Having now access to the derived probability space representations of the crucial thermodynamic quantities at any point during the cyclic operation, their impact on the heat engine's efficiency and power output can be analyzed and compared across a multitude of parametric regimes.

Net work, efficiency, and power output
The cyclic operation of the QH engine according to a specific finite-time coupling/de-coupling and driving protocol can be evaluated with respect to its characteristic thermodynamic quantities. Figure 3(a) displays the strong coupling dependence of the net work W d +W I . It turns from negative (net work output) to positive, thus highlighting the coupling work W I as an essential contribution in the work balance. The SLN approach allows to reveal two distinct components, i.e. = + W W W I I I ,var ,corr . The first one, determined by á ñ q 2 , also exists at high temperatures and close to the adiabatic limit (equilibration), while the second one, depending on á + ñ qp pq -correlations, becomes especially relevant at finite cycle times and in the deep quantum regime. Since á + ñ qp pq only contains a 2 and (a † ) 2 when expressed in annihilation and creation operators of corresponding Fock states (while á ñ q 2 also contains occupations Table 1. Work and heat expressed in terms of the unitary evolution of system and reservoir and in terms of SLN propagation of a periodic steady state. , thus counteracting W I,var (b). In turn, control of qp-quantum correlations opens ways to tune the impact of W I on the energy balance as will be discussed in section 8. Heat Q h , see (c), follows a non-monotonous behavior with γ, also a genuine quantum effect that cannot be captured by standard weak coupling approaches. Its decrease beyond amaximum can be traced back to enhanced momentum fluctuations due to damping.

Hamiltonian dynamics SLN dynamics
We are now in a position to discuss the ratio In regimes where η is nominally negative, the system is not a QH, but merely a dissipator in the sense that driving work adds to the energy flow from hot to cold reservoir. The theory of the adiabatic Otto cycle and its extension using an adiabaticity parameter predicts some regimes of pure dissipation (see equation (6) of [10]), however, without recognizing the coupling work W I as an essential ingredient. As seen above, its detrimental impact can be soothed by quantum correlations, see (d).
The combined dependence on γ and thermalization adiabaticity parameter w t I 0 yields a phase diagram pointing out QH phase (η > 0) and a dissipator phase (η 0) over a broad range of thermal couplings up to γ/ω 0 ∼1, figure 4. A QH is only realized if τ I exceeds a certain threshold which grows with increasing mediumreservoir coupling. To make this more quantitative, progress is achieved for small compression ratios to estimate W d and W I as we will show in the next section.
As expected, in figure 4 values obtained for η are always below the Curzon-Ahlborn and the Carnot efficiencies, but yet, even beyond weak coupling, they do exceed η∼0.2. The coupling work appears as an essential ingredient also to predict for the power output the optimal cycle time and correct peak height, figure 5(a). If it is ignored, misleading data are obtained. Beyond the harmonic case, i.e. for stiffer anharmonic potentials, dynamical features discussed in figure 2(d) reduce the efficiency, figure 5(b). They have a similar impact as enhanced thermal couplings, both having the tendency to suppress (increase) fluctuations in position (momentum).

Analytic estimates for driving and coupling work
In order to achieve a better insight into the numerical results at least in certain regimes, we start from the formulation of the driving and the coupling work in the SLN context as specified in table 1. The driving work  provides contributions along the isentropic strokes of length τ d with w The coupling work provides contributions within time intervals of length τ I at the beginning (l > a 0  ) and the end (l < a 0  ) of the isochoric strokes. To estimate the above integrals we now assume the following: during τ d and l = t const. ( )  during τ I , (ii) the contribution of the variances dominates in W I and (iii) within τ d the variances change linearly for Δω/ω 0 =1. We then obtain with t X corresponding to point X in the cycle. The coupling work can be approximated by where we considered á ñ » q const. 2 during switching the coupling off or on for sufficiently short τ I . Due to the general relations á ñ > á ñ á ñ > á ñ q q q q contact to cold bath , contact to hot bath 21 for γ>0, one first concludes that W d <0 while apparently W I >0. Further, one finds for the total net work  . Now, to qualify for a heat engine, the condition W d +W I <0 needs to be fulfilled which implies w t w g w t w g . Due to (21) one always has 0R1. The above relation can be easily solved for τ I to read (ii) In the quantum regime (lower temperatures) in quasi-equilibrium one has R>R cl .
(iii) For low temperatures and sufficiently large gt R to allow for states close to thermal equilibrium, the variances in position depend only weakly on γ for g b p <  2 . Accordingly, the γ-dependence of the threshold of τ I is predominantly given by γ/Δω.
(iv) For ÿβ c ω c ?1 while ÿβ h ω h =1 and γτ R >1 to allow for approximate equilibration with etc, one finds R≈ÿβ h ω 0 . For the parameters in figure 4 of the main text we then have ω 0 τ I >12γ/ω 0 which, for the given value of τ R , describes the minimal τ I sufficiently accurate for all γ/ω 0 >0.1.

Role and control of the coupling work
As we have seen above, the overall impact of the coupling work is detrimental to the efficiency of QH engines. However, it can be reduced by qp-correlations if they are properly controlled. To see this in detail, we realize that according to ,corr : the first part which survives the adiabatic limit is determined by the variance á ñ q ; 2 the second part, breaking time reversal symmetry, vanishes for adiabatic driving and is determined by the x á ñ q average and the contribution of the correlations á + ñ qp pq . The contribution W I,corr is dominated by the part depending on the qp-correlations. It is this part which we want to consider in the following in more detail:  , where the sign depends on the sign (phase) of the qp-correlations at the beginning of the medium-reservoir coupling segments. Its dominant contribution is provided by á + ñ qp pq C , i.e. after expansion and before coupling to the cold reservoir. For the parameters chosen in figure 3, one always has á + ñ > qp pq 0 to lead to a reduced W I . An instructive example is shown in figure 7. By extending the unitary time evolution after expansion/ compression such that the system evolves for about half a period π/ω α at constant frequency, the phase of the qp-correlations at C turns from positive to negative: coupling work is further enhanced and the efficiency further suppressed. By applying more advanced techniques, e.g. from optimal control, one could be able to shape the impact of the coupling work properly.
We note that experimentally a direct time dependent medium-reservoir coupling as considered here, is often difficult to implement. In recent developments for superconducting devices [59], for example, one instead uses cavities interfacing medium and reservoir with strong thermal contact to the respective broadband reservoir and tunable resonance frequency. The effective spectral distribution seen by the medium is thus Lorentzian which allows to continuously vary between weak and strong spectral overlap and thus between no and maximal thermal coupling. Qualitatively, this scenario is captured by our model.

Summary and outlook
In conclusion, by simulating non-perturbatively and within a systematic formulation the dynamics of quantum thermal machines with single degrees of freedom as work medium, we have obtained a complete characterization of their properties. The medium-reservoir boundary appears as an internal feature of the model so that full control over the medium as well as its thermal contact to reservoirs is possible. The example of the four stroke Otto engine demonstrates the decisive role of the coupling work that must be considered as an integral part of the total energy balance. Its overall impact is detrimental to the efficiency of QH engines, however, can be reduced by qp-correlations if they are properly controlled. This sensitivity of QH engines to changes in the driving protocol can be exploited by optimal control techniques in future devices. The presented approach provides the required tools to follow theoretically these activities.