Three-dimensional harmonic oscillator as a quantum Otto engine

A quantum Otto engine based on a three-dimensional harmonic oscillator is proposed. One of the modes of this oscillator functions as the working fluid, while the other two play the role of baths. The coupling between the working fluid and the baths is controlled using an external central potential. All four strokes of the engine are simulated numerically, exploring the nonadiabatic effects in the compression and expansion phases, as well as the energy transfer during the working fluid's contact with the baths. The efficiency and power of several realizations of the proposed engine are also computed with the former agreeing well with the theoretical predictions for the quantum Otto cycle.


I. INTRODUCTION
The purpose of a quantum engine is the same as that of a classical one: converting energy into work.The main distinction is that quantum components play the role of the working fluid.[1,2] Examples of quantum systems that can act as the fluid include two-level systems [3,4], one [5][6][7][8][9] or multiple [10] harmonic oscillators, or photons [11].The second difference between classical and quantum varieties is their energy sources.While classical versions get their energy from hot reservoirs, quantum engines can be powered by measurements, as described in several recent studies.[12][13][14][15][16][17][18][19][20][21] Despite the variety of quantum engines, perhaps the most common one is based on the classical Otto cycle, where the working fluid is a gas.This cycle involves four processes: adiabatic compression, constant-volume heat transfer to the gas, adiabatic expansion, and constantvolume heat transfer away from the gas.The two heat exchange processes indicate that the engine requires two heat reservoirs for its operation.In the quantum variety of the Otto cycle, the working fluid is often chosen to be a harmonic oscillator [5][6][7][8][9][10] because the compression and expansion can be achieved by adjusting the oscillator's force constant.
This work proposes an implementation of a quantum Otto engine, consisting of a three-dimensional harmonic oscillator with anisotropic force constants in the three orthogonal directions.The three corresponding modes act as a hot reservoir, a cold reservoir, and the "gas" with an adjustable force constant.The coupling between the gas and the baths can be controlled using an external central potential, as shown in Ref. [22].Although single "bath" modes do not function as true thermodynamic reservoirs, they can be regarded as ancillae that deliver energy to or extract it from the gas mode.While a particular bath mode is disconnected from the gas, it can be either cooled or heated to return it to its pre-contact state.The cooling and heating mechanisms for the ancillae may involve optics (such as laser cooling) or coupling the mode to a thermodynamic reservoir at a particular temperature.
For the sake of brevity, the two ancillary modes will be referred to as "baths" in the text while keeping in mind that there are additional external reservoirs.
A reasonable question is why the auxiliary bath modes are necessary instead of directly cooling and heating the gas mode using the same means.The reason is complexity: cooling and heating of a single mode means one must be able to swap the heating and the cooling mechanisms.If there are dedicated ancillary modes, however, each one of them needs to be connected to, at most, a single external system.Moreover, if the "resetting" process of the bath modes is sufficiently slower than their rate of energy exchange with the gas mode, the external system can remain connected to the bath mode throughout the engine's operation.The idea of always keeping the bath connection on for different implementations of the quantum Otto cycle was explored earlier in Refs.[23,24].
The efficiency of an ideal quantum Otto cycle is obtained by assuming that the working fluid is compressed and expanded adiabatically, and that it reaches the bath temperature.[1] For an engine to deliver a nonvanishing power, they must be able to complete their cycles in a finite amount of time.Therefore, the compression and expansion phases are likely to give rise to nonadiabatic effects.Additionally, for the engine proposed here, one does not expect the gas mode to reach the temperature of the bath mode.In fact, the state of the gas mode is likely not to be thermal generally.Nevertheless, this work demonstrates that, despite these deviations from the ideal cycle, the resultant efficiency can be close to the predicted value.
To show that the engine proposed does not need to adhere to the ideal configuration to deliver such efficiency, it is useful to first focus on individual strokes of the Otto cycle.Thus, after presenting the theoretical model describing the engine in Sec.II, this work dedicates Secs.III and IV to the compression/expansion phases and energy exchange with the baths, respectively.The entire cycle of the engine operation is presented in Sec.V. Summary and conclusions are found in Sec.VI.
All computations are performed using the julia programming language.[25] The plots are made using Makie.jlpackage [26] using the color scheme designed for colorblind readers.[27] The scripts used for computing and plotting can be found at https://github.com/rodinphysics/quantum-oscillator-engine.

II. MODEL
As discussed in the introduction, the engine consists of a three-dimensional harmonic oscillator, where two of the dimensions function as baths, and the remaining one operates as the "gas" in the Otto cycle.The Hamiltonian for such a system is given by where the c, h, and g correspond to cold, hot, and gas, respectively.The last term of the first line allows the gas to be compressed and expanded, as required by the engine operation.The second line contains the interaction terms between the gas and the baths.
The most natural way to describe the system is using the Fock basis |j⟩ c ⊗ |k⟩ h ⊗ |l⟩ g ≡ |j, k, l⟩, where j, k, and l are the energy levels of the three modes.Therefore, one may be tempted to write the portions of Eq. ( 1) corresponding to the baths as ℏΩ b ( b † b + 1/2).Although entirely valid, this choice makes the subsequent computations messier.Consequently, it is easier first to rewrite Eq. (1) as Here, it was assumed that the gas force constant takes a range of values during the compression and expansion phases.Writing the first line of Eq. ( 2) using the ladder operators gives the following Hamiltonian where Ω = k/m and α(t) = κ(t)/mΩ 2 = κ(t)/k.Because the frequency of the hot mode is given by Ω h = (k + κ)/m = Ω √ 1 + α max , one gets α max = ω 2 − 1 with ω = Ω h /Ω giving the ratio of compressed and uncompressed frequencies.It is also convenient to express all the energies in terms of ℏΩ, time as t = 2πτ /Ω, and lengths in terms of the corresponding quantum oscillator length to get where Φ = U/ℏΩ.It is worth noting that only some terms must be included from Eq. ( 4) during the engine operation.For example, during the compression and expansion phases, the interaction terms are switched off as α(τ ) changes in time.On the other hand, only one of the interaction terms is nonzero during the bath coupling phases, and the uncoupled bath can be ignored as it evolves independently.

III. COMPRESSION AND EXPANSION
When the gas undergoes compression or expansion, it is disconnected from the baths and the relevant normalordered portion of Eq. ( 4) is The process is guided by the time evolution operator, given by the solution to the time-dependent Schrödinger equation where the factor of 2π arises from the definition of τ .In this work, Eq. ( 6) is solved using the fifth order Runge-Kutta method with the numerical benchmarks provided in Appendix A.
Although it is common to treat the compression/expansion phases of the quantum Otto cycle as adiabatic, the finite duration of the strokes is likely to give rise to nonadiabatic effects.To quantify the amount of nonadiabaticity, it is convenient to proceed as follows.Let a system, described by a Hamiltonian Ĥi start in some state ρi .If the system is compressed or expanded adiabatically, the final energy will be E A f = tr Ĥf ρi , where Ĥf is the final form of the Hamiltonian.For a non-adiabatic modification, E f = tr Ĥf Û ρi Û † .The level of nonadiabaticity is obtained by taking the ratio 1 illustrates the effects of nonadiabatic tuning of thermal states at two temperatures which, as expected, are greatest for higher rates of change, corresponding to small τ α , and large ω.A system initialized at a 41level thermal state at temperature ωT = kBT /ℏΩ with α = 0 (top row) or α = ω 2 − 1 (bottom row).Next, the state is either compressed (until it arrives at α = ω 2 − 1) or expanded (until α = 0) uniformly over time τα.The final energy of the system is then divided by what would be the energy of the system if the modification were adiabatic.The nonadiabatic effects are most prevalent when the rate of the confinement change is high (corresponding to large ω and small τα).Hightemperature compression demonstrates the effects of the finite Fock space, as discussed in the main text.
In addition to illustrating the consequences of dα(τ )/dτ ̸ = 0, Fig. 1 highlights another aspect of this stroke that needs to be kept in mind when performing numerical simulations.One can see that the difference between the finite-time and adiabatic results for the hightemperature compression phase, shown in Fig. 1(c), is greater than for all other cases.This is a spurious effect arising from the finite size of the Fock space used here.For any modification of the oscillator frequency, there is a finite probability of the oscillator transitioning to higher energy levels, especially during a compression.At low temperatures, only the lowest-energy states are occupied at τ = 0 so fewer Fock states are needed to correctly represent the compression.As the temperature is increased, more levels become occupied, requiring a larger Fock space to allow the required transitions.All panels in Fig. 1 use 41 states and, for ω T = 5, this is not sufficient.Halving the number of states aggravates the problem specifically for panel (c) while keeping the other panels relatively intact.When demonstrating the engine operation in Sec.V, the hot bath's temperature will be set to ω T = 5, while the cold one will be kept at ω T = 1/10.The fact that panels (a) and (d) do not change when reducing the number states indicates that this basis size is sufficient to avoid the finite-basis effects for these temperatures with ω ≤ 3 and τ α ≥ 1.

IV. BATH COUPLING
The purpose of the bath modes is to either add energy to or remove it from the gas mode.The state of the baths does not have to take a particular form as long as the energy flow occurs in the right direction.It is, however, convenient to consider thermal baths as illustrative examples so that their unnormalized density operators are ρb = e − Ĥb /ω T .
Coupling the gas mode to a thermal bath will not, generally, set the gas to a thermal state.However, repeating the process multiple times by introducing new identical thermal baths is expected to bring the gas to the bath temperature eventually.One can use this physical intuition to check that the energy exchange between the modes functions appropriately before focusing on the engine itself.
When the gas couples to a bath, α remains fixed, leading to Assuming that the coupling between the gas and the baths is switched on and off quickly, it is reasonable to suppress the time argument inside the interaction term so that the matrix elements in the Fock space become where Ψ n (x) = ⟨x|n⟩ are the harmonic oscillator wavefunctions.
As discussed in the introduction, the coupling between the modes is controlled using a central potential that can be switched on and off.Equation (8) shows that if the coupling term is even in x b and x g (as is the case for a central potential aligned with the origin), only the states with the same parity in each of the modes can couple.One can introduce coupling between more modes by positioning the extremum of the central potential away from the symmetric (0, 0) point so that the interaction term takes a general form Φ(r − r 0 ).For illustration, it is convenient to use a Gaussian interaction Φ 0 exp −|r − r 0 | 2 /2σ 2 = Φ 0 e −(xg−xg,0) 2 /2σ 2 e −(x b −x b,0 ) 2 /2σ 2 The advantages of this interaction are twofold.First, its amplitude and extent are easily tunable.Second, because the term is separable in x g and x b , the integrals in Eq. ( 8) can be easily computed numerically.
If the extremum of the central potential is positioned diagonally from the origin (x g,0 = x b,0 = x 0 ), the procedure of computing the interaction matrix becomes particularly simple.First, one obtains the matrix Φ single with elements ⟨j|e −(x−x0) 2 /2σ 2 |k⟩ for all the Fock states in the single-oscillator basis.From this, the full interaction matrix becomes Φ = Φ 0 Φ single ⊗ Φ single .
A set of four simulations is performed to demonstrate the energy flow between oscillator modes.To this end, four thermal oscillator states using 41 energy levels are generated defined by (α, ω T ) with α ∈ {0, 8} and ω T ∈ {ω cold T = 1, ω hot T = 5}.The interaction width σ and x 0 are set to 1.The coupling strength will be allowed to vary, as described below.For a given α, one of the states is designated as the gas, while the other acts as the bath with the gas (bath) state denoted by ρg (ρ b ).
The key idea is that by having the gas interact with a series of baths, the state of the gas should approach the bath state.The trace distance is used to measure how close the gas is to ρb .
For this demonstration, the gas will interact with twelve baths.For the first five of them, the coupling strength Φ 0 = 1, for the next four Φ 0 = 1/5, and for the final three Φ 0 = 1/20.The reason for ramping down Φ 0 is the reduction of energy stored in the interaction between the two oscillators: as the gas state approaches ρb , most of the energy needs to be stored in the oscillator energy, not in the coupling term between the gas and the bath.
At τ = 0, as the gas is brought into contact with the first bath, the trace distance is calculated.At this moment, the full state of the system is ρ0 total = ρg ⊗ ρb .Because the Hamiltonian does not vary in time, Û(δτ ) = e −2πi Ĥδτ can be calculated exactly for Φ 0 = 1 and δτ = 5 so that the state of the composite system at later times is given by ρtotal (n × δτ ) = Ûn ρ0 total ( Û † ) n .For 1 ≤ n ≤ 10, the partial trace of the gas with respect to the bath tr b [ρ total (n × δτ )] is computed and its trace distance to ρb is calculated to show how the reduced density operator evolves while the gas is coupled to the bath.At n = 10, the bath is traced out and the gas becomes coupled to a new bath in state ρb so that the total state is ρtotal = tr b [ Û10 ρg ⊗ ρb ( Û † ) 10 ] ⊗ ρb .The process is repeated for each of the twelve bath using the Û for the appropriate Φ 0 .
The computed trace distance as a function of time is plotted in Fig. 2(a).One can see that the trace distance is less than 10 −2 at the end of the simulation for all four configurations.The final distance is smaller for the hot bath because, for a given Φ 0 , the interaction energy is proportionally smaller compared to the average oscillator energy at large temperatures.One can reduce the distance for small ω T by further lowering Φ 0 .
Figure 2(a) shows that the gas state approaches the bath.However, it is not evident that the gas is not also approaching some temperature different the bath's ω T .One can verify that the gas is indeed not drifting towards some wrong temperature by computing the trace distance between the final gas state and thermal states for a range of ω T with the corresponding α.The results, shown in Fig. 2(b) confirm that the distance is smallest when ω T equals the bath temperature.Thus, even though the gas state is not quite thermal [because of the nonzero trace distance in panel (a)], the thermal state that it is closest to has the correct temperature.
It is important to note that the thermal equilibration protocol is somewhat ad hoc with the aim to demonstrate the general process, not make the energy transfer fast or efficient.The contact times, σ, x 0 , as well as the values of Φ 0 were not optimized.A more carefully designed procedure can speed up the process and bring the gas closer to ρb .The most important result of this section is that the coupling model used here results in the correct energy transfer between oscillator modes, validating its use in the engine description.

V. ENGINE OPERATION
Having demonstrated the compression and expansion strokes of the engine, as well as the correct energy exchange with the baths, it is now possible to turn to the engine operation.The main goal for this section to is confirm that the engine states are cyclical and that work can be extracted from the gas mode.
Without the loss of generality, it is convenient to have the cycle begin with the gas compression.Thus, if the state of the gas in the beginning of the nth cycle is given by ρn , the state of the gas in the beginning of the following cycle is Equation ( 9) should be read from the inside outward to follow the transformations that the gas undergoes.First, sandwiching ρn between Û and its conjugate compresses the gas.Next, the compressed gas is coupled to the hot bath in thermal state ρh , as shown by the tensor product.
After that, the composite system is allowed to evolve in time by applying operators Ĥ and Ĥ † , followed by a decoupling represented by the partial trace operator tr b .
Then, the gas is expanded, as can be seen by the reversed application of Û † and Û, and connected to a cold bath ρc .Following the evolution guided by Ĉ, the gas is finally separated from the bath by the partial trace operator, completing the cycle.For a particular engine configuration, Ĉ, Ĥ, Û, and ρc/h need to be computed only once.The first two operators are calculated by multiplying the total Hamiltonian by −2πiτ contact and exponentiating the result, where τ contact is the gas-bath interaction time.Û is computed using the same Runge-Kutta scheme as was employed in Sec.III for a linear compression.
Four configurations are chosen for demonstration with two different ω's (2 and 3) and two different sets of times used for expansion/compression and bath contact.For the first pair, the duration of all strokes of the cycle is equal to 1.For the second pair, the expansion/compression time is set to 4, while the bath contact time is 10.Hence, there are two "fast" and two "slow" cycles.The bath temperatures are the same for all realizations with ω cold T = 1/10 and ω hot T = 5 and the basis consists of 41 states.
For each engine, the gas is initialized in the thermal state at ω cold T .Next, it is taken 50 times through the cycle described by Eq. ( 9).At the end of each stroke, the energy of the gas is calculated by taking the trace of the product of the gas density operator and the Hamiltonian in Eq. ( 5) with α = 0 or α = ω 2 − 1, as appropriate.
The computed energies for the four realizations are given in Fig. 3(a)-(d).The x-coordinate labels the cycle and, for each cycle, the order of the points is "expanded cold" → "compressed cold" → "compressed hot" → "expanded hot," after which one moves to "expanded cold" of the next cycle.
First, one can observe that, for the slow engine, the energies stabilize within a few cycles.The fast engine, on the other hand, requires a substantially longer time.Aside from this effect, the largest difference between For the fast engine (top row), the time of each stroke is 1.For the slow engine (bottom row), the time of expansion/compression is 4, while the bath contact time is 10.
The system is evolved following Eq.( 9) with ω hot T = 5 and ω cold T = 1/10.The interaction between the baths and the gas is the same as in Fig. 2 with Φ0 = 1.For the compression and expansion phases, Û is calculated for a 41-state basis.(e) Efficiency for each of the engines as a function of the cycle number, computed by dividing the total work output by heat input.(f) The power of the engines, obtained by dividing the total work output by the duration of a cycle.the two speeds is the amount of energy transferred to and from the baths, as can be seen from the difference between "compressed cold" → "compressed hot" and "expanded hot" → "expanded cold" transitions.Naturally, the ω = 3 configuration demonstrates larger energy changes during the compression and expansion phases, as expected.
To quantify the efficiency of the engine, one first calculates the work output by the engine, given by − [(E exp hot − E comp hot ) + (E comp cold − E exp cold )].Dividing this value by the heat delivered by the hot bath E comp hot − E comp cold yields the efficiency, plotted in Fig. 3(e).The plot demonstrates that, for the slow engine, the efficiency is ≈ 1/2 for ω = 2 and ≈ 2/3 for ω = 3.These values agree well with the expected efficiency of an Otto cycle that operates in a fully adiabatic regime with thermal baths, where the efficiency is given by 1 − ω −1 .[1] Curiously, the efficiency of the fast engine is slightly higher, but still close to these values.It is worth noting that both engines operate close to the adiabatic regime, as can be seen from Fig. 1 for the values of τ α employed.For the heat exchange, on the other hand, even the slow engine is not expected to reach the temperature of the bath, as confirmed by Fig. 2.
As is clear from panels (a)-(d), the fast engine yields much less work per cycle.However, given that its cycle is seven times slower, it is more appropriate to compare the power of the two setups, as is done in Fig. 3(f) by dividing the work output by the cycle period.This figure shows that despite a higher efficiency, the fast engine delivers less power by about a factor of two.
The most important message of this section is that the engine reaches a stable cycle and its efficiency agrees well with the predicted value both in a slow and fast operation regimes.Additionally, for the realizations here, the factor that reduces the work output of the engine is not the nonadiabatic effects associated with fast compressing/expansion, but a shorter contact with the bath, limiting the amount of energy transferred.

VI. SUMMARY
This work has introduced and simulated a realization of a quantum Otto engine comprising of a single threedimensional harmonic oscillator.One of the modes of the oscillator functions as a compressible working fluid, while the others act as hot and cold reservoirs.The coupling between the baths and the working fluid is controlled by a nonlinear external potential.Individual finite-time strokes of the engine were simulated numerically to explore the role of adiabaticity during the compression and expansion phases, as well as the energy flow during thermal contact with the baths.It has been shown that even for a ninefold increase of the working fluid's force constant, performing the compression and expansion over a few oscillator periods essentially eliminates the nonadiabatic effects.Additionally, it has been confirmed that having the working fluid interact with a series of baths eventually brings the working fluid to the bath temperature, as expected.Finally, the study demonstrates that the working fluid reaches a stable state as it goes through multiple engine cycles.The efficiency obtained here agrees well with the theoretically predicted value for the quantum Otto cycle operating in the adiabatic regime with thermal reservoirs.
There are several research directions that this study opens up.In this work, it was assumed that the baths manage to reach thermal states while they are decou- A system is initialized in a 101-dimensional thermal mixed state with ωT = 5 and α = 0 at τ = 0.It is then evolved numerically until τ = 5 using the fifth order Runge-Kutta method.For one set of simulations, α is set to a finite value at τ = 0 + and kept constant (labeled α fixed in the legend).For the second set, α increases linearly until it reaches its maximum value at τ = 5 (labeled as αmax).
The final state is calculated for 3 ≤ ε/δτ ≤ 16 and the trace distance is obtained between the smallest time step and all others for each configuration.The trace distance is then taken as the error.
pled from the working fluid.It is worth investigating the importance of the state being thermal and how the operation of the engine changes if it is not.More importantly, because the bath modes become "reset" while decoupled, they act as ancillae to connect the working fluid to thermodynamic reservoirs.If the bath modes do not have to be in thermal states and simply need to be able to exchange energy with the working fluid, it is interesting to explore the possibility of the confining potentials acting as the energy reservoirs.For example, the cold mode could be laser cooled to remove excess energy, while the hot mode experiences an external driving force that contributes energy to it.As an extension, it is useful to explore the possibility of keeping the bath modes connected to their respective reservoirs throughout the engine operation.
Finally, the operation protocol presented here was not optimized for power or efficiency.It would be useful to determine which parameters enhance the energy and power output of the engine.

ACKNOWLEDGMENTS
The author acknowledges the National Research Foundation, Prime Minister Office, Singapore, under its Medium Sized Centre Programme and the support by Yale-NUS College (through Start-up Grant).The author is grateful to Keian Noori and Silvia Lara for their input and discussion.To ensure numerical stability and physical realism, there are several guidelines that need to be followed when solving Eq. ( 6).Most obviously, the Fock space has to be represented by a finite number of states.When working with thermal states, the maximum level should be chosen so that n max ≫ ω T to guarantee the correct level occupancy, where ω T = k B T /ℏΩ is the thermal frequency.With n max fixed, the largest elements in the Hamiltonian are ≈ [1 + α(τ )/2]n max .Consequently, the time step δτ has to be substantially smaller than the period associated with this frequency: 2π[1 + α(τ )/2]n max ≪ 1/δτ .Hence, one needs to guarantee that ω T ≪ n max ≪ [2πδτ (1 + α max /2)] −1 .It is useful to define ε = [2πn max (1 + α max /2)] −1 so the time step requirement becomes δτ ≪ ε.
As the first step, it is important to demonstrate the accuracy of the Runge-Kutta approach by performing two benchmark procedures.For the first one, α(0) ̸ = 0 is kept fixed and the system is initialized in a thermal mixed state with the density operator ρ0 = e − Ĥ0/ωT /tr e − Ĥ0/ωT for Ĥ0 = ĝ † ĝ + 1/2.To illustrate the dependence of the numerical error on the step size, Eq. ( 6) is solved for several values of δτ and the final density operator is computed using these Û's.By taking the trace distance between each of the results and the density operator obtained for the smallest δτ , the convergence of the numerical results is observed.One should note that, because the Hamiltonian does not change in time, the analytical form of Û(τ, 0) = e −2πiτ Ĥ so that ρ(τ ) = e −2πiτ Ĥ ρ0 e 2πiτ Ĥ .The reason for not comparing the numerical results to the known analytical form has to do with the second benchmarking procedure where α(τ ) starts at zero and increases to some maximum value.In this case, the analytical result is generally not known and the best one can do is show the convergence of the numerical calculations.Hence, the same procedure is employed even in the constant-α case for the sake of consistency.The results for these two checks are given in Fig. 4, showing that one can achieve errors less than one part per million without making the time step drastically smaller than ε.As a balance between accuracy and speed, the time step for RK solutions in the main text is taken δτ = ε/5.

15 FIG. 1 .
FIG. 1. Nonadiabatic effects.A system initialized at a 41level thermal state at temperature ωT = kBT /ℏΩ with α = 0 (top row) or α = ω 2 − 1 (bottom row).Next, the state is either compressed (until it arrives at α = ω 2 − 1) or expanded (until α = 0) uniformly over time τα.The final energy of the system is then divided by what would be the energy of the system if the modification were adiabatic.The nonadiabatic effects are most prevalent when the rate of the confinement change is high (corresponding to large ω and small τα).Hightemperature compression demonstrates the effects of the finite Fock space, as discussed in the main text.

FIG. 2 .
FIG.2.Interaction with a series of baths.(a) Trace distance between the gas and bath modes as a function of time.The density operators include 41 Fock states."Expanded" ("compressed") corresponds to α = 0 (α = 8)."Heat" means that the oscillator starts at ωT = 1 and the bath is at ωT = 5; for "cool" the temperatures are switched.The interaction between the modes is Gaussian with σ = 1 and offset x0 = 1.The interaction strength starts with 1, later being reduced to 1/5 and, finally, to 1/20.The vertical lines mark the times when the bath is switched our for a new one.For dashed lines, Φ0 is unchanged; solid lines show where the interaction strength is reduced.(b) Trace distance between the final states from panel (a) and thermal states for a range of ωT .Vertical lines indicate the temperatures of the cold and hot baths.These lines coincide with the minima of the distance curves, indicating that the final gas state approaches the correct thermal state.

= 3 FIG. 3 .
FIG. 3. Engine operation.(a)-(d) Energies of the gas inan Otto engine at each phase of the cycle as a function of time.For the fast engine (top row), the time of each stroke is 1.For the slow engine (bottom row), the time of expansion/compression is 4, while the bath contact time is 10.The system is evolved following Eq.(9) with ω hot

25 FIG. 4 .
FIG.4.Numerical benchmarks.A system is initialized in a 101-dimensional thermal mixed state with ωT = 5 and α = 0 at τ = 0.It is then evolved numerically until τ = 5 using the fifth order Runge-Kutta method.For one set of simulations, α is set to a finite value at τ = 0 + and kept constant (labeled α fixed in the legend).For the second set, α increases linearly until it reaches its maximum value at τ = 5 (labeled as αmax).The final state is calculated for 3 ≤ ε/δτ ≤ 16 and the trace distance is obtained between the smallest time step and all others for each configuration.The trace distance is then taken as the error.