Quantum and classical dynamics of a three-mode absorption refrigerator

We study the quantum and classical evolution of a system of three harmonic modes interacting via a trilinear Hamiltonian. With the modes prepared in thermal states of different temperatures, this model describes the working principle of an absorption refrigerator that transfers energy from a cold to a hot environment at the expense of free energy provided by a high-temperature work reservoir. Inspired by a recent experimental realization with trapped ions, we elucidate key features of the coupling Hamiltonian that are relevant for the refrigerator performance. The coherent system dynamics exhibits rapid effective equilibration of the mode energies and correlations, as well as a transient enhancement of the cooling performance at short times. We find that these features can be fully reproduced in a classical framework.


Introduction
Historically, thermodynamics started with the goal of trying to understand how to convert heat into useful mechanical motion. For that purpose, steam engines have been developed which revolutionized our lives. Other useful thermal machines, such as refrigerators and heat pumps, followed. Naturally, these thermal machines are all large macroscopic entities that one uses classical laws to describe them and study their performances.
In the context of refrigerators, two quantum approaches stand out: the smallest possible refrigerator consists of a qutrit [29] or three interacting qubits [30], but a quantum absorption refrigerator can also be realised with three interacting harmonic modes [31,32]. This latter interaction, which has been demonstrated with trapped ions [33], is the object of this study. We focus on the unitary dynamics of the refrigerator itself and highlight phenomena like effective equilibration, as well as the challenge of identifying genuine quantum features. Elements on the process of refrigeration, i.e. the transfer of energy from a cold to a hot bath, mediated by the machine, are given in the Appendix.
Section 2 introduces the model Hamiltonian, the initial states that will be considered (independent thermal states of different temperatures) and the dynamical variables whose evolution is studied (occupation number of each mode). We discuss why we focus on the unitary dynamics. This dynamics is then solved and discussed in Section 3, first in general, then in the context of refrigeration. The following two sections are devoted to the two main observed features of the dynamics of the occupation numbers. In the short-term regime, there is a cooling enhancement which would be absent if the interaction were incoherent (Section 4). In the long-term regime, one observes effective equilibration of the occupation numbers even if there is no dissipative dynamics (Section 5). Finally, in Section 6 we study the classical model obtained by replacing the mode operators with conjugate variables in the Hamiltonian. We show that the dynamics of the occupation numbers exhibits the same features observed in the quantum formalism.

Model
The model we study here is a system of three interacting harmonic oscillators, with the free Hamiltonian and the interaction Hamiltonian It is convenient to work in the interaction picture whereĤ with ∆ = ω h − ω w − ω c and g the coupling constant. We will focus on the resonant case where ∆ = 0 in what follows. This Hamiltonian describes a wide range of physical processes: parametric amplification, frequency conversion, and second harmonic generation [34][35][36]. This work is based on considering the three interacting modes as an absorption refrigerator [31][32][33].
The refrigeration process works as follows. Each of the three modes is in contact with a thermal bath: the cold bath at T c , the hot bath at T h > T c , and the work mode at T w . We shall correspondingly refer toâ h ,â w andâ c as the hot mode, work mode, and cold mode, respectively. The dynamical variables of interest are the occupation numbersn i (t) = tr{ρ(t)â † iâ i } with i = h, w or c. Absorption refrigeration occurs if the interaction can induce a stationary heat flow from the cold to the hot mode (ṅ c (t) <ṅ c (0),ṅ h (t) >ṅ h (0)), the work mode providing a sufficient amount of free energy.
A full description of the refrigeration process takes into account both the unitary interaction generated by the trilinear Hamiltonian and the dissipative dynamics due to the interaction of each mode with its bath. This is a numerically involved task, further complicated by the subtle issue of formulating a master equation for composite systems coupled to multiple baths that is consistent with the second law of thermodynamics [37][38][39]. Yet, we are interested in the regime of weak coupling to the baths such that all the interesting aspects of the dynamics of refrigeration are captured by restricting the analysis to the unitary dynamics. This is confirmed in the Appendix, where we provide a comparison of the dynamics with and without the dissipative coupling. In the rest of this paper, we shall thus focus on the unitary dynamics, which has been experimentally implemented with trapped ions [33].
The initial state of the system consists in the three modes being prepared in uncorrelated thermal states, This corresponds to the situation in which each modes have been kept in contact with their respective baths for long time, before turning on the interaction. From here, focusing on the unitary dynamics implies that the modes are effectively decoupled from their baths during the evolution; this will capture the full refrigeration dynamics in the limit of slow thermalization rate.
Before proceeding, let us notice that some states of this family are stationary with respect to the unitary dynamics, and therefore also to the complete dynam- 3 Unitary dynamics

Methods and general features
We will mostly focus on temperatures T i that correspond to comparably small initial average occupation numbers,n i (0) = tr{ρ(t = 0)â † iâ i } with i = h, w or c. For later comparison with the classical framework, we will plot the initial average energy of each mode in the diagrams, i (0) =hω i (n i (0) + 1/2), instead of the occupation number.
At time t = 0, the interaction HamiltonianĤ int is switched on and the system evolves unitarily according to ρ(t) =Û ρ 0Û † withÛ = exp(−iĤ int t). This coherently transfers populations between the three modes and changes the average energies i (t). However, we note that even the closed system is not amenable to an exact analytical solution; previous studies focused on either short-time behavior, resorted to using products of coherent or Fock states as the initial state, or considered limiting cases of average occupation numbers much larger than one [34][35][36]40].
For an efficient simulation of the system dynamics, we take advantage of the fact that the interaction Hamiltonian couples only Fock states of the form |n h , n w , n c = |n, N − n, M − n with fixed integers N and M and 0 ≤ n ≤ min{N, M }. That is,Ĥ int is block-diagonal with respect to finite-dimensional subspaces characterized by the two conserved quantities N and M and the dimension d = min{N, M } + 1. The unitary evolution can then be efficiently computed by diagonalizing the Hamiltonian in each of the contributing subspaces, up to a cutoff for both N and M that truncates the originally infinite dimensional Hilbert space. For the simulations presented in this paper, the cutoff is chosen to ignore populations in the initial density matrix smaller than 10 −4 . Fig. 1 shows the typical time evolution of the modes' energy under the unitary dynamics for a thermal state of the form (4). In particular, we have chosen initial average occupations ofn h (0) = 0.5, n w (0) = 2.5, andn c (0) = 2.0, so that the order of magnitude of the contributing N and M are 10 0 or 10 1 . It is readily seen thatN =n h (t) +n w (t), and M =n h (t) +n c (t) are conserved.
The curves in Fig. 1 exhibit two noteworthy features. Firstly, the evolution starts with a transient oscillatory stage, where the largest fluctuation away from the initial value is found. We shall come back to this feature in Section 4. Secondly, the system energies seem to approach some apparent equilibrium values (dashed lines), around which only small residual oscillations persist. One can obtain those values from the infinite time-averaged state The observation that the expectation value of certain dynamical variables approach their long-time average, in spite of the fact that the system is not converging to a steady state because the dynamics is unitary, is a phenomenon known as effective equilibration [41]. We will get back to it in Section 5. Conservation of N and M implies that any increase in the average energy of the hot mode corresponds to an equal decrease in the work mode and cold mode energies. Hence it is sufficient to focus on just one of the modes. In Fig. 2, we plot the change of the average cold mode energy (or occupation number) between the long-time equilibrium state σ and the initial state for different initial energies of the work and the hot mode. The initial occupation number in mode c is fixed atn c (0) = 2.0. One sees that, depending on the values ofn h (0) andn w (0), the equilibrium energy of mode c can either increase (red), decrease (blue), or stay unchanged (thick line).
Let us remark again that, by preparing the ini- tial system in a product of thermal states and then bringing them into interaction, we mimic the thermodynamic situation where these three subsystems each have thermalized with their own baths and are then allowed to exchange heat among each other.
In Appendix A, we explicitly compare this two-step treatment with a simultaneous interaction-dissipation model based on a heuristic master equation. The latter contains the interaction Hamiltonian and independent thermal dissipators for each mode. When the dissipation rates are small compared to the coupling rate g, we find that the time evolution of the system energies is similar to that of the purely unitary dynamics. Specifically, the initial transient features are approximately the same, and the only relevant difference is a slight deviation of the steady-state energies from the long-time averages predicted by the unitary model. Strong dissipation rates, on the other hand, would only thwart the three-mode interaction dynamics and freeze the system state close to the initial thermal state.

Analogies and differences with refrigeration processes
After presenting the general features of the closed system dynamics governed by the trilinear Hamiltonian (3), we now focus on the thermodynamic aspects. As already mentioned, the interaction is capable of driving a quantum absorption refrigerator, in close analogy to the earlier proposed three-qubit fridge [30] with its interaction HamiltonianĤ int = g(|010 101| + |101 010|). The cooling performance of the latter has been extensively studied, and we will point out similarities and crucial differences for the present three-oscillator case. The initial hot and cold mode energies are fixed at h = 1.5hω h and c = 2.5hωc, while the work mode energy starts from w = 1.5hωw in the top panel and is increased in steps ofhωw to the bottom. Note that this still admits a lower temperature in the cold than in the hot mode, as ωc < ω h . The middle panel corresponds to a stationary configuration.
As proposed in theory [31] and backed by observations with normal modes of trapped ytterbium ions [33], the average energy of the c-mode can decrease or increase, depending on the initial work mode valuē n w (0). In this sense, we can speak of cooling or heating, as demonstrated in Fig. 3 for an exemplary choice of initial values,n h (0) = 1.0,n c (0) = 2.0, andn w (0) varying from 1.0 to 5.0 in unit steps. In fact, the contour plot in Fig. 2 can be viewed as a phase diagram, where the thick solid line separates the heating (red) from the cooling (blue) regime. For the latter, the initial occupation numbers must satisfȳ Since the initial states are thermal, the inequality translates into a relation between the temperatures, The cooling regime (7) is consistent with the general virtual qubit framework developed recently [42]. In this picture, one thinks of the hot mode and work mode forming a set of virtual qubits by the levels |n h , n w and |n h − 1, n w + 1 at an effective virtual temperature The cold mode at initial temperature T c then interacts with the virtual mode at T v . In the regime 0 ≤ T v < T c , cooling occurs since net energy is transferred from the cold mode to the virtual mode at a lower virtual Figure 4: Thermal Fano factor (9) of the infinite timeaveraged state reduced to the cold mode for varied initial conditions as in Fig. 2. A value of zero indicates a thermal energy distribution, i.e. that the energy variance is that of a Gibbs state. We find that the variance is consistently above (below) thermal in the cooling (heating) regime.
temperature. At T v = T c , which corresponds to the equilibrium condition (5), no energy exchange takes place between the two modes. Two important remarks are now in order: Firstly and unsurprisingly, the asymptotic three-mode state (6), which determines the long-time average energies after interaction, is in general not a tensor product of thermal states σ = ρ th h ⊗ ρ th w ⊗ ρ th c . Secondly, the reduced single-mode states tr j =i {σ} = ρ th i do not exhibit a thermal energy distribution either. This can be easily illustrated by looking at the second moments of the local distributions. The deviation of, say, the cold mode energy variance in the time-averaged state σ from that of a Gibbs distribution can be captured in the Fano factor Fig. 4 shows a contour plot of this quantity as a function of initial average energies as in Fig. 2. One sees that the variance is consistently higher (q > 0) than that of a thermal state of the same energy in the cooling regime, and lower in the heating regime. At the same time, we found that the single-mode entropy is always lower than that of the thermal state, suggesting that the time-averaged energy distribution is more biased.
Hence, and even though the reduced single-mode states remain diagonal at all times, we cannot assign temperatures to the asymptotic state σ that the trilinear system approaches in the long-time average. This is in contrast to the three-qubit scenario, where one can formally associate a temperature to any diagonal single-qubit state [43,44]. Here, the three modes are correlated and explicitly driven out of thermal equi-librium by the trilinear interaction (which remains the case in the presence of simultaneous weak bath couplings). This highlights the difference between genuine thermalization in open systems undergoing thermodynamic cycles and effective equilibration at an athermal state in quantum systems with non-trivial Hamiltonians [41,45,46].

Enhanced cooling in the single-shot regime
One important signature of the cooling dynamics that can be seen from Fig. 3 is that the cold mode always overshoots to a lower-energy state at transient time scales, before it approaches the long-time average value. Hence, if one can control the interaction and halt the dynamics at the appropriate time, the cooling performance can be enhanced. A similar, measurement-induced transient cooling of a single qubit was found in Ref. [47]. In the three-qubit scenario, this feature termed single-shot cooling was linked to the presence of quantum coherence [43,44]. Here, a similar effect can be observed as well [33]. The difference is that, in the three-qubit case, the system passes through many transient oscillations before reaching an equilibrium with help of the simultaneous thermalization with three independent baths. Here, the closed system alone approaches an effective equilibrium rather rapidly after the first overshooting oscillation. Irregular oscillations around the asymptotic state σ then prevail, but at comparably small amplitudes. The higher the initial temperatures, the more negligible these oscillations become.
Let us now examine the role of coherence in this transient effect. Following the qubit analogy [44], we compare the coherent energy exchange driven by the trilinear Hamiltonian to incoherent implementations of the same exchange process. The simplest incoherent model is to assume the excitation of the hot modê L =â † hâ wâc and the de-excitationL † occur spontaneously at the same rate γ. The corresponding master equation in the interaction frame reads aṡ which replaces the von Neumann equation with the trilinear Hamiltonian (2) of the coherent model. One can easily check that the master equation still conserves both N and M , and that it does not build up coherence between the combined Fock basis vectors |n, N − n, M − n ; initially diagonal states remain diagonal at all times. The latter facilitates an efficient numerical implementation. Alternatively, one can conceive a random unitary model where the coherent energy exchange driven bŷ H int ∝L +L † is switched on and off at random times  Fig. 3). The solid line represents the coherent evolution, which exhibits a characteristic transient overshoot to below the long-time average for gt < 1. This behavior is not predicted by the incoherent models (10) and (11) for the same energy exchange between the modes, based on spontaneous excitation processes (dashed) or on phase-incoherent exchange (dotted), respectively. The latter yields the correct long-time average. We assume exchange rates equal to the coherent coupling constant g.
with an average exchange rate γ, The master equation describes dephasing in the eigenbasis of the interaction Hamiltonian. The final equilibrium state is given by the fully dephased, or infinite time-averaged, initial state (6). Note however that the final state σ does contain nondiagonal coherence terms when represented in terms of the |n, N − n, M − n basis. At the same time, neither the coherent nor the two incoherent models generate coherence locally, i.e. nondiagonal elements in the Fock state representation of the reduced single-mode states. An exemplary comparison of the models in terms of the time evolution of the average cold mode energy is shown in Fig. 5, assuming γ = g for simplicity. As expected, neither the fully incoherent model (green dashed line) nor the dephasing model (orange dotted) reproduce the enhanced cooling feature of the coherent model (blue solid). Instead, both incoherent models lead to an exponential decay towards their equilibrium values; the one of the dephasing model (11) is the value that the coherent evolution oscillates around and approaches in the long time average.
The time evolution and the equilibrium states differ slightly here, whereas in the previously studied threequbit scenario there is no difference between the two incoherent models. This can be understood if we restrict to initial states ρ(0) from the two-dimensional subspace N = M = 1. The trilinear interaction is then equivalent to that of three qubits, and the two master equations (10) and (11) are identical.
Finally, we note that the comparison between the coherent and incoherent exchange process might suggest a quantum advantage in the transient cooling performance over "classical" implementations. This argument is based on the notion of "classical" states as a subset of quantum states without coherence in the relevant basis representation. In Section 6, we provide an alternative view by comparing the coherent quantum model to its counterpart in the fully classical framework. It will turn out that the cooling enhancement is not a genuine quantum effect. Indeed, the fact that the reduced single-mode quantum states always remain diagonal in the single-mode Fock basis suggests that the trilinear Hamiltonian only establishes coherence between the modes. Such a notion of inter-mode coherence, in the sense of a finite interference capability and well-defined phase relation over a given coherence time, can also be understood in classical optics [48].

Effective Equilibration
A striking feature of the studied three-mode interaction, which we have observed and exploited in the previous sections, is the fast effective equilibration of the average single-mode energies: even though the time evolution is strictly reversible, they appear to converge towards the stationary values of the timeaveraged, fully dephased state, with only minor oscillations around those values for times t 1/g. The phenomenon of effective equilibration has been explored in finite systems of high dimension [41,[49][50][51][52][53][54][55][56][57][58][59][60][61], and is observed here for the case of initially thermal states. The effect is most pronounced at high initial temperatures, when the thermal occupation of the modes extends to high-dimensional subspaces associated to large values of the conserved quantities N and M . A proper equilibration in the presence of a bath would lead to the same values, as long as the three-body coupling g remains the dominant rate (see Appendix A). This underpins the thermodynamical assessment based on the time-averaged steady state of the three modes as a model absorption refrigerator system in Sect. 3.2.
Quantitatively, one may speak of effective equilibration with respect to a set of observables if the cumulative time averages of their expectation values over the coherent time evolution quickly converge close to the values associated to a dephased stationary state [51,53,62,63]. The characteristic convergence time scale and the residual difference should typically decrease with the (finite or effective) system dimension [55,57]. This implies that recurrences of the initial values, if present, must be short-lived and rare under unitary evolution.
We observe such behaviour in the present, formally infinite-dimensional, scenario where the observables of interest are the single-mode energies. Our numerical studies show that both the characteristic equilibration time and the residual difference to the expectation values associated to the dephased state (6) decrease with increasing initial occupation numbersn i . We could not observe recurrences in our simulations, some of which extended over time intervals a few orders of magnitude greater than 1/g, see e.g. Appendix A. In Fig. 6(a), we plot the time at which the cold mode energy assumes its extremal overshoot value, which can serve as an estimate for the characteristic equilibration time. The times are extracted from unitary time evolutions with time increments of 0.01/g for varying initial conditions as in Fig. 2. We observe that the values decrease with growing mode temperatures and effective Hilbert space dimension covered by the initial state, except at thermal equilibrium (black line).
We note that for pure initial Fock product states |n, N − n, M − n , one rather observes persistent strong and irregular oscillations of the average energies instead of the described equilibration. The latter is a consequence of averaging over the incommensurate sub-spectra of the trilinear Hamiltonian associated to its different orthogonal subspaces. Nevertheless, signatures of effective equilibration also exist at the level of quantum states in high-dimensional subspaces, e.g. by looking at entangling dynamics [64].
As an example, we show in Fig. 6(b) how the reduced single-mode purity decays as a function of time for an initial Fock product state with n = 100 at varying subspace dimension. The blue, red, and green lines (top to bottom) in the double-logarithmic plot correspond to N = M = 300, 3000, and 30000, respectively. The dashed lines represent the reduced single-mode purities of the respective dephased states. As expected, the trilinear interaction rapidly entangles the three modes, and we observe a fast decay of single-mode purity close to the dephased values. Both these values and the characteristic half-life time decrease roughly in proportion to the inverse subspace dimension d = min{N, M } + 1.
On the other hand, rapid equilibration in finite high-dimensional (e.g. interacting many-body) systems has been linked to the non-Poissonian nature and spread of the underlying energy gap spectrum [41,[65][66][67], and it is assumed to appear in nonintegrable systems whose classical counterpart may exhibit signatures of chaos [68][69][70].
The energy gap spectrum of the quantum model varies strongly with the choice of the conserved quantities (N, M ), if restricted to a single invariant subspace. At large subspace dimension, non-Poissonian energy gap distributions are observed, in agreement with the observation of rapid equilibration. Three exemplary histograms of energy gaps (i.e. differences between subsequent energy levels) are shown in Fig. 6(c). They correspond to increasing values of M at fixed N = 10000 and subspace dimension d = N + 1, using a bin width of 2hg. In all cases, the probability does not drop with the gap size in a Poisson-like fashion, as associated to nonequilibrating Hamiltonians.
Taken together, these intricate equilibration features distinguish the resonant three-mode Hamiltonian (3) from other simple exchange coupling models. The resonant two-mode exchange coupling, for instance, can be diagonalized explicitly and leads merely to a normal mode splitting with a single fixed energy gap given by the coupling frequency. The same reversible dynamics can be seen in the twodimensional subspaces (for N = 1 or M = 1) of the present system. At high initial excitations, the threeoscillator interaction model presented here could provide a minimal, accessible, and practically relevant testbed for studies of effective equilibration.

Classical analysis
The present quantum refrigerator model of three trilinearly coupled harmonic oscillators admits direct benchmarking against the predictions of classical physics. We arrive at the classical version of the system straightforwardly by replacing the mode operatorsâ i andâ † i with canonical complex variables α i and α * i . The Hamilton function becomes The corresponding Hamilton equations of motion can be solved analytically after switching to the action-angle representation [71]. For this, the complex amplitudes are expressed in terms of magnitude and phase, where the ι i represents the energy contained in mode i in units ofhω i . As in the quantum case, one finds that the sums I 1 = ι h + ι w and I 2 = ι h + ι c are conserved, as well as the total phase difference Φ = φ w + φ c − φ h . Given the initial ι i (0) and φ i (0) and the additional constant of motion L = ι h ι w ι c cos 2 (Φ), this leads to the solution Here, sn(u | m) denotes the Jacobi elliptic function 1 and a > b > c are the three ordered solutions of the equations: a + b + c = I 1 + I 2 , ab + bc + ca = I 1 I 2 , and abc = L. The correct sign in the argument is given by that ofι h (0) ∝ − sin Φ(0). The solution (13) is periodic with period 2K(m), where K denotes a complete elliptic integral of the first kind. The time average of the energy contained in each mode is obtained by integrating (13) over one period, yielding a compact expression in terms of elliptic integrals of first and second kind This is the classical counterpart of the quantum long time average, which is the mean energy associated to a completely dephased state.  (12), averaging over random samples of initial conditions. The two curves each correspond to 10 7 trajectories drawn from a thermal distribution using either the same average mode energies (I, green) or the same initial temperatures (II, orange) as in the quantum case. The dotted lines indicate the infinite time averages.
For a comparison to the quantum simulation with thermal initial states, we evaluate the classical time evolution in a Monte-Carlo simulation by drawing initial angles φ i (0) according to a uniform distribution and ι i (0) according to an exponential distribution with mean valuesῑ i = k B T i /hω i . However, the different energy statistics leaves room for a certain ambiguity in the matching of quantum and classical initial conditions. One can either match the initial temperatures T i , arguing that the quantum and classical predictions shall be based on the same physical boundary conditions. Then the initial mean energies hω iῑi and i =hω i (n i + 1/2) will differ slightly. Or one matchesῑ i =n i + 1/2, which implies a better matching of the predicted energy trajectories at the cost of slightly different reservoir temperatures.
In Fig. 7, we compare both classical options to the quantum prediction for an exemplary initial configuration in the cooling regime. The classical simulations were averaged over 10 7 trajectories each. The relative difference to the quantum result decreases at higher initial temperatures where the continuous classical distribution of energies approximates well the discrete quantum statistics.
As expected, the classical result with matching initial energies (I) starts at the same point and remains closer to quantum case than the result with matching temperatures (II). Most importantly, the classical model predicts the same dynamical features as the quantum model, namely the short-time cooling enhancement and the effective equilibration at long times. The latter can again be attributed to thermal averaging: Even though single classical trajectories (13) describe periodic orbits, their periods depend on the initial energies and relative phases between the modes, and so they vary broadly over the thermal distribution of initial conditions. The averaged time evolution is not periodic and does not exhibit recurrences even for very long interaction times gt 1. In Section 4, we have seen that the transient enhancement known as single-shot cooling is linked to the presence of coherence between eigenstates of the interaction Hamiltonian in the quantum framework. Now the quantum-classical comparison shows that this feature is not relying on genuinely quantum coherence, but is rather a generic feature of the particular three-mode exchange interaction. Indeed, it appears as well in the classical framework where coherence between harmonic modes exists in the form of a fixed phase relation between their oscillations. The transient cooling enhancement can be suppressed by subjecting the classical phase coordinate Φ(t) to external noise, just as the feature will disappear in the quantum case if dephasing is present.

Conclusion
In this article, we have provided an in-depth study of a quantum absorption refrigerator model based on three resonantly coupled harmonic oscillator modes, as recently realized in an ion-trap experiment [33]. A closed-system analysis already exhibits the key features of the earlier studied analoguous three-qubit refrigerator model, in particular, the description of heating and cooling regimes in terms of virtual temperatures, and the enhanced cooling performance at transient time scales.
On the other hand, we have found important differences and new insights in the more complex threemode system. Even in the absence of simultaneous bath coupling, the initially thermal energy distributions of the cold, the hot, and the work modes equilibrate effectively around athermal distributions corresponding to a state that is fully dephased in the eigenbasis of the interaction Hamiltonian. This state exhibits residual coherences in the combined three-mode Fock representation, and it differs from the steady state of a fully incoherent implementation based on spontaneous excitation and de-excitation.
Most strikingly, we found that all essential features of the model can be explained by entirely classical means, with only minor quantitative differences that become irrelevant with growing temperatures. This includes, in particular, the enhanced transient cooling feature that had previously been attributed to quantum coherence in the three-qubit scenario. With the possibility of a one-to-one comparison between quantum and classical predictions at hand, the three-mode system may in fact serve as an ideal testbed to eluci-date the role of quantum resources in thermodynamics. Centres of Excellence Programme and through the Competitive Research Programme (Award No. NRF-CRP12-2013-03).

A Influence of simultaneous bath coupling
Here we investigate the performance of the refrigerator when the three harmonic modes are simultaneously interacting via the trilinear Hamiltonian and each coupled to their individual thermal reservoirs. For that, we employ the master equation where L j is the superoperator describing the effect of the bath coupled with mode j. We explicitly assume that each mode couples to its own bath locally and separately, i.e. we employ the usual Lindblad dissipators, with D[L]ρ =LρL † − 1 2 {L †L , ρ}. Here κ j is the coupling strength of the mode j with its bath andn j is the mean occupation number at the bath temperature T j . In the following, we assume that the bath couplings of the three modes are equal, κ h = κ w = κ c = κ.
It was pointed out that using local dissipators (A.2) in a system of inter-coupled modes might not always yield consistent predictions for transient quantities such as the energy flow between the subsystems [37][38][39]. However, since we are concerned here with the qualitative behavior of local observables (namely, single-mode energies), local dissipators may be used as an approximation.
In order to observe the refrigeration effect, the bath coupling κ should not be too strong compared to the trilinear coupling g. Otherwise the modes will stay thermalized with their baths and cannot evolve away from the initial state. The system dynamics will be frozen.
For κ g on the other hand, the dynamics is not expected to deviate much from the unitary evolution studied in the main text. Indeed, this is confirmed in Figure 8: Time evolution of the average cold mode energy for different coupling rates κ between the modes to their thermal baths. The open-system dynamics is described by the master equation (A.1), which was simulated numerically using the QuTiP package, whereas κ = 0 denotes the closed unitary evolution (solid line). We start from thermal initial states at h (0) =hω h , w (0) = 3hωw, and c(0) = 2.5hωc. Fig. 8, which compares the evolution of the average cold mode energy for various thermalization rates to the unitary case. In particular, the transient oscillations at short times are well reproduced for κ g (green dashed and blue dotted line) and the steadystate value approached in the long-time limit does not differ significantly from the infinite time average of the unitary case (black solid line). We ran simulations up to gt = 10 4 to confirm that deviations as significant as the initial oscillations do not recur.
Increasing the bath coupling κ to the level of g and beyond thwarts the coherent energy exchange between the modes; the transient oscillations disappear for κ ≥ g (red dash-dotted and purple fine-dotted line) and the mode energy quickly approaches equilibrium close to the initial value.
Concluding, the unitary dynamics studied in the main text approximates the transient refrigerator dynamics well in the relevant regime of operation, κ g, and it gives reasonably good account of the longtime cooling performance.