Autonomous quantum heat engine based on non-Markovian dynamics of an optomechanical Hamiltonian

We propose a recipe for demonstrating an autonomous quantum heat engine where the working fluid consists of a harmonic oscillator, the frequency of which is tuned by a driving mode. The working fluid is coupled two heat reservoirs each exhibiting a peaked power spectrum, a hot reservoir peaked at a higher frequency than the cold reservoir. Provided that the driving mode is initialized in a coherent state with a high enough amplitude and the parameters of the utilized optomechanical Hamiltonian and the reservoirs are appropriate, the driving mode induces an approximate Otto cycle for the working fluid and consequently its oscillation amplitude begins to increase in time. We build both an analytical and a non-Markovian quasiclassical model for this quantum heat engine and show that reasonably powerful coherent fields can be generated as the output of the quantum heat engine. This general theoretical proposal heralds the in-depth studies of quantum heat engines in the non-Markovian regime. Further, it paves the way for specific physical realizations, such as those in optomechanical systems, and for the subsequent experimental realization of an autonomous quantum heat engine.


Introduction
The global societal impact of heat engines has been enormous through history.Various types of heat engines such as those powering pumps and motors have been extensively Consitently, the heat obtained from the hot reservoir by the system Q h exceeds the heat released from the system to the cold reservoir Q c .
In the spirit novel quantum technologies and inspired by the simple QHO model, we propose a way of utilizing systems realizing an optomechanical Hamiltonian [26,27] as an autonomous quantum heat engine.In a prototypical optomechanical system, an optical cavity mode with frequency ω a is coupled to a slower mechanical mode with frequency ω b so that the frequency of the optical mode is modulated by the mechanical displacement x, approximately as ω a (x) ≈ ω 0 a + x ∂ωa ∂x x=0 .Let the optical cavity be further coupled to two heat reservoirs with peaked spectral densities such that the center angular frequencies of the hot reservoir ω h and cold reservoir the ω c fulfill ω h > ω a (x) > ω c .Assigning temperatures T h and T c to the hot and cold reservoirs, respectively, such that T h > T c , leads to a device where the state of the optical cavity may undergo the quantum Otto cycle.The cavity mode at ω a will traverse between the characteristic frequencies of the heat reservoirs driven by the mechanical mode.The optical mode will therefore pump photons from the high-frequency high-temperature environment to the low-frequency low-temperature environment, leading to energy boost to the mechanical mode owing to energy conservation.In this paper we aim to show, that the thermal energy released over the cycle is transformed, at least partially, into the amplitude of the coherent state of the mechanical mode, thus available for useful work in the realm of quantum devices.
For the sake of simplicity, we refer to the modes used in our proposal as the optical and the mechanical mode, but our theoretical proposal is hardware agnostic.If there are two modes and reservoirs that fulfill the discovered conditions, such a physical system is prone to be used as an autonomous quantum heat engine.
The proposed device is a heat engine in the most fundamental sense-a heat engine is, by definition, a device that transforms heat into deterministic mechanical motion [28,29].It should be noted, however, that this is not the first time an optomechanical system is used as the basis for a QHE.References [30][31][32] extensively analyze the possibility of realizing a coherently driven quantum Otto cycle in an optomechanical system.The above-mentioned references [19][20][21] utilize optomechanical systems, but rely on periodic incoherent thermal drives.Consequently, the coupling to the heat bath needs to be controlled by some external method, possibly consuming significant power.In contrast to such models, we propose a method which is fully free of temporally varied external controls and aim for a device that can work fully autonomously drawing power from constanttemperature heat reservoirs with static coupling spectra.Further, the references [22][23][24] study theoretical models of autonomous QHEs based on optomechanical systems on the level of Markovian master equations.In general, it may be difficult to justify the Markov approximation if the heat reservoirs are spectrally separated and strongly structured in the relevant frequency range of the system, as illustrated in Fig. 2. Thus, in contrast to the previous works, we do not invoke the Markov approximation, but rather consider the non-Markovian dynamics.Consequently, we rigorously take into account the spectral shape arising in all bosonic thermal sources due to quantum effects.As a result, our QHE exhibits autonomous cyclic dynamics in contrast to the traditionally classified completely continuous QHEs and externally controlled reciprocating QHEs.
This paper is organized as follows.
In section 1, we analyze the the proposed quantum heat engine by building an approximate analytical model for it.Subsequently, we find reasonable estimates for the net output power and efficiency of the device along with some physical intuition.Next, in section 2, we develop a more detailed quasiclassical model based on Heisenberg-Langevin equations to describe the intricate dynamics of the QHE.We further outline how to numerically solve the resulting coupled non-linear stochastic integro-differential equations.In section 3, we investigate the temporal evolution of the device according to the quasiclassical model.We obtain estimates for the power and efficiency, and analyze the power fluctuations inherently present in any QHE.Finally, in section 4, we summarize our findings, draw conclusions, and provide our insight for the fruitful future research on this autonomous quantum heat engine.

Analytical Model of the Quantum Heat Engine
In order to gain physical intuition of the working mechanism of our device, we estimate the efficiency of the proposed QHE by an approximate analytical model.We begin with the usual model of the adiabatic quantum Otto cycle [25] and proceed by approximately taking into account the finite widths of the reservoir spectra.First, let us assume that the heat reservoirs have a step-function-like spectra rather than some experimentally encountered smooth peaked line shape, as depicted in Fig. 2. In this case the coupling to each reservoir is sharply turned on and off without any residual coupling when the frequency of the optical mode moves into and out of the corresponding spectrum, respectively.The step function spectra are centered around ω h and ω c , respectively, and are assumed to have equal widths l h = l c .Our optical mode frequency ω a (t) oscillates about ω 0 a , which is assumed to be centered in between ω h and ω c .We assume that the optical mode frequency ω a (t) traverses from one reservoir to the other periodically in time with the period τ b = 2π/ω b .Further, we define the interaction time τ , which we assume to be equal for both reservoirs since the frequency ω a (t) spends an equal amount of time per cycle in each frequency range defined by the step function of each reservoir.(See Fig. 2 for an illustration and Supplementary materials for theoretical details.)We further define the average angular frequency over the interaction period ωh/c a , which is merely the time average of ω a (t) over τ , respectively for each reservoir.We use this average value as the effective angular frequency of the optical mode in the corresponding region.Note that we allow the possibility of ω a (t) reaching values beyond the extents of the step function spectra.
To estimate the net work output of this QHE, we need to compute the mean number of photons transferred from the hot to the cold reservoir per cycle.Let the optical mode be coupled to each of the reservoirs so that the dissipation rates are Γ h and Γ c , respectively.The mean number of thermal photons in the optical mode under the steadystate operation after the interaction period τ with the hot and cold reservoir is given, respectively by ) where − 1 is the mean photon number of the optical mode if it is fully thermalized with the hot or the cold reservoir, respectively, ℏ is the reduced Planck constant, and k B is the Boltzmann constant.This set of equations can be solved for the steady state photon numbers given by ) Note that we have used an equal interaction time τ for the reservoirs, without loss of generality since the only the product Γ h/c τ affects the end result, and hence different interaction times may be compensated by the dissipation rates.
In this paper, we follow the traditional division of work and heat in quantum systems.Work is considered to be the energy related to changing the frequency of the optical mode operating here as the working fluid and heat is considered to be the energy related to changing the population of the fluid.We assume that all the net energy released in a cycle from the reservoirs is transferred to the mechanical mode and that the photon number of the optical mode only changes during the interactions with the reservoirs.In this case, we may compute the average energy output per cycle and the average output power as where ∆ω a = ωh a − ωc a and ∆n is the difference of the mean photon number in the optical mode during a cycle.In the limit of ideal adiabatic operation, the efficiency of this device is defined in the usual way, which we find to coincide with the Otto efficiency [25] given by the effective frequencies as η eff = ∆W/∆Q h = E cyc /(ℏω h a ∆n) = 1 − ωc a /ω h a .This efficiency is smaller than the often-used efficiency given by the peak-to-peak compression ratio κ = (ω 0 a + ∆ω a /2)/(ω 0 a − ∆ω a /2): η eff ≤ η max = 1 − κ −1 .This model already has a considerable number of adjustable parameters.We aim to illuminate the meaning and implications of the most important parameters here.In Fig. 3, we show some derivative parameters of interest as a function of more fundamental, system parameters.In Fig. 4 we study the power in Eq. (1.6) and various limits as a function of the most interesting parameters.These graphs illustrate the dynamics and parameter ranges of the device quite well.Note that in realized optomechanical systems, the frequency of the mechanical mode is typically much lower than what is used here, but such a case will merely lead to closely adiabatic dynamics with respect to the thermalization and is theoretically easier and less interesting to analyze than the case studied here.
First, we observe from Fig. 4(a) that the output power increases up to a limit with increasing ω b .This limit is a result of the fact that there is a trade-off between increasing speed of extracting E cyc and decreasing the interaction time τ resulting in lower ∆n and hence lower E cyc as shown in Fig. 3(c).One of the most important fundamental parameters here is how much we vary ω a (t) over the cycle.We define ∆ω a as the peakto-peak modulation amplitude of ω a (t), i.e., ω a (t) will oscillate between ω 0 a − ∆ω a /2 and ω 0 a + ∆ω a /2.In Fig. 4(b), we show the output power as a function of ∆ω a .We clearly observe the limits of too low and too high modulation amplitudes, i.e., not reaching the reservoir at all or overshooting.This phenomenon is also observed in Figs.3(a) and 3(b).In Fig. 4(c) we see compatible behaviour as a function of the difference between the reservoir frequencies.Finally, in panel (d) of Fig. 3 we plot the efficiency of the heat engine.We note that utilizing the whole width of the step function reservoir in frequency gives the best efficiency, as expected.Increasing the peak-to-peak modulation amplitude and the energy gap between the reservoirs would increase the compression ratio, therefore allowing the increasing of the efficiency closer to the Carnot efficiency η C = 1 − T c /T h ≈ 80%.and (d) indicate the threshold where the amplitude is just high enough for the optical-mode frequency to enter the non-vanishing regions of the step-function-like power spectra of the reservoirs (left, dark red) and the threshold where the frequency leaves at its extrema the non-vanishing spectra.The parameters used here are the following: To finalize our examination of the analytical model, we estimate the dissipation rate and the resulting quality factor of the mechanical mode required for steady-state operation of the QHE.Note that in practice, this effective dissipation of the mechanical mode is a result of using the power flowing into the mechanical mode in a desired use case.
We assumed above that all the net thermal energy extracted over a cycle is transferred to the coherent motion of the mechanical mode.Under this assumption, which we justify below with a numerical model, we can match the power calculated above and the dissipated power from the mechanical mode in order to estimate the required quality factor of the mechanical mode.The average dissipated power of the mechanical mode Assuming that the mechanical mode modulates the optical mode frequency by ∆ω a ≈ 2 √ n b g 0 , where g 0 is the optomechanical coupling constant [26] discussed below, one can solve for the unknown parameter n b .By setting the dissipated power equal to the power obtained in Eq. (1.6) and using Eq.(1.5), we obtain Γ b = 2∆ω a ∆ng 2 0 /(π∆ω 2 a ).We read the optimal values for ∆ω a , ∆ω a , and ∆n from Fig. 4 and set the optomechanical coupling to a reasonable value, g 0 = 0.01 × ω 0 a , to retrieve a numerical value of Γ b ≈ 2.6 × 10 −5 × ω 0 a , which translates to a quality factor of about Q b = ω b /(2Γ b ) = 960 for the mechanical mode.Note this estimation does not depend on how the excitations are removed from the mechanical mode.We merely state that one can extract energy at this rate from the the mechanical mode.The extracted energy is interpreted as the work output of our heat engine.Note that, even though in reality there is also incoherent thermal occupation in the mechanical mode, we have not considered it here, nor does the analytical model capture the possibility.This will be discussed more in section 3.
The analytical model presented here provides us intuition about the physics of the device and reasonable estimates on the output power and other parameters.It, however, assumes step-like power spectral densities of the baths and that the net thermal energy extracted from the reservoirs is fully transferred into the coherent motion of the mechanical mode.The analytical model does not allow us to investigate the dynamics of the whole system from arbitrary initial conditions.Of course, we also made a number of other rather crude approximations in order to simplify the calculations to the bare minimum.Below, we develop a more detailed theoretical description of the device allowing the investigation of the temporal evolution given arbitrary initial states, describing the optomechanical coupling more carefully and taking into account experimentally feasible spectral shapes of the reservoirs.

Heisenberg-Langevin Equations for the Quantum Heat Engine 2.1 Derivation of the Equations of Motion
As discussed above, we consider a system realizing an optomechanical Hamiltonian Ĥopt [26,27] of an optical mode with bare angular frequency ω 0 a and a mechanical mode with angular frequency ω b .The optical mode is coupled to two non-Markovian heat reservoirs [33] with peaked spectral line shapes at different central frequencies.The total Hamiltonian of the optomechanical system and its reservoirs is written as where where â and b are the annihilation operators of the optical and the reservoirs mode, respectively, g 0 is the optomechanical coupling constant [26], the operators ĥk and ĉk are the annihilation operators of the k:th modes at angular frequencies ω h,k and ω c,k in the hot and cold reservoirs, respectively, and the corresponding coupling constants between the reservoir modes and the optical mode are λ h,k and λ c,k , respectively.Thus, Ĥh and Ĥc describe the reservoirs as infinite sums of bosonic modes.In addition, Ĥh int and Ĥc int describe the linear interaction between the reservoir modes and the optical mode.The Heisenberg equations of motion for the optical and mechanical modes, as well as for the reservoir degrees of freedom are given by Let us integrate out the reservoir modes from the above set of equations.To this end, we solve Eqs.(2.2c) and (2.2d) by explicit integration as and substitute the solutions into equation (2.2a).Consequently, we arrive at the Heisenberg-Langevin equation (HLE), where we define the memory kernels with J h (ω) and J c (ω) being the spectral functions of the hot and cold reservoir, respec-tively, and the noise operators are given by At this point, the Markovian approximation [33][34][35] is typically invoked such that the noise is replaced by white noise, characterized by ξ † h/c (t) ξh/c (t ′ ) = δ(t − t ′ ), where the average is taken over the reservoir degrees of freedom, and consequently the integrals in the above equations are reduced to linear dissipation [26,27,36].However, as described above, the spectral density of the environment is here not frequency-independent, but rather, has a peaked line shape.Therefore, we consider the full non-Markovian model [37].The noise operators defined by Eqs.(2.6) have non-local temporal correlation relations, and together with the memory kernels, they characterize the non-Markovian features of the reservoirs and of the resulting dynamics of the system.
The coupled Eqs.(2.2b) and (2.4) determine the dynamics of the optical and mechanical modes driven by the noise terms with damping arising from the memory kernel integrals.Unfortunately, these equations are challenging to solve exactly, even numerically, since they are non-linear operator equations of high dimension [26].However, we can proceed by decomposing the operators into classical and quantum components, â = α + δâ and b = β + δ b, and linearize in terms of the quantum components [26,38].Namely, we assume the quantum components to provide a small correction to the classical dynamics.Consequently, we obtain ) where G = g 0 α(t), we have transformed the noise operators ξh/c into stochastic complexvalued variables ξ h/c , and ].These equations are exact up to the linearization and describe the dynamics of the system under the influence of the reservoirs.The quasiclassical equations (2.7a) and (2.7b) retain the inherent non-linearity of the system, whereas the quantum equations (2.7c) and (2.7d) yield the quantum corrections.
The classical equations can be understood as a mean-field-like approximation to the full quantum equations (2.4) and (2.2b).Thus, they may provide a reasonable approximate solution in many physical situations.For instance, in reference [21], where a QHE based on optomechanical system pumped by periodic incoherent drive is studied, the difference between quantum and classical models is found to be nighly negligible.

Considerations for Numerical Solution of the Equations
To analyze the system dynamics based on the approximate Heisenberg-Langevin equations derived in previous section, we resort to numerics.In order to relax the requirements of the computational resources, let us investigate a quasiclassical [39] model based on Eqs.(2.7a) and (2.7b).Further, we are mostly interested in the non-linear dynamics of the system, captured by the classical equation, and have already made the assumption that the classical component of the dynamics is dominating.Using noise arising from thermal quantum fluctuations, the quantum statistics of the reservoirs is taken into consideration and a quasiclassical model for the system dynamics is achieved.Although we neglect the quantum fluctuations of the system itself [40], we consider the model to likely be accurate enough especially for average quantities since quasiclassical models similar to this have been shown to produce valuable insight into quantum systems interacting with heat sources [39].
Here, we numerically solve the coupled equations 2.7a and 2.7b in the time domain.Strictly speaking, stochastic differential equations (SDEs), and more importantly the numerical solution methods, are only defined for Wiener processes, i.e. white noise [41,42].Physicist and engineers have however stretched the limits of mathematical definitions, and considered equations with coloured noise obtaining physically meaningful results [41,43,44].Therefore, we follow these footsteps and proceed with our analysis assuming that in the sense of small enough time-steps our noise is locally white.We interpret our SDE in the sense of Stratonovich stochastic calculus and apply a second-order predictor-corrector method [41,42,45] to solve the above stochastic integro-differential equation [46].The coloured noise is generated from the power spectral densities, given below by Eq. 2.8, by filtering zero-mean Gaussian white noise with a filter, the frequency response function FRF h/c (ω) of which satisfies FRF h/c (ω) 2 = S h/c (ω).Note also that since we are dealing with a stochastic system, single realizations of the dynamics exhibit random behaviour with little conclusions to be drawn.Thus, we simulate the system with a large number of repetitions and compute ensemble averages of the quantities of interest.
Let us specify the reservoir spectra and the temperature dependence of the heat reservoirs.First, we consider the following power spectral density of the thermal fluctu-ations [47,48] in the reservoir: where F [f (t)] (ω) denotes the Fourier transform of function f (t) and where g h/c is a dimensionless constant scaling the coupling strength between the optical mode and the corresponding reservoir, γ h/c equals the full width of the spectrum at half maximum, ω h/c is the renormalized centre frequency, and T h/c is the temperature of the reservoir as defined above.The hyperbolic cotangent term in the above expression arises by directly calculating the variance of noise operators following the usual analysis of generalized bosonic heat reservoirs [47,49,50].The line shape given by J h/c (ω) is chosen merely to consider spectra which are peaked in a reasonable way without an attempt to optimize them for output power.
Let us make a minor remark about the parameters.We introduce two new parameter here, namely g h/c and γ h/c .Unfortunately, these are not directly quantitatively comparable to their counterparts in the analytical model.Where Γ h/c in the analytical model is the relaxation time constant of the optical mode coupled to the reservoir, g h/c describes the coupling strength with the reservoir.These are closely related, but not the identical, and can vary from physical realization to another.The full width at half maximum of the spectrum in Eq. (2.9) is given by γ h/c .This can be compared to the width of the step function spectrum, l h/c , more or less directly, keeping in mind that the spectrum here is smooth with polynomially decreasing tails, in stark contrast with the step function.(See Fig. 2 for comparison.)

Dynamics and Performance of the Quantum Heat Engine 3.1 Temporal Evolution and Output Power
In Fig. 5, we show the temporal evolution of the amplitudes α and β and excitation numbers n a = ⟨|α| 2 ⟩ and n b = ⟨|β| 2 ⟩ of the optical and the mechanical mode, respectively, assuming negligible dissipation on the mechanical mode.Here, the average is an ensemble average with respect to many independent realizations of the noise trajectories ξ h (t) and ξ c (t).This theoretical scenario allows us to gain insight into the operation principle of the quantum heat engine.The parameters are chosen to roughly optimize the power generation in the mechanical mode with a constant temperature difference between the reservoirs.Full systematic optimization of the parameters, being somewhat computationally intensive, is left for future work.
At time instant t = 0, the modes are initialized into coherent states corresponding to n a = 0.5 and n b = 39.The mean values are found by averaging over 1000 realizations of the noise.
Initially, the modes are set to coherent states with mean excitation numbers of n a = 0.5 and n b = 39, correspondingly.Firstly and most importantly, we observe from Fig. 5(b) that the initially prepared coherent state tends to be preserved in the evolution and even increasing its amplitude rather than winding towards a thermal state.The increasing amplitude of the oscillation is in agreement with Fig. 5(a) where the average excitation number of the mechanical mode steadily increases in time.These key observations justify the assumption we introduced for the analytical model that the net thermal energy extracted from the heat baths tends to convert into the coherent motion of the mechanical mode.
From the results of Fig. 5(c), we observe for the optical mode that instead of thermalizing to a steady state with a constant average photon number between the thermal occupation numbers N th h and N th c given by the Bose-Einstein distribution at the reservoir temperatures, the average photon number of the optical mode exhibits clear temporal oscillations.These oscillations take place between the boundaries set by the thermal photon numbers with a period given by the frequency of the mechanical mode ω b .This behaviour arises from driving the frequency of the optical mode by the mechanical mode between the peaked spectra of the hot and cold reservoirs that inject and absorb optical photons, respectively.In stark contrast to the traditional externally driven quantum heat engines however, driving is here purely caused by the internal dynamics of the device.In a qualitative agreement with our analytical model and Fig. 1, the average excitation number of the mechanical mode exhibits oscillations phase shifted by π/2 from those of the optical mode, i.e., when the optical mode has a peak in its average photon number, the mechanical mode is gaining energy at maximum speed from the optical mode.
Since the dissipation of the mechanical mode is neglected and the mechanical mode has a reasonably low amplitude, we do not observe the mechanical mode to stabilize into a steady state.However, we can estimate the upper limit of power, at which the mechanical mode accumulates energy.In Fig. 5(a), we employ a linear fit on the average excitation number of the mechanical and find that based on the slope of the fitted line, this device yields an average power of P = 8.63 × 10 −5 ℏ(ω 0 a ) 2 , which is clearly higher than the power achieved according to the analytical model in section 1.This is partially attributed to the slightly different parameter values and approximations used in the two models, but the main reason for observing a higher power output is, however, expected to be the initial thermalization transient experienced by the mechanical mode.From Figs. 5(c) and 5(d), it is clear that the optical mode reaches a state dominated by thermal noise very quickly, but the saturation of the mechanical mode may take significantly more time due to its remarkably weaker coupling to the thermal sources.Of course, we do not even wish the mechanical mode to reach a thermal state, but to retain its coherence as well as possible.
A reasonable estimate for the average thermal occupation in the mechanical mode after a long time is given by the weighted average is the mechanical-mode occupation fully thermalized at the reservoir temperature T h/c given by the Bose-Einstein distribution.This further supports the observation that increasing occupation well above n b = 40 must be at least partially coherent generation, as we would expect from Figs. 5(b) and 5(d).Therefore, it is reasonable to assume that the increasing occupation is partially due to coherent generation and partially due to thermalization.Let us next study the steady states of the mechanical mode with finite dissipation and after a long evolution in order to shed some light on the matter.In Fig. 6(a), we show the temporal evolution of the average excitation number of the mechanical mode for different levels of excess dissipation, i.e., non-vanishing quality factors.As expected in this case, the excitation number exponentially reaches a steady state where we constantly extract work from the mechanical mode.Note, that the work here is simply dissipated away, but in practice it can be released to a waveguide and measured.From Fig. 6(a), we note that the quality factors of the mechanical mode are on a similar level to those found with the analytical model above.
Let us estimate the dissipated power based on the dissipation rate by the method used above in the analytical model.However, in order to separate between work and heat in the mechanical mode, we make the reasonable assumption that after a long time the mechanical mode stabilizes into a thermal coherent state [51], and use the fact that the average occupation number of the thermal coherent state is given as a sum of the coherent and thermal occupations.We treat the average occupation number resulting from our simulation, n b , as the total occupation number of the thermal coherent state, and use the above-obtained estimate for the thermal contribution N th b .Thus the coherent contribution is given by n coh b = n b − N th b , which allows us to interpret the coherent dissipated power as work and compute the corresponding power output as P = Γ b ℏn coh b ω b .We show the obtained power in Fig. 6(b) for different levels of excess dissipation.We point out that for the highest efficiency with the mechanical-mode quality factor Q b = 1125, the thermal contribution in the mechanical mode occupation is only about N th b /n b = 14% of the total occupation number.Furthermore, the power output obtained here is less than half of the upper-bound estimate obtained above.This observation supports the hypotheses of the thermalization transient.
In addition, we find in Fig. 6(b) that we can vary the dissipation level quite significantly without much disturbing the power generation.As long as the dissipation is on a reasonable level compared with the power generation rate found above, the systems seems to be able to find a steady state.Bear in mind that this is a dynamic process involving non-linear back and forth coupling between the mechanical and the optical mode as a function of the occupation number n b .As opposed to the analytical model, this quasiclassical model takes into account the intricate dynamics owing to the optomechanical coupling.

Efficiency
From the above dynamics of our quantum heat engine, it is somewhat challenging to rigorously define and determine certain quantities traditionally studied in quantum heat engines, stemming from the lack of the external driving field.Typically, the work done on the working fluid by the external drive plays an important role in defining these quantities.As already mentioned in section 1, work is the energy related to changing the Hamiltonian, or frequency, of the working fluid whereas heat is the energy related to changing the populations.The issue arising from the lacking external drive is somewhat more pronounced in the quasiclassical model than in the analytical model since for the quasiclassical model we are truly solving the dynamics of the whole optomechanical system, whereas in the analytical model, we still treat them as separate systems.
Let us, nevertheless, try to estimate the efficiency of our quantum heat engine by using the method we utilized for the analytical model.We consider the optical mode to be the working fluid and assume that the thermal occupation of the mechanical mode does not significantly vary over a cycle once the steady state is reached.We further approximate that the oscillation in the occupation of the optical mode is almost solely due to the interactions with the thermal reservoirs.This is a very good approximation, since the mechanical mode varies the occupation of the optical mode very little over a cycle due to the relatively low optomechanical coupling, g 0 ≪ ω 0 a , and excitation level of the mechanical mode.(See Supplementary Fig. S(1).)In this scenario, we find the heat absorbed by the device from the hot reservoir per cycle, ∆Q h , by finding the average variation of the photon number in the optical mode over cycle.This can be computed from the maxima and minima of the photon number, n a , over many cycles, as depicted by the red and blue crosses in Fig. 5(c).The average work per cycle is simply given by the power: ∆W = P τ b .With this method we estimate the efficiency, η, shown in Fig. 6(b) as a function of excess dissipation levels of the mechanical mode.We note that the efficiency is qualitatively remarkably close to the one achieved in the analytical model, which would suggest that the analytical model captures the essence of the physics relatively well.As one may expect, the efficiency is far below the Carnot efficiency, η C ≈ 80%, and also significantly below the Otto efficiency given by the compression ratio η max ≈ 16%.
Finally, we highlight an important detail, i.e., we cannot claim that the heating and cooling processes of the working fluid are perfectly isochoric, nor can we say that the compression and decompression are fully adiabatic.If we take into account the slowly vanishing tails of the reservoir spectra, it becomes evident that there is constant heat exchange between the optical mode and the reservoirs over the cycle.These considerations inevitably lead to the fact that we cannot, strictly speaking, call this cycle an Otto cycle in its ideal form.The cycle of our device displays more continuous dynamics, not perfectly split into well defined phases.Moreover, we point out that the above discussion considers only the internal efficiency of the device.In an experimental realization, there is also direct heat leakage from the hot reservoir to the cold reservoir through spurious parallel channels to the quantum heat engine, decreasing the total efficiency in practice.

Stability and Power Fluctuations
In this section, we briefly analyze the ensemble fluctuations of the output power and the stability of the device with respect to some parameters.In Fig. 7, we compare the net output power obtained from the quasiclassical model with the analytical model of section 1 as a function of a few of the most important parameters.Interestingly, there is a relatively narrow range of optimal mechanical-mode frequencies in the quasiclassical model depicted in Fig. 7(a).In the analytical model, we do not restrict the mechanicalmode frequency in any way, and seems that it can be increased indefinitely, whereas in the quasiclassical model we clearly find a sweet spot for the frequency ω b .This discrepancy may be a result of the fixed peak-to-peak modulation amplitude ∆ω a assumed in the  We use cubic spline method to interpolate (dashed cyan color) the numerical data (cyan color dots) of the Heisenberg-Langevin equation model.Note that because of the differences in the models and their parameters, we do not expect the models to accurately quantitatively agree.
analytical model, whereas the amplitude in the quasiclassical model is a result of the temporal evolution determined by the other parameters.In Fig. 7(c), we find that the two models behave qualitatively similarly as a function of the difference between the centre angular frequencies of the reservoir power spectra.Note that the curve of the quasiclassical model does not display any sharp cusp as observed for the analytical model, owing to the continuous spectral functions in the quasiclassical model.The behaviour of the output power as a function of the hot-reservoir temperature is found in Fig. 7(d) qualitatively quite similar to the analytical model, as expected.However, the analytical model predicts a steeper increase in the power as function of increasing hot-reservoir temperature than the HLE model, which is attributed to the accordingly rising mechanical-mode thermal occupation in the HLE model, which does not contribute to the work.It is interesting to investigate the relative occupation number fluctuations of the mechanical mode, because it gives insight to the excitation counting statistics one would observe when measuring the mechanical mode.Further, because the output power is directly related on the occupation number, as defined above through dissipated power, the relative output power fluctuations of the QHE [52,53] are consequently directly related to the occupation number fluctuations.In Fig. 8(a) we illustrate the occupation number fluctuations by showing hundred different realization of the temporal evolution of the excitation number of the mechanical mode under excess dissipation characterized by the quality factor Q b = 1125.In Fig. 8(b), we display the probability density for the occupation number after reaching a steady state and compare it to the Poisson distribution.We use the Fano factor [54], which quantifies the deviation from a Poisson distribution, as a measure to compare the occupation number distribution to the Poissonian.The Fano factor is defined as F = σ 2 P / ⟨P ⟩, where σ 2 P and ⟨P ⟩ are the variance and the expectation value of the occupation number, respectively.For the data in Fig. 8(b), we find F = 0.27, which signifies far sub-Poissonian (F < 1) statistics of the occupation number fluctuations.This suggest a smaller output power fluctuations for the QHE than for a Poisson process.

Conclusions and Discussion
We have proposed and analyzed an autonomous quantum heat engine in the spirit of an Otto cycle, where one harmonic mode, the working fluid, traverses the cycle driven by another harmonic mode.We utilized an optomechanical Hamiltonian for this realization and consequently referred to the working fluid as the optical mode and the driving mode as the mechanical mode.However, any physical system realizing the optomechanical Hamiltonian in the appropriate parameter range will suffice our heat engine.
The feasible parameter range for the quantum heat engine is straightforward to find thanks to the analytical model we established by assuming that the net thermal energy extracted from the thermal reservoirs is fully converted into the coherent motion of the driving mode.By deriving a dynamical non-Markovian quasiclassical model and utilizing it in demonstrating the operation of the quantum heat engine, we justified the assumption made in the analytical model.We achieved output powers and efficiencies of the heat engine in a very good agreement between the two models.The efficiency obtained was roughly 10% and the output power was found to be a reasonable fraction of the energy quantum of the working fluid times the frequency of the driving mode.Importantly, the output power was found to monotonically increase with the temperature of the hot reservoir providing prospects for possible applications.
Our initial analysis showed that the output power of the quantum heat engine significantly fluctuates from noise realization to another.Although this fluctuation seems suppressed in comparison to Poisson processes, it calls for future studies on the noise spectrum of the coherent-field output of the heat engine, and possible ways to decrease the observed phase noise and power fluctuations.
Thanks to the versatility of the proposal, we consider a future experimental realization of our quantum-heat-engine proposal likely feasible, but requiring prior specific theoretical analysis for each experimental realization.Although we leave such specialized analysis for future work, we note that our example calculations show that the quality factor the the working-fluid mode may be of the order of hundred and that of the driving mode roughly thousand for the optomechanical coupling constant of roughly one percent of the angular frequency of the working fluid.In typical optomechanical systems, for example, such quality factors are straightforward to achieve, but the optomechanical coupling constant may be smaller, which calls for an optimization of the separation of the reservoir spectra [26,27,36].In addition to optomechanical systems, it is interesting to consider a realization of the quantum heat engine in superconducting circuits [14,16,21,55] which may offer more flexibility in engineering the optomechanical coupling constant [27,56,57] and angular frequencies of the two modes [57,58].For high-power operation, it may be interesting to consider engineered environments [18,59] which have recently been demonstrated to provide fast preparation of thermal states [60].
Even with the possible loss of quantum character due to our quasiclassical methods of solution, we deemed retaining the non-linearity and non-Markovianity as the important novelty of our approach, and although investigating the full quantum character of the device is left for future work, the results obtained here remain significant and interesting proof of concept.In the future, it is interesting to study several different possible experimental realizations of the proposed quantum heat engine and their advantages and disadvantages with respect to each other.In addition, the mitigation of the phase and amplitude noise in the coherent output field of the heat engine call for further studies together with possible use cases of this work obtained from the quantum heat engine.The interaction time τ of the analytical model for the hot reservoir in the most general case can be calculated, as depicted in Fig. S (5), from the inequality ω h − l h /2 ≤ ω 0 a + ∆ω a /2 sin(ω b t) ≤ ω h + l h /2.
Solving the two equations arising at the limit of equality, we find the conditions

Figure 1 :
Figure 1: Quantum Otto cycle considered in this work.The angular frequency ω a of the harmonic oscillator depicted as the quadratic potential traverses between values ω c and ω h at average photon numbers n c and n h , respectively.The work obtained from the cycle W out exceeds the work done on the oscillator W in .Consitently, the heat obtained from the hot reservoir by the system Q h exceeds the heat released from the system to the cold reservoir Q c .

Figure 2 :
Figure 2: Spectral densities of the reservoirs.(a) Experimentally realizable Lorentzian line shapes of the cold (blue color) and hot (red color) reservoirs.The vertical dotted lines denote the center angular frequencies of the cold and hot reservoirs, ω c and ω h respectively, and the unbiased angular frequency of the optical mode ω 0 a .(b) Effective step function spectra denoted by the dashed lines.The modulation of the optical-mode frequency ω a (t) in time is shown by the sinusoidal purple line.The shaded areas represent the interaction regions of the optical mode with the heat reservoirs.The interaction time τ is indicated in the figure and assumed equal for the two reservoirs.

Figure 3 :
Figure 3: Derivative parameters and efficiency of quantum heat engine according to the analytical model.(a) Difference of average frequencies over the interaction periods of the hot and cold reservoirs, ∆ω a = ωh a − ωc a , as a function of the peak-to-peak modulation amplitude ∆ω a of the optical-mode frequency ω a (t) = ω 0 a + ∆ω a /2 sin(ω b t).(b, c) Difference of the mean photon number of the optical mode over the cycle as a function of (b) ∆ω a and (c) the mechanical mode frequency ω b .(d) Efficiency of the device η eff = 1 − ωc a /ω h a as a function of ∆ω a .The dashed vertical lines in (a), (b), and (d) indicate the threshold where the amplitude is just high enough for the optical-mode frequency to enter the non-vanishing regions of the step-function-like power spectra of the reservoirs (left, dark red) and the threshold where the frequency leaves at its extrema the non-vanishing spectra.The parameters used here are the following: ∆ω a = 0.139 × ω 0 a , ω b = 0.05 × ω 0 a , k B T h = 0.56 × ℏω 0 a , k B T c = 0.11 × ℏω 0 a , ω h = 1.03 × ω 0 a , ω c = 0.97 × ω 0 a , Γ h = Γ c = 0.022 × ω 0 a , and l h = l c = 0.04 × ω a

Figure 4 :
Figure 4: Output power of the quantum heat engine according to the analytical model.(a-d) Output power as a function of (a) the mechanical mode angular frequency ω b , (b) the peak-to-peak modulation amplitude of the opticalmode angular frequency ∆ω a , (c) the difference between centre angular frequencies of the power spectra of the heat reservoirs ω h − ω c , and (d) the temperature of the hot reservoir T h .The dashed vertical lines indicate the cases of zero distance (left, dark red), ideal distance (centre, turquoise), and too large distance (right, yellow) of the reservoir spectra with respect to the modulation amplitude ∆ω a .We use identical parameters values to those in Fig. 3.

Figure 5 :
Figure 5: Temporal evolution of the quantum heat engine according to the quasiclassical Heisenberg-Langevin equations (2.7a) and (2.7b).(a) Mean occupation number of the mechanical mode as a function of time.The dashed line shows a linear fit to the data in order to find out the average rate of occupation generation.(b) Single realization of the real and imaginary parts of the complex-valued amplitude of the mechanical mode as functions of time.(c) Mean occupation number of the optical mode as a function of time.We show the thermal photon numbers for the optical mode, N th h (red dashed line) and N th c (blue dashed line), given by the Bose-Einstein distribution at the reservoir temperatures.The red and blue crosses indicate the maxima and minima of the optical-mode occupation in each cycle after the transient.(d) Single realization of the complex-valued amplitude of the optical mode.The parameters used are ω b= 0.048 × ω 0 a , ω h = 1.04 × ω 0 a , ω c = 0.964 × ω 0 a , k B T h = 0.56 × ℏω 0 a , k B T c = 0.11 × ℏω 0 a , g 0 = 0.012 × ω 0 a , g h ≈ 0.007, g c ≈ 0.0082, γ h ≈ 0.031 × ω 0 a , γ c ≈ 0.025 × ω 0 a .At time instant t = 0, the modes are initialized into coherent states corresponding to n a = 0.5 and n b = 39.The mean values are found by averaging over 1000 realizations of the noise.

Figure 6 :
Figure 6: Operation of the quantum heat engine with finite extracted power.(a) Temporal evolution of the average excitation number of the mechanical mode under different levels of dissipation, as indicated by the given quality factors.The dashed black lines show the fitted exponential evolution.(b) Output power (left vertical axis) and efficiency (right vertical axis) of the QHE as a function different levels of excess dissipation.The dots are the calculated values based on the simulation results extracted at the different levels of dissipation found in panel (a).The dots are connected by a solid line based on cubic spline interpolation.Details of the calculation can be found in the main text.The parameters and initial state used in the calculations are identical to those in Fig. 5.

Figure 7 :
Figure 7: Output power of the quantum heat engine according to the quasiclassical model compared with the analytical model.(a-d) Output power of the quasiclassical (cyan color) and the analytical (purple color) model as a function of (a) the mechanical mode angular frequency ω b , (b) the optomechanical coupling g 0 , (c) the difference between the centre angular frequencies of the power spectra of the heat reservoirs ω h − ω c , and (d) the temperature of the hot reservoir T h .For the quasiclassical model, we use identical parameters values to those in Fig. 6 for Q b = 1250.The data for the analytical model is identical to those in Fig. 4. The inset in (a) shows the output power in the full range of Fig. 4(a).We use cubic spline method to interpolate (dashed cyan color) the numerical data (cyan color dots) of the Heisenberg-Langevin equation model.Note that because of the differences in the models and their parameters, we do not expect the models to accurately quantitatively agree.

Figure 8 :
Figure 8: Mechanical mode occupation fluctuations of the quantum heat engine (QHE).(a) Hundred realizations of the temporal evolution of the average excitation number of the mechanical mode under dissipation characterized by Q b = 1125.The black dashed line shows the average exponential behaviour as found in Fig. 6(a).(b) Probability distribution of the mechanical mode occupation number after reaching a steady state together with a Poisson distribution for reference.The parameters and initial state used in the calculations are identical to those in Fig. 5.
R e [ α ] , R e [ β ] −6 −4 −2 0 2 4 6 8 I m [ α ] , I m [ β ] Dynamics of the amplitudes of the optical and mechanical modes.Phase space evolution of the complexvalued amplitudes of optical (purple color) and the mechanical (orange color) mode illustrated in two plus one dimension.These data is identical to those shown in Figs.5(b) and 5(d).
Using these definitions, the interaction time can be calculated from τ = (t out − t in ) − (t ′ in − t ′ out ).The average angular frequency over the interaction period ωh a in the analytical model is found by time averaging over the interaction period, b t) [Θ(t ′ out − t) + Θ(t − t ′ in )] dt,where Θ(t) is the Heaviside step function.Analogous calculations apply for the cold reservoir.