Dynamics of a strongly coupled quantum heat engine -- computing bath observables from the hierarchy of pure states

We present a fully quantum dynamical treatment of a quantum heat engine and its baths based on the Hierarchy of Pure States (HOPS), an exact and general method for open quantum system dynamics. We show how the change of the bath energy and the interaction energy can be determined within HOPS, for arbitrary coupling strength and smooth time dependence of the modulation protocol. The dynamics of all energetic contributions during the operation can be carefully examined both, in its initial transient phase and also later, in its periodic steady state. A quantum Otto engine with a qubit as inherently nonlinear work medium is studied in a regime where the energy associated with the interaction Hamiltonian plays an important role for the global energy balance and, thus, must not be neglected when calculating its power and efficiency. We confirm that the work required to drive the coupling with the baths depends sensitively on the speed of the modulation protocol. Remarkably, departing from the conventional scheme of well-separated phases by allowing for temporal overlap, we discover that one can even gain energy from the modulation of the bath interactions. We visualize these various work contributions using the analogue of state change diagrams of thermodynamic cycles. We offer a concise, full presentation of HOPS with its extension to bath observables, as it serves as a universal tool for the numerically exact description of general quantum dynamical (thermodynamic) scenarios far from the weak-coupling limit.


I. INTRODUCTION
Quantum thermodynamics tries to extend the powerful classical theory to the realm where the dynamics of the involved systems cannot be described by classical physics 1 .First models and results have been investigated way back in the last century 2,3 .Through the enormous experimental advances in preparing, controlling, and measuring small individual quantum systems during the last decades, however, the subject has gained significant attention more recently.Quantum thermodynamic setups are no longer limited to thought experiments, but can actually be implemented on a broad variety of platforms such as trapped ions 4 , superconducting circuits 5,6 , semiconductor quantum dots 7,8 , or nuclear spins [9][10][11][12][13][14][15][16] .
Among the various quantum thermodynamic scenarios that have been proposed and investigated, quantum heat engines have always been of particular interest.It is the conversion of heat into work which made classical thermodynamics a practically very relevant theory and an extension to the quantum realm seems natural.However, it turns out that this task is not as straightforward as one might expect.The concepts of work and heat have proved to be elusive in a quantum setting, and there is an ongoing debate on a meaningful (measurable) definition of these quantities in the quantum domain [17][18][19][20][21][22][23][24] .Still, a vast amount of proposals for quantum heat engines and refrigerators have been investigated in recent years [25][26][27][28][29][30][31][32][33][34][35] .Importantly, all of those obey the Carnot bound and are thus constrained by the same limits we know from classical thermodynamics.Examples of quantum heat engines which seem to break the Carnot limit rely on some sort of non-thermal energy reservoir 1,36,37 .
Even though quantum heat engines will most likely not be able to outperform classical ones, there are microscopic setups where quantum effects play an important role for the description of the thermodynamic behavior.For example, it has been shown that coherences in the work medium or non-classical correlations can have a great influence on the performance of quantum engines, increasing power or efficiency in some scenarios [38][39][40][41][42][43][44] , while resulting in a quantum disadvantage in others 45,46 .
As a complex many-body quantum system, a first principle theoretical treatment of a quantum heat engine based on a full system-bath model is notoriously difficult.Regarding the reduced dynamics of the work medium, a variety of exact approaches has matured over the last years, steadily extending the accessible parameter regime [47][48][49][50][51][52][53] .Nonetheless, to address the relevant properties such as efficiency and power of a quantum heat engine, a handle on the bath dynamics and its observables is needed in addition to the sole system dynamics.In particular, the change of the bath energy and the energy attributed to the interaction Hamiltonian are of relevance.
Advanced methods for open system quantum dynamics allow for access to both quantities.Most notably, bath related observables can be obtained numerically from the auxiliary density operators of the Hierarchical Equations of Motion (HEOM) formalism of open system dynamics 33,34,[54][55][56] .Alternatively, a stochastic Liouville-von Neumann equation approach to open system dynamics 57 can be adapted for time dependent modulations and the determination of interaction and change of bath energy.It was used by Wiedmann et al. in 31 to study the dynamics, power and efficiency of an Otto engine with an oscillator work medium, serving as a guidance for our treatment of the Otto cycle here, based on HOPS.It is also known that knowledge of appropriate system multi-time correlation functions allows for the determination of relevant bath properties 58 .In practice, this is numerically hard too, as it requires the evaluation of the full influence functional or process tensor 59 , nowadays obtained from tensor network contraction methods 53,60 .
In this work, we show how to employ the stochastic Hierarchy of Pure States (HOPS) approach to open quantum system dynamics 48,61,62 for a full dynamical treatment of a quantum heat engine and its baths.Notably, HOPS reflects a stochastic representation of the full state of system and environment and thus, it comes as no surprise that collective bath observables are accessible, too.We show that they can be calculated in terms of the auxiliary states of the HOPS, complementing similar findings in HEOM 33,55,56,63 .Note that our HOPS formalism does not impose any constraints on the spectral density of the baths or the nature of the work medium, unlike the approach in Ref. 31 , which is based on the infinite high energy cutoff limit.Here, we choose a bath frequency cutoff comparable to system energy scales.
Special features of the stochastic approach, e.g. the unitary stochastic representation of finite temperature, are shown to be of great benefit, as the hierarchy itself is determined from zero temperature properties only.Furthermore, all the advantages of HOPS are preserved.Foremost among them is the improved memory efficiency that results from working with pure state trajectories of dimension D of the system Hilbert space instead of density matrices of dimension D 2 .Importantly, in this article we work out how to apply the non-linear variant of HOPS to the study of system-bath interaction and bath observables, which allows to efficiently treat the regime of strong system-bath coupling.
As a central part of this work, we use HOPS to investigate the power and efficiency of a quantum heat engine from first principles, similar to the studies performed in Refs. 31,34.A general model of a heat engine consists of a work medium S with Hamiltonian H S which is coupled to N baths with bare Hamiltonians H (n) B .The interaction is mediated by Hamiltonians H (n) I .Note that in generaland in particular for engines that work in cycles -both the system Hamiltonian and the system-bath interaction are explicitly time-dependent.Thus, the general global Hamiltonian that we consider in this study has the form with bosonic heat baths.On this global level, all energy contributions or their changes are well-defined and acces- FIG. 1.The Otto cycle exemplified with a single two-level system as work medium.In general, the unitary "compression" (I)/"expansion" (III) of the work medium increases/decreases the level spacing while changing the orientation of the qubit, too.
sible by the HOPS.By taking a global unitary approach, difficulties often encountered in reduced descriptions, such as the ambiguity in the splitting of the interaction energy in a system and a bath contribution, do not occur.We will rely on the well-defined expectation values of the respective parts of the Hamiltonian: ⟨H S ⟩, ⟨H I ⟩, and ∂ t ⟨H B ⟩ (see also Refs. 31,34).No ad hoc separation of work-like and heat-like energy is called for.
We will illustrate the strength of such a global approach by studying a common model of a quantum heat engine using HOPS: the Otto cycle.A thorough and related analysis with an oscillator instead of our qubit work medium was performed by Wiedmann et al. 31 based on a stochastic Liouville-von Neumann equation approach.Simplified investigations of the Otto engine in terms of master equations can be found in Refs. 3,38,64.An Otto cycle consists of four separate phases, which are depicted in Fig. 1 and listed in the following (consider the system initialized in equilibrium with the cold bath).
the cold bath is decoupled from the system, that is, H cold I (t 0 ) = 0. Now, the system Hamiltonian changes from H S (t 0 ) = H 0 → H S (t 1 ) = H 1 and the work medium evolves unitarily.Depending on these two system Hamiltonians, the machine will work either as a heat engine or as a refrigerator.To operate in the heat-engine regime, the level-spacing of H 1 should be greater than that of H 0 .This ensures that energy can be extracted from the system through inverting this step later in the cycle.The term "compression" is chosen as the level spacing increases.This is analogous to the spectrum of a particle in a box whose volume is decreased.II t 1 → t 2 isochoric heating: The system Hamiltonian is kept constant.The interaction H hot I with the hot bath is turned on and the system equilibrates with this reservoir.
III t 2 → t 3 isentropic "expansion": The hot bath is decoupled, H hot I (t 2 ) = 0, and the system Hamiltonian changes back to its initial value, i.e.H S (t 2 ) = H 1 → H S (t 3 ) = H 0 .The energy of the entire system decreases.Ideally, the energy decrease is larger than the energy increase in phase I.
IV t 3 → t 0 isochoric cooling: The system Hamiltonian is kept constant.The interaction with the cold bath H cold I is turned on, the system equilibrates with the cold bath and, finally, the interaction is turned off.The cycle has closed and can be repeated.
The clear-cut separation of the four phases is an idealization.In a real Otto engine, be it classical or quantummechanical, all operations are continuous in time.Moreover, the actions that are performed in the four phases may overlap in time.In fact, this additional freedom will be exploited later to optimize the heat engine.
Crucially, cyclic modulations will in general lead to nonequilibrium dynamics in the system.Still, in (classical) macroscopic engines -due to the many degrees of freedom and their fast internal thermalization -the dynamics of the work medium is well described by quasistatic state changes.Heat engines based on current quantum technologies, by contrast, use work media that are very small -most notably qubits or single bosonic modes 4,5,8,14,16,65 .For such microscopic systems, a quasistatic description is no longer appropriate.Foremost, the assumption of an instantaneous internal equilibration does not hold, in particular for finite work media as considered here.Second, even if the cycle is slow enough such that the system can always be considered in a steady state with respect to the instantaneous external parameters, this steady state is generally not given by a Gibbs state with respect to the instantaneous local Hamiltonian.Such a description is valid in an ultra-weak coupling regime only 66,67 .However, in a generic scenario with a finite (or even large) coupling strength between the system and the baths, for the steady state of the system the interaction Hamiltonian has to be taken into account 31,68,69 .
Thus, in order to fully describe a quantum heat engine it is necessary to consider transient effects, dynamics well beyond the quasistatic regime, as well as finite coupling strengths with the baths.As the HOPS method can satisfy all these demands, we are able to determine the dynamics of the engine for an arbitrary driving of the work medium through a time-dependent system Hamiltonian, as well as arbitrary (continuous) protocols for coupling/decoupling to/from the heat baths through the timedependent interaction Hamiltonians.Moreover, there is no restriction on the dimension of the work medium nor on a specific form of the spectral density.
In this article, we go beyond previous considerations of the HOPS approach to open system dynamics by exploiting that the change of the bath energy and the interaction energy can be obtained from within that formalism naturally.We will provide illustrative examples with incomplete thermalization and different coupling strengths that confirm the impact of smoothly time-dependent couplings to the cold and hot baths, as shown in 31 .We stress again that the work contributions associated with modulation of the interaction Hamiltonians, i.e. switching the coupling to the baths on and off, significantly contribute to the overall work and energy balance.It turns out that these contributions are sensitive to slight changes of the modulation protocol, as are power and efficiency (see also 31 ).Thus, a correct characterization of a quantum heat engine requires a careful assessment regarding all details of the cycle.Furthermore, we extend these considerations by exploring overlapping cycle phases with an eye on optimizing the performance of the engine.Remarkably, we find that the coupling/decoupling of the baths can contribute positively to the power-output of the machine.
The remainder of this article is structured as follows.In Sec.II we review natural definitions of power and efficiency of a heat engine in terms of the global energy balance of system, interaction, and bath.Sec.III introduces a qubit (two-level) heat engine with an Otto-like cycle (OLC) modulation protocol, on which all results of this work are based.We discuss the dynamics of the work medium and the flow of the various energetic contributions.We then turn to the parameter dependence of power and efficiency of the engine.In Sec.IV we analyze modulations beyond the OLC by allowing overlapping cycle phases and continue in Sec.V with an investigation on how cycle length and coupling strength affect the performance of the machine.In Sec.VI, we present details of the employed HOPS method, focusing on how relevant bath observables can be computed.We close with conclusions and an outlook in Sec.VII.

II. ENERGY BALANCE, POWER, AND EFFICIENCY
In classical thermodynamics, the performance of a heat engine is quantified in terms of the heat that is transferred from and to the system and the work extracted during a cycle.For quantum heat engines, work and heat are not easily defined in such a setting.Fortunately, in our global approach we have access to the total energy balance of system, interaction and bath.Thus, we do not need to rely on ambiguous definitions of work and heat for the subsystems.Instead, power and efficiency can be defined solely through the well-defined expectation values of the involved Hamiltonians 31,34,56 and are introduced below to make this work self-contained.
It has been demonstrated by Wiedmann et al. 31 that for small quantum systems, the energy associated with the system-bath interaction Hamiltonian needs to be studied carefully.While in an ultra-weak coupling scenario one simply neglects this contribution to the total energy balance, for the finite-time cyclic heat engine studied here, it turns out to play an important role.

A. Power
The global dynamics of system and bath evolves unitarily, i.e. without any entropy production.Thus, a global energy change can always be considered work.Accordingly, the total power is the rate of change of the global energy, where the last equality follows from the time independence of the bath Hamiltonians (see Eq. ( 1)).The total power in Eq. ( 2) consists of individual contributions associated with respective Hamiltonians, namely the system power P S (t) and, likewise, power contributions P (n) I (t) from the coupling/decoupling of the baths.Note that a usable power output corresponds to a positive value for the power, as the global system loses energy.
For the performance of a heat engine, one is usually interested in the cycle-averaged long-time behavior.The average power P and the work W performed over one such limit cycle are defined through where Θ denotes the length of the cycle.

B. Efficiency
In conventional thermodynamics, the efficiency is defined as a benefit to cost ratio.In this case, the benefit is the work output W .The cost is the energy consumed from the hot bath.Recall that due to the infinite number of bath degrees of freedom, the bath energy is infinite.Nonetheless, the change of the bath energy ∆ H hot B is finite and is accessible by the global HOPS method.Note that this energy change is not necessarily only heat-like.In Ref. 31 it is shown that in the periodic steady state this quantity can be decomposed into terms that stem from the change in the system energy due to the interaction and changes due to the modulation of the interaction.Note however, that care has to be taken with such an interpretation as these contributions are not entirely independent through the coupled dynamics of the entire, global system.
As we consider a general scenario with time-dependent and strong couplings, the modulation of the interaction Hamiltonian will also perform work on the baths 70 .We could not care less: clearly, the efficiency of an engine is determined by the net energy loss of the hot bath over one cycle that is converted into work, accounting for all energetic processes in the bath, be it heat or work 63 .Consequently, the efficiency of the engine is where ∆ denotes the change of a quantity over one limit cycle of length Θ.We note in passing that the limit cycle condition ensures ∆ ⟨H S ⟩ = 0 and ∆ ⟨H I ⟩ = 0, such that holds.It can be shown that Eq. ( 4) is bounded by the Carnot efficiency η C = 1−T cold /T hot under the assumption that a limit cycle has been reached.This is generally true whenever all system quantities are time-periodic as for our model (see 56 , the argument given there can be easily extended to time dependent couplings).For an OLC the efficiency is also bounded by the efficiency of an ideal Otto cycle using a qubit as work medium η Otto = 1 − ω cold /ω hot , where ω i is the energy gap of the qubit during coupling to the hot and cold bath, respectively -a result arising from rather general considerations 71 .

III. QUBIT HEAT ENGINE MODEL AND ITS DYNAMICS
In this section, we let a single qubit take the role of the work medium for the heat engine, which is coupled to an infinite bosonic hot and cold bath along the cycle in a spin-boson like fashion.The total Hamiltonian of the Otto engine thus reads where i = hot, cold.With bosonic annihilation and creation operators a, a † , the bath Hamiltonians are = λi ω λi a † λi a λi , and the coupling agents read B i = λi g λi a λi .The explicit time dependence of both, system and system-bath interaction Hamiltonian are encoded in the real-valued smooth and periodic (period Θ) functions ⃗ s(t), s 0 (t) and h i (t).Their choice determines the specific protocol for system modulation and bath coupling/decoupling.The time-dependent identity contribution to the system Hamiltonian is introduced to control the ground state energy of the work medium (fixed to zero, in most cases).A similar scheme with, however, an instantaneous instead of a smooth switching on/off of the interaction Hamiltonian has been studied in Ref. 63 .
In order to model a somewhat generic scenario for the qubit work medium, we choose , and s y = 0, where 0 ≤ f (t) ≤ 1 is a smooth function modulating the qubit energy gap periodically in line with the phases of the cycle (see also Fig. 2).In terms of the classical Otto cycle, this system modulation corresponds to expansion and compression of the work medium.In the expanded phase at f (t) = 0 the (smaller) level spacing is ε 0 = Ω 1 + (s x /Ω) 2 whereas in the compressed phase at f (t) = 1, we have (the larger) Whenever s x ̸ = 0, the time dependence of f (t) results in a time-dependent orientation ⃗ s(t) of the system Hamiltonian.This implies that the instantaneous eigenstates of the system Hamiltonian become time-dependent -the thermodynamic relevance of which will be discussed later.For s x = 0, by contrast, the eigenstates of σ z are the time-independent eigenstates of H S (t) throughout the cycle.
With B i (t) denoting the above coupling agent in the interaction picture with respecfunctiont to its H (i) B , the microscopic structure of the two baths and their coupling to the work medium is captured in the zero-temperature bosonic bath correlation function 72 We choose an Ohmic spectral density J i (ω) for each bath with identical exponential cutoff (see Methods Sec.VI for details) but possibly different coupling strengths ηi , The cutoff frequency ω c is set to the rather small value ω c = Ω of the order of the qubit level spacing.Thus, the memory time of the environment is of the same order as the system timescale, and the dynamics is expected to be non-Markovian.Unless otherwise stated, we set the temperature of the cold bath to T cold = Ω/2 and that of the hot bath to T hot = 4Ω.The latter is such that a reasonable occupation of p ↑ ≈ 0.38 of the upper level can be expected from the corresponding Gibbs state of the qubit (see also lower panels of Fig. 2).The lower temperature (up state occupation of p ↑ ≈ 0.12) is small enough to yield a significant system energy difference during one cycle.Note that we refer to the Gibbs state here as a rough estimate for the choice of temperaturesthe true quasi-equilibrium state while coupled to a bath during the cycle will in general not be the local Gibbs state due to the strong interaction.
In order to achieve a comparable effect of both baths on the system dynamics, we choose the coupling strengths ηi such that the value of the effective thermal spectral densities agree at the approximate level spacing, i.e.
and inverse temperature β = (kT ) −1 .This choice is motivated from the ultra-weak coupling regime where, according to Fermi's Golden Rule, the environmental influence is mainly determined by the value of the thermal spectral density at the resonance frequencies 66 .The specific value serves as global measure for the coupling strength between the qubit and the baths.Unless otherwise stated, we choose δ = 0.7.
For the following calculations, the initial condition is fixed to the product state where |↓⟩ denotes the eigenstate of σ z with eigenvalue −1 and ρ β the Gibbs state of the corresponding bath Hamiltonian.As the system starts off in a pure state, we expect some transient dynamics before the heat engine reaches a periodic limit cycle (the third, for our choice of parameters).

A. Otto-like cycle modulation protocol for the qubit engine
We begin our investigation with a modulation protocol inspired by Ref. 31 which aims to be close to the ideal Otto cycle (see upper row in Fig. 2), i.e.: I Compression phase: Change the system Hamiltonian from expanded f = 0 to compressed f = 1 (see Eq. ( 5) and Eq. ( 6)), while the interaction with the baths is off (h i (t) = 0).
II Hot Phase: Turn on, hold and turn off the interaction with the hot bath using a smooth, non-zero 0 ≤ h hot (t) ≤ 1.

III Expansion phase:
Undo the change of the system Hamiltonian (suitable f (t)).
IV Cold phase: Turn on, hold and turn off the interaction with the cold bath (through 0 ≤ h cold (t) ≤ 1).

Restart the cycle.
We call such a protocol Otto-like cycle (OLC) to emphasize that the four phases are non-overlapping, as in the Otto cycle, while any modulation of the global Hamiltonian is modelled by a smooth transition.In particular, we utilize smooth sigmoid-like functions that are twice continuously differentiable.The time interval to perform a single transition will be denoted by τ tr and is initially set to τ tr = 0.06Θ.The cycle time Θ is chosen to be Θ = 60/Ω, which allows the system to relax to an effective steady state during the phase where it is in contact with either of the baths (given that the coupling strength is no less than roughly δ ∼ 0.5, see Eq. ( 10)).This also means that we consider a regime where the correlation time of the bath (∼ 1/Ω, see above) is short compared to the cycle length.These parameters are chosen to be in a regime of positive power output without too much fine-tuning.
Before focusing on power and efficiency, we first discuss the dynamics of the work medium for such an OLC, visualized in Fig. 2. In particular, the evolution of the Bloch vector is shown for two different choices of the modulation of the system Hamiltonian in the middle and lower panel of that figure.
First, the middle panel shows the more general case where the scale and the orientation (s x ̸ = 0) of the qubit Hamiltonian is modulated.During the initial compression phase I (grey background), the qubit state rotates on the surface of the Bloch sphere due to the σ x component of the system Hamiltonian.This rotation is visible as initial coherent oscillations of the x and y Bloch vector components.At the beginning of phase II, the interaction with the hot bath is smoothly turned on (light red background).As expected 73 , the qubit loses purity, as is indicated by the decrease of the absolute value of the Bloch vector (black solid line, middle panel of Fig. 2).Once the interaction is fully turned on (red background) the qubit relaxes towards an effective steady state.In case of the hot bath, this steady state agrees fairly well with the local Gibbs state ∼ e −β hot H S (f ) | f =1 of the qubit at the hot bath temperature (thick red horizontal lines).Notably, during the decoupling from the hot bath (light red background at the end of phase II) the qubit experiences a slight change towards a purer state.This is closely related to a small decrease of the system energy during the decoupling stage (further discussed in the next paragraph).The isolated system modulation in phase III (grey background) produces unitary dynamics (constant length of the Bloch vector).Since the orientation of the system Hamiltonian changes, its instantaneous eigenstates change as well, as does the qubit state.As soon as the interaction with the cold bath becomes active (blue background) the qubit state becomes purer which indicates its cooling.The qubit then reaches an effective steady state while in contact with the cold bath.Yet, that state deviates markedly from the corresponding local Gibbs state (thick horizontal blue lines).As our choice of parameters leads beyond the usual weak coupling regime, this is to be expected and shows the influence of the interaction Hamiltonian on the steady state.Again, turning off the interaction with the bath results in a slight increase of purity.Now the cycle is repeated, and we see that the initial transient dynamics -which depended on the chosen initial qubit state -are already strongly suppressed at the start of the second cycle.There is no noticeable difference between the second and third cycle.Consequently, we can regard the third cycle as a reasonable approximation of the limit cycle.
Now we consider the second, somewhat simpler scenario, where we set s x = 0 in Eq. ( 6).Here, the modulation keeps the system Hamiltonian diagonal in the σ z -basis for all times.In that case, it holds generally true that the spin-boson Hamiltonian acts on two distinct subspaces of the global Hilbert space and, thus, the Schrödinger equation does not induce any transitions between them (parity symmetry).As a consequence, our choice of initial condition results in dynamics confined to the z-axis of the Bloch sphere.This behavior can be seen in the computed simulations shown in the bottom panel of Fig. 2. We note, however, that the individual pure state trajectories of the employed stochastic HOPS method (details in Sec.VI) do not share that property.This is why very small deviations of r x and r y from zero show up in the simulations.These deviations can in principle be made arbitrarily small by increasing the number of stochastic trajectories used to calculate the reduced state.We thus conclude that in the lower panel for s x = 0, the initial transient dynamics is very different from the middle, and moreover, r x and r y remain zero for all times.Still, the overall limit cycle behavior, featuring the relaxation towards an effective steady state with higher energy when in contact with the hot bath and lower energy when in contact with the cold bath, is very similar in both cases.Therefore, in the following discussion on power and efficiency, we will focus on that simpler case without the static σ x contribution.
So far we have discussed the exact dynamics of the work medium, i.e. of the reduced qubit state.Considering the explicit time dependence, the strong coupling to the environment, and not least a bath correlation time which is of the order of the system timescale, this is challenging on its own.Additionally, the HOPS approach allows access to the dynamics of the bath degrees of freedom, to which we turn next.In particular, changes of bath energies as well as the energetic contribution due to the interaction Hamiltonian (details in App.VI C and App.VI D) can be computed.Thus, we are in a position to calculate power and efficiency of the quantum heat engine from first principles.
First, we show in Fig. 3 the dynamics of the individual energy contributions while operating the engine.Using HOPS, the detailed access to the dynamics of all energy contributions, including the interaction energy, allows us to present results beyond the earlier treatment in Ref. 31 .The system energy ⟨H S ⟩ (black) and both interaction energies ⟨H c./h.I ⟩ (dashed blue, dashed red) are shown in terms of their absolute value.For the infinite baths, it only makes sense to consider energy changes.We choose ∆⟨H c./h. B ⟩ (solid blue, solid red) to denote the change of the bath energy from the value at their respective initial conditions.energy, all quantities are zero.The pure compression during the first system modulation (gray shading) only lifts the upper qubit level, which leaves the system energy (solid black) at zero (referring also to the dynamics shown in the bottom panel of Fig. 2).

Initial system compression (I):
Hot phase (II): When turning on the interaction with the hot bath (II↑, light red shading), the energy of the system rises while the contributions of the interaction and the bath decrease.Summing up the three parts shows that the change in total energy ∆⟨H⟩ (purple) decreases, too.This means that while switching on the coupling to the hot bath, work is extracted (in contrast to the standard thermodynamic Otto-cycle).Once the interaction with the hot bath is fully established (red shading) the total energy remains constant as the global Hamiltonian is now time-independent.Powered by the energy flow from the hot bath, the system heats up further as the qubit relaxes towards its temporary (hot) steady state.Notice, however, that part of that missing hot bath energy flows into the interaction energy, which increases slightly after the initial sharp drop.Before the interaction with the hot bath is switched off, the dynamics of the energy contributions have relaxed and reach constant values.The process of switching off the hot bath (II↓, light red shading) requires energy since the interaction energy goes from its temporarily negative value to zero.That energy is only partly provided by the system (and, to a lesser degree, by the hot bath), meaning that the total energy rises: switching off the hot bath requires work.

System expansion (III):
The decrease of the system energy upon expansion of the hot qubit directly carries over to a decrease of total energy, i.e. further work is extracted.This is clearly the main mechanism for work extraction (as for the standard Otto-cycle that motivated the design of the model).
3. Dynamics of the various energy contributions: work medium (system) energy (black), change of the cold bath energy (blue), cold bath interaction energy (dashed blue), change of the hot bath energy (red), hot bath interaction energy (dashed red) and change of the global energy (purple -not constant due to the explicit time dependence of the global Hamiltonian).The colored dotted lines show the overall trend for the respective quantity after a full cycle (same parameters as in Fig. 2 but with sx = 0 only).
Cold phase (IV): Switching on the interaction lowers the total energy and some work is extracted (IV↑, light blue shading).This means that in absolute values, the energy decrease of the system and the interaction is larger than the increase of the cold bath energy.When fully coupled (blue shading) the cold bath takes energy from the system.The qubit is cooled as intended.At the same time -but to a lesser degree -the cold bath also takes up energy from the interaction.Again, the process of decoupling the bath requires work (IV↓, light blue shading).The increase of the interaction energy is mainly powered by that external work, but receives small contributions from the system and the cold bath, too.Overall, the full cold phase requires work.
2 nd system compression (I.2):In contrast to the initial compression, as the qubit has not completely relaxed to its ground state, the compression of the system now requires some work, too.After that, the general behavior of the following cycles agrees with the first one.
Following the overall change of the bath energies over consecutive cycles (dashed lines in Fig. 3) reveals that the energy flow into the cold bath is far less than the energy lost from the hot bath.Therefore, part of that energy lowers the total energy which has been extracted as work due to the various time dependent modulations of system and interactions.In other words, we have demonstrated that the OLC based on the spin-boson model with Ohmic environments works as a heat engine.Note that once the system has reached its limit cycle, the change of system energy over one cycle is zero.By construction, this holds true for the interaction energy, too.
In numbers, during the third cycle, the energy change is ∆⟨H hot B ⟩ = −0.489Ωfor the hot bath and ∆⟨H⟩ = −0.147Ωfor the total energy.This amounts to an average cycle power of P = 2.45 • 10 −3 Ω 2 and an efficiency of η = 30.0%.The first protocol in Fig. 2, with a σ xcomponent of s x = 0.3Ω/2 of the system Hamiltonian, achieves ≈ 93% of the above mentioned power output and efficiency.Increasing s x to 0.5Ω/2 (not shown in Fig. 2) decreases the work output and efficiency further to ≈ 80% of the model with s x = 0.The primary work extraction mechanism in the present model is the modulation of the Hamiltonian of the work-medium.If the modulated part of the system Hamiltonian does not commute with the constant part, coherences are generated.Thus, the Bloch vector of the work medium moves away from the axis along which energy is measured, leading to decreased work output for a similar change in hot bath energy.Therefore, the efficiency and power both change by the same relative amount.The decrease in power output due to coherences is known as "quantum friction" 1,74 .
Next, we address the sensitivity of these operationally relevant quantities with respect to the modulation protocol.

IV. INFLUENCE OF THE MODULATION PROTOCOL, OVERLAPPING CYCLE PHASES
It is clear that the timing and modulation of the bath interactions will influence power and efficiency.For instance, a refrigerator is realized by exchanging the hot and cold bath in the sequence of the cycle.A shift in time of the bath couplings may result in a temporal overlap with the expansion/compression phase: it turns out that the efficiency of the engine can be significantly enhanced by carefully tailoring such an overlap.In the following FIG. 4. The different modulation schemes discussed in the text.
In contrast to Ref. 31 , the bath modulation phases are shifted relative to the system modulation, and relative to each other.
From the bottom up, the modulation of the OLC is shown.This cycle features no shifts and all transitions happen within the time τtr = 0.06Θ.Above, the bath-coupling modulation is shifted by the amount τ .Next, the modulation shift is the same as below, but with a slower transition time τtr = 0.12Θ.
In the top row only the cold bath coupling modulation is shifted (slowly as below).The shading of the background corresponds to the modulation of the system shown in the lowest row.
we investigate the influence of the modulation protocol by considering three cases: (i) a simultaneous shift in time of the cold and hot bath interaction, (ii) that same simultaneous shift but with slower transitions, and (iii) the case where only the cold bath is shifted, with those slower transitions (see Fig. 4).This third case allows for an overlap of the hot and cold phases, with even further room for optimization.We also tried to modulate the various parts of the Hamiltonian simultaneously, yet didn't see further improvement of the performance.
A common temporal shift of both couplings amounts to the replacement h h./c.(t) → h h./c.(t − τ ) with the shift parameter τ .In Fig. 5, upper row, the power and efficiency of the heat engine as a function of the coupling shift τ relative to the modulation of the work medium is shown.It turns out that a positive shift -a delay of the hot/cold phases -enhances the power up to a noticeable factor of ≈ 1.7 at a shift of τ ≈ 0.2Θ (orange line).For larger shifts, the power drops and eventually becomes negative, which signifies the regime of a refrigerator.As shown in the bottom row of Fig. 5, left, the primary source for the boost in power is indeed the contribution of the interaction energy.This underlines the importance of considering that contribution to the power output in agreement with Ref. 31 , where under their choice of modulations, switching on the interaction always requires work.Indeed, for the original OLC protocol from above (lower left panel at τ = 0) we observe the same behavior (see Fig. 6 for details).Now, however, we see that this work can be reduced to zero, or work can even be extracted for a newly designed protocol with shifted modulation (overlapping phases).Once the shift is too large, the compressed (expanded) state is no longer heated (cooled) and the work extraction due to the system modulation quickly drops.This leads to a dramatic decrease in average power output and efficiency, as holds true for negative shifts, too.Although the power shows a plateau over roughly 1/6 of the cycle length, the efficiency starts to decline earlier, for shifts τ ≳ 0.2Θ, as displayed in Fig. 5, upper row, right.Note, however, that power and efficiency can be improved by the shift in the modulation timing at the same time.
In addition to the timing of the modulation, we show in Fig. 5 (green) that the speed of the modulation has a significant impact, too, as was found in Ref. 31 .As an example, we consider a protocol where all modulation steps are performed at half speed (shifted, slow), while keeping the cycle length constant.Fig. 4 as well as the background shades at the bottom of Fig. 6 depict the phases of that protocol.The parameter τ tr is used to specify the duration of a single transition ("normal speed": τ tr = 0.06Θ, half speed: τ tr = 0.12Θ referred to as slow modulation).The lower middle panel of Fig. 5 shows that the system contribution to the average power is nearly unaffected, while the interaction contribution is raised compared to the case with original speed.Notably, for slow modulation, the average power contribution from the interaction becomes positive for a wide range of positive shifts with a maximum at τ ∼ 0.25Θ.At τ = 0.18Θ, the total power has increased by more than a factor of 2 compared to the original OLC.At the same time, slowing down the modulation results in a larger efficiency for the heat engine of more than 40% compared to the value of ∼ 30% for the original speed (Fig. 5, upper row, right).For a detailed analysis of these findings, see Fig. 6, where the evolution of system and interaction work contributions for the different modulation protocols are displayed.
As an intuitive explanation for the observed improvement in the interaction power contribution, consider the following: In the expanded phase, the spectral density of the baths (recall that we have chosen the sames spectral density for both baths) has its maximum at the system level-spacing (ε 0 = Ω).This resonance leads to a larger magnitude of the negative interaction energy.In the sudden limit, where the interaction is turned-off instantaneously, this energy would be entirely compensated by the interaction work.In the adiabatic limit, the opposite would be true, and the negative interaction energy would be entirely supplied by changes in system and bath energy.Our regime lies between those two extremes, and thus it requires work to decouple the bath.This work input grows with the magnitude of the interaction energy.If this energy is reduced by any means other than external modulation, less external work is needed to turn off the  5. Upper row: Power and efficiency as a function of the shift of the timing of the coupling.A shift of zero (gray vertical lines) corresponds to the original OLC protocol.We consider three cases: Hot and cold phases are both shifted (orange lines), they are both shifted yet the switching on to maximum coupling is slower (green lines), and the cold bath only is shifted with slower transition to maximum coupling (red lines).Efficiency is shown for scenarios with positive work only.Bottom row: The individual contributions to the total average power (black lines) due to change in the system (magenta lines) and in the interaction Hamiltonian (brown lines) are shown.The three diagrams correspond to the different cases of the upper row (parameters: sx = 0, δ = 0.7, ωc = Ω, T cold = Ω/2, T hot = 4Ω, Θ = 60/Ω, τtr = 0.06Θ, slow modulation τtr = 0.12Θ).
interaction.This is accomplished by detuning the system from the bath (compression) before decoupling.The effect is also present for the hot bath.However, the effect should be reversed compared to the cold bath.Recall, in the OLC, the system is decoupled from the hot bath when it is compressed (ε 1 = 2Ω) and thus off-resonant with the (zero-temperature) spectral density.Shifting the decoupling of the hot bath into the expanded phase actually brings the system in resonance.This doesn't seem to damage the performance of the cycle, which might be due to the higher temperature of the hot bath which obscures the zero-temperature structure of the bath through thermal fluctuations (see Section VI C 2).According to the above, not shifting the coupling to the hot bath should improve the performance of the engine even further.Thus, we consider a set of protocols where the modulation of the interaction with the cold bath is shifted only (in combination with slow modulations).The results shown in Fig. 5 and Fig. 6 (red lines) indicate that the mean interaction power is even larger, compared to a simultaneous shifting, but the system power decreases slightly.In total, a very similar behavior of the average power is observed for values of the shift that are not too large.However, for the protocol where the cold bath timing is changed only, the efficiency for the shift with the largest average power is increased from 40.3% to 43.4%.
A more detailed analysis shows that the improvement in efficiency originates from the reduced energy flow out of the hot bath.If the hot bath is decoupled after the work medium is expanded, additional energy can flow into the system from the hot bath during and after the expansion.This is because we disrupt the equilibrium state by modulating the system.As the system energy gap is now smaller and more on resonance with the hot bath, additional energy can be transferred to the system before reaching a new equilibrium.Although the work output of the interaction modulation is improved, the system work output suffers slightly.This is due to the slight reduction in system energy when uncoupling the hot bath, as can be observed for the OLC in Fig. 3. Overall, we get similar work output, but improved efficiency.This is due to the additional energy taken from the hot bath instead of the interaction work.We conclude that considering the structure of the environment is indeed important when optimizing the cycle performance.
Using a different representation as in Fig. 6, it is instructive to discuss and compare the overall energy changes in terms of thermodynamic work diagrams (Fig. 7, see also Ref. 34 ).We stress that even though time no longer appears explicitly in these diagrams, the curves are simply FIG. 6. Work performed by the system (solid lines) and by the interaction modulation (dashed lines), shown as a function of time during a limit cycle and slightly beyond (same parameters as in Fig. 5) for the best-performing cycles (maximum power/efficiency) with τ = 0.18Θ.As reference, we show again the work for the original OLC modulation (upper shading of the phases) with blue lines.The blue arrow indicates that work needs to be done in order to drive the interaction with the baths.
For the shifted modulation (orange lines) the compression/expansion of the qubit takes place when the system is in full contact with the cold/hot bath (see middle shading of the phases).In this case, the difference of interaction work after a full cycle is almost zero.For the shifted cycle and a slower modulation (green lines, lower shading of the phases) one actually gains work from the driving of the interaction, as emphasized by the green arrow.When only shifting the cold phase (red lines, no shades for that protocol shown but see Fig. 4) the work gain due to the bath interactions is even larger, although the system work is smaller.Note that we compare the work curves of the cycles where both modulations are shifted just after decoupling from the cold bath (dash-dotted vertical lines).
derived from the time-dependent calculations as displayed in Figs. 3 and 6 and thus by no means we imply quasistatic state changes.
For the work diagrams, we regard the change of the Hamiltonian as a function of the externally controlled modulation functions f and h i . From we see that, e.g., the work due to the system modulation amounts to the area enclosed by the f -⟨∂ f H⟩-diagram.
In Fig. 7 we compare state changes using the original OLC protocol with the improved protocol where the cold bath is shifted by τ = 0.18Θ in combination with slower transitions, i.e. τ tr = 0.12Θ.The plots offer various insights.First, it is evident that the work contribution due to the modulation of the interaction Hamiltonian must not be neglected in the overall energy balance, as the trajectories in the corresponding state space extend over similar regions.Note that the individual subplots have been scaled such that work identified with the visual appearance of some area is indeed the same for all plots.Second, although the shape of the system work diagram changes when the phases overlap (lower panel in Fig. 7) compared to the OLC (upper panel), the area and, thus, the work remains almost unaffected.In contrast, the work due to the interaction terms changes significantly.Most importantly, it increases for the hot as well as the cold bath contributions.For the OLC the positive work output due to the interaction with the hot bath is easily overtaken by the consumption due to the cold bath modulation (blue arrow in Fig. 6).On the other hand, the well-chosen overlap of the phases can maintain overall positive work output (green arrow in Fig. 6), leading to higher power and efficiency as explained earlier.
These results suggest that a systematic optimization of the coupling modulation speeds and timings should prove fruitful -Bayesian optimization appears to be the method of choice, as the individual simulations are somewhat costly in computational resources.Another interesting question for future study would be how the work extraction from the interaction modulation behaves if coherences in the work medium are present with the aim to make coherences less detrimental to the engine performance.

V. INFLUENCE OF CYCLE LENGTH AND COUPLING STRENGTH
In the previous section, the variation of the modulation protocol has been investigated.Now we focus on how efficiency and power are affected by the overall length of the modulation cycle, and the coupling strength.Our In correspondence with standard thermodynamic pressure-volume diagrams, the conjugate variables ⟨∂X H⟩ vs. X are plotted for the limit cycle (see also Eq. ( 13)).Here, X is either the system modulation f (left column), or the interaction modulation h with the hot (middle) and cold bath (right).From such a representation (negative) work amounts to the enclosed area.The coloring of the area is chosen according to the scale presented in Fig. 8. Green shades correspond to positive (= usable) power output, whereas red shades indicate negative power.The upper row pertains to the OLC whereas the lower row displays the work diagrams for the best performing slow cycle with τtr = 0.12Θ and with the cold bath modulation shifted by τ = 0.18Θ (see the rightmost lower panel in Fig. 5).The plots show that the slower modulation and the time-shift in the cold modulation reduces the energy cost of the hot phase and boosts the energy output from decoupling the hot bath.Some diagrams feature "spikes" (upper left, lower right).These appear when the expectation value ⟨∂X H⟩ overshoots the final value at the end of a phase.This is a token of the departure from standard thermodynamics.Take the left-most column as an example.When the hot bath is decoupled at the end of the hot phase (II), the system energy is lowered slightly.This does not usually happen in standard thermodynamics, where the thermalization process is complete and monotonic.In addition, the "kink" seen in the lower middle panel is due to the onset of the cold phase before the hot phase has finished.
results are in qualtitative agreement with Ref. 31 , where a similar investigation has been performed with a larger parameter range for an oscillator work medium.
While cycle duration can easily be adjusted from a practical perspective, tuning the coupling strength, especially beyond the weak coupling regime, is significantly more experimentally challenging 75 .The results of our simulation for a rough scan over these two parameters are shown in Fig. 8.In order to keep things as transparent as possible, the simulation was done for the OLC modulation protocol without a σ x contribution in the system Hamiltonian (see also Fig. 2).
As seen in the upper left panel of Fig. 8, the efficiency of the engine improves with weaker coupling and longer cycles, consistent with the findings of Ref˙3 1 .This agrees with the intuition that a weaker coupling amounts to a shallower slope of the coupling and decoupling from the environments.The same holds true when increasing the cycle length, eventually approaching the quasistatic limit.However, the power output (upper right panel in Fig. 8) is not trivially correlated with efficiency.Instead, we find an optimal coupling strength of δ ≈ 0.3.The decrease for weak couplings is evident.In the limit of zero coupling, no work can be extracted.The decline for larger coupling strength is due to the larger losses associated with switching on and off the environments (bottom left panel of Fig. 8) while the system yield saturates (bottom right panel).This behavior agrees with the findings of Ref. 31 , where a maximum in the power output with respect to the coupling strength was observed, too.
For a fixed coupling strength, our results indicate that there is also an optimal cycle length with respect to the average power output.It is evident that in the very slow driving regime (large Θ, quasistatic regime) the average cycle power should be small and vanishes in the limit Θ → ∞.The behavior for small Θ, i.e. fast driving, is FIG.8.The average total power P and the efficiency η are shown in the upper row.The interaction and the system contribution which make up the total power are shown in the lower row.At the gray dots, we have exact data from our simulation.The values in between were obtained by linear interpolation.As they serve as the basis for the contour lines, the plots should be considered as a rough estimate of the landscape while, still, showing the major features.These figures can be compared to Figs. 3, 4 and 5 of Ref. 31 where qualitatively similar results have been found.not as obvious.Referring to the bottom right panel of Fig. 8, we see that the average system power increases while shortening the cycle.This means that even though there is only little time for the system to be influenced by the baths, i.e. the system picks up/gives off only a small amount of energy, due to the shortening of the cycle, the power still increases.It is sensible that for very short cycles, the work done by the system modulation scales linearly with the cycle length, thus leading to a finite power in the limit of Θ → 0. On the other hand, a faster change of the system-bath interaction results in larger energy loss due to the coupling to the baths being switched on and off.Summing up these two contributions yields a net power decrease when shortening the cycle.

VI. METHODS
Our approach to determine the dynamics and to access bath observables of the microscopic open quantum system model Eq. ( 1) relies on the stochastic description of the dynamics by means of the Non-Markovian Quantum State Diffusion (NMQSD) equation 76,77 , and its numerical treatment using the HOPS 48,61,62 .For a comprehensive reading, we briefly review these methods before turning to the computation of bath observables in Section VI C and following sections.

A. Non-Markovian quantum state diffusion (NMQSD)
In the usual model for open quantum systems 78 each bath consists of an infinite set of harmonic oscillators and an interaction Hamiltonian that is linear in the ladder operators, with The a λ are bosonic annihilation operators acting on the nth bath Hilbert space, g (n) λ denotes the individual coupling strengths, and L n (t) is the system coupling agent.For NMQSD they can be considered arbitrary, time dependent, not necessarily Hermitian operators acting on the system Hilbert space, mediating the influence of the nth bath.The bare system dynamics is governed by the (possibly time dependent) system Hamiltonian H S (t).Quite general baths can be mapped to this model in the thermodynamic limit, when the bath consists of an infinite number of non-interacting sub-systems 79 .Despite the simple harmonic structure of the bath, Eq. ( 14) is generally very hard to solve in a non-perturbative manner 66 .The NMQSD approach tackles that problem by recasting Eq. ( 14) into a stochastic differential equation on the level of the system Hilbert space, in which the bath degrees of freedom are accounted for by Gaussian stochastic processes and a convolution term.Note that in contrast to many of the perturbative approaches, the formalism is agnostic to a possible time dependence of the system Hamiltonian H S (t) and/or the coupling operators L n (t).
For the sake of simplicity, we review the NMQSD (and the HOPS) method for a single bath -the generalization to several baths is straightforward (see A). Furthermore, it is sufficient to treat the zero temperature initial condition only, i.e. a bath which is initialized in its ground state |0⟩.The effect of finite temperature can be accounted for correctly by an additional stochastic Hermitian contribution to the system Hamiltonian 48,62 .
To derive the NMQSD equation, the bath degrees of freedom are written in a basis of non-normalized Bargmann coherent states 80 .These are given by |z λ ⟩ = exp za † λ |0⟩ for a single bath mode.The total systembath state, thus, takes the form where z is the vector of all coherent state labels z λ , one label for each environmental oscillator.The properties of Bargmann coherent states make the relative system states |ψ(t, z * )⟩ a holomorphic function of z * .After transforming Eq. ( 14) into the interaction picture with respect to H B and using the properties of the coherent states ( ⟨z λ |a λ |ψ⟩ → ∂ z * λ ⟨z λ |ψ⟩, ⟨z λ |a † λ |ψ⟩ → z * λ ⟨z λ |ψ⟩) we deduce from the Schrödinger equation of the total state an equation for |ψ(z * , t)⟩ The key to eliminate the infinite sum is to shift perspective and regard |ψ(t, z * )⟩ to be a functional of the Gaussian stochastic process For a product initial condition |ψ S (0)⟩ |0⟩ B , this leads to the NMQSD equation 76,77 When relating these pure states to the total system-bath state, we can read the integral in Eq. ( 16) in a Monte-Carlo sense.From the Gaussian distribution of the coherent state labels z λ it follows that η t is a Gaussian stochastic process which obeys This shows that in the NMQSD formalism the microscopic structure of the bath, i.e. the values of the individual coupling strengths g λ , is fully accounted for in terms of the zero temperature bath correlation function (BCF) as introduced earlier in (7).In the continuous limit (infinitely many bath oscillators), the spectral density becomes a continuous function of ω and the BCF decays to zero for τ → ∞.Note that Eq. ( 19) reveals that the propagation from 0 to t requires knowledge of the bath correlation function from 0 to t, only.
Tracing out the bath degrees of freedom amounts to averaging the dyadic product of stochastic pure states for many realizations of the stochastic process and, thus, the reduced state reads Notably, Eq. ( 19) does not preserve the norm of the state vector, leading to flawed stochastic convergence of Eq. ( 22).This can be overcome within the nonlinear NMQSD formalism 76 From the solutions of the nonlinear NMQSD Eq. ( 26), the reduced system state is recovered through so that all trajectories contribute with "equal weight".
B. Hierarchy of stochastic pure states (HOPS) In equation ( 19) the bath degrees of freedom have been removed from explicit consideration, replacing them with a Gaussian stochastic process and a rather complicated term containing a convolution integral with a functional derivative In special cases, there are analytical and perturbative expressions for this term 76 , but we keep the approach as general as possible and instead choose a numerical avenue.
The key to proceed is to cast the complicated term containing the functional derivative into a set of auxiliary states 48,61,62,81 .As is done in the HEOM approach, we expand the BCF into exponentials α(τ and kµ (30)   we can define the kth hierarchy state The origin of the normalization chosen in Eq. ( 31) is the desire to give all hierarchy states the same unit and to formulate the final HOPS equations unified into one equation in an enlarged Hilbert space as is done in Ref. 82 .This is a very convenient choice, but other choices of normalization may be made.For a hierarchy state, the following equation of motion can be derived 61,62 where k = (k 1 , k 2 , . . ., k M ) with k µ ≥ 0 is a multi index and (e µ ) ν = δ µ,ν .The norm of the multi index k is denoted |k| = µ k µ .We call this norm the hierarchy level of ψ k .The zero-level hierarchy state corresponds to the pure state trajectory obeying Eq. ( 19), i.e, |ψ(η * t , t)⟩ ≡ ψ 0 .Since the HOPS does no longer involve the functional derivative, we omit the explicit dependency on the stochastic process in the notation for the elements |ψ⟩ k of the composite hierarchy state (. . ., |ψ⟩ k , . . .).Eq. ( 32) is the stochastic Hierarchy of Pure States (HOPS) -each state couples to states of neighboring level.This is similar to the HEOM approach 56 , but with the advantage of reducing the dimensionality of the auxiliary objects from dim H S 2 to dim H S by treating pure states instead of density matrices.
It has been shown that the HOPS for the nonlinear NMQSD equation is obtained by, again, replacing η → η and L † → L † − L † t 48,61,62 .When truncating the hierarchy, the resulting finite system of nonlinear differential equations with time dependent coefficients can be solved numerically using standard ODE-solvers.By increasing the cutoff level, the numerical error can be made arbitrarily small 62 .There are various truncation schemes available 62,81,83 .As we are not resource limited in this work, a simple simplex truncation condition was used, i.e. |k| ≤ k max .We show in the next section that they can be of great use beyond the mere calculation of the root stochastic pure state |ψ(η * t , t)⟩ = ψ 0 .The HOPS formalism presented here relies on an exponential representation of the BCF.Finding such multiexponential representations is a common problem in the field of open quantum systems 51,84,85 .For an Ohmic spectral density (SD) with exponential cutoff (as used in this work, see Eq. ( 8)) J(ω) ∼ ω e −ω/ωc (33)   the resulting BCF decays algebraically.Thus, any finite sum of exponential terms will only ever serve as an approximation.However, by numerically fitting the exponential expression to the BCF, the algebraic decay is well re-produced over several orders of magnitude using only a few exponential terms 48,81 .This allows us to faithfully propagate the open quantum system efficiently.
C. Bath energy change

Preliminaries
After reviewing the fundamentals of the NMQSD equation and the HOPS approach in Sections VI A and VI B, we are now in a position to calculate the change of the bath energies -a novel theoretical aspect within the NMQSD/HOPS framework.
As elucidated before, the NMQSD equation determines the dynamics of the stochastic trajectory |ψ(η * t , t)⟩, which is a projection of the global state of system and bath onto a coherent bath state.In this way, the complete information about the global state is encoded in the ensemble of trajectories.Accessing individual properties of a single bath oscillator would be challenging, since the influence of the bath as whole is mediated and encoded through the collective stochastic process.Still, we can expect to be able to access collective bath quantities, such as λ , i.e. a part of the interaction Hamiltonian of the nth bath in Eq. (14).
In order to demonstrate and elaborate further this general idea, we determine the change of the bath energy The choice of the negative sign ensures J is attributed with the flow of energy out of the bath.For the time being, we stick to the simplest version of the general model Eq. ( 14) with a single bath at zero temperature and no explicit time dependence of the Hamiltonian While the bath energy expectation value ⟨H B ⟩ is generally infinite, ∂ t ⟨H B ⟩ will remain finite.
From the Ehrenfest theorem, it follows with where we switched to the interaction picture with respect to H B to also connect to the NMQSD formalism.Thus, the energy flow out of the bath is We will remain in the interaction picture and drop the subscript in the following.Notably, the expectation value L † ∂ t B(t) can be connected to the first-level auxiliary states of the HOPS.As for the derivation of the NMQSD equation (Section VI A), using non-normalized Bargmann coherent states |z⟩ we write and find which yields an expression that only depends on the stochastic trajectory and the BCF.
Within the HOPS formalism (Section VI B), where we have assumed α(t) = µ G µ e −Wµt , the energy flow can thus be expressed in terms of an ensemble average involving the first-level hierarchy states |ψ eµ (t)⟩,

Finite temperature
In the general NMQSD/HOPS formalism, finite temperature can be included in many different ways 76,77,[86][87][88] .Using HOPS, it is most convenient to include finite temperature through a stochastic (Hermitian) contribution to the system Hamiltonian 48,62 : for an arbitrary observable O, not necessarily of the system, its expectation value at time t can be written as where U (t) denotes the global unitary time evolution operator and ρ S ⊗ ρ B is the factorized system-bath initial condition.For a Gibbs state ρ B = e −βHB /Z in its Glauber-Sudarshan P representation we may write where the coherent states have already been expressed in terms of the multimode shift operator Further, n λ denotes the Bose-Einstein distribution n λ = n(βω λ ) = (e βω λ −1) −1 .Using the unitarity of D(y), cyclic permutations under the trace and reading the y λintegrals in a Monte-Carlo sense as an average over the Gaussian y-coherent state labels results in This expression involves the propagation of an initial product state where the bath is now in its ground state and can, thus, be calculated using the zero-temperature-HOPS introduced earlier.Note, however, that the time evolution is determined by the stochastically shifted Hamiltonian H(y, t) = D † (y)H(t)D(y).In our case, referring to the NMQSD formalism (open system Hamiltonian ( 14) in the interaction picture with respect to H B ), the transformation simply adds to the original Hamiltonian H(t) a stochastic system contribution The function depends on the values of y, which are complex valued Gaussian random variables with M(y λ ) = 0 = M(y λ y λ ′ ) and M(y λ y * λ ′ ) = n λ δ λ,λ ′ .Thus, ξ(t) can be seen as a Gaussian stochastic process completely specified through its moments M(ξ(t)) = 0 = M(ξ(t)ξ(s)) and closely related to the thermal part of the usual bath correlation function.Specifically, for the change of the bath energy we have O = L † ∂ t B(t).Its transformed expression becomes It follows that for a non-zero temperature initial condition we obtain where the expectation value on the right-hand side has to be taken with respect to states that have been propagated with a zero-temperature initial condition and a system Hamiltonian with additional stochastic shift H shift S (ξ, t).Finally, within the HOPS formalism the change of the bath energy takes the form The appearance of the derivative ξ(t) of the stochastic process may look troublesome.However, if the autocorrelation function M(ξ(t)ξ * (s)) is twice differentiable, as in our case, the sample trajectories are smooth.When sampling such trajectories using Fourier-methods 81,89 , their time derivative can be calculated equally well.

N baths in the nonlinear formalism
Referring to the general Hamiltonian (14) with N baths, the NMQSD/HOPS formalism can be generalized with no major difficulties (see A). From the N BCFs α n (τ ) one deduces N stochastic processes η * n (t) and their thermal correspondent ξ n (t).Further, the hierarchy index vector k generalizes to N such vectors k = (k 1 , . . .k N ).Employing the nonlinear NMQSD equation to benefit from the superior stochastic convergence, we arrive at the following expression for the change of the energy of the nth bath The average M is taken with respect to the 2N stochastic processes (η * 1 , . . .η * N , ξ 1 , . . .ξ N ).The HOPS vector (. . ., |ψ⟩ k , . . . ) is propagated using the nonlinear HOPS (see Sections VI A and VI B).

D. Interaction energy
The interaction energy between the system and the nth bath is directly accessible by the HOPS formalism, too.From the expression and the derivation shown in Section VI C, one immediately finds where, again, the average M is taken with respect to (η * 1 , . . .η * N , ξ 1 , . . .ξ N ) and nonlinear HOPS determines the evolution of ψ k (t) .

E. Total power
Since the derivations in the previous sections did not include any explicit time derivatives of the system Hamiltonian H S and the coupling operators L n , we can simply substitute H S → H S (t) and L → L(t) in any of the expressions.In case of such explicit time-dependencies of the open system Hamiltonian (1), the total power will in general be non-zero.It can be calculated using the HOPS by simply replacing L(t) with L(t) in Eq. ( 55).

F. Other collective bath observables
Finally, for the record, we note that the dynamics of any observable that can be expanded in terms of the collective bath operators B n and B † n of Eq. ( 14) can be calculated using the HOPS formalism.It turns out that the expectation value of an operator O with normal ordering in B n and B † n , i.e.
where F acts on the system only, can be expressed in terms of auxiliary states of level l and m.Details are given in B.
The expectation value of the interaction energy discussed above is a particular example of this kind.

VII. CONCLUSION AND OUTLOOK
In this article, we expand the numerically exact HOPS method and apply it to the dynamics of a quantum heat engine.Both the work medium and bath degrees of freedom are explicitly included.Since the approach allows for arbitrary time-dependent couplings and system Hamiltonians, we are able to tackle realistic quantum thermodynamic cycles.The HOPS has originally been developed for the numerically exact description of the reduced dynamics of open quantum systems.In this paper, we extend the method to also access the dynamics of relevant bath observables.In particular, we are able to follow the various energy contributions in time, i.e. the system, interaction, and bath energies.On such a global scope it is then straightforward to define power and efficiency based on well-defined energy expectation values, circumventing ambiguities that may arise when working on a reduced level (e.g. with master equations).We find that the energetics of a quantum heat engine can be highly sensitive to the details of the modulations during the cycle.Even slight changes in the protocol can significantly affect the power and efficiency of the machine.For instance, overlapping phases can lead to increased power output and efficiency.
We illustrate our approach to quantum heat engines by taking a closer look at the Otto cycle with a qubit as work medium and two Ohmic baths.Special emphasis is placed on the per-phase relaxation dynamics of the qubit, showing once again that the steady state does not necessarily correspond to the Gibbs ensemble of the instantaneous system Hamiltonian.Furthermore, a detailed discussion of the behavior of the individual energy contributions of the global Hamiltonian with a fine resolution on the phases of the cycle is given.The loss in overall energy quantifies the extracted work and yields the cycle averaged power.Viewing the energy loss from the hot bath as the energy needed to drive the cycle, the efficiency is calculated rigorously.The dependence of power and efficiency on the parameters of the model is addressed carefully, investigating a broad range.We find that when keeping the overall coupling strength and cycle length fixed, altering the speed and timing of the modulation protocol can yield significant enhancements.Notably, for a particular overlap between the hot and the cold phases, we discover that the work contribution due to the modulation of the interaction Hamiltonian becomes positive.Thus, the process of coupling to and decoupling from the baths itself contributes positively to the amount of work done by the engine.Using thermodynamic work diagrams, a helpful visualization of these findings is possible.Finally, we discuss in detail how the coupling strength and the cycle length influence power and efficiency.As expected, the efficiency increases when moving towards the quasistatic case.For the total power, which is composed of competing contributions from the system and the interaction part, we find that there is an optimal choice for both parameters.In general, our thorough analysis, enabled by the HOPS method, reveals that thermodynamic properties of a quantum heat engine are very sensitive to the details of a fully quantum dynamical treatment, which are often lost by a too simplistic reduced description.Finally, we want to emphasize that the HOPS is equally applicable to other work media, non-Ohmic spectral densities, arbitrarily low temperatures, and any time-dependent system Hamiltonian and coupling, making it a universal method for describing general quantum dynamical (thermodynamic) scenarios.
For concrete applications, the modulation protocol with the highest efficiency and the largest power has yet to be found, depending on specific experimental parameters.Our investigations, however, show that such an optimized protocol will most likely include overlapping phases.Also, instead of introducing a semi-classical external modulation, autonomous heat engines with an explicit quantum energy storage would be of interest in our HOPS framework 1 .Such a more explicit inclusion of the work extraction process would also be necessary to answer the question of how the work generated by such a single-qubit heat engine can actually be used or measured.All thermodynamic quantities presented in this article have to be understood in the usual sense of quantum expectation values.It is still an open question, and the subject of ongoing research, to what extent the work extracted from single-shot experiments is measurable and deterministically usable at all 21,[90][91][92][93][94] .On a more fundamental level, we are confident that the possibility to access bath observables will allow for further insights in the connection between strongly coupled open quantum system dynamics and thermodynamics.With the help of the extended HOPS method, recent proposals for entropic quantities and generalizations of the thermodynamic laws at strong coupling can be tested for specific, realistic model systems far from the weak coupling limit [95][96][97][98][99][100][101][102][103] .
Further improvements of the method are desirable and are subject of ongoing research.In particular, depending on the details of the model, reaching a limit cycle in the HOPS simulations can be computationally expensive and the initial, transient dynamics is of little interest with respect to the thermodynamic properties of the machine.Therefore, a more direct access to the steady-state behavior of the heat engine would be valuable.The modulation protocol in this work was still slow enough so that the work medium could reach a steady state in most of the phases.It remains to conduct a more detailed study for faster protocols with off-diagonal system modulations to explore the effects of coherences in the work medium, especially with respect to the shifted (overlapping) protocols.
Finally, it would be fascinating to investigate manybody work media using the methods suggested in this work.

FIG. 2 .
FIG. 2. The modulation protocol31 (upper panel) of the Otto-like cycle (OLC) shows the separated cold (blue line and background) and hot (red line and background) phases.The system is modulated (black line) during the phases where the interaction with both baths is off.For that protocol, the reduced dynamics of the qubit in terms of the Bloch vector is shown in the middle (sx = 0.3Ω/2) and the lower panel (sx = 0).The length of the Bloch vector (indicating the purity of the state) is shown in black.The horizontal thick lines correspond to the Bloch vector components of the Gibbs state of the qubit with respect to the current system Hamiltonian and the temperature of the corresponding bath.The vanishing x and y components of the Gibbs state in the lower panel are not shown (parameters: δ = 0.7, ωc = Ω, T cold = Ω/2, T hot = 4Ω, Θ = 60/Ω, τtr = 0.06Θ).
FIG.7.In correspondence with standard thermodynamic pressure-volume diagrams, the conjugate variables ⟨∂X H⟩ vs. X are plotted for the limit cycle (see also Eq. (13)).Here, X is either the system modulation f (left column), or the interaction modulation h with the hot (middle) and cold bath (right).From such a representation (negative) work amounts to the enclosed area.The coloring of the area is chosen according to the scale presented in Fig.8.Green shades correspond to positive (= usable) power output, whereas red shades indicate negative power.The upper row pertains to the OLC whereas the lower row displays the work diagrams for the best performing slow cycle with τtr = 0.12Θ and with the cold bath modulation shifted by τ = 0.18Θ (see the rightmost lower panel in Fig.5).The plots show that the slower modulation and the time-shift in the cold modulation reduces the energy cost of the hot phase and boosts the energy output from decoupling the hot bath.Some diagrams feature "spikes" (upper left, lower right).These appear when the expectation value ⟨∂X H⟩ overshoots the final value at the end of a phase.This is a token of the departure from standard thermodynamics.Take the left-most column as an example.When the hot bath is decoupled at the end of the hot phase (II), the system energy is lowered slightly.This does not usually happen in standard thermodynamics, where the thermalization process is complete and monotonic.In addition, the "kink" seen in the lower middle panel is due to the onset of the cold phase before the hot phase has finished.
which naturally implements importance sampling of the trajectories.The formal expression s.