A quantum Otto engine with finite heat baths: energy, correlations, and degradation

We study a driven harmonic oscillator operating an Otto cycle between two thermal baths of finite size. By making extensive use of the tools of Gaussian quantum mechanics, we directly simulate the dynamics of the engine as a whole, without the need to make any approximations. This allows us to understand the non-equilibrium thermodynamics of the engine not only from the perspective of the working medium, but also as it is seen from the thermal baths' standpoint. For sufficiently large baths, our engine is capable of running a number of ideal cycles, delivering finite power while operating very close to maximal efficiency. Thereafter, having traversed the baths, the perturbations created by the interaction abruptly deteriorate the engine's performance. We additionally study the correlations generated in the system, and relate the buildup of working medium-baths and bath-bath correlations to the degradation of the engine's performance over the course of many cycles.


I. INTRODUCTION
The second law of thermodynamics prohibits extracting mechanical work from systems in thermal equilibrium. Therefore, in order to obtain work, one has to have access to systems out of thermal equilibrium. The theoretically simplest out-of-equilibrium system is one composed by two subsystems that are each at individual equilibrium and at different temperatures. This is the traditional setup for a heat engine: a working medium (WM) reciprocating between two thermal baths, pumps heat from the hotter bath (at temperature T h ) to the colder one (at temperature T c ) and outputs work as a result. The ideal engine converts the internal energy of the hot bath into work with an efficiency given by Carnot's formula, η C = 1−T c /T h . The idealizations needed for the machine to operate at such an efficiency are that (i) the baths interact with the working medium weakly [1,2], (ii) the cycle is a quasiequilibrium process and hence it takes infinite time to complete [1,3,4], and (iii) the baths are infinitely large [1,[5][6][7][8]. It has to be noted, however, that the size of the working medium itself is of no relevance -it can be anything from a two-level quantum system [9][10][11] to a giant steam engine [12].
Strictly speaking, conditions (i) and (ii) can never be satisfied: any interaction has finite strength and any process that can be observed takes finite time. In the generic setup where the bath is a many-body system with shortrange interactions and the WM couples to it locally, the breakdown of (ii) entails the failure of (iii) even if the bath is infinitely large [13]. Indeed, in such systems, the Lieb-Robinson bounds [14,15] imply that, roughly speaking, the correlations spread with finite velocity. This means that, in finite time, the WM can have access to only a finite region of the bath (see Ref. [16], where this idea was brought to use for the first time). However, it should be emphasized that the said finite region gets re-thermalized by the rest of the bath, so this scenario is not entirely equivalent to a finite bath.
The engine runs a strong-coupling adaptation of the Otto cycle [12]: the two "isochoric" thermalizations are intermediated by two "adiabatic" changes of the WM's Hamiltonian (see Sec. IV for the precise description). For the first cycle, the WM starts uncoupled from the baths and at equilibrium with the cold bath. This makes the initial state of the overall system a Gaussian state. Given that the total Hamiltonian is quadratic at any moment of time, the dynamics of the system can be described within the formalism of Gaussian quantum mechanics (GQM) [46]. The latter maps the intractable Schrödinger equation in the infinite-dimensional Hilbert space of the overall system onto a linear evolution of the finite-dimensional phase space. This allows us to perform a comprehensive analysis of the machine's operation without the need to adhere to any of the many approximations usually made when dealing with quantum open-system dynamics [44,45]. Moreover, by directly simulating the overall system's evolution, we gain access to the states of the baths at any moment of time, which allows us to reveal the physical mechanisms governing the degradation and eventual exhaustion of the initial disequilibrium provided by the baths in the finite-size, finite-time, and strong-coupling regime. With our ap-proach, we can easily work with baths of size up to 300 times the size of the WM with just a standard table-top computer.
The paper is organized as follows. First, in Sec. II, we give a short account on the notions from GQM that will be needed throughout the rest of the paper. This section is intended as an introduction and can be safely omitted by those familiar with GQM. In Sec. III, we describe the interaction of the WM with a single bath. In Sec. IV, we explore the physics of the Otto cycle, focusing first on the performance of the cycle (Sec. IV A), and then on the dynamics and the role of correlations (Sec. IV B). Finally, we summarize our conclusions in Sec. V. The codes, both in Matlab and Python, of all the numerical computations performed in this work are available in Ref. [47].

II. REVIEW OF GAUSSIAN QUANTUM MECHANICS
In this section, we review the formalism of Gaussian quantum mechanics, focusing on the aspects necessary for our study. For a much broader introduction to the topic, the reader is referred to Ref. [46]. Note that throughout this paper all expressions are given in natural units, i.e., we assume = k B = 1.
The primary computational advantage of this formalism is that it allows us to study interacting systems via a direct system-plus-bath perspective, without having to resort to perturbation theory [48] or other open-systems techniques. This provides access to the exact evolution of the bath in addition to the system, a fact we take great advantage of in this work.
Consider one or more quantum systems ascribed with bosonic canonical quadrature operators, satisfying the canonical commutation relations (CCRs), [q i , p j ] = i δ ij , where the indices label the systems (henceforth referred to as oscillators or modes). If one were to think about a harmonic oscillator with Hamiltonian P 2 2µ + µω 2 Q 2 2 , then a convenient choice of quadratures would be q = Q √ ωµ and p = P √ ωµ . In terms of the creation and annihilation operators, the quadratures are expressed through For a system of N modes, the quadratures form a phase space that we represent as the vector of operators Due to the CCRs, the phase space is a symplectic space, endowed with the structure [x a , x b ] = i Ω ab . Ω ab are the components of the so-called symplectic form, given by In GQM one works with Gaussian states. A state of an N -mode system is Gaussian if and only if it is an exponent of a quadratic form in {x a } 2N a=1 . Importantly, thermal states of quadratic Hamiltonians fall within this class. The defining feature of Gaussian states is that they are fully described by the first and second moments of their quadratures, i.e., their mean position and their variances in phase space. The mean quadratures of all the states we consider in this work will be zero, and so the formalism further simplifies. We thus characterize the state of our system via the 2N × 2N covariance matrix σ, the entries of which are given by An important aspect of GQM is that creating ensembles and performing partial traces is trivial. This is due to working in phase space rather than in a Hilbert space, where partitions are represented as a direct sum rather than as a tensor product. Thus, any combined state of two systems A and B takes the form where σ A and σ B are the reduced states of systems A and B respectively, and the matrix γ AB specifies the correlations between the systems. The superscript T denotes the operation of transposition. A fact crucial for GQM is that any unitary evolution generated by a time-dependent Hamiltonian that is quadratic at any moment of time will preserve the Gaussianity of a state [49]. Any such unitary, U , on the Hilbert space corresponds to a linear symplectic transformation on the phase space of quadratures: x → U † xU = Sx, with S satisfying The symplecticity of S, expressed by Eq. (5), ensures that the CCRs are preserved throughout the change of basis. On the level of the covariance matrix, it is easy to see that this transformation acts as A. Energy, Evolution, and Thermality Another convenient aspect of GQM is that it allows us to compute average energies, evolve the system over time according to some time-dependent quadratic Hamiltonian, and diagonalize the system into its normal mode basis without ever referencing a Hilbert space object.
The average energy of a state represented by the covariance matrix σ, with respect to a purely quadratic Hamiltonian H = x T F x, is given by The symplectic (i.e., unitary in the Hilbert space) evolution matrix S(t) generated by this (in general, timedependent) Hamiltonian obeys a Schrödinger-like equation: where F s = F + F T . For a constant Hamiltonian the solution trivially takes the form S(t) = exp(ΩF s t), and for general driven systems, the equation can be straightforwardly integrated by standard numerical techniques. When speaking of a "free" system we mean that we are working in the basis that diagonalizes the system's Hamiltonian. This is called the normal mode basis, in which the Hamiltonian takes the form where, in the second equality, we have ignored the (constant) zero-point energy.
The corresponding phase-space matrix is diagonal in this basis: . By definition, the normal modes do not interact with each other. This means that any thermal state on the entire system is given by the tensor product (in phase space, the direct sum) of the individual normal modes' thermal states.
In general, the system may have couplings between pairs of modes (for example, between nearest neighbours), which give non-diagonal elements to the matrix F . The normal-mode basis can be obtained by symplectically diagonalizing this matrix: SF S T = F free , where S is a symplectic matrix, and F free is diagonal as above.
In the normal-mode basis the covariance matrices of the system's thermal states are given by where ω i are the normal frequencies. We can thus find the thermal covariance matrix of any interacting system by first identifying the normal basis, specifying the covariance matrix σ as above, and then applying the inverse transformation to this matrix to put it back into the physical-mode basis. The values ν (th) i in Eq. (10) are referred to as the thermal state's symplectic eigenvalues. In general, every Gaussian state of N modes has N symplectic eigenvalues ν i , which are obtained by symplectically diagonalizing the covariance matrix: there always exists a symplectic matrix S such that The symplectic eigenvalues can be directly computed by taking the regular eigenvalues of the matrix iΩσ, which come in ±ν i pairs.

B. Entropy and Correlations
Consider a two-party state of the form of Eq. (4). The off-diagonal matrix γ AB contains the correlation functions between the two systems, and these systems are uncorrelated if and only if γ AB = 0. As a measure of correlations we use the mutual information, defined as Here, S(σ) is the von Neumann entropy of the state with covariance matrix σ, given by where This shows that the symplectic eigenvalues of a statewhich are invariant under symplectic transformationsgive a measure of mixedness for that state. For example, the entropy is zero, i.e., a Gaussian state is pure, if and only if all its symplectic eigenvalues are equal to one. Note that no state can have eigenvalues smaller than one (this is a statement of the uncertainty principle). We are thus able to very easily compute the mutual information across any partition in our system, independent of how many modes each partition contains.
Note that the entanglement is also computable, but, for most situations, it is considerably more difficult. In the particular case of two modes it is nevertheless easy [50], and we discuss some findings in that regard in the next sections. However, due to the thermality of our system, quantum correlations are hard to maintain, and we have found that generally entanglement does not play a significant role in the scenarios we consider below. Interestingly, this aspect is in accord with (yet by no means logically necessitated by) the fact that, although capable of manifesting many interesting quantum features, GQM is an essentially classical, noncontextual sector of quantum mechanics in that it can be described by a local hidden variable model [51].

III. GAUSSIAN INTERACTION WITH A SINGLE BATH
Before performing the analysis of the Otto cycle, let us study some relevant features of the isochoric interaction of the WM with a single thermal bath. We will thereby introduce the specific Hamiltonians that describe the components of the Otto engine in the next sections.
Throughout this work, we model thermal baths as collections of harmonic oscillators arranged in onedimensional, translation-invariant rings with nearestneighbour interactions.
We consider only positionposition couplings so that the free Hamiltonian of a bath is given by where N is the number of oscillators in the bath, ω b is the bare frequency of each of them, and α controls the coupling strength. Note that, because of the periodic boundary conditions, q N +1 = q 1 . The phase-space matrix corresponding to this Hamiltonian is where 0 is the 2 × 2 matrix of zeros, and At the beginning of the process, the bath is initialized in a thermal state Due to the interactions, the covariance matrix will not be given by a simple direct sum as in Eq. (10). Rather, we must first identify the normal mode basis that symplectically diagonalizes the Hamiltonian matrix ω bath /2 = SF bath S T , where ω bath is the diagonal matrix composed of the normal mode frequencies of H bath . We then identify the thermal state as per Eq. (10) and finally transform back to the physical basis to find the thermal state of the ring, σ bath = S −1 σ T (S −1 ) T (this calculation is included in the function Initialize in Ref. [47]).
As a WM we employ yet another harmonic oscillator, with bare frequency ω m . Its coupling to the bath is described by where {int} is the set of the bath's nodes which the WM interacts with. For most of our analysis, we choose this set to contain just one and always the same node of the bath, which we label the first node q 1 . However, also a situation where the interaction set has more than one element is discussed in Sec IV A.
The function λ(t) is a switching function that modulates the interaction in time. In particular, we choose the following compactly-supported, smooth switching function where τ ≥ 2δ is the total duration of the interaction with the bath and δ is the time that takes to fully switch on (and fully switch off) the interaction. We will refer to δ as the ramp-up time of the isochoric interaction. The phase-space matrix of the overall Hamiltonian will matrix containing all zeros except for the first entry γ 11 = λ(t)γ, corresponding to the q m q 1 interaction we are imposing. In order to construct the symplectic evolution matrix S(t) of the overall system, which is generated by F tot , we numerically integrate Eq. (8). The covariance matrix at the moment t will then be simply given by σ is the initial state of the WM (its values being given by Eq. (10)) [47].
During the evolution of the overall (WM plus bath) system, we have found that the state of the WM remains very close to being thermal. That is, at any moment of time, its covariance matrix is very close to that given by Eq. (10) for some ν (th) (for a detailed discussion and an explicit characterization of the distance of the actual state of the WM to a thermal state, see Appendix A). Given this, we are able to assign a meaningful effective temperature to the WM by computing the temperature associated with its symplectic eigenvalue. An example of this is shown by the green, solid line of Fig. 1.
Importantly, we notice that at t ≈ 93 the temperature of the WM becomes equal to that of the bath. We note that this moment is not the thermalization time in the proper sense because the interaction is still on. Rather, the exact thermalization time is τ th ≈ 98.5. It turns out that, in order to achieve precise thermalization, one needs to match the frequencies, so that ω m is also the frequency of the individual bath oscillators (the ω b in Eq. (15)), and the couplings, so that the WM-baths interaction strength is equal to the intra-ring coupling strength, i.e., γ = α. Intuitively, this matching ensures that the rate of information transfer between the WM and the bath is the same as between oscillators within the bath. Whenever ω m (resp. γ) is outside a small neighbourhood of ω b (resp. α), the WM does not thermalize with the bath at all (similar frequency filtering phenomenon in the classical setting was reported in [52]). In view of this, we from now on set ω m = ω b and γ = α.
With such a configuration of the parameters, and fixed δ/τ , the thermalization time τ th scales as for α 1. In fact, the above scaling is, to a good approximation, preserved also for large α. A remark is in order here. As τ th , we choose the smallest τ that achieves thermalization. Since the bath is finite and the WM couples to it strongly, the temperature of the WM will not approach the bath's temperature T b in a monotonic way: with passing time, the WM's final temperature will first go slightly above T b -the maximum being T b +O(α 2 ) [53] -then go below T b , and continue an oscillatory behavior as that depicted in Fig. 1 for t > 100.
Another important aspect of our thermalization process is that, due to the finite duration and the finite strength of the WM-bath interaction, it has a non-zero work cost. More specifically, the extracted work, as quantified by the difference between the initial and final average energies of the total system, is not zero. However, despite the strong non-equilibrium character of the process, this amount is small (compared to, e.g., the energy exchanged between the WM and the bath). In fact, for small α, this work cost, W i (where the subscript i stands for isochoric), scales as and is almost independent of the ramp-up time, δ. Taking, for example, together with the fact of exact thermalization discussed above, means that the fine-tuning of the frequencies and couplings provides us with an example of almost workfree thermalization in finite time, resulting from a strong interaction between the WM and the bath. A similar example, where the structure of the bath is known and the Hamiltonian of the WM is finely-tuned, was constructed in [54]. This is not a standard, exponential relaxation behaviour [45], and it can be argued that such behaviour cannot occur for general baths of unknown structure [24]. We also note that the fact of almost zero work justifies the usage of the term "isochoric" for this process. Indeed, in this case, most of the energy exchange is heat transfer, which is the characteristic of isochoric processes [12]. In the strong coupling regime, strictly isochoric (or, equivalently, constant-Hamiltonian processes) cannot exist as any non-zero coupling will change the system Hamiltonian, and therefore the term needs to be adapted.
Furthermore, subtle processes such as the evolution of the correlations between the WM and the bath or information exchange between the WM and the bath can be examined in very great detail within the framework of GQM. As an illustration, in Fig. 1 we examine the evolution of the correlations between various partitions during the interaction.
We compare the correlations, as measured by the mutual information, between the WM and the whole bath (dashed blue line), between the WM and the node in the bath it interacts with (black dotted line), and between the latter and the rest of the bath (red dot-dashed line). This provides us with a number of insights into the nonperturbative interaction of the WM and the bath. First, during the phase of switching on the interaction, the WM and the bath quickly build up strong correlations, which decay later on. This decay is caused by the fact that the bath nodes to which the WM is coupled also interact with the rest of the bath, and the bath, due to its tendency to thermalize, forces these correlations to decay. We can see this process in more detail by examining the other (dashed blue) between the WM and the bath as a whole, (black dotted) between the WM and the specific bath oscillator with which it interacts, and (red dot-dashed) between this oscillator and the rest of the bath. The bath contains N = 30 oscillators, all with frequency ω b = 2, and is initialized in a thermal state at temperature T b = 4. The WM has also frequency ωm = ω b , but is initialized in a thermal state at temperature Tm = 0.5. The total time of interaction is τ = 245, the ramp-up time is δ = 0.1τ , and the interaction strengths are α = γ = 0.1. For these parameters, the exact thermalization time as defined in the main text is τ th = 98.5. Note that the red curve is not initially zero (and it should not be, because there is initial correlation from the ring couplings). However, since we are working with a relatively hot bath, these correlations are very small (of the order of 10 −3 ), and its magnitude cannot be appreciated in full detail in the figure.
two lines. The correlation between the machine and the interacting node similarly rises and then falls, and the decay occurs exactly as this node becomes significantly correlated with the rest of the bath.
This gives us an important intuitive picture. The interaction between the WM and the bath generates correlations between the two (specifically, between the WM and the interacting node). Due to the intra-bath couplings, the WM also becomes correlated with other ring nodes in an outwards-propagating manner. However, these couplings also mean that the correlation between the WM and bath will, over time, be swapped to correlations between different bath nodes, as we see, for example, in the red dot-dashed line in Fig. 1. Over the course of many interaction sessions, the bath nodes therefore become more and more intercorrelated, which will eventually result in a halt of the machine. We elaborate on this process in the next section, where we discuss the performance of a WM operating cyclically between two finite-sized baths.
One does not need to move to the two-bath scenario to observe the effects of having finite-sized baths, though. In fact, one only needs to interact with the bath for a time that is long enough. We do so in Fig. 2, where we compute the effective temperature of the machine after interacting with a bath composed of N nodes during a time τ , for different values of N and τ , and all the other parameters being the same as those used for Fig. 1.
In Fig. 2, we observe two very distinct behaviors that are clearly separated. For τ < c·N , where c indicates the slope of the "causal cone", the temperature of the WM is insensitive to the size of the bath. Indeed, the interaction time in this case is short enough so as to allow the interaction to finish before the perturbations that propagate through the bath return to the region which interacts with the WM (i.e., the interacting node of the ring). Therefore, there is no difference between the temperature that the WM achieves in this case and the temperature that it would achieve from interacting with an infinite bath. The opposite occurs for τ > c · N : in this case, the interaction time is long enough so as to permit the perturbations generated by the interaction with the WM to return to the interacting node. These perturbations modify the local state of the interacting node, which in turn translates into a response in the WM that diverges from that expected for infinite baths.
It is also worth noting that, for short interaction times, the WM does not have enough time to fully thermalize with the bath. We observe that the effective temperature of the WM increases with the interaction time until the point where thermalization is achieved. After this point, increasing the interaction time further has no major influence on the WM's temperature until, of course, it is long enough for the perturbations to go around the bath. FIG. 3. Visualization of the Otto cycle. The standard sequence of isochoric thermalisations and adiabatic compressions/expansions defining the Otto cycle in phenomenological thermodynamics are, in our case, implemented as a sequence of q − q interactions (as described in Sec. III) and sudden changes of the WM's Hamiltonian. More specifically, the cycle consists of the following steps: (i) the WM interacts with the hot bath (the red harmonic chain) by a coupling that is smoothly switched on, kept constant, and smoothly switched off, (ii) the Hamiltonian of the WM is suddenly changed so that the frequency matches the individual frequencies in the cold bath, (iii) the WM is brought into contact with the cold bath (the blue harmonic chain), with the same pattern of interaction as in step (i), and (iv) the Hamiltonian of the WM is suddenly changed back to its original value.

IV. THE GAUSSIAN OTTO CYCLE
We now study the performance of the WM running an Otto cycle between two thermal baths at temperatures T h (hot) and T c (cold), as depicted in Fig. 3.
The Otto cycle we consider is composed of two isochoric interactions between the machine and each of the baths (as described in Sec. III) separated by two sudden changes of the WM's Hamiltonian. Specifically, between subsequent interactions, we instantaneously swap the WM's Hamiltonian, so that the WM's state remains unchanged. The fact that the WM is detached from the baths during the swap ensures that the process is adiabatic, i.e., thermally isolated, in the thermodynamic sense [55]. It is important to note that the change in Eq. (22) is not equivalent to simply quenching the frequency of the oscillator. Rather, it requires simultaneously changing both the mass and the frequency: µ → µ ωc ω h and ω c → ω h . Here, ω c and ω h are the frequencies of the WM used during the interactions with the cold and hot baths, respectively. For a discussion of the case when only the frequency is quenched, see Appendix B. As mentioned in Sec. III, ω c and ω h are chosen to coincide with the frequencies of the nodes of, respectively, the cold and the hot baths. Moreover, we also match the interaction strength with the ring coupling strengths, i.e., γ = α c = α h (we choose α c = α h for simplicity only, without losing generality).
The total work extracted during a cycle is given by the sum of works extracted during each of the four parts of the cycle. As we showed in Sec. III, the work contributions from WM-bath interactions are small, hence most of the work is generated during the adiabats. The work produced by a sudden change in Hamiltonian in Eq. (22) is given by a particularly simple expression. Indeed, since the baths remain intact during the adiabat, the work is given by the energy change of the WM only. For example, in the adiabat after an interaction with the hot bath, the energy of the WM is decreased by W h→c = (ω h − ω c ) Tr σ m . If we choose ω c < ω h , not only will W h→c > 0, but also the net work will be positive. Note also that if ω c > ω h we would be running a refrigerator. Lastly, note that the engine cycles are not cyclic in the standard thermodynamic sense. Indeed, since the baths are finite and the interaction with WM perturbs them non-negligibly, at the end of each cycle, the state of the WM (and, of course, the state of the baths) will be different from that at the beginning. Nevertheless, the deviation from cyclicity is small during the period of "perfect" cycles that we describe below.

A. Cycle performance
We begin the cycle by the interaction with the hot bath, so the starting Hamiltonian of the WM is ω h a † m a m and its state is thermal, at temperature T c with respect to ω c a † m a m . Due to the finite size of the baths and the strong perturbations that the interactions with the WM causes in them, we expect that the performance of the engine will drop over time. This intuition is confirmed in Fig. 4, where we plot the work output and efficiency of the engine as a function of the number of cycles of operation.
In Fig. 4a, each bar represents the work output during an adiabat, and the red line represents the total work output in each cycle (the sum of the works in the adiabats plus the sum of the works in the isochores), as described above. The heat Q is defined as the energy the hot bath loses per cycle. We define the energy of the bath with respect to the Hamiltonian in Eq. (15), and, for the n-th cycle, the heat is given by where σ h (t) is the covariance matrix describing the state of the hot bath as a function of time. The efficiency of the engine is defined as usual: η = W/Q.
In Fig. 4, we see that the engine's performance has two regimes. First, the work output and absorbed heat are approximately constant, decreasing very slowly, for the first 15 complete cycles. We call these "perfect" cycles. The degradation of the engine's performance during these cycles is due to the residual perturbations near the interaction site that the outward-propagating perturbations created by the WM-bath interaction leave behind. As the cycles proceed, these small deviations from the interaction site's equilibrium state accumulate, causing the gradual decrease in work and heat.
What is more, during the perfect cycles, the perturbations, created by the WM-bath interaction, propagate through the baths in the same way as they would do were the baths infinite. Therefore, any given perfect cycle is unaffected by the further increase in the size of the baths. This implies that the engine's degradation in the course of perfect operation is not a finite-size effect, and hence occurs also when the baths are infinite. Moreover, the difference between the work outputs in, say, the first and second cycles, does not vanish when the coupling is taken to zero. Therefore, the perfectregime degradation is not a strong-coupling effect either. Rather, it strongly depends on the ramp-up time and can be decreased noticeably by increasing δ/τ . However, going to high values of δ/τ prevents the WM from thermalizing with the baths, thereby impairing the functioning of the engine. Hence, the degradation cannot be eliminated completely in our model. We explicitly compute the correction to the optimal figures of merit due to this degradation in Appendix C. We note that, while the dependence of single-cycle characteristics on the ramp-up time is in line with the general intuition that non-commutative, time-dependent interactions generate excitations that cause thermodynamic friction (see, e.g., [26,56]), the important fact of the cycle-to-cycle accumulation of the imperfections caused by finite switching time is a separate phenomenon.
We furthermore observe that the number of perfect cycles, N p , increases asymptotically linearly with N ≡ N c = N h , the number of nodes in the baths, as is to be expected given the constant, finite speed of propagation of the perturbations in the bath [57]. However, when the interaction time τ is close to τ th , N p does not depend on α for small α. Indeed, although the thermalization time increases with decreasing α (see Eq. (20)) and this requires longer interaction times τ with the baths, the propagation of perturbations within the baths also slows down, and the two effects almost exactly compensate each other. Along with the fact that the degradation is slow, the linear dependence of N p on N makes the perfect regime relevant for practical engines with large baths.
Differences from the perfect-cycle behavior begin to appear only when the perturbations return to the region of the bath that directly interacts with the WM. This is the point at which the finite-size effects take relevance, and it is marked by the drastic, discontinuous drop in the work output in Fig. 4. The performance of the engine becomes unreliable due to large variations that heat and work undergo both in magnitude and sign. The above discontinuous behaviour of the engine's figures of merit is contrasted with the conventional gradual degradation of the performance of an engine operating between finite reservoirs (see, e.g., [40]). The contrast is further sharpened by the observation that, as is also the case in the said conventional picture, the baths diverge from their initial states in a gradual, continuous manner. This is illustrated in Fig. 4b, where the distance, as measured by relative entropy [58], of the hot bath's state at the beginning of the i-th cycle, ρ (i) h , to the bath's initial state, ρ pendix D for a more detailed discussion), is plotted as a function of i. In Fig. 4b it can be seen how, during the "perfect" cycles, the characteristics of the engine stay almost constant despite the fact that the bath's state changes at a steady rate.
An important implication of Fig. 4 is that, during the perfect cycles, the efficiency of the engine η is very close to η O = 1 − ω c /ω h . The latter is the theoretical maximum for an oscillator running an idealized Otto cycle between two infinite thermal baths to which it is coupled weakly enough for the standard Markovian open quantum system techniques [45] to be applicable [56,59,60]. Such idealized engines are known to obey the so-called power-efficiency trade-off, which states that the power output of the engine has to approach to zero whenever the efficiency comes close to the reversible maximum (see, e.g., [4,24,60]). Our model respects the power-efficiency trade-off for the Otto cycle in the following manner: for α 1, the efficiency approaches η O from below as while for the work output of a perfect cycle we have . We refer the reader to Appendix C for a more detailed discussion on these quantities. Taking into account Eq. (20), this leads us to Here P is the power output of the engine: P = W/τ cycle , where τ cycle = 2τ is the duration of the cycle. We note that, although the setting of our problem is different from that in Ref. [38], the scalings in Eqs. (24) and (25) agree with (and saturate) the optimal scalings derived there. One can also consider coupling the WM to more than one, evenly spaced ring sites. It turns out that adding more interacting sites reduces the amount of perfect cycles, which matches the intuitive picture described earlier. Indeed, the reduced distance between the sites leads to shorter time needed for the perturbations created by the interaction to reach the nearest site of interaction. Interestingly, the work output of a single perfect cycle is insensitive to the cardinality of the set {int} (as long as the perturbations generated in one interacting site do not have time to arrive to any other), but of course the total work output of the engine over several cycles does get reduced by increasing the number of interaction points. On the other hand, the more sites the WM interacts with, the smaller is the time necessary for it to thermalize. This leads to an increased power output for the initial perfect cycles, albeit at the cost of decreasing the number of such cycles.

B. Propagation of correlations
As noted before, the formalism presented in Sec. II allows for an easy way of identifying whether two systems are correlated. In this subsection, we use this property to study how correlations distribute along the baths and the WM. We consider this as one of the (probably many) paths to obtain a better understanding of the phenomenology presented above.
In Fig. 5, we show the strength of the correlations between the WM and each of the oscillators in each bath and how these correlations evolve in time for five consecutive cycles. Throughout this subsection, each bath is composed of N = 30 oscillators, with all other parameters being the same as in Sec. IV A. The vertical lines denote the instants of time at which the machine stops interacting with one bath and, after the corresponding adiabat, begins interacting with the other.  Fig. 3. Note how during the first three cycles the machine observes no differences from the interaction with infinite baths, and after this point the perturbations in the chains arrive back to the interacting oscillator, modifying its local state.
One feature we immediately observe is the explicit propagation of the perturbations in the form of localized wavepackets at finite speed, in full agreement with the Lieb-Robinson bound [14,15]. Although it is hard to see in Fig. 5, while propagating, these wavepackets leave residual perturbations behind. The latter are small and, during the first three cycles of operation (t ∈ [0, 600]), the WM appears to interact with almost unperturbed baths. These are the perfect cycles described above. The time t = 600 is when the perturbations that were generated during the first three cycles manage to intercept the WM as it is currently interacting with the bath, leading to the sudden drop of the work output that has been discussed in Sec. IV A.
We also observe that the propagating correlations quickly fade. This is unsurprising, and carries the same explanation as that given for Fig. 1. Our computations show that, to a surprisingly good approximation, during an interaction, the WM becomes correlated with just a single non-local mode in the bath-the mode that propagates outwards-as can be appreciated in Fig. 5. However, both the WM and this propagating mode are interacting with the rest of the bath as well, and thus this correlation is quickly lost and distributed among bath modes. This also explains why the decay occurs much faster in the hot bath than in the cold bath. Indeed, the hotter the bath, the larger the thermal noise that will break the correlations.
It is also instructive to observe how the correlations are built up and distributed along the baths. This is illustrated in Fig. 6, in which the mutual information between each pair of oscillators in each bath is shown at various times. One immediately notices the outwardpropagating nature of these correlations. An important insight into the process of the engine's degradation is gained by looking at the bath-bath correlations instead. Indeed, given that the WM acts as a carrier of both energy and correlations between the baths, and that the baths gradually evolve away from their initial states, one would expect that, over time, the baths get more and more correlated and end up reaching a global passive state (see Ref. [61] for the characterization of passivity within GQM). We explore this intuition in Fig. 7, in which we show that, surprisingly, the mutual information between the two baths remains close to zero during the ideal cycles, and starts abruptly increasing after the last ideal cycle is complete. This can be explained by noticing in Fig. 5 that, during the perfect cycles, the WM is virtually uncorrelated with the baths both at the beginning and at the end (but not in the middle) of each interaction session, which means that the WM does not transmit correlations during these cycles. This picture obviously breaks down once the perturbations reach the interaction site. This thereby establishes a clear quantitative link between the correlations among the baths and the optimal performance of the engine.
It is worth noticing that, despite the fact that the mutual information between different elements of the system can be substantially large, for our choice of parameters, none of the correlations built involve entanglement. In- deed, it is well known that entanglement in quantum fields decays very rapidly with temperature, reaching zero at a finite value [50]. However, for a sufficiently cold bath, one could still expect some entanglement to be present, although it is not clear whether it will play a significant role in the engine's performance.

V. SUMMARY AND CONCLUSIONS
Using the formalism of Gaussian quantum mechanics we have been able to circumvent the standard assumptions of weak coupling, slow driving, and infinite size of the baths, usually employed in studying thermodynamic phenomena. The focus of our study was on a single, driven harmonic oscillator undergoing an Otto cycle between two finite harmonic thermal reservoirs. Despite the attention the physics beyond these assumptions has received in recent years, to the best of the authors' knowledge, this is the first work where none of these assumptions is made.
We first study the interaction of a machine with a single bath, modeled as an N -mode translationally-invariant harmonic ring. GQM allows to observe not only how the machine thermalizes, but also how the interaction creates correlations between the WM and the region of the bath that directly interacts with it, and how these correlations later on propagate across the bath.
In our study of the quantum Otto cycle, we conclude that the crucial element that determines the performance of the engine is the propagation of the perturbations created by the WM-bath interaction. During the first cy-cles of operation, the WM interacts with the baths in such a way that an infinite-sized-baths behaviour is observed. We call these cycles "perfect" and find that the figures of merit of the engine during these cycles remain almost unchanged with the efficiency being very close to its optimal value, while maintaining finite power and respecting the power-efficiency trade-off. We furthermore observe that the perturbations generated during these interactions propagate through the baths as wavepackets moving with constant velocity. These wavepackets leave residual perturbations behind, which causes a slow, but gradually-accumulating degradation of the engine's performance -an effect that persists even for infinitely large baths. After enough time, the wavepackets return to the interaction region and start disrupting the thermalization of the WM, thereby drastically affecting the work output and efficiency of the machine. We expect this picture to also hold beyond the Gaussian regime, provided the speed of sound within the baths is finite [14,15].
We have also explored the interplay between the degradation of our engine and the creation of correlations within the overall system. As discussed, the process of running the engine inevitably creates an increasing number of correlations with the baths, within the baths, and among the baths. Remarkably, the bath-bath correlations remain very close to zero during the perfect cycles, and start increasing abruptly right after. This represents the overall system's gradual evolution to a more and more passive state. We believe that further study into this interesting dynamics is warranted.
By its own example, this work demonstrates the capabilities of Gaussian quantum mechanics as a workhorse for assessing finiteness effects in a field that has historically relied on infinite (time, size, and subtlety of the interactions) idealizations. Not only can GQM address fundamental questions in quantum thermodynamics, but it also provides us with more tractable numerical computations. This we believe may be of great use for the community of quantum thermodynamics.
Note added -After having submitted this manuscript, we became aware of related work [62], which studies a finite-size Gaussian engine consisting of a single-oscillator working medium operating an Otto cycle between two heat baths composed of single oscillators each. Due to the strong coupling to the bath during the isochoric interaction, the state of the WM will in general acquire non-diagonal terms, leading to the state being, in general, not a thermal state. Indeed, the thermal state is a function of the Hamiltonian and hence cannot have non-diagonal terms in the energy eigenbasis. In this appendix, we introduce a measure of athermality for Gaussian states of an oscillator and show that for the isochoric interactions we consider in the main text (Secs. III and IV) the athermality is negligible, especially at the end of the interaction. This additionally justifies the view that the thermal bath thermalizes the WM.
Let a single-oscillator Gaussian state be described by a covariance matrix σ m = ν 1 κ κ ν 2 . If the state were thermal, then its covariance matrix would simply bê σ m = ν 0 0 ν , i.e., ν 1 = ν 2 ≡ ν and κ = 0. Then, following Eq. (10), its temperature would be (A1) This will not be, however, the case of our WM for every instant during the interaction with a bath. In order to prescribe a temperature to a general state of the WM we take its covariance matrix, σ m , symplectically diagonalize it toσ m = ν 0 0ν as prescribed in Sec. II, and compute its temperature via Eq. (A1). The obtained temperature will be the effective temperature of the original state described by σ m . This is the definition we use in Fig. 1. Note that, for thermal states, this effective temperature coincides with the real temperature of the system. In order to define the athermality of the state given by σ m , ρ(σ m ), we compute the Uhlmann fidelity [58] between ρ(σ m ) and the thermal state at the effective temperature, ρ(σ m ), given by F [σ m ,σ m ] is equal to 1 if and only if σ m =σ m , and is < 1 otherwise.
For Gaussian states the Uhlmann fidelity can be directly expressed in terms of their covariance matrices. For purely quadratic states (which recall is the case of thermal states) the fidelity is given by [63] where the quantities A and B are given by Now we can look at how much the state of the WM, as given by the covariance matrix σ m (t), differs from the thermal state at temperature T ef f -given by the covariance matrixσ m (t)-at any moment t during the isochoric interactions described in Secs. III and IV. To that end, we define the athermality of the state at a moment t as Due to the properties of the fidelity, we thus have that 0 ≤ A(t) ≤ 1 and A(t) = 0 iff σ m (t) =σ m (t), i.e., if the state of the system is really thermal.
In Fig. 8 we plot the athermality of the WM for an example isochoric interaction with a bath consisting of N = 30 oscillators. We clearly see that, except for a period in the beginning of the process, the WM is rather close to being thermal. It is especially interesting that by the end of the process, the WM is almost exactly thermal. This information, in addition to the final temperature of the WM shown in Fig. 1, indicates that the WM is thermalized by the bath.
where Q m and P m are the canonical position and momentum of the WM, and µ is its mass. In this case, the quadratures and the creation-annihilation operators change as well: In these terms, the change in Eq. (B1) takes the form Now, if we perform this change instantaneously, the initial thermal state of the WM, ρ ∝ e −Hm/Tm , will remain unchanged. Since [H m , H m ] = 0, this will result in the state having coherences in the new energy eigenbasis [60], thereby significantly reducing the efficiency of the engine (this can be straightforwardly deduced from the analysis presented in Appendix C). If, on the other hand, as a result of Eq. (B1), the state were also changed to its covariance matrix, as defined by Eq. (3) with x m instead of x m , would remain unchanged. If now the interaction Hamiltonian, H int (see Eq. (18)), would couple to the bath degrees of freedom with the new quadrature, q m (instead of q m ), the dynamics of the overall covariance matrix would be the same as that presented in the main text. Put in other words, were we to change to the new quadratures in all formulas, the dynamics on the level of covariance matrices and Hamiltonian matrices would remain unchanged. Hence, we would have exactly the same results for such an engine that those shown in the main text.
For such a modification to work, we need to have ρ evolving into ρ as a result of Eq. (B3). Since the instantaneous change in Hamiltonian cannot produce any changes in the state, let us allow the change to take some non-zero time τ ad . The problem is now the following: is it possible to devise a quadratic Hamiltonian path connecting H m to H m that is capable of evolving the state in time τ ad so that the covariance matrix remains the same?
To answer the above question, let us notice that the covariance matrix remains unchanged as a result of quantum adiabatic evolution. The latter is defined by the unitary evolution operator where |n Hm and |n H m are the n-th eigenvalues of, respectively, H m and H m (see, for instance, Ref. [64]). Indeed, we have that σ ab = Tr(U ρU † (x a x b + x b x a )) = Tr(ρ (x a x b + x b x a )), where we used the fact that Now, as is shown in Ref. [64], for a single oscillator, a shortcut to adiabaticity can be constructed for implementing the unitary in (B5) by adding the timedependent term to the Hamiltonian of the oscillator for the period of the adiabat τ ad . Here, the function ω m (t) is arbitrary provided that it satisfies ω m (0) = ω m and ω m (τ ad ) = ω m . In order for this new cycle to coincide with the one in the main text, we need τ ad to approach to zero. This would require very fast generation of the term (B8) and very quick driving of the frequency. Although we leave the question of the experimental accessibility of such quick driving open, we note that the codes provided in the computational appendix [47] straightforwardly allow for simulating an Otto cycle with any non-zero τ ad .

Appendix C: Full-cycle energetics
In this appendix, we perform a detailed analysis of the work and heat involved in the perfect cycles of operation of the WM. As in the main text, we take as a starting point the moment of the cycle where the WM is still in a state ∝ e −ωca † m am/T (0) c but its Hamiltonian has already been changed to ω h a † m a m , and the isochoric interaction with the hot bath is about to start. We label the initial temperature of the bath with superscript (0) to indicate that it was its temperature before the first cycle started. At this moment, the total Hamiltonian is just the sum of the individual Hamiltonians, and the energy before starting the isochore is thus where we have defined After the isochore, the total Hamiltonian returns to being the sum of the individual Hamiltonians and the state of the WM is again thermal (albeit now correlated with the bath). Hence, the energy of the system right after the isochore is where T (1) h is the temperature of the WM after the interaction with the hot bath. Recall that, as discussed in Section IV A, although the parameters can be chosen so that T (1) h = T h exactly, this temperature does not need to be exactly the temperature of the bath T h . In the example in the main text, namely, when N c = N h = 30, ω c = 1, ω h = 2, T c = 0.5, T h = 4 and τ = 100, T (1) h is slightly greater than T h : T (1) h − T h ≈ 1.3 × 10 −4 . This small difference does not affect our analysis because such deviations, if not appearing in the first cycle, do appear in the subsequent ones. The work extracted during the hot isochore is then where H h is the Hamiltonian of the hot bath, Z h = Tr e −H h /T h , and [58] S(ρ||σ) = Tr [ρ(ln ρ − ln σ)] .
The relative entropy has several features desirable for a distance measure. In particular, those that are of interest for the case studied are that S(ρ||σ) ≥ 0, (D3) i.e., that it is a positive quantity, and that S(ρ||σ) = 0 iff ρ = σ. (D4) Although the relative entropy does not satisfy the triangle inequality and is not symmetric, which means it is not a distance measure in the proper sense, it decreases monotonically under completely positive trace preserving operations [58], which makes it a distinguishability measure of choice in many situations [58].
Whenever the second argument in S(•||•) is a Gibbs state, as, e.g., is the case in Eq. (D1), where S is the von Neumann entropy, E h is the energy of the hot bath at the beginning of the i-th cycle, and F = E − T S is the free energy. Eq. (D6) means that relative entropy measures the "thermodynamic" distance between ρ (i) h and ρ (1) h , which additionally motivates our choice of relative entropy as a distance quantifier.
The quantities in Eq. (D5) can be readily calculated directly in the covariance-matrix picture via Eqs. (7), (13), and (14). An example of such a calculation is presented in Fig. 4b, where it can be seen that the relative entropy distance of the state of of the bath from the initial state of the bath increases at a steady rate as the "perfect" cycles proceed. This is contrasted by how the efficiency and work per cycle remain almost constant during these cycles.