Out-of-equilibrium operation of a quantum heat engine

Realistic quantum heat engines are operated at finite cycle times and thus represent open quantum systems in non-equilibrium. Here, their full time evolution including control of both the work medium and thermal contact to reservoirs is simulated non-perturbatively starting from a microscopic modeling. This gives access to explore the full dynamical range of specific models and identify domains, where they qualify for quantum thermal machines. As an example a finite-time generalization of the four stroke Otto engine is considered to reveal how microscopic dynamics, correlations between medium and reservoirs, and thermodynamic properties are related. The approach may provide accurate predictions and design concepts for nanoscale devices.

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; the motion of the relevant degrees of freedom is then captured by Langevin equations or related approaches [1][2][3].
At atomic scales and low temperatures, quantum mechanics takes over, and concepts of classical thermodynamics may need to be modified [4][5][6][7][8]. 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 [9][10][11][12]. While the first nanoscale engines implemented with trapped ions [13,14] or solid state circuits [15] still operated in the classical regime, more recent experiments entered the quantum domain [16][17][18][19]. In the extreme limit, the work medium may even consist of only a single quantum object [20].
One is thus faced with the theoretical challenge of finding quantum analogues of Langevin dynamics. In open quantum dynamics, particular challenges are nonnegligible correlations between work medium and thermal reservoirs, time dependence of system-reservoir couplings and non-perturbative external driving according to the engine's protocol [21][22][23][24][25][26][27]. In particular, the coupling/de-coupling processes to/from the reservoirs require work input which is non-negligible in microscopic contexts. Standard perturbative approaches to open quantum systems have provided some insight in this context [20,[28][29][30][31][32][33][34][35][36][37][38] but are limited in their ability to explore the full dynamics of these phenomena throughout the relevant parameter ranges of finite-time operation at low temperatures of quantum machines [39].
In this Letter we use a non-perturbative simulation platform and apply it to a finite-time generalization of the Otto cycle as paradigmatic example, see which can be rigorously derived [40] from either the exact Feynman-Vernon path integral formulation [41] or quantum Langevin equations [42]. Thermal reservoirs are thus not abstracted out, invoking approximations or tentative thermodynamic arguments, but are explicitly included in the modeling. In particular, accommodating timecontrolled thermal contact between medium and reservoirs is straightforward. Simulation results of periodic steady states can be identified with the observables of a full microscopic model of a heat engine. Work media with either a single harmonic or anharmonic degree of freedom are discussed to make contact with current experiments. Since the approach can easily be extended to few continuous or discrete degrees of freedom, it may provide detailed predictions for actual implementations using, e.g., trapped ions, NV-centers in diamond, or superconducting circuits.
Modeling-A quantum thermodynamic device with cyclic operation involving external work and two thermal reservoirs is described by the generic Hamiltonian where H m , H c/h denote the Hamiltonians of the work medium and the cold/hot reservoirs, respectively, with interactions H I,c/h . Not only the working medium is subject to external control, but also the couplings-this is required when using quantum dynamics to describe operation of the full compound according to specific engine protocols. As a minimal model we consider a particle in a one-dimensional potential, H m (t) = p 2 /(2m) + V (q, t), also motivated by recent ion-trap experiments [13,14]. Reservoirs are characterized not only by their temperatures T c < T h , but also by a coupling-weighted spectral density [41,43]. Assuming the free fluctuations of each reservoir to be Gaussian, these can be modeled as a large collection of independent effective bosonic modes with coupling terms of the form The coefficients µ c/h are conventionally chosen such that only the dynamical impact of the medium-reservoir coupling matters. In a quasi-continuum limit the reservoirs become infinite in size; thermal initial conditions are therefore sufficient to ascertain their roles as heat baths. This setting will be explored with the goal to determine its dynamics 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 in cyclic operation, either as quantum heat engine (QH) or refrigerator (QR), is not determined or known a priori. However, non-thermal states are to be expected, except in the combined limit of adiabatically slow driving/weak coupling.
In order to tackle this formidable task, we use a dynamical approach that originates from the exact Feynman-Vernon path integral approach [41,44], which provides a formal expression for the reduced density operator ρ m (t) = Tr R {ρ tot (t)} of the working medium (distinguished system) [40,45]. The quantum correlation func- are memory kernels of a nonlocal action functional, representing the influence of the reservoir dynamics on the distinguished system as a retarded self-interaction. We re-formulate this approach by constructing classical Gaussian processes which exactly match the correlation functions L c/h (t−t ), thereby mapping the dynamics of the distinguished system to a noisedriven propagation in its own Liouville space, governed by a Stochastic Liouville-von Neumann equation (SLN) [40]. Cases of strong or fast driving and low tempera-tures, which forbid the usual Markovian approximation, have been successfully treated using this method [46][47][48][49][50].
In the present context, we also consider timedependent control of the system-reservoir couplings (see [51]). For a reservoir with Drude-like spectral characteristics with a large UV cutoff ω cut [45], this leads to a SLN of the forṁ with reservoir-related terms Averaging over samples of the operator-valued process ρ ξ (t) yields the physical reduced density ρ m (t) = E[ρ ξ (t)]. The independent noise sources ξ h (t), ξ c (t) are related to the reservoir correlation function through is represented in the deterministic terms in (2). It will be important in the following that mixed system-reservoir correlation functions involving X c/h are faithfully represented by stochastic correlations involving ξ c/h and system observables; this is in contrast to deterministic approaches to reduced-system dynamics [30,35].
Equations (1) and (2) represent the first main result of this paper: It provides a rather general and nonperturbative platform to efficiently capture the quantum dynamics of set-ups with working media consisting of single or few continuous or discrete (spin) degrees of freedom; it also offers the flexibility to run different protocols such as for QHs or QRs, with explicit, unambiguous identification of per-cycle energy transfers to work or heat reservoirs. Next, we will demonstrate this for a finite-time generalization of the Otto cycle.
Engine cycle-For this purpose, steering of both the time-dependent potential V (q, t) and time-dependent couplings λ c/h (t) in an alternate mode is implemented, see Fig. 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 a protocol, where ω(t) varies around a center frequency ω 0 between ω h and ω c 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) (cf. Fig. 1). The isochore strokes are divided into an initial phase raising the coupling parameter λ c/h from zero to one with duration τ I , a relaxation phase of duration τ R , and a final phase with λ c/h → 0, also of duration τ I . The cycle period is thus T = 4τ I + 2τ d + 2τ R , as indicated in Fig. 1. In conventional treatments, τ I is often neglected, assuming the modulation of the interaction to be instantaneous and without effect on the energy balance of the cycle [52]. However, in a fully dynamical setting, such effects play a crucial role as discussed below. The total simulation time covers a sufficiently large number of cycles to approach a periodic steady state (PSS) with ρ m (t) = ρ m (t + T ).
Cyclic dynamics-A typical example is shown in Fig. 2 for a purely harmonic system, for which analytical results have been derived in limiting cases recently [13,22,53].
Here, we use it as a starting point to refer to the situation in ion trap experiments [14] and also to later identify 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 ω c/h ; the time to reach a PSS typically exceeds a single period. An alternative representation of the covariances is their parameterization through a squeezing amplitude r and a squeezing angle ϕ [54], shown over one cycle in (b, c). The phase space ellipsoid makes a full turn (2ϕ changes by 2π) with its width oscillating around equilibrium values for hot/cold baths. Further dynamical properties are depicted in (d, e). There, we take the oscillator at its mean frequency ω 0 as a reference and employ the corresponding Fock state basis to monitor populations p n (t) = n|ρ(t)|n and coherences ρ nm (t) = n|ρ(t)|m , n = m. 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 in the medium [35,55].
Work and heat-The key thermodynamic quantities of a QH are work and heat per cycle. Note that even though we operate the model at finite times with a medium far from equilibrium, these thermodynamic 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 systemreservoir correlations [25,56].
During the isochores heat is exchanged between medium and reservoirs with the heat released from a reservoir per cycle. Within the SLN the integrand can be transformed into an expression of the form "force×momentum", consisting of separate terms associated with reservoir noise and medium back action [51].
Work associated with expansion and compression is measured as injected power, i.e., while work associated with the coupling/de-coupling processes between medium and respective baths W I,c/h (coupling work) during the isochores is obtained in the same way from H I,c/h . While ambiguities in the definition of the internal energies of system and reservoirs can arise from finite terms H I,c/h [25,35,56,57], these terms do not enter heat and work of a complete cycle. In fact, careful analysis indicates [51] that the coupling work W I is completely dissipated [58]. Figure 3 displays W d , W I , and Q c/h versus the maximal medium-reservoir coupling γ at different temperatures. While W d saturates with increasing γ, Q h follows a non-monotonous behavior: Increasing thermal contact first facilitates heat transfer to the medium but suppresses it for larger γ. This phenomenon can be traced back to the growing coupling work W I . Lowering the temperature of the hot bath (keeping T c fixed) reduces the amount of both work and heat. The friction dependence of W d associated with unitary time evolution illustrates global features of the PSS. In Fig. 3d the ratio is the efficiency of a QH when work and heat have the appropriate signs. The coupling work generally obeys the inequality W I ≥ 0, thus diminishing the efficiency. In regimes where η is nominally negative, the system is not a QH, but merely a dissipator. The theory of the adiabatic Otto cycle and its extension using an adiabaticity parameter predicts some regimes of pure dissipation [59], however, without recognizing the crucial role of the coupling strength γ. The combined dependence on γ and thermalization adiabaticity parameter ω 0 τ I yields a phase diagram pointing out QH phase (η > 0) and a dissipator phase (η ≤ 0), Fig. 4. A QH is only realized if the coupling times τ I exceed a certain threshold which grows with increasing medium-reservoir coupling, due to increasing coupling work. Appreciable values of η typically require |W d | |W I,c/h |. In obtaining these data, we noted that the duration of transient dynamics differs between data points and verified a PSS in each case. The values of η we observe are always below the Curzon-Ahlborn efficiency [60] η CA = 1 − T c /T h and the Carnot efficiency η C = 1 − T c /T h [51].
Slowing the cycle tends to increase efficiency, as shown in Fig. 5. Output power (Fig. 5a) shows maximum power peaks which are significantly lowered and shifted when coupling work is no longer ignored. Indeed, in the QH regime the PSS substantially deviates from a mere sequence of equilibrium states as also illustrated by the von Neumann entropy S vN in Fig. 5. Even with a relatively long contact time compared to 1/γ, the entropy values alone indicate incomplete thermalization with the colder reservoir, even when the non-thermal nature of squeezing is disregarded.
We now turn to the impact of anharmonicities. While the overall picture resembles the harmonic case, details are quite different: Dynamical features such as covariances display smoother traces due to non-equidistant energy level spacings; this broken symmetry also reduces the efficiency as seen in Fig. 6. One can relate this reduction also to W d which due to the parametric drive is determined by a reduced q 2 for stiffer potentials. In terms of efficiency, anharmonicities seem then to play a similar role as enhanced thermal couplings, both having the tendency to localize the oscillator degree of freedom in position. The versatility of the setting is also eluci- dated when tuning parameters from the operation as a QH into that of a QR (ω h /ω c > β c /β h ). Thermodynamic quantities can be obtained from the microscopic dynamics as for the QH. In essence, we find that the QR regime defined by Q c /(W d + W I ) > 0 is only approached if medium-reservoir couplings are not too strong and temperatures of the cold baths are not too low, a signature of the third law in this dynamical setting (details see [51]).
In conclusion, by simulating non-perturbatively the dynamics of quantum thermal machines with single (or few) degrees of freedom as work medium we have obtained a complete characterization of their properties, from microscopic details to thermodynamic quantities. 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 specific example of the four stroke Otto engine demonstrates that also at finite cycle times the dynamical system operates as QH/QR in certain domains and with reasonable efficiency and output power. In nanoscale realizations it may thus be beneficial to exploit the relatively fast time scales on which control is possible, e.g. with complex pulse shapes obtained via optimal control techniques. The approach provides the required tools to follow theoretically experimental activities to design quantum thermal devices.
We thank R. Kosloff We indicate how to adapt the SLN simulation method to the case of reservoirs with time-dependent coupling. This results in a complete, exact dynamical simulation of the quantum Otto engine in terms of the full, interacting dynamics of the work medium and its two thermal and harmonic reservoirs. We also provide expressions for work and heat, relating expectation values of the full system-plus-reservoir dynamics to SLN-based expressions. Some more detailed numerical results are presented and discussed.

STOCHASTIC MAPPING OF OPEN SYSTEM DYNAMICS WITH TWO RESERVOIRS AND DRIVING
We consider a distinguished quantum system H m (t) = p 2 2m + V (q, t) which is interacting with two thermal, harmonic reservoirs H c/h under influence of time-dependent bilinear coupling Starting from factorizing initial conditions for the global density matrix ρ tot and reservoirs with ohmic characteristics, i.e. spectral densities of the form J(ω) = mγω/(1 + ω 2 /ω 2 cut ) 2 up to a high frequency cutoff ω cut (significantly larger than any other frequency of the problem, including β c/h ) and potential renormalization µ = 2 π ∞ 0 dω J(ω) ω , the Feynman-Vernon path integral formulation for the reduced density operator of the medium can be converted into the Stochastic Liouville-von Neumann equation with dissipation (SLN) of the form [1] This procedure involves mapping the reduced system evolution to a stochastic propagation in probability space that is governed by two Gaussian noise sources ξ h (t), ξ c (t) that are constructed to match the real part of the respective reservoir correlation function ξ α (t)ξ α (t ) = L α (t − t ) − 2mγα βα δ(t − t ), α = c, h. The subtracted white-noise term is treated separately, leading to the last term in (2). In the case of ohmic reservoirs, the limit of large ω cut also makes the dynamics response of the reservoir near instantaneous, allowing its representation by a Dirac delta term of the form L α (t − t ) = mγα 2 d dt δ(t − t ), α = c, h. An integration by parts transforms L α (τ ) into a delta function, with a boundary term removing the q 2 term in the coupling, and introducing terms dependent on momentum p and the time derivativeλ. A single stochastic sample obtained by propagating the density with an individual trajectory possesses no direct physical meaning, except in the classical limit → 0. The physical density matrix is obtained by averaging over a sufficiently large number of samples ρ m (t) = E[ρ ξ (t)].
Taking the spectral density of the reservoirs to be a smooth function of frequency implies the limit of an infinite number of environmental modes, with a correspondingly infinitesimal coupling of each individual mode. This procedure also gives the reservoirs infinite heat capacity, allowing them to be treated as thermodynamic reservoirs even if only their initial state is specified as thermal.

HEAT ENGINE CYCLE
We will now consider a four stroke quantum heat engine with an harmonic or an anharmonic oscillator as working medium. Our description is based on a non-perturbative propagation of the reduced density matrix within the dissipative SLN probabilistic framework. The model comprises three control parameters. The hot and cold baths are controlled with the time-dependent coupling parameters λ c/h (t), the oscillator potential is frequency modulated arXiv:1903.11368v1 [quant-ph] 27 Mar 2019

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 Eq. (1) are due to work terms, with separate terms indicating the different modes of performing work associated with the parameter ω(t) and λ c/h (t).
Heat is identified as energy transferred into the reservoirs over a cycle, At the end of a period of cyclic operation, the microstate of the medium reverts its initial state, ρ m (t + T ) = ρ m (t). Moreover, the collective response function of the reservoirs decays in time sufficiently fast that the collective reservoir coordinate entering H I,c/h shows periodic behavior at long times. We thus have Defining per-cycle work terms through it is easily verified that these quantities obey the first law of thermodynamics It is instructive to note that Eq. (6) can be used to give alternative expression for the heat terms per period, The first term of (9) is an energy flow from the system due to the coupling; the second is the work performed through changes of λ c/h . Viewing this equation as a continuity equation, we conclude that the work described by the second term is completely dissipated.

Work and heat in the probabilistic SLN context
The work terms W d and W I as well as the expression (9) for the heat transfer Q c/h have equivalent expressions in the SLN dynamics (2), even for those terms involving H I,c/h . In order to obtain this stochastic equivalent, the equal-time correlations qX α and pX α , (α = c, h), involving system coordinate/momentum and the reservoir operator X α = k c k,α (b † k,α + b k,α ) are needed. In simpler SLN approaches [2], the noise variable ξ can serve as a direct substitute for X [3]. Here we need to treat the terms which have been contracted to delta terms in the more complicated SLN equation (2) separately. A careful consideration of short-time dynamics on timescales of order 1/ω cut before taking the limit of large ω cut leads to the results Hamiltonian dynamics SLN dynamics where α = c, h. Table I 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 Eq. (9), thus avoiding expressions involving momenta or velocities of the reservoir.

Details about the numerical implementation
The numerical solution of the dissipative SLN leads to one single trajectory of the reduced system density. This is realized by moving to position representation using symmetric and antisymmetric coordinates, i.e. ρ ξ (r, y) = r − y 2 |ρ ξ |r + y 2 . This representation allows an efficient split-operator technique. In addition to the commonly used FFT method (alternating between diagonal potential and kinetic terms) we employ a third step related to the operator [q, {p, ·}], which can be understood to be the generator of a re-scaling operation. For typical parameters, e.g., parameters ω 0 β h = 0.25, ω 0 β c = 3, ω cut /ω 0 = 30, γ/ω 0 = 0.25 and characteristic cycle time scales ω 0 τ I = 10, ω 0 τ d = ω 0 τ R = 5, ω 0 T = 60 up to the first three cycles of a periodic steady-state takes approximately 72 CPU core hours on a Intel Xeon CPU (Sandy Bridge architecture). A typical number of samples n samp ≈ 500; the resulting statistical errors are less than the line width or symbol size used in our figures.

Additional results
Here we show additional results that are not included in the main text but provide further indications on the features and versatility of our simulation platform. Fig. 1a gives an alternative rendition of the PSS's operating regime as a QH, depending on the maximal medium-reservoir coupling and the coupling time as parameters. Numbers on top indicate the number of driving periods before the PSS cycle is reached. Fig. 1b compares the efficiency of the finite-time harmonic engine to the Carnot and the Curzon-Ahlborn efficiencies. The efficiencies lie always below the analytical values of η CA and η C . For decreasing damping strength, the efficiencies tend to decrease as a result of diminished heat transfer and driving work. Stronger damping meanwhile fosters energy losses due to irreversible coupling work which also leads to decreased efficiency. The cycle's efficiency can be increased by extension of period T (see main text). Fig. 2 shows details of the dynamics of the medium in a QR setting. These include (a) moments which fully characterize a non-thermal Gaussian state as well as work, heat, and efficiency for a refrigerator setting (b-d). Fig. 3 shows variances in position and momentum in thermal equilibrium for anharmonic oscillators with varying anharmonicity parameter, illustrating the fact that the oscillator "stiffness" becomes effectively temperature dependent. The resulting deformations of equilibrium position and momentum variances have significant impact on heat, work and efficiency properties for increasing κ (see discussion in the main text). Fig. 4 shows the driving work  Anharmonic oscillator in thermal equilibrium with one single thermal reservoir: equilibrium position q 2 th and momentum p 2 th dispersion for γ/ω0 = 0.1, ω0 β th = 3 (red) and γ/ω0 = 0.5, ω0 β th = 5 (blue) versus the anharmonicity κ ∈ [0, 0.9]; with increasing κ, the oscillator potential gets stiffer and leads to decreasing second moments in position while those of the momentum are broadened; the horizontal lines correspond to analytical equilibrium values in the harmonic case.