A self-contained quantum harmonic engine
Brendan Reid Simon Pigeon Mauro Antezza Gabriele De Chiara




 
 
1 Centre  for  Theoretical  Atomic,  Molecular  and  Optical  Physics, Queen's  University  Belfast,  Belfast  BT7 1NN,  United  Kingdom

2 Laboratoire Kastler Brossel, Sorbonne Université, CNRS, ENS-PSL Research University, Collège de France, 4 place Jussieu Case 74, F-75005 Paris, France

3 Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Universit de Montpellier, F- 34095 Montpellier, France

4 Institut Universitaire de France, 1 rue Descartes, F-75231 Paris, France
B. Reid 

05.30.-dQuantum statistical mechanics
05.70.-aThermodynamics
07.20.PeHeat engines




We propose a system made of three quantum harmonic oscillators as a compact quantum engine for producing mechanical work. The three oscillators play respectively the role of the hot bath, the working medium and the cold bath. The working medium performs an Otto cycle during which its frequency is changed and it is sequentially coupled to each of the two other oscillators. As the two environments are finite, the lifetime of the machine is finite and after a number of cycles it stops working and needs to be reset. 
Remarkably, we show that this machine can extract more than 90% of the available energy during 70 cycles. 
Differently from usually investigated infinite-reservoir configurations, this machine allows the protection of induced quantum correlations and we analyse the entanglement and quantum discord generated during the strokes. Interestingly, we show that high work generation is always
accompanied by large quantum correlations.
Our predictions can be useful for energy management at the nanoscale, and can be relevant for experiments with trapped ions and experiments with light in integrated optical circuits.



A self-contained quantum harmonic engine
    B. Reid1 S. Pigeon1,2 M. Antezza3,4 G. De Chiara 1,3
    December 30, 2023
========================================================






§ INTRODUCTION
The recent renewed interest in quantum aspects of out-of-equilibrium thermodynamics <cit.> has triggered an intensive research program to design and realise thermal engines whose working substance is a quantum system <cit.>. Most of the proposals employ a few qubits or quantum harmonic oscillators as the working substance but a few works explore the possibility of using quantum many-body systems <cit.>. The interest in quantum thermal engines has received a huge thrust thanks to the experimental realisation of prototypes with trapped ions <cit.> and solid state devices <cit.>. A remarkable open problem is whether quantum effects, such as coherence, entanglement or more general quantum correlations like discord, improve the engines' performance compared to their classical counterparts. 
A unified response to the question is still lacking as many of the results found in the literature are model dependent.

In all these proposals, the hot and cold environments required for the functioning of the machine are assumed to be infinite although not necessarily in equilibrium. In this work, we challenge this paradigm and consider a compact quantum engine made entirely of three quantum harmonic oscillators (see Fig. <ref>(a)). One of them, the working substance, interacts sequentially with each of the other oscillators, acting as the hot and cold reservoirs, and is subject to a compression and expansion stage, thus realising the Otto cycle. Contrary to a recent proposal employing a few qubits subject to dephasing for the hot and cold reservoir <cit.>, the time evolution of our three-oscillator engine is always unitary.



Thermalisation of quantum systems in contact to finite structured environments has been discussed in the past
<cit.>. The modelling of open quantum system dynamics in the presence of a small environment used as a calorimeter has been studied in Ref. <cit.>.

In our model, due to the finiteness and harmonicity of the environments, the working substance does not thermalise in the isochoric strokes and does not evolve in a perfect cyclic fashion. Nevertheless, we show in the following that our harmonic engine is capable of performing many cycles and producing more than 90% of the maximum extractable work. An important consequence of the finiteness of the reservoirs is that the working substance becomes  and remains quantum correlated with the reservoirs as we show using entanglement negativity and quantum discord, which we analyse using initial thermal and squeezed states.

Our engine is self-contained and can be embedded in a larger quantum system without the need of controlled external reservoirs. Thanks to the simplicity of the model, which can be solved exactly, we believe that our engine could be implemented in any experiment with controlled quantum harmonic oscillators as, for example, trapped ions or optical modes (see Fig. <ref>(b)). Detailed discussion of these experimental proposals follow.







§ MODEL
 
We consider three quantum harmonic oscillators of equal mass m and frequency ω_i subject to the local Hamiltonians H_ii=mω_i^2 x_i^2/2+p_i^2/2m. The position and momentum operators for oscillator i are given by x_i=√(ħ/2mω_i)(a_i+a^†_i) and p_i=-i√(ħ mω_i/2)(a_i-a^†_i), respectively, with a_i (a_i^†) the annihilation (creation) operators satisfying [a_i,a_j^†]=δ_ij corresponding to the vacuum state |0⟩_i. We have set ħ=1 everywhere.
The interaction Hamiltonian between oscillators i and j is 
H_ij=α_ij(a_i^† a_j+a_ia_j^†) characterised by the coupling strength α_ij. If the oscillators are resonant (ω_i=ω_j), it follows that [H_ii+H_jj,H_ij]=0 such that turning the interaction on and off does not change the energy of the system or require external work. The oscillators are initially uncorrelated and each prepared either in a thermal state with density matrix ρ_i=e^-β_i H_ii/Z_i where β_i=1/(k_BT_i) is the inverse temperature, k_B the Boltzmann constant (we set k_B=1 everywhere) and Z_i=(e^-β_i H_ii) the corresponding partition function, or in a squeezed vacuum state |r_i⟩=exp[r_i( a_i^2-a_i^† 2)/2]|0⟩_i with real squeezing parameter r_i.

The total Hamiltonian is a quadratic form of the oscillators' positions and momenta: H=∑_ij=1^3 H_ij=∑_α,β=1^6 h_αβ R_α R_β, with  R={x_1,x_2,x_3,p_1,p_2,p_3}. As a consequence, the state of the three oscillators is always Gaussian and can be fully characterised by the covariance matrix σ, with elements: σ_αβ=1/2⟨ R_α R_β+R_β R_α⟩-⟨ R_α⟩⟨ R_β⟩. The
 time evolution of the covariance matrix is governed by the Lyapunov equation:

    σ̇(t)=Aσ(t)+σ(t)A^T,
 
  where A=Ω h and Ω_αβ=-i[R_α,R_β].



§ THE CYCLE
 
We study the Otto cycle, described by two isentropic and two isochoric processes, or strokes. In our finite environment model oscillators 1, 2 and 3 play the role of the hot environment, the working medium and the cold environment respectively. The frequencies of oscillators 1 and 3 are constant in time ω_1>ω_3, while the frequency of the second oscillator is changed in time during the compression/expansion strokes between ω_1 and ω_3. Initially all couplings are zero, α_ij=0.

Using the notation E_2i,f^(s) to denote the energy of the second oscillator at either the initial or final time of the stroke s, and noticing that E_2f^(s)≡ E_2i^(s+1), the Otto cycle can be described as follows (see Fig. <ref>(a)):



  * Compression: The frequency of the second oscillator is changed from ω_3 to ω_1 in a time τ_ comp according to the schedule: ω_2^2(t) = ω_3^2 + (ω_1^2-ω_3^2) t /τ_ comp. The work done on the system during this stroke is defined as W_1=E_2i^(2)-E_2i^(1).
 

  * Heating: The first and second oscillators are suddenly coupled for a time τ_H with coupling constant α_12≠ 0, transferring energy between them. The energy transferred in this stroke is Q_1=E_2i^(2)-E_2i^(3).

  * Expansion: The frequency of the second oscillator is reduced back to ω_3: ω_2^2(t) = ω_3^2 + (ω_3^2-ω_1^2) (t-τ_  comp-τ_H)/τ_ exp in a time τ_ exp=τ_ comp. The work done is W_2=E_2i^(4)-E_2i^(3). 


  * Cooling: The second and third oscillators are suddenly coupled for a time τ_C with coupling constant α_23≠ 0, transferring energy Q_2=E^(4)_2i-E^(4)_2f.


The evolution of the covariance matrix during each stroke can be calculated analytically. For the heating and cooling strokes, when two oscillators are coupled, the evolution operators can be computed by performing a Bogoliubov transformation from the local coupled modes to the normal uncoupled modes. For the schedule we are considering the evolution operators for the compression and expansion strokes can be expressed in terms of Airy functions <cit.>. Details are provided in the Supplementary Material. 

Two extreme cases, the instantaneous change of frequency and the very slow ramp, are particularly simple and we will discuss them in detail in the next section. 
According to our definition negative work means extracted/produced and negative heat means absorbed by the working medium.

Contrary to conventional engines, at the end of the Otto cycle it is not guaranteed that the second oscillator has returned to its initial state, so
the energy balance dictated by the first law reads:
W_1+W_2=Q_1+Q_2+Δ U,
where Δ U=E_2f^(4)-E_2i^(1) is the variation of the internal energy of the second oscillator in each cycle. 
Although at each cycle the amount of work generated and heat exchanged might be different, we define a cycle efficiency as the ratio of the work  𝒲=W_1+W_2 produced in one cycle plus the variation of the internal energy Δ U and the heat absorbed from the two environmental oscillators:

    η = max(0, -𝒲+Δ U)/∑_Q_i<0|Q_i|

Thus, the machine has a positive efficiency if it produces work or if it accumulates energy in the working medium (like a battery) to be used in a subsequent cycle. Thanks to the energy balance 0≤η≤ 1 and η=1 when both Q_1 and Q_2 are negative, i.e. no heat is released. In the standard case with infinite reservoirs Q_1<0 and Q_2>0 we obtain the more common formula: η = 1-Q_2/|Q_1|.

Looking at the whole system of three quantum oscillators as a driven system initially in a non passive state (the product of local thermal passive states is not necessarily passive), we can assess the performance of the process by looking at the ratio of the work done, if any, and the ergotropy ε, which we compute numerically <cit.>. Erogotropy is described as the maximum amount of work that can be extracted from a quantum system through a unitary operation. Briefly, consider the initial state of a system ρ_i=∑_jq_j|q_j⟩⟨q_j| with Hamiltonian H=∑_jϵ_j|ϵ_j⟩⟨ϵ_j| where we assume the ordering q_j≥ q_j+1 and ϵ_j≤ϵ_j+1. With this ordering, the least energetic state that can be obtained with a unitary operation can be written as ρ_f=∑_jq_j|ϵ_j⟩⟨ϵ_j|. The ergotropy can then be defined as the work extracted in the process:
    ε=∑_ijq_jϵ_i(|⟨q_j|ϵ_i|⟩|^2-δ_ij).





§ ONE CYCLE
 
We start our analysis with the case in which the interactions between the oscillators are very small thus mimicking the weak coupling regime of open quantum systems. In this regime we can expand all our expressions up to second order in α_ijτ_H,C.
We start with the case in which the compression/expansion strokes are instantaneous. We also assume that the second and third oscillator are initially thermalised (β_2=β_3) for thermal states or have the same initial squeezing parameter (r_2=r_3) for squeezed states. Under these assumptions the total work obtained after one cycle for thermal states is given by:






    𝒲_th=τ_H^2(ω_1^2-ω_3^2)/4ω_1ω_3[ω_1(α_12^2+ω_1^2-ω_3^2)c_3-α_12^2ω_3 c_1],

where for ease of notation we have written c_i=(β_iω_i/2). The equivalent expression for initially squeezed states is given by:

    𝒲_sq=τ_H^2(ω_1^2-ω_3^2)/4ω_1ω_3×
       [ω_1(α_12^2+ω_1^2-e^4r_3ω_3^2)e^-2r_3
    -α_12^2ω_3 e^2r_1].

In order for work to be extracted, we require Eqs. (<ref>, <ref>) to be negative. In the limit of zero initial temperature (similarly zero squeezing) in oscillators two and three, the special case we will consider here, the following inequality must be satisfied:





    X>ω_1/ω_3(1+ω_1^2-ω_3^2/α_12^2),

where X=(β_1ω_1/2) for initially thermal states and X=e^2r_1 for initially squeezed states. This inequality holds for sufficiently large temperatures T_1 or for large initial squeezing r_1. Since X represents the excitation energy of the first oscillator, inequality (<ref>) is a requirement on the initial available energy of the hot oscillator. 

For these two cases we address the entanglement generated during the heating stroke between the first and second oscillator. This is the largest entanglement created during the process since the following strokes deteriorate or destroy it. For Gaussian states it is convenient to use the logarithmic negativity as a measure of entanglement <cit.>. For two oscillators i and j, this is defined as E_ij=max(0,-log|2ν_ij|) where ν_ij is the smallest symplectic eigenvalue of the partially transposed covariance matrix: σ^T_j=Pσ_ijP, where σ_ij is the restriction of the full covariance matrix to the oscillators i and j and P=(1,1,1,-1). Notice that entanglement is nonzero only if ν_ij<1/2.




When the oscillators are initially prepared in thermal states we obtain:




    ν_12=
       c_3[(1+A)ω_1ω_3c_1^2-A(ω_1^2+ω_3^2)c_1c_3-(1-A)ω_1ω_3c_3^2/2ω_1ω_3(c_1^2-c_3^2)]

with A=2(α_12τ_H)^2.
If we consider T_3=0 such that c_3→1, the condition for entanglement generation (ν_12<1/2) becomes (β_1ω_1)<(ω_1^2+ω_3^2)/(2ω_1ω_3) which happens for low enough temperatures. However, due to Eq. (<ref>), there is no work produced. Thus for initial thermal states, work extraction and entanglement are incompatible. As we will see later, however, other forms of quantum correlations might arise in the process.

















If the first oscillator is instead in an initial squeezed state the symplectic eigenvalue ν_12 takes the form:



    ν_12=1/2-√(A)e^-(r_1+r_3)|Φ|/2√(2ω_1ω_3)+Ae^-2(r_1+r_3)|Φ|^2/8ω_1ω_3,
 where Φ=(ω_1-e^2(r_1+r_3)ω_3). 


For quasi-static compression/expansion, using the asymptotic formulae presented in the Supplementary Material, we find that work is always extracted after the first cycle:




    𝒲_th=-1/2α_12^2τ_H^2(ω_1-ω_3)(c_1-c_3),


if we have c_1>c_3, when the oscillators are initially prepared in thermal states. For oscillators prepared in initially squeezed states, the expression for work extraction after the first cycle is given by:

    𝒲_sq=-α_12^2τ_H^2(ω_1-ω_3)/4(e^2r_1-e^2r_3)(1-e^-2(r_1+r_3)),
which is negative for r_1>r_3.



The symplectic eigenvalue between the first and second oscillator for the quasi-static case is given by, for thermal states, 

    ν_12=(β_3ω_3/2){1/2-α_12^2τ_H^2[sinh(β_1ω_1-β_3ω_3/2)/sinh(β_1ω_1+β_3ω_3/2)]}.
 Taking T_3=0, for which we would expect a higher value of entanglement between the first and second oscillators, the limit of this function is strictly greater than 1/2, so no entanglement is generated.





For initially squeezed states the expression for the lowest symplectic eigenvalue is too lengthy to report here. However, preparing the second and third oscillator with zero initial squeezing, entanglement generation is guaranteed between the first and second oscillator: ν_12 = 1/2 -α_12τ_H sinhr_1



for non-zero r_1. To summarise, to get entanglement created in the initial heating stroke it is sufficient, but not necessary, to prepare the first oscillator in a squeezed state.

In all the cases discussed under the weak coupling approximation, we find that Q_2=0 and thus, if 𝒲<0 , η=1. However it is important to stress that since these analytical calculations refer to a single cycle, the work extracted to ergotropy ratio is very small. In the following we discuss our numerical results in a more general scenario. 





§ MANY CYCLES
 
To fully characterise the performance of the engine over many cycles we use numerical optimisation to enhance work extraction in our machine.
We set β_1=10^-2ω_1^-1, β_2,3→∞ and ω_1=1.  For values of 0<ω_3<ω_1, the variables to be optimised are the coupling values α_ij and corresponding coupling times τ_H,C as well as the compression/expansion time τ_ comp. To this end, we employ the built-in Mathematica function NMinimize <cit.>. To avoid large interaction times and large coupling values, which would result in a short-lived engine, each parameter was restricted to a finite interval: 10^-4≤α_ijω_1^-1≤0.05, 10^-3≤ω_1τ_H,C≤1, 1≤ω_1τ_R≤100. 



For ω_3=0.1ω_1, the result of the optimisation is shown in Fig. <ref>(a) where the total work extracted to ergotropy ratio is shown. The machine extracts work until the 70th cycle and then reverses itself absorbing work and dumping heat into the first oscillator, returning almost completely to its initial state. At the end of the 70th cycle the total work extracted from the engine is 𝒲_T≃ -89.5ω_1,  approximately 91% of the ergotropy of the initial state. It should be noted that the parameters reported in the caption of Fig. <ref> are not unique and other sets of parameters can provide a similar amount of work extraction.

Due to the small value of α_23 we find that the heat exchanged to the third oscillator is negligible: Q_2∼ 0, thus the thermal efficiency of the machine is practically 1. For an optimised engine such as the one presented here, Q_2 will always be small as most of the energy will be extracted during the expansion stroke of the cycle. This feature is specific to an optimised engine: there exist other sets of parameters, leading to a functioning engine, albeit with a smaller value of extracted work, that have larger values of α_23 giving Q_2 small but non-zero. It is interesting that we obtain the same result for the total work extracted, within our working precision, by eliminating the third oscillator altogether. This is a consequence of the non-cyclic nature of the engine.

The choice of preparing the first and second oscillators in their vacuum states is not accidental, we found this case allows for the highest work to ergotropy ratio. However, non-zero initial values of T_2,3 also allows for the extraction of work. In fact, we found that for T_2=T_3>0, and disallowing instances where the second oscillator can “ignore” the third with small coupling values and/or times (such as the case presented above) optimal work extraction decreases linearly with T_3. Setting T_3=0 and T_2>0 we found that maximum work extraction is identical to that presented above: 𝒲_T≃-89.5ω_1. However, the work to ergotropy ratio is lower.

 As shown in the previous section entanglement between initially thermal states is not present in a regime where work can be extracted, and so we evaluate the non-classicality of this process via the Gaussian quantum discord <cit.>. We provide a brief, mathematical definition of Gaussian quantum discord in the Supplementary Material. Figure  <ref>(b) shows the discord dynamics between the first and second oscillator as a function of time. It is symmetric about the point of reversal for the engine and is the only non-classical correlation present in the engine with 𝒟_23, 𝒟_13=0.  Detailed information on the internal energies of the three oscillators performing many cycles can be found in the Supplementary Material.
 
 







Figure <ref> shows optimised work extraction for changing values of ω_3. The maximum amount of extractable work, as a proportion of the ergotropy of the initial state, decreases linearly with increasing ω_3, until it disappears for ω_3∼ω_1. While work extraction is possible for any value of ω_3<ω_1, clearly our machine operates better in the regime ω_3≪ω_1. 











§ QUANTUM CORRELATIONS
 

The optimised case presented above provides an insight on the performance of the engine for maximum work production. While this is arguably the most important aspect of any machine, it is unclear on what connection, if any, exists between work extraction and non-classical correlations in the system. In order to investigate any link between thermodynamic quantities (work, efficiency) and non-classical correlations (entanglement, discord) we simulate the engine working with random parameters. 

We consider two cases: initially thermal states and initially squeezed states. We fix the frequency of the first oscillator (ω_1=1) and the initial temperatures (squeezing parameters) for thermal (squeezed) states. The variables to be randomized are ω_3, the coupling strengths α_ij and times τ_H,C and τ_ comp. The engine is prepared with these random variables and performs Otto cycles until the criteria 𝒲<0 is no longer satisfied. This process is repeated ten thousand times. For thermal states, we have chosen variables β_1=10^-2ω_1^-1 and β_2,3→∞. For squeezed states, we fix r_2,3=0 and r_1=1/2arccosh[(β_1ω_1/2)] so the initial energy of the first oscillator is equal in both cases.

For initial thermal states, there is no significant generation of entanglement and we assess the amount of quantum correlations between any two oscillators using quantum discord. For initial squeezed states, instead we show entanglement measured by the logarithmic negativity defined above.


For thermal states, Fig. <ref>(a-b-c) show the scatter plots of the maximum discord generated between any two oscillators at the end of a cycle and the total work produced with a given set of parameters. The discord 𝒟_12 between the first two oscillators show a clear pattern: high values of work are necessarily accompanied by large quantum discord. In order to produce work, there must be a coherent energy transfer between the first and second oscillator that can only happen if the two oscillators become correlated. This is also reminiscent of synchronisation phenomena in coupled harmonic oscillators (see for example <cit.>). In contrast, the quantum discords of the other pairs of oscillators do not show such a connection, and large values of work are not necessarily followed by large values of discord. Notice also that the discord 𝒟_13 created between the first and third oscillator is mediated through the central oscillator. 



For initial squeezed states, the results for the maximum negativity created are shown in Fig. <ref>(d-e-f). As in the case of initial thermal states, a strong entanglement between the first two oscillators is necessary for extracting work. There also seems to be a maximum value of entanglement that can be created  connected to the initial amount of squeezing of the first oscillator. 
Entanglement between the second and third oscillator and between the first and third oscillator tends to be low and correlated to small values of total work. This suggests a minimal involvement of the third oscillator in the work extraction.





§ IMPLEMENTATION
 Our self-contained engine can be implemented in any experiment with controlled quantum harmonic oscillators. Two platforms fit especially well our proposal: trapped ions and integrated optical circuits. With trapped ions, the role of each oscillator is played by the ion displacement from equilibrium. The trapping frequency and the ions distance, determining their interactions, can be accurately controlled <cit.>. Thanks to the advance in state engineering,  it is possible to prepare the ions in thermal and squeezed states <cit.>. 


Integrated optical circuits are another interesting candidate. The time evolution of a quantum oscillator can be mapped to the propagation of an electromagnetic field in a single mode nanofibre. Preparation of thermal and squeezed states has been well mastered in quantum optics <cit.>. The circuit proposal in Fig. <ref>(b) allows optical states to undergo frequency shifting and coupling via evanescent fields which realise H_ij <cit.>. The main difficulty with this system is to shift the frequency of the oscillators. For frequency shifting based on acoustic-optic effect only a modulation of a few percent can be reached with high efficiency <cit.>. Frequency conversion is an active field of research and new methods based on doped fibres have been recently developed, allowing for a larger frequency modulation but still with smaller efficiency <cit.>.   




§ CONCLUSIONS
  In this work we have designed the simplest quantum engine made of three quantum harmonic oscillators and performing the Otto cycle. This machine does not require the use of external infinite thermal reservoirs but it operates thanks to an initial non-passive state. As such it can be integrated in a larger quantum machine that needs mechanical work in a very localised region of space. Since its resources are finite so is its lifetime. However we have shown that the optimised engine works for about 70 cycles before stopping, extracting approximately 91% of the available energy. Our investigation confirms the common belief, although not conclusively proven, that quantum correlations enhance work extraction as evidenced numerically in our Fig. <ref>.

Our research opens a new avenue in the study of quantum thermodynamic devices with finite reservoirs. Extending the lifetime of the engine could be achieved by coupling more oscillators onto the first and third oscillators, thus increasing the size of the energy reservoirs <cit.>. It has also been shown that systems coupled to a finite reservoir can mimic the action of a system connected to an infinite reservoir <cit.>. Finally, it would be interesting to extend a program such as the one we presented here to spin systems and many-body lattice systems.



We acknowledge illuminating discussions with O. Abah, A. Ferraro, A. Hewgill, V. Parigi, J. P. Paz and A. J. Roncaglia. GDC thanks the CNRS and the group Theory
of Light-Matter and Quantum Phenomena of the Laboratoire
Charles Coulomb for hospitality during his stay in Montpellier. BR would like to thank S. Campbell for illuminating discussions.
This work was supported by the EU Collaborative project TherMiQ (grant agreement 618074), the John Templeton Foundation (grant ID 43467) the French ANR (ACHN C- Flight), the UK EPSRC and the Professor Caldwell Travel Scholarship. 







eplbib





Supplementary Material





§ EVOLUTION OPERATORS FOR THE VARIOUS STROKES

In this section we find analytically the evolution operators needed to evolve the covariance matrix in time. To this end we need to solve the Lyapunov equation σ̇(t)=Aσ(t)+σ(t)A^T. For a time independent Hamiltonian, i.e. the matrix A is constant, the solution is simply:

    σ(t) = exp(At) σ(0) exp(A^T t)

This is the case of the heating and cooling strokes when two resonant oscillators are coupled. 
When the first and second oscillators, both at frequency ω_1, are coupled α_12≠ 0 and α_23=0. The matrix A reads:

    A=([        0        0        0        1 α_12/ω_1        0;        0        0        0 α_12/ω_1        1        0;        0        0        0        0        0        1;   -ω_1^2 -α_12ω_1        0        0        0        0; -α_12ω_1   -ω_1^2        0        0        0        0;        0        0   -ω_3^2        0        0        0 ])

and we obtain:

    exp(At)= 
       ([        cosα_12tcosω_1t       -sinα_12tsinω_1t                      0 1/ω_1 cosα_12tsinω_1 t 1/ω_1 sinα_12tcosω_1 t                      0;       -sinα_12tsinω_1t        cosα_12tcosω_1t                      0 1/ω_1 sinα_12tcosω_1 t 1/ω_1 cosα_12tsinω_1 t                      0;                      0                      0                cosω_3t                      0                      0         1/ω_3 sinω_3 t;    -ω_1cosα_12tsinω_1t    -ω_1sinα_12tcosω_1t                      0        cosα_12tcosω_1t      -sinα_12tsinω_1 t                      0;    -ω_1sinα_12tcosω_1t    -ω_1cosα_12tsinω_1t                      0      -sinα_12tsinω_1 t        cosα_12tcosω_1t                      0;                      0                      0            -ω_3sinω_3t                      0                      0                cosω_3t ])

A similar expression for A and exp A t is obtained as well when α_23≠ 0, α_12=0 and ω_2=ω_3. 

Now let us consider the compression/expansion strokes. In this case the oscillators are uncoupled and the second oscillator is subject to a time-dependent frequency. The reordered Hamiltonian matrix is:

    A(t) =([     0_3     1_3; -ω^2(t)     0_3 ])

where 0_3 and 1_3 are the 3× 3 zero and identity matrices, respectively, while ω^2(t)=(ω_1^2,ω_2^2(t),ω_3^2).
The solution to the Lyapunov equation (1) is thus σ(t)=Uσ(0)U^T, where:

    U = ([ a b; c d ])

and 

    a   =   (cosω_1t, y(t), cosω_3 t)
    
    
    b   =   (1/ω_1sinω_1t, x(t), 1/ω_3sinω_3 t)
    
    
    c   =   ȧ
        d   =   ḃ

The functions x(t) and y(t) fulfil the equations:

    ẍ(t)   =   -ω_2^2 x(t),    x(0)=0,  ẋ(0)=1
    
    ÿ(t)   =   -ω_2^2 y(t),    y(0)=1,  ẏ(0)=0

These equations can be solved analytically for the schedule 
ω_2^2(t) = ω_ in^2 + (ω_ fin^2-ω_ in^2) t /τ. Their solutions are given in terms of Airy functions:

    x(t)   =    -π(τ/ω_ in^2-ω_ fin^2 )^1/3{[z(t)](w)-(w)[z(t)]}
    
    y(t)   =    -π{[z(t)]'(w)-[z(t)]'(w)}

where 

    w   =   -ω_ in^2 [τ/ω_ in^2-ω_ fin^2]^2/3
    
    z(t)   =   -ω_2^2(t)[τ/ω_ in^2-ω_ fin^2]^2/3


In our calculations we specialise to two extreme cases, one in the limit τ→ 0 in which case we simply have U=1_6 and the other for very slow ramp: τ≫max(ω_ in^-1,ω_ fin^-1). In this case we can use the asymptotic expansion of the Airy functions for large arguments. Under this approximation the evolution operator takes the form:

    U=([             cosω_1τ                   0                   0        1/ω_1sinω_1τ                   0                   0;                   0 √(ω_ in/ω_ fin)cosϕ                   0                   0 sinϕ/√(ω_ inω_ fin)                   0;                   0                   0             cosω_3τ                   0                   0        1/ω_3sinω_3τ;         -ω_1sinω_1τ                   0                   0             cosω_1τ                   0                   0;                   0 -√(ω_ inω_ fin)sinϕ                   0                   0 √(ω_ fin/ω_ in)cosϕ                   0;                   0                   0         -ω_3sinω_3τ                   0                   0             cosω_3τ ])

where we have set

    ϕ=2/3τω_ fin^2+ω_ inω_ fin+ω_ fin^2/ω_ fin+ω_ in.

Equation (<ref>) can be interpreted as the evolution operator for the free evolution of oscillators 1 and 3 with their respective frequencies for a time τ times a squeezing operator dependent on the ratio ω_ in/ω_ fin combined with a phase space rotation with angle ϕ.



§ DEFINITION OF GAUSSIAN DISCORD


In this section we will give a brief mathematical definition of Gaussian discord <cit.>. For convenience, we use a different ordered basis from the main text, R̃=(x_1, p_1, x_2, p_2), and recall that the elements of a covariance matrix of two oscillators are defined as σ_αβ=1/2⟨R̃_αR̃_β +R̃_αR̃_β⟩ - ⟨R̃_α⟩⟨R̃_β⟩. For a bipartite covariance matrix in this basis, 




    σ=([   A   C; C^T   B ]),

where A, B and C are 2×2 matrices, the matrix A (B) relates only to the mean values in subsystem 1 (2) and the matrix C relates to correlations in these mean values. The Gaussian quantum discord of this state can then written as 

    𝒟=h(√(I_2))- h(d_-)-h(d_+)+h(√(E_ min))

where we have 


    h(x)=   (x+1/2)log(x+1/2)-(x-1/2)log(x-1/2),
    
    d_±^2=   1/2(λ±√(λ^2-4I_4)),
    λ=   I_1+I_2+2I_3

for I_1=(A), I_2=(B), I_3=(C) and I_4=(σ). The quantity E_ min is given by the following formula:




    E_ min=
    (|I_3|+√(I_3^2-(I_1-4I_4)(I_2-1/4))/2(I_2-1/4))^2    if(I_1 I_2-I_4)^2≤(I_1+4I_4)(I_2+1/4)I_3^2
    
    1/I_2(I_1I_2-I_3^2+I_4-√((I_1I_2+I_4-I_3^2)^2-4I_1I_2I_4))    otherwise.








§ INTERNAL ENERGY AND QUANTUM CORRELATIONS FOR THE OPTIMISED MANY-CYCLE ENGINE

In this section, we present the internal energy dynamics of each oscillator for the optimised many cycle case presented in the main text. We also analyse the discord dynamics between each of the oscillators. The parameters used are ω_3=0.1ω_1, τ_ comp=85ω_1^-1, τ_H=0.59ω_1^-1, τ_C=0.9996ω_1^-1,α_12=0.038ω_1,α_23=10^-4ω_1. Figure <ref> shows the internal energy dynamics E_i=⟨ H_i ⟩ for each oscillator plotted against time in units of the Otto cycle time τ=τ_H+τ_C+2τ_ comp. The system performs 70 Otto cycles satisfying 𝒲<0 and the plots also show the subsequent Otto cycles where the engine reverses. At the end of the 140th Otto cycle each oscillator has almost returned to its initial state. As the coupling strength between the second and third oscillator is extremely weak, the internal energy of the third oscillator does not change. The efficiency of this process is unity, as Q_2=0.




The number of cycles shown in Fig. <ref> makes the fine details of the dynamics difficult to parse. Figure <ref> shows a closer view of the internal energies of the first and second oscillators from the 40th to the 44th stroke. The internal energy of the first oscillator is shown in Fig. <ref>(a). As it is only interacting with the second oscillator during the heating stroke its energy is constant for most of the cycle. Due to the disparity of the time scales of the compression, expansion stages and the heating, cooling stages the latter appear instantaneous on this scale. In Fig. <ref>(b) the internal energy of the second oscillator is shown. The large slopes of each arch correspond to the energy increase and decrease caused by the compression and expansion stages. The small peak at the top of each arch corresponds to the heating stroke, where the height of each peak corresponds directly to the height of each energy reduction in panel <ref>(a). 












