Reconciliation of quantum local master equations with thermodynamics

The description and manipulation of energy transfer at the quantum scale is a problem of fundamental and technological importance, receiving increasing attention in recent years. On one side, the study of energy transport through mesoscopic or atomic systems has found applications in disparate systems: from solid state devices [1, 2] to light harvesting complexes [3–6]; from ultra-cold atomic systems [7–9] to trapped ions [10–12]. On the other side, coherent manipulation of energy is at the core of quantum thermodynamics, whose aims include the understanding of the emergence of thermodynamic laws from quantum mechanics [13–21] and the design of thermal machines made with quantum devices [22–41].

The starting point for addressing these phenomena is the theory of open quantum systems [42–45]. In this framework, the reduced dynamics of the state of the system under scrutiny, when in contact with an environment, is cast in the form of a master equation after a series of approximations. Under this category fall very common models studied in the literature, such as the spin-boson or the Caldeira–Leggett models [46–48].

The situation becomes more involved when the system S is made up of several interacting subsystems S₁, ..., S_n and each subsystem interacts with a local environment E_i, possibly held at different temperatures T_i (see figure 1(a)). This type of scenario is the basis for the description of transport properties in quantum systems. Microscopic derivations in this case usually lead to global master equations (GMEs), in which the environment introduces jump operators that allow for incoherent excitation transfers between the different subsystems [49–53]. However, these derivations are in general quite involved since they require knowledge of the full set of eigenvalues and eigenvectors of the system's Hamiltonian, something which quickly becomes prohibitive when the number of subsystems increases. Moreover, depending on the approximations employed, one may also arrive at equations which do not generate completely positive maps (the so-called Redfield equations [50]), or equations which contain unphysical heat currents [54]. For these reasons, microscopic derivations of master equations for systems connected to multiple environments still continues, nowadays, to be a topic of great interest.

Zoom In Zoom Out Reset image size

An alternative, more heuristic, approach consists in deriving a master equation for the individual subsystems, neglecting the interaction with the remaining subsystems. The resulting master equation will then contain only local jump operators describing exchanges between the environment E_i and its corresponding subsystem S_i. Such equations, which we shall henceforth refer to as LMEs (also frequently called boundary-driven master equations), are typically accurate when the dissipation rates are larger than the interaction between subsystems. Due to their computational simplicity, they have been extensively employed over the last two decades in the study of transport in non-equilibrium quantum systems [55–73].

It turns out, however, that the nonlocal terms neglected in the LME may still lead to non-thermal steady-states [74] and play a significant role if the heat exchanges are small, even for weakly interacting parts. As a consequence, it has been found that LMEs may lead to apparent thermodynamic inconsistencies, as pointed out recently by Levy and Kosloff [75]. They have shown that the LME for two coupled quantum harmonic oscillators (QHO) may predict currents from a cold to a hot thermal reservoir, or the existence of currents even in the absence of a temperature gradient. Moreover, both may occur even in the limit of weakly coupled oscillators. The origin of these effects, as shown in [75], lies in the fact that the heat fluxes become of the same order of magnitude as the neglected nonlocal terms They thus suggest to use the GME to fix such thermodynamic anomalies. These findings generated a considerable stream of investigations on the comparison between global and LMEs [49, 52, 53, 61, 76–82] as well as clever possible alternatives [83, 84].

Recent results indicate, however, that it is possible to construct a consistent thermodynamic framework for LMEs, resolving these seeming contradictions. First and foremost, it is important to mention that, unlike Redfield equations, LMEs are in Lindblad form and therefore generate completely positive trace preserving (CPTP) maps. Second, LMEs, being CPTP maps, can be microscopically derived in a controlled way using the idea of repeated interactions (collisional models) [55, 58, 85–90]. In this approach, each environment is divided into an ensemble of identically prepared auxiliary systems, called units, which interact sequentially with an individual subsystem for a short time τ (see figure 1(b)). In the limit τ → 0 this leads to a LME, irrespective of the internal system interactions. Thirdly, Barra recently showed that since the interaction between the auxiliary units and the system is time dependent, there is an inherent external work required for generating the dissipative evolution described by the LME [88]. Hence, it is possible to put forth a consistent thermodynamic framework for the repeated interactions scheme, including definitions of heat, work and entropy production [88, 89].

The goal of this paper is to advance further on this reconciliation between LMEs and thermodynamics. Using the techniques of eigenoperators, we show that while LMEs satisfy local detailed balance, in general global detailed balance is broken. Moreover, we find that there is a fundamental work cost associated with this breaking of global detailed balance. By taking this exact work cost into account, we then show that LMEs become fully reconciled with thermodynamics. Quite surprisingly, the discrete nature of the repeated interactions method allows us to establish a direct mapping between non-equilibrium steady-states and limit-cycles of quantum heat engines. Hence, we find that the anomalous situation in which heat flows from cold to hot baths corresponds simply to a refrigerator, requiring positive injection of work to operate.

As our main application, we focus on the case of a system of two coupled QHO, which has been the subject of considerable interest due to potential applications in trapped ions [10, 11], opto-mechanical systems [8, 91] and ultra-cold atoms [8, 9]. Several interesting features of this model are studied in detail. First we show how to engineer the system to behave as a refrigerator, a heat engine or an accelerator. We study the efficiency and/or coefficient of performance (COP) for this model and show that they reach the optimal value of the Otto cycle in the case where the oscillators interact without counter-rotating terms in the Hamiltonian: including counter-rotating terms only degrades the operation of the system. We also discuss how to reconcile these results with a theorem by Martinez and Paz [92], which states that no autonomous refrigeration is possible with equilibrium harmonic oscillators linearly coupled.

Let us mention that, given a physical system, the actual choice of whether to describe it using an LME or a GME approach depends ultimately on the system itself. Here we stress that the LME approach per se does not give rise to thermodynamic inconsistencies, provided the mechanism for its emergence is fully taken into account.

In order to clearly present the motivation behind our study and the main features of our approach, we begin in section 2 with the simplest example of a LME, consisting of two coupled QHOs. We first review the seeming thermodynamic inconsistencies that may stem from this model and then go on to show how the results obtained using the method of repeated interactions can be used to completely remove them. Then, in section 3 we carry out our main theoretical development, showing how the eigenoperator technique may be used to directly access thermodynamic quantities and describe the possible breaking of global detailed balance. Finally, in section 4 we return to the oscillator model and show how the full machinery of our approach may be employed to study the interplay between different types of interactions between the two oscillators. Finally, in section 5 we summarize our results and conclude. Our paper also contains three appendices with various details of the derivations of our main results. In appendix A we review the technique of Lyapunov equations for dealing with Gaussian continuous variable systems. In appendix B we provide additional mathematical details on the developments of section 3. In appendix C we compute the entropy production showing that it is always positive.

In order to motivate our study and provide an intuitive summary of our main results, we begin by discussing the model studied by Levy and Kosloff [75], consisting of two coupled harmonic oscillators, S₁ and S₂, each coupled to its own local heat bath via a LME (see figure 1(a)). The Hamiltonian is taken to be

$$\begin{eqnarray}&&{H}_{S}={H}_{1}+{H}_{2}+{H}_{I},\end{eqnarray} \tag{ 1 }$$

where

$$\begin{eqnarray}&&{H}_{i}=\displaystyle \frac{1}{2}({p}_{i}^{2}+{\omega }_{i}^{2}{x}_{i}^{2})={\omega }_{i}({a}_{i}^{\dagger }{a}_{i}+1/2),\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&{H}_{I}=\epsilon ({a}_{1}^{\dagger }{a}_{2}+{a}_{1}{a}_{2}^{\dagger }).\end{eqnarray} \tag{ 3 }$$

Here a_i are the annihilation operators of each harmonic oscillator and $({x}_{i},{p}_{i})$ the corresponding position and momentum (ℏ = k_B = 1). Both oscillators are assumed to have unit mass, but different frequencies ω_i. The interaction (3) is taken to be of the typical tight-binding form, which conserves the number of quanta in the case of equal frequencies, ω₁ = ω₂. More general interactions will be discussed in section 4.

The system is also subject to two baths E₁ and E₂, which we assume can be modeled by LMEs acting only on S₁ and S₂. Hence, we take the system to evolve according to the Lindblad master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{S}}{{\rm{d}}t}=-{\rm{i}}[{H}_{S},{\rho }_{S}]+{D}_{1}({\rho }_{S})+{D}_{2}({\rho }_{S}),\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{D}_{i}({\rho }_{S})=\gamma ({n}_{i}+1){ \mathcal L }[{a}_{i},\rho ]+\gamma {n}_{i}{ \mathcal L }[{a}_{i}^{\dagger },\rho ].\end{eqnarray} \tag{ 5 }$$

Here γ is the dissipation rate, n_i is the average excitation number for the Bose–Einstein distribution

$$\begin{eqnarray}&&{n}_{i}=\displaystyle \frac{1}{{{\rm{e}}}^{{\beta }_{i}{\omega }_{i}}-1},\end{eqnarray} \tag{ 6 }$$

with inverse temperature ${\beta }_{i}={({T}_{i})}^{-1}$ and

$$\begin{eqnarray}&&{ \mathcal L }[A,\rho ]=A\rho {A}^{\dagger }-\displaystyle \frac{1}{2}\{{A}^{\dagger }A,\rho \},\end{eqnarray} \tag{ 7 }$$

is the dissipator in Lindblad form. Henceforth, we always assume that T₁ < T₂ (see figure 1(a)). Below, it will also be useful to keep in mind that n_i is a monotonically increasing function of T_i/ω_i.

2.1. Predictions for the steady-state

Due to the form of the Hamiltonian and the jump operators, the steady-state of the master equation (4) will be Gaussian and can therefore be easily found using the standard Lyapunov equation technique. A review of this method is given for completeness in appendix A. Here we only focus on the main results, which can be summarized by the following expectation values:

$$\begin{eqnarray}&&\langle {a}_{1}^{\dagger }{a}_{1}\rangle ={n}_{1}+\displaystyle \frac{2{\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}({n}_{2}-{n}_{1}),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&\langle {a}_{2}^{\dagger }{a}_{2}\rangle ={n}_{2}-\displaystyle \frac{2{\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}({n}_{2}-{n}_{1}),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\langle {a}_{2}^{\dagger }{a}_{1}\rangle =\displaystyle \frac{\epsilon }{{{\rm{\Delta }}}^{2}}({\omega }_{1}-{\omega }_{2}+{\rm{i}}\gamma )({n}_{1}-{n}_{2}),\end{eqnarray} \tag{ 10 }$$

where ${{\rm{\Delta }}}^{2}={\gamma }^{2}+4{\epsilon }^{2}+{({\omega }_{1}-{\omega }_{2})}^{2}$. Using equation (2) we then find that the energy of oscillator 1 will be

$$\begin{eqnarray}&&\langle {H}_{1}\rangle =\langle {H}_{1}{\rangle }_{\mathrm{th}}+\displaystyle \frac{2{\omega }_{1}{\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}({n}_{2}-{n}_{1}),\end{eqnarray} \tag{ 11 }$$

where $\langle {H}_{1}{\rangle }_{\mathrm{th}}={\omega }_{1}({n}_{1}+1/2)$ is the energy the oscillator would have if it were in thermal equilibrium with its corresponding bath. An equivalent expression can be obtained for the second oscillator replacing $1\leftrightarrow 2$. Similarly, the heat flux from the first oscillator to the second will be

$$\begin{eqnarray}&&{\dot{Q}}_{1\to 2}=-{\rm{i}}\langle [{H}_{S},{H}_{1}]\rangle =-{\rm{i}}{\omega }_{1}\epsilon \langle {a}_{1}^{\dagger }{a}_{2}-{a}_{1}{a}_{2}^{\dagger }\rangle =\displaystyle \frac{2\gamma {\epsilon }^{2}{\omega }_{1}}{{{\rm{\Delta }}}^{2}}({n}_{2}-{n}_{1}),\end{eqnarray} \tag{ 12 }$$

which is consistent with the expression obtained in [75].

The relevant point to highlight from these results is that the sign of the ²-terms in equations (11) and (12) will depend on the relative magnitudes of n₁ and n₂ and not of T₁ and T₂. This simple fact, when combined with the definition of n_i in equation (6), introduces several properties which, at first sight, seem to violate the second law of thermodynamics. First, it implies that oscillator 1 can be cooled to a temperature smaller than the temperature T₁ of the coldest reservoir, in the sense that it could reach a state with lower energy than that of the thermal state it would reach if it were in contact only with the coldest reservoir. Second, one may obtain a non-zero heat current even if both baths are at the same temperature, T₁ = T₂. And third, if T₁ < T₂ but we happen to satisfy

$$\begin{eqnarray}&&{n}_{1}\gt {n}_{2}\quad \leftrightarrow \quad \displaystyle \frac{{T}_{1}}{{\omega }_{1}}\gt \displaystyle \frac{{T}_{2}}{{\omega }_{2}},\end{eqnarray} \tag{ 13 }$$

then there will be a current flowing from the cold to the hot bath (${\dot{Q}}_{1\to 2}\lt 0$).

As pointed out in [75], these results seem to violate the second law of thermodynamics, a fact which the authors attribute to the inadequacy of LMEs to describe transport processes. The result of the cooling of the first oscillator below the temperature of the coldest environment also seems to suggest that this master equation could be used to design Gaussian (quadratic) absorption refrigerators, which would contradict a general theorem proved by Martinez and Paz [92].

2.2. Resolving the thermodynamic inconsistencies

Using equation (4) one finds that the total rate of change of the system Hamiltonian (1) is given by

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}\langle {H}_{S}\rangle }{{\rm{d}}t} & = & \gamma ({n}_{1}-\langle {a}_{1}^{\dagger }{a}_{1}\rangle )+\gamma ({n}_{2}-\langle {a}_{2}^{\dagger }{a}_{2}\rangle )\\ & & -\gamma \epsilon \langle {a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1}\rangle .\end{array}\end{eqnarray} \tag{ 14 }$$

However, this change in energy cannot be associated only with heat flowing to the reservoirs [88]. Instead, there is an associated work cost which, as we show below, is related to the breaking of global detailed balance. Quoting the results that will be derived in section 3, we now discuss how one can remove all the inconsistencies by including this work cost.

In particular, we find that the correct splitting of equation (14) is of the form

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {H}_{S}\rangle }{{\rm{d}}t}={\dot{Q}}_{1}+{\dot{Q}}_{2}+{\dot{W}}_{\mathrm{ext}},\end{eqnarray} \tag{ 15 }$$

where

$$\begin{eqnarray}&&{\dot{Q}}_{i}=\mathrm{tr}[{D}_{i}({\rho }_{S}){H}_{i}]=\gamma {\omega }_{i}({n}_{i}-\langle {a}_{i}^{\dagger }{a}_{i}\rangle ),\end{eqnarray} \tag{ 16 }$$

is the heat rate flowing from each reservoir to the system (see equation (41) for the general expression) and

$$\begin{eqnarray}&&{\dot{W}}_{\mathrm{ext}}=\mathrm{tr}\{[{D}_{1}({\rho }_{S})+{D}_{2}({\rho }_{S})]{H}_{I}\}=-\gamma \epsilon \langle {a}_{1}^{\dagger }{a}_{2}+{a}_{2}^{\dagger }{a}_{1}\rangle ,\end{eqnarray} \tag{ 17 }$$

is the required work rate (see equation (43) for the general expression). From equations (8)–(10) we then find that in the steady-state

$$\begin{eqnarray}&&{\dot{Q}}_{1}=\displaystyle \frac{2\gamma {\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}{\omega }_{1}({n}_{1}-{n}_{2}),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{\dot{Q}}_{2}=-\displaystyle \frac{2\gamma {\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}{\omega }_{2}({n}_{1}-{n}_{2}),\end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray}&&{\dot{W}}_{\mathrm{ext}}=-\displaystyle \frac{2\gamma {\epsilon }^{2}}{{{\rm{\Delta }}}^{2}}({\omega }_{1}-{\omega }_{2})({n}_{1}-{n}_{2}).\end{eqnarray} \tag{ 20 }$$

These are very important results, which we shall now discuss in depth. But before doing so, some consistency checks are in order. First, ${\dot{W}}_{\mathrm{ext}}$ tends to zero as  → 0, in which case global detailed balance is recovered. Second, the work required is proportional to γ². As discussed in [52], the LME can be shown to be correct for small interactions up to the same order (see comment in section 2.4 in that paper). Third, in the steady-state ${\rm{d}}\langle {H}_{S}\rangle /{\rm{d}}t\,=\,0$ so that we should have

$$\begin{eqnarray}&&{\dot{W}}_{\mathrm{ext}}=-{\dot{Q}}_{1}-{\dot{Q}}_{2},\end{eqnarray} \tag{ 21 }$$

as can indeed be easily verified.

We now show that, once this work cost is correctly attributed, the two-oscillator system at steady-state functions precisely as a heat engine, with two reservoirs and a work source. Indeed, we find from equations (18)–(20) that the system can be tuned to function as a refrigerator, a thermal engine or an accelerator (oven). The description of each mode of operation and what they mean in terms of ${\dot{Q}}_{1}$, ${\dot{Q}}_{2}$ and $\dot{W}$ is presented in table 1. We also show in figure 2 a plot of these 3 quantities illustrating the typical operation regimes.

Zoom In Zoom Out Reset image size

Table 1.  The three modes of operation of the steady-state of the two oscillators, always assuming E₁ is the cold reservoir (T₁ < T₂). In this paper positive heat or work always means energy entering the system. So ${\dot{Q}}_{i}\gt 0$ means energy entered S through E_i and ${\dot{W}}_{\mathrm{ext}}\gt 0$ means work was performed on the system by an external agent.

Mode of operation	${\dot{Q}}_{1}$	${\dot{Q}}_{2}$	$\dot{W}$	Operation regime
Refrigerator	${\dot{Q}}_{1}\gt 0$	${\dot{Q}}_{2}\lt 0$	$\dot{W}\gt 0$	ω1 < ω2T1/T2
Consumes work to make heat flow from cold to hot.
Engine	${\dot{Q}}_{1}\lt 0$	${\dot{Q}}_{2}\gt 0$	$\dot{W}\lt 0$	ω2T1/T2 < ω1 < ω2
Uses hot bath to produce useful work, dumping the remainder in the cold bath.
Accelerator (oven)	${\dot{Q}}_{1}\lt 0$	${\dot{Q}}_{2}\gt 0$	$\dot{W}\gt 0$	ω1 > ω2
Consumes work to heat the cold reservoir more than it would with spontaneous thermal conduction.

We begin with the case of a refrigerator, ${\dot{Q}}_{1}\gt 0$. According to equation (18), this requires n₁ > n₂ which is tantamount to the condition (13). But we now see that when this condition is satisfied, one also has ${\dot{W}}_{\mathrm{ext}}\gt 0$. Thus, for the two oscillators to act as a refrigerator, they must absorb work from an external agent. Hence, we see that there are no violations of the second law. Instead, by taking into account this work cost, one recovers a fully consistent thermodynamic description. The COP of the refrigerator is the ratio of the heat extracted from the cold environment divided by the work invested:

$$\begin{eqnarray}&&\mathrm{COP}=\displaystyle \frac{{\dot{Q}}_{1}}{{\dot{W}}_{\mathrm{ext}}}=\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}-{\omega }_{1}}\end{eqnarray} \tag{ 22 }$$

which coincides with the COP of the Otto cycle. As will become clear in section 3, this connection between a steady-state and the limit-cycle of a heat engine, can actually be traced back to the basic idea of the repeated interactions method, which models the steady-state as a sequence of (extremely short) strokes, each having an associated work rate and heat rate. Notwithstanding, we should also mention that an Otto COP is a special consequence of this type of interaction. As will be shown in section 4, the introduction of counter-rotating terms only reduces the COP below the Otto value.

Increasing ω₁ in the interval

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{1}}{{T}_{2}}{\omega }_{2}\lt {\omega }_{1}\lt {\omega }_{2},\end{eqnarray} \tag{ 23 }$$

we find that the two-oscillator system in its steady-state functions as a thermal machine, extracting heat from the hot bath, producing work and dumping heat in the cold bath. In this regime the machine's efficiency will be:

$$\begin{eqnarray}&&\eta =\displaystyle \frac{| {\dot{W}}_{\mathrm{ext}}| }{{\dot{Q}}_{2}}=1-\displaystyle \frac{{\omega }_{1}}{{\omega }_{2}}\end{eqnarray} \tag{ 24 }$$

which is again the Otto efficiency for a thermal machine. Finally, for ω₂ < ω₁ we necessarily have n₁ < n₂ and therefore ${\dot{Q}}_{1}\lt 0$, ${\dot{Q}}_{2}\gt 0$ and ${\dot{W}}_{\mathrm{ext}}\gt 0$. Thus the device transports heat from the hot to the cold environment and at the same time transforms work into heat which is damped in the cold bath: $| {\dot{Q}}_{1}| ={\dot{W}}_{\mathrm{ext}}+{\dot{Q}}_{2}\gt {\dot{Q}}_{2}$. Hence, it functions as an accelerator (i.e. an oven) [93], heating the cold reservoir faster than it would with spontaneous thermal conduction.

Notice that at the special value ω₁/ω₂ = T₁/T₂, the system operates as a thermal machine with the Carnot efficiency and zero power.

We conclude this section by emphasizing that within the repeated interactions framework one is able to resolve all of the seeming thermodynamic inconsistencies of LMEs. The zeroth law does not apply since the system is externally driven. The first law is satisfied, as we accounted for the energy balance of the whole system plus environment, including the required work cost. The second law is satisfied, as refrigeration is not spontaneous, but accompanied by external work. And finally, the LMEs do not violate the theorem by Martinez and Paz [92], which states that no autonomous refrigeration is possible with equilibrium harmonic oscillators. In fact, in our framework, there is also work associated with the environment–system interaction so the machine is not autonomous.

We now move on to the mathematical development of LMEs using eigenoperators and the repeated interactions.

We now turn to the microscopic derivation of LMEs using the method of repeated interactions [55, 58, 85–90, 94–102], which will form the basis for our thermodynamic description. The basic idea behind the method of repeated interactions is to model the environment as a set of units, each prepared in a thermal state. At time t = 0 the system S is allowed to interact for a certain time τ with one such unit, which we generically refer to as E. After this time the unit is discarded and a fresh new one is introduced, again in the same initial state. The process is then repeated sequentially. For τ → 0 this method generates a continuous time description which can be modeled by a Lindblad master equation that is, in general, local. We emphasize that a larger class of non-Markovian master equations including the GME can be engineered in a similar way [87] and realized experimentally as proposed in [103].

We consider here this method implemented in the general scenario of figure 1(b) and derive the results of the previous section as particular cases. Our system S is therefore assumed to be composed of N subsystems S₁, ..., S_N, where each S_i interacts with its own environment E_i. The Hamiltonian of the system is taken to be

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \sum _{i=1}^{N}{H}_{{S}_{i}}+{H}_{I},\end{eqnarray} \tag{ 25 }$$

where ${H}_{{S}_{i}}$ is the Hamiltonian of subsystem S_i and H_I summarizes all interactions between subsystems. Moreover, the total Hamiltonian comprising the system, the environments and all interactions, is taken as

$$\begin{eqnarray}&&{H}_{\mathrm{tot}}={H}_{S}+\displaystyle \sum _{i=1}^{N}{H}_{{E}_{i}}+\displaystyle \frac{1}{\sqrt{\tau }}\displaystyle \sum _{i=1}^{N}{V}_{i},\end{eqnarray} \tag{ 26 }$$

where ${H}_{{E}_{i}}$ is the Hamiltonian of environment E_i and V_i is the interaction between S_i and E_i. We also follow the customary approach of scaling the V_i by the interaction time τ, which is convenient, although not necessary, for taking the continuous time limit [58, 88, 89]. Following the usual repeated interactions approach, the reduced density matrix of the system ρ_S(t) will then evolve according to the map

$$\begin{eqnarray}&&{\rho }_{S}((n+1)\tau )={\mathrm{tr}}_{E}\{{{\rm{e}}}^{-{\rm{i}}\tau {H}_{\mathrm{tot}}}\ {\rho }_{S}(n\tau ){\rho }_{E}\ {{\rm{e}}}^{{\rm{i}}\tau {H}_{\mathrm{tot}}}\},\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&{\rho }_{E}=\underset{i=1}{\overset{N}{\bigotimes }}\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{i}{H}_{{E}_{i}}}}{\mathrm{tr}{{\rm{e}}}^{-{\beta }_{i}{H}_{{E}_{i}}}},\end{eqnarray} \tag{ 28 }$$

is the thermal state of the environments.

Expanding equation (27) in a power series in τ one then finds, up to first order, the LME

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{S}}{{\rm{d}}t}=-{\rm{i}}[{H}_{S},{\rho }_{S}]+\displaystyle \sum _{i=1}^{N}{D}_{i}({\rho }_{S}),\end{eqnarray} \tag{ 29 }$$

where H_S is given in (25) and

$$\begin{eqnarray}&&{D}_{i}({\rho }_{S})=-\displaystyle \frac{1}{2}{\mathrm{tr}}_{{E}_{i}}[{V}_{i},[{V}_{i},{\rho }_{S}{\rho }_{{E}_{i}}]].\end{eqnarray} \tag{ 30 }$$

The interesting aspect of this approach is that the structure of the dissipative terms are completely independent on the choice of H_S. Or, more specifically, on the interaction terms H_I. This is a consequence of the short interaction times between the system and environment which do not allow for the information of each E_i to scramble towards different S_j [104, 105].

3.1. Local versus global detailed balance

Up to now the discussion makes no reference to the structure of the system–environment interactions V_i. A particularly interesting situation, which is in practice the most widely studied case, is when the V_i satisfy local detailed balance with respect to each subsystem, but do not necessarily satisfy global detailed balance due the system–system interactions H_I. To make this argument more precise we introduce the idea of eigenoperators.

By using the method of eigenoperators, we can cast the conditions of detailed balance in terms solely of the algebra of operators. This therefore allows us to introduce the notion of local versus global detailed balance, as referring to the local or global Hamiltonian of the system. We note that here this terminology is being used in a slightly different context as, for instance, [88]. There the term is employed in the standard sense of 'detailed balance' [106].

Let H be an arbitrary Hamiltonian. An operator A is called an eigenoperator of H if it satisfies [H, A] = −ωA, for some ω > 0 (usually referred to as a Bohr, or transition, frequency). Due to this algebra, A and A^† function as lowering and raising operators for the H, causing transitions between energy levels E and E' separated by an energy ω.

Returning to our problem, we shall assume that the V_i have the form

$$\begin{eqnarray}&&{V}_{i}=\displaystyle \sum _{k}{g}_{i,k}({L}_{i,k}^{\dagger }{A}_{i,k}+{L}_{i,k}{A}_{i,k}^{\dagger }),\end{eqnarray} \tag{ 31 }$$

where g_i,k are constants, while L_i,k and A_i,k are eigenoperators of ${H}_{{S}_{i}}$ and ${H}_{{E}_{i}}$, respectively. That is:

$$\begin{eqnarray}&&[{H}_{{S}_{i}},{L}_{i,k}]=-{\omega }_{i,k}{L}_{i,k},\end{eqnarray} \tag{ 32 }$$

$$\begin{eqnarray}&&[{H}_{{E}_{i}},{A}_{i,k}]=-{\omega }_{i,k}{A}_{i,k},\end{eqnarray} \tag{ 33 }$$

where, for each S_i, {ω_i,k} represents the set of transition frequencies of the Hamiltonian ${H}_{{S}_{i}}$. The rationale behind these expressions is the following. For each subsystem S_i, the set of Bohr frequencies {ω_i,k} can be viewed as the set of transitions which are activated in that subsystem, due to the contact with its local bath. For instance, in the case of a harmonic oscillator, a dissipator such as (5) is characterized by an eigenoperator a, which only induces transitions between neighboring levels. Hence, the only Bohr frequency would be ω (the natural frequency of the oscillator). Similarly, a system–environment interaction containing a² would induce transitions with Bohr frequencies 2ω.

In an interaction such as (31), however, there is also the additional assumption that whenever there is a transition of +ω_i,k in the system, the corresponding transition in E_i occurs with energy −ω_i,k Hence, it follows that

$$\begin{eqnarray}&&[{H}_{{S}_{i}}+{H}_{{E}_{i}},{V}_{i}]=0,\end{eqnarray} \tag{ 34 }$$

meaning that all the energy that leaves system S_i enters bath E_i. This equation, together with the assumption that the bath is thermal, implies local detailed balance.

However, a similar relation does not hold for the total system Hamiltonian H_S, due to the interaction H_I between subsystems. The reason is that, in general, L_i,k are not eigenoperators of H_I. Hence

$$\begin{eqnarray}&&[{H}_{S}+{H}_{E},V]=[{H}_{I},V]\ne 0,\end{eqnarray} \tag{ 35 }$$

where ${H}_{E}={\sum }_{i=1}^{N}{H}_{{E}_{i}}$ and $V={\sum }_{i=1}^{N}{V}_{i}$. This means that even though detailed balance may hold locally, it is violated globally. As will be shown below, taking into account this fact is essential for providing a consistent thermodynamic description of LMEs.

Substituting equation (31) in equation (30) one may obtain a more explicit formula for the LME dissipators. To do so one notices that since the environments are in equilibrium, it follows that $\langle {A}_{i,k}{A}_{i,q}\rangle =0$ and $\langle {A}_{i,k}^{\dagger }{A}_{i,q}\rangle \propto {\delta }_{k,q}$. Hence, equation (30) reduces to

$$\begin{eqnarray}&&{D}_{i}({\rho }_{S})=\displaystyle \sum _{k}{\gamma }_{i,k}^{-}\ { \mathcal L }[{L}_{i,k},{\rho }_{S}]+{\gamma }_{i,k}^{+}\ { \mathcal L }[{L}_{i,k}^{\dagger },{\rho }_{S}],\end{eqnarray} \tag{ 36 }$$

where ${ \mathcal L }[A,\rho ]$ is the Lindblad dissipator defined in equation (7), whereas

$$\begin{eqnarray}\begin{array}{rcl}{\gamma }_{i,k}^{-} & = & {g}_{i,k}^{2}\ \langle {A}_{i,k}{A}_{i,k}^{\dagger }\rangle ,\\ {\gamma }_{i,k}^{+} & = & {g}_{i,k}^{2}\ \langle {A}_{i,k}^{\dagger }{A}_{i,k}\rangle .\end{array}\end{eqnarray} \tag{ 37 }$$

Since the A_i,k are eigenoperators of the bath which is in equilibrium at inverse temperature β_i, it follows that

$$\begin{eqnarray}&&\displaystyle \frac{{\gamma }_{i,k}^{+}}{{\gamma }_{i,k}^{-}}={{\rm{e}}}^{-{\beta }_{i}{\omega }_{i,k}},\end{eqnarray} \tag{ 38 }$$

which is another manifestation of local detailed balance.

3.2. Example: two harmonic oscillators

To show how these results may be applied, we consider again the two harmonic oscillators a₁ and a₂ treated in section 2. We assume that the environmental units E₁ and E₂ are themselves harmonic oscillators with operators b₁ and b₂⁸ , and we take the total Hamiltonian to be

$$\begin{eqnarray}&&{H}_{\mathrm{tot}}={H}_{S}+\displaystyle \sum _{i=1,2}\{{\omega }_{i}{b}_{i}^{\dagger }{b}_{i}+g({a}_{i}^{\dagger }{b}_{i}+{b}_{i}^{\dagger }{a}_{i})\},\end{eqnarray} \tag{ 39 }$$

where H_S is given in equation (1). Since the frequency ω₁ of S₁ and E₁ are the same, the operators a₁ and b₁ are eigenoperators of ${H}_{1}={\omega }_{1}{a}_{1}^{\dagger }{a}_{1}$ and ${H}_{{E}_{1}}={\omega }_{1}{b}_{1}^{\dagger }{b}_{1}$. Hence, local detailed balance, equation (34), is satisfied. However, a₁ is not an eigenoperator of H_S due to the interaction term ${H}_{I}=\epsilon ({a}_{1}^{\dagger }{a}_{2}+{a}_{1}{a}_{2}^{\dagger })$ (equation (3)) so that global detailed balance is in general broken.

Applying equation (36) with ${L}_{i,k}={a}_{i}$ and A_i,k = b_i we then immediately find the master equation (4) with the transition rates ${\gamma }_{i,k}^{\pm }$ given by

$$\begin{eqnarray}&&{\gamma }_{i}^{+}={g}^{2}\langle {b}_{i}^{\dagger }{b}_{i}\rangle := \gamma {n}_{i},\qquad {\gamma }_{i}^{-}={g}^{2}\langle {b}_{i}{b}_{i}^{\dagger }\rangle =\gamma ({n}_{i}+1),\end{eqnarray} \tag{ 40 }$$

which are the coefficients in equation (5).

3.3. Thermodynamics of the repeated interactions method

We now turn to a description of the thermodynamics of the repeated interactions method. Since the global S + E interaction is unitary, we may study the changes in the energy of the system and environment individually for a given stroke. The heat δQ_i exchanged between the system and environment i in one stroke must then be simply $\langle {H}_{{E}_{i}}{\rangle }_{0}-\langle {H}_{{E}_{i}}{\rangle }_{\tau }$. Dividing by τ and taking the limit τ → 0 will give the heat current ${\dot{Q}}_{i}$ to environment E_i. The explicit calculation is postponed to appendix B. As a result, we find the surprisingly simple formula:

$$\begin{eqnarray}&&{\dot{Q}}_{i}=\displaystyle \sum _{k}{\omega }_{i,k}\{{\gamma }_{i,k}^{+}\langle {L}_{i,k}{L}_{i,k}^{\dagger }\rangle -{\gamma }_{i,k}^{-}\langle {L}_{i,k}^{\dagger }{L}_{i,k}\rangle \}.\end{eqnarray} \tag{ 41 }$$

For example, in the case of the two harmonic oscillators, equation (41) gives precisely the heat rate formula (16) that was used in section 2.

In addition to the heat rates, however, there will in general also be a work contribution. The simplest way of seeing this is to note that for each interaction stroke one must turn the system–environment coupling on and off, so that the total Hamiltonian (26) must actually be time dependent. Since the global S − E system is isolated, the work in an individual stroke between times nτ and (n + 1)τ may be unambiguously defined as

$$\begin{eqnarray}&&\delta {W}_{\mathrm{ext}}={\int }_{n\tau }^{(n+1)\tau }\left\langle \displaystyle \frac{\partial {H}_{\mathrm{tot}}}{\partial t}\right\rangle {\rm{d}}t.\end{eqnarray} \tag{ 42 }$$

Again, dividing by τ and taking the limit τ → 0 yields the work rate ${\dot{W}}_{\mathrm{ext}}$. As shown in appendix B, carrying out this computation yields

$$\begin{eqnarray}&&{\dot{W}}_{\mathrm{ext}}=\displaystyle \frac{1}{2}\displaystyle \sum _{i,k}\{{\gamma }_{i,k}^{-}\langle {L}_{i,k}^{\dagger }{F}_{i,k}\rangle -{\gamma }_{i,k}^{+}\langle {F}_{i,k}{L}_{i,k}^{\dagger }\rangle \}+{\rm{c}}.{\rm{c}}.,\end{eqnarray} \tag{ 43 }$$

where F_i,k = [H_I, L_i,k] and c.c. stands for complex conjugate. Hence, we see that the work rate is related precisely to the non-commutativity of the jump operators L_i,k with the system–system interaction H_I. If global detailed balance, equation (35), is recovered, then the work cost is identically zero. Otherwise, whenever global detailed balance is broken, there will be an associated work cost. As an example, in the case of the two QHOs, [H_I, a₁] = −a₂, so that a direct application of equation (43) leads to expression (17).

From a practical point of view, equations (41) and (43) are our main results, as they offer general expressions that may be applied over a broad range of situations. These results generalize Barra's findings [88] to arbitrary bath and system's structure. In appendix C, we calculate the entropy production showing that it is indeed always positive in agreement with the second law. We now return to the problem of coupled harmonic oscillators and show how they can be applied to the study of more complex situations.

We consider once again the problem studied in section 2, but we now make the following generalizations. First, we assume to have a one-dimensional chain of N, instead of 2, oscillators, with the first and last coupled to local baths. Second, we consider a more general type of interaction between them, so that we take the system Hamiltonian to be

$$\begin{eqnarray}&&{H}_{S}=\displaystyle \sum _{i=1}^{N}{\omega }_{i}{a}_{i}^{\dagger }{a}_{i}+\displaystyle \sum _{i=1}^{N-1}\{\epsilon ({a}_{i}^{\dagger }{a}_{i+1}+{a}_{i+1}^{\dagger }{a}_{i})+\eta ({a}_{i}^{\dagger }{a}_{i+1}^{\dagger }+{a}_{i+1}{a}_{i})\}.\end{eqnarray} \tag{ 44 }$$

When η = 0 we recover the excitation-conserving tight-binding model. Conversely, η =  corresponds to a position–position coupling x_ix_i+1. The total system evolves according to the master equation

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}{\rho }_{S}}{{\rm{d}}t}=-{\rm{i}}[{H}_{S},{\rho }_{S}]+{D}_{1}({\rho }_{S})+{D}_{N}({\rho }_{S}),\end{eqnarray} \tag{ 45 }$$

where D_i(ρ_S) is given in Eq (5).

Using equation (41) we find that the heat rate to the first and last environments will be

$$\begin{eqnarray}&&{\dot{Q}}_{i}=\gamma {\omega }_{i}({n}_{i}-\langle {a}_{i}^{\dagger }{a}_{i}\rangle ),\qquad i=1,N.\end{eqnarray} \tag{ 46 }$$

The work rate, on the other hand, is found using equation (43) and reads

$$\begin{eqnarray}\begin{array}{rcl}\dot{W} & = & \displaystyle \frac{\gamma \epsilon }{2}\langle {a}_{1}^{\dagger }{a}_{2}+{a}_{1}{a}_{2}^{\dagger }+{a}_{N}^{\dagger }{a}_{N-1}+{a}_{N}{a}_{N-1}^{\dagger }\rangle \\ & & +\displaystyle \frac{\gamma \eta }{2}\langle {a}_{1}{a}_{2}+{a}_{1}{a}_{2}+{a}_{N}{a}_{N-1}+{a}_{N}{a}_{N-1}\rangle .\end{array}\end{eqnarray} \tag{ 47 }$$

If N = 2 and η = 0 this reduces to equation (17).

4.1. Degradation of the engine's operation due to counter-rotating terms

As a first application, we return to the case of N = 2 oscillators but generalize the results of section 2 to include the counter-rotating terms η. We then show that, interpreting our steady-state as a heat engine, as in figure 2, the presence of this term only degrades the machine's operation. We once again use the Lyapunov equation technique of appendix A. Although the steady-state may be trivially found numerically, it turns out that the analytical expression for the covariance matrix becomes extremely cumbersome when $\eta \ne 0$. For the purpose of illustration, we therefore present here results that are valid when the bath coupling γ is much larger than all other energy scales in the problem. In this case we find the following simple expressions for the heat rates (46) and the work rate (47):

$$\begin{eqnarray}&&{\dot{Q}}_{1}=\displaystyle \frac{2{\omega }_{1}}{\gamma }[{\epsilon }^{2}({n}_{1}-{n}_{2})-{\eta }^{2}({n}_{1}+{n}_{2}+1)],\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{\dot{Q}}_{2}=-\displaystyle \frac{2{\omega }_{2}}{\gamma }[{\epsilon }^{2}({n}_{1}-{n}_{2})+{\eta }^{2}({n}_{1}+{n}_{2}+1)],\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\dot{W}}_{\mathrm{ext}} & = & -\displaystyle \frac{2}{\gamma }[{\epsilon }^{2}({n}_{1}-{n}_{2})({\omega }_{1}-{\omega }_{2})+\\ & & \ -{\eta }^{2}({\omega }_{1}+{\omega }_{2})({n}_{1}+{n}_{2}+1)].\end{array}\end{eqnarray} \tag{ 50 }$$

If η = 0 we recover equations (18)–(20), provided we also approximate Δ ≃ γ.

The influence of η in these thermodynamic quantities is shown in figure 3, which compares the results with those of figure 2 corresponding to η = 0. As can be seen, the region in which the machine operates as a refrigerator or an engine is severely reduced by the presence of η. Moreover, as can be easily verified from equations (48)–(50), in the case η =  of a position–position coupling all of these operation regimes become inaccessible. This type of behavior has already been observed experimentally in the case of opto-mechanical systems, in which the interplay between η and can be tuned by appropriately choosing the pump's sideband. And indeed, what is observed, is that since the η interaction tends to entangle the two oscillators (since it leads to an evolution described by a two-mode squeezing operator), it has the general tendency of heating them up [91, 107, 108], which is precisely what is observed in figure 3.

Zoom In Zoom Out Reset image size

4.2. Analysis for a chain of $N$ oscillators

We next carry on a numerical analysis for the case of a chain of N oscillators, assuming for simplicity that η = 0. The steady-state in this case is obtained by solving the Lyapunov equation numerically (appendix A). A set of results is presented in figure 4 for the case of a linear function profile

$$\begin{eqnarray}&&{\omega }_{i}=\displaystyle \frac{(N-i){\omega }_{1}+(i-1){\omega }_{N}}{N-1}.\end{eqnarray} \tag{ 51 }$$

Zoom In Zoom Out Reset image size

As can be seen, in general the thermodynamic quantities depend on the system size. However, the ratios between heat and work are independent of size and only depend on the frequencies ω₁ and ω_N of the initial and final oscillators. This can be proven rigorously after writing the continuity equation for the number of excitations of each oscillator using techniques similar to those employed in [57, 58]:

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {a}_{1}^{\dagger }{a}_{1}\rangle }{{\rm{d}}{t}}=-{\rm{i}}\epsilon \langle {a}_{1}^{\dagger }{a}_{2}-{a}_{1}{a}_{2}^{\dagger }\rangle +\gamma ({n}_{1}-\langle {a}_{1}^{\dagger }{a}_{1}\rangle ),\end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {a}_{{\rm{i}}}^{\dagger }{a}_{{\rm{i}}}\rangle }{{\rm{d}}{t}}=-{\rm{i}}\epsilon \langle {a}_{{\rm{i}}}^{\dagger }{a}_{{\rm{i}}+1}-{a}_{{\rm{i}}}{a}_{{\rm{i}}+1}^{\dagger }\rangle ,\quad {i}=2,\,\ldots ,\,N-1\end{eqnarray} \tag{ 53 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}\langle {a}_{N}^{\dagger }{a}_{N}\rangle }{{\rm{d}}{t}}={\rm{i}}\epsilon \langle {a}_{N-1}^{\dagger }{a}_{N}-{a}_{N-1}{a}_{N}^{\dagger }\rangle +\gamma ({n}_{N}-\langle {a}_{N}^{\dagger }{a}_{N}\rangle ).\end{eqnarray} \tag{ 54 }$$

Hence, at steady-state we have:

$$\begin{eqnarray}&&\gamma ({n}_{1}-\langle {a}_{1}^{\dagger }{a}_{1}\rangle )=-\gamma ({n}_{N}-\langle {a}_{N}^{\dagger }{a}_{N}\rangle )\end{eqnarray} \tag{ 55 }$$

and replacing this in the expressions for the heat flows (46) we obtain:

$$\begin{eqnarray}&&{\dot{Q}}_{1}=\gamma {\omega }_{1}({n}_{1}-\langle {a}_{1}^{\dagger }{a}_{1}\rangle )=-\displaystyle \frac{{\omega }_{1}}{{\omega }_{N}}{\dot{Q}}_{N}.\end{eqnarray} \tag{ 56 }$$

Finally, using the first law of thermodynamics: ${\dot{W}}_{\mathrm{ext}}+{\dot{Q}}_{1}+{\dot{Q}}_{N}=0$ we obtain:

$$\begin{eqnarray}&&\displaystyle \frac{{\dot{Q}}_{1}}{{\dot{W}}_{{\rm{ext}}}}=\displaystyle \frac{{\omega }_{1}}{{\omega }_{N}-{\omega }_{1}};\quad \displaystyle \frac{{\dot{W}}_{{\rm{ext}}}}{{\dot{Q}}_{N}}=-1+\displaystyle \frac{{\omega }_{1}}{{\omega }_{N}}.\end{eqnarray} \tag{ 57 }$$

We remark that these ratios coincide with those of Otto thermal machines and refrigerators operating with a single harmonic oscillator driven between frequencies ω₁ and ω_N. Moreover, the universal ratios (57) would hold if we replace the linear chain with a harmonic lattice of arbitrary geometry or replace the QHO with coupled qubits for which a continuity equation equivalent to equations (52)–(54).

In this paper we have systematically analyzed the external work required for LMEs based on the repeated interaction model. Specifically, we have provided explicit expressions, equations (41) and (43), for the relevant thermodynamics quantities as a function of the reservoir jump operators for a generic master equation in Lindblad form (with time independent rates). We have found that, at steady-state, this external work is directly related to the breaking of global detailed balance stemming from the internal system–system interactions. Our analysis is general and can be applied to any array of interacting d-level systems.

For the purpose of illustration, we have chosen to provide a detailed analysis of a chain of harmonic oscillators described by LMEs, a problem which has attracted considerable attention from the quantum thermodynamics and open quantum systems communities. We find that the heat transport through the chain is dramatically affected by the boundary driving and can even be inverted, leading to a flow of heat from the cold to the hot bath. However, this does not violate any thermodynamic law as this regime requires introducing positive work into the system, the latter effectively realizing a quantum refrigerator.

Our findings therefore provide clear evidence that in order to be consistent with thermodynamics, one must take into consideration the microscopic model associated with a particular dynamical equation. Moreover, for a given process, complete positivity ensures that a quantum microscopic model can always be found. Hence, by taking into account all sources of energy exchange, we prove that the thermodynamic laws will always be satisfied.

We thank D Alonso, F Barra, M Campisi, J P Paz for discussions. GTL acknowledges the support from the São Paulo Research Foundation, under grant number 2016/08721-7. GDC and MA acknowledge the KITP program 'Thermodynamics of quantum systems: Measurement, engines, and control' 2018, where part of this work has been done. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958. AJR acknowledges support from CONICET, ANPCyT and UBACyT. AF acknowledges support from the EPSRC project EP/P00282X/1.

We discuss here how to compute the steady-state of Gaussian preserving master equations. For concreteness, we focus on the problem discussed in section 2, more specifically equation (4). The extension to multiple oscillators (section 4) is straightforward. Since the master equation is Gaussian preserving, the steady-state must necessarily be Gaussian. Hence, it is completely determined by the second moments of the oscillators' positions and momenta, arranged for convenience in the vector:

$$\begin{eqnarray}&&{Y}^{T}=({x}_{1},{p}_{1},{x}_{2},{p}_{2}).\end{eqnarray} \tag{ A1 }$$

We introduce the drift matrix containing the couplings and decay constants:

$$\begin{eqnarray}K=\left(\begin{array}{cccc}-\gamma /2 & 1 & 0 & \lambda \\ -{\omega }_{1}^{2} & -\gamma /2 & -\mu & 0\\ 0 & \lambda & -\gamma /2 & 1\\ -\mu & 0 & -{\omega }_{2}^{2} & -\gamma /2\end{array}\right),\end{eqnarray} \tag{ A2 }$$

where:

$$\begin{eqnarray}&&\lambda =\displaystyle \frac{2\epsilon }{\sqrt{{\omega }_{1}{\omega }_{2}}},\quad \quad \mu =2\epsilon \sqrt{{\omega }_{1}{\omega }_{2}}\end{eqnarray} \tag{ A3 }$$

and the system's covariance matrix V with entries:

$$\begin{eqnarray}&&V(i,j)=\displaystyle \frac{1}{2}\langle \ \{Y(i),Y(j)\}\ \rangle \end{eqnarray} \tag{ A4 }$$

the average being done on the state of the system. By calculating the equation of motion for all the elements of the covariance matrix, we obtain the following Lyapunov equation for the covariance matrix:

$$\begin{eqnarray}&&\dot{V}={KV}+{{VK}}^{T}+D,\end{eqnarray} \tag{ A5 }$$

where

$$\begin{eqnarray*}\begin{array}{rcl}D & = & \mathrm{diag}\,\left[\displaystyle \frac{\gamma (2{n}_{1}+1)}{2{\omega }_{1}},\displaystyle \frac{\gamma {\omega }_{1}(2{n}_{1}+1)}{2},\right.\\ & & \times \ \left.\displaystyle \frac{\gamma (2{n}_{2}+1)}{2{\omega }_{2}},\displaystyle \frac{\gamma {\omega }_{2}(2{n}_{2}+1)}{2}\right].\end{array}\end{eqnarray*}$$

The steady-state is then determined by the algebraic equation

$$\begin{eqnarray}&&{KV}+{{VK}}^{T}+D=0.\end{eqnarray} \tag{ A6 }$$

In absence of coupling between the first and second oscillator,  = 0, the stationary state of each oscillator is thermal at its respective temperature:

$$\begin{eqnarray}V=\left(\begin{array}{cccc}\displaystyle \frac{2{n}_{1}+1}{2{\omega }_{1}} & 0 & 0 & 0\\ 0 & \displaystyle \frac{(2{n}_{1}+1){\omega }_{1}}{2} & 0 & 0\\ 0 & 0 & \displaystyle \frac{2{n}_{2}+1}{2{\omega }_{2}} & 0\\ 0 & 0 & 0 & \displaystyle \frac{(2{n}_{2}+1){\omega }_{2}}{2}\end{array}\right).\end{eqnarray} \tag{ A7 }$$

When $\epsilon \ne 0$, one finds instead the following expression for the entries of V (recall that V_ij = V_ji):

$$\begin{eqnarray}&&{V}_{11}=\displaystyle \frac{1}{2{\omega }_{1}}\left(\displaystyle \frac{4{\epsilon }^{2}({n}_{2}-{n}_{1})}{{{\rm{\Delta }}}^{2}}+2{n}_{1}+1\right),\end{eqnarray} \tag{ A8 }$$

$$\begin{eqnarray}&&{V}_{22}=\displaystyle \frac{1}{2}{\omega }_{1}\left(\displaystyle \frac{4{\epsilon }^{2}({n}_{2}-{n}_{1})}{{{\rm{\Delta }}}^{2}}+2{n}_{1}+1\right),\end{eqnarray} \tag{ A9 }$$

$$\begin{eqnarray}&&{V}_{33}=\displaystyle \frac{1}{2{\omega }_{2}}\left(\displaystyle \frac{4{\epsilon }^{2}({n}_{1}-{n}_{2})}{{{\rm{\Delta }}}^{2}}+2{n}_{2}+1\right),\end{eqnarray} \tag{ A10 }$$

$$\begin{eqnarray}&&{V}_{44}=\displaystyle \frac{1}{2}{\omega }_{2}\left(\displaystyle \frac{4{\epsilon }^{2}({n}_{1}-{n}_{2})}{{{\rm{\Delta }}}^{2}}+2{n}_{2}+1\right),\end{eqnarray} \tag{ A11 }$$

$$\begin{eqnarray}&&{V}_{13}=\displaystyle \frac{\epsilon ({n}_{1}-{n}_{2})({\omega }_{1}-{\omega }_{2})}{\sqrt{{\omega }_{1}{\omega }_{2}}{{\rm{\Delta }}}^{2}},\end{eqnarray} \tag{ A12 }$$

$$\begin{eqnarray}&&{V}_{14}=\displaystyle \frac{\gamma {\omega }_{2}\epsilon ({n}_{2}-{n}_{1})}{\sqrt{{\omega }_{1}{\omega }_{2}}{{\rm{\Delta }}}^{2}},\end{eqnarray} \tag{ A13 }$$

$$\begin{eqnarray}&&{V}_{23}=\displaystyle \frac{\gamma {\omega }_{1}\epsilon ({n}_{1}-{n}_{2})}{\sqrt{{\omega }_{1}{\omega }_{2}}{{\rm{\Delta }}}^{2}},\end{eqnarray} \tag{ A14 }$$

$$\begin{eqnarray}&&{V}_{24}=\displaystyle \frac{\epsilon ({n}_{1}-{n}_{2})({\omega }_{1}-{\omega }_{2})\sqrt{{\omega }_{1}{\omega }_{2}}}{{{\rm{\Delta }}}^{2}},\end{eqnarray} \tag{ A15 }$$

$$\begin{eqnarray}&&{V}_{12}={V}_{34}=0,\end{eqnarray} \tag{ A16 }$$

which are tantamount to equations (8)–(9).

In this section we provide additional details on how to establish the thermodynamic quantities within the method of repeated interactions. A single interaction stroke between system and environment is described by the unitary map

$$\begin{eqnarray*}&&{\rho }_{{SE}}^{\prime} ={{\rm{e}}}^{-{\rm{i}}\tau {H}_{\mathrm{tot}}}\ {\rho }_{S}{\rho }_{E}\ {{\rm{e}}}^{{\rm{i}}\tau {H}_{\mathrm{tot}}},\end{eqnarray*}$$

where H_tot is given in equation (26) (with ${V}_{i}\to {V}_{i}/\sqrt{\tau }$, as discussed in the main text). Hence, the evolution of any observable ${ \mathcal O }$ (from either the system or from the environments) due to this interaction is given by

$$\begin{eqnarray}&&\langle { \mathcal O }\rangle ^{\prime} =\langle {{\rm{e}}}^{{\rm{i}}\tau {H}_{{\rm{tot}}}}\ { \mathcal O }{{\rm{e}}}^{-{\rm{i}}\tau {H}_{{\rm{tot}}}}\rangle ,\end{eqnarray} \tag{ B1 }$$

where the unprimed average in the right-hand side is computed with respect to the initial state before the evolution. Using this idea we find that the change in the energy of the system H_S and of each environment ${H}_{{E}_{i}}$ is up to order τ

$$\begin{eqnarray}&&{\rm{\Delta }}{H}_{S}=-\displaystyle \frac{\tau }{2}\displaystyle \sum _{i=1}^{N}\langle [{V}_{i},[{V}_{i},{H}_{S}]]\rangle ,\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&{\rm{\Delta }}{H}_{{E}_{i}}=-\displaystyle \frac{\tau }{2}\langle [{V}_{i},[{V}_{i},{H}_{{E}_{i}}]]\rangle .\end{eqnarray} \tag{ B3 }$$

Carrying out the computations, using equation (31), we find

$$\begin{eqnarray*}\begin{array}{rcl}{\rm{\Delta }}{H}_{S} & = & \displaystyle \frac{\tau }{2}\displaystyle \sum _{i,k}{g}_{i,k}^{2}\{\langle {A}_{i,k}^{\dagger }{A}_{i,k}\rangle \langle {L}_{i,k}[{H}_{S},{L}_{i,k}^{\dagger }]+[{L}_{i,k},{H}_{S}]{L}_{i,k}^{\dagger }\rangle \\ & & +\langle {A}_{i,k}{A}_{i,k}^{\dagger }\rangle \langle {L}_{i,k}^{\dagger }[{H}_{S},{L}_{i,k}]+[{L}_{i,k}^{\dagger },{H}_{S}]{L}_{i,k}\rangle \}\\ {\rm{\Delta }}{H}_{{E}_{i}} & = & -\tau \displaystyle \sum _{k}{g}_{i,k}^{2}{\omega }_{i,k}\{\langle {A}_{i,k}^{\dagger }{A}_{i,k}\rangle \langle {L}_{i,k}{L}_{i,k}^{\dagger }\rangle -\langle {A}_{i,k}{A}_{i,k}^{\dagger }\rangle \langle {L}_{i,k}^{\dagger }{L}_{i,k}\rangle \}.\end{array}\end{eqnarray*}$$

Here the averages over the environment operators are always computed with respect to the same state ρ_E, whereas the average over system operators are computed over the instantaneous state ρ_S(t).

Thus, identifying the transition rates ${\gamma }_{i,k}^{\pm }$ in equation (37), we may then divide both sides by τ and take the limit τ → 0, which yields

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{{\rm{d}}\langle {H}_{S}\rangle }{{\rm{d}}t} & = & \displaystyle \frac{1}{2}\displaystyle \sum _{i,k}\{{\gamma }_{i,k}^{+}\langle {L}_{i,k}[{H}_{S},{L}_{i,k}^{\dagger }]+[{L}_{i,k},{H}_{S}]{L}_{i,k}^{\dagger }\rangle \\ & & +{\gamma }_{i,k}^{-}\langle {L}_{i,k}^{\dagger }[{H}_{S},{L}_{i,k}]+[{L}_{i,k}^{\dagger },{H}_{S}]{L}_{i,k}\rangle \}.\end{array}\end{eqnarray} \tag{ B4 }$$

Moreover, with respect to the change in energy of the environments, we define the heat rates as ${\dot{Q}}_{i}=-{\rm{\Delta }}{H}_{{E}_{i}}/\tau$, which then yields precisely equation (41).

The heat rate may now be found in two ways. First, using equation (25) and noting that

$$\begin{eqnarray}&&[{H}_{S},{L}_{i,k}]=-{\omega }_{i,k}{L}_{i,k}+[{H}_{I},{L}_{i,k}],\end{eqnarray} \tag{ B5 }$$

we see that ${\rm{d}}\langle {H}_{S}\rangle /{\rm{d}}t$ can be split into two terms, one of which is precisely ${\sum }_{i}{\dot{Q}}_{i}$. Hence, the remainder must be attributed to the work that has to be performed, which yields precisely equation (43).

Instead, another way of defining the work rate is by noting that since one must decouple the system from the environments at each stroke, the total Hamiltonian is actually time dependent. That is, if we focus on just a single interaction stroke, then instead of equation (26), the Hamiltonian for this interaction is more appropriately defined as

$$\begin{eqnarray}&&{H}_{\mathrm{tot}}={H}_{S}+{H}_{E}+\displaystyle \frac{\lambda (t)}{\sqrt{\tau }}\displaystyle \sum _{i=1}^{N}{V}_{i},\end{eqnarray} \tag{ B6 }$$

where λ(t) has the value 1 when t ∈ [nτ, (n + 1)τ] and zero otherwise. Since the global S + E interaction is unitary, work can be unambiguously defined as in equation (42). Carrying out the integration one then finds

$$\begin{eqnarray}&&\delta W=\displaystyle \sum _{i=1}^{N}\displaystyle \frac{\langle {V}_{i}\rangle -\langle {V}_{i}\rangle ^{\prime} }{\sqrt{\tau }},\end{eqnarray} \tag{ B7 }$$

where, once again, the primed expectation value refers to the state after the interaction. Using once again equation (B1) for evaluating this expectation value, we then find

$$\begin{eqnarray}&&\delta W=\displaystyle \frac{\tau }{2}\displaystyle \sum _{i=1}^{N}\langle [{V}_{i},[{V}_{i},{H}_{S}+{H}_{{E}_{i}}]]\rangle ={\rm{\Delta }}{H}_{S}+\displaystyle \sum _{i=1}^{N}{\rm{\Delta }}{H}_{{E}_{i}}.\end{eqnarray} \tag{ B8 }$$

Hence, taking again the limit τ → 0, one identifies the work rate as $\dot{W}=\tfrac{{\rm{d}}\langle {H}_{S}\rangle }{{\rm{d}}t}-{\sum }_{i}{\dot{Q}}_{i}$.

In this appendix we discuss how to express the second law of thermodynamics within the repeated interactions scheme. Since the interaction strokes only last for a small time τ, one may neglect any potential bath–bath correlations that may appear when a system is interacting with multiple environments. Hence, following [89], we can define the entropy production in a single stroke as

$$\begin{eqnarray}&&{\rm{\Sigma }}:= { \mathcal I }(S:E),+S({\rho }_{E}(\tau )| | {\rho }_{E}^{\mathrm{th}})\geqslant 0,\end{eqnarray} \tag{ C1 }$$

where

$$\begin{eqnarray}&&{ \mathcal I }(S:E)={\rm{\Delta }}{S}_{S}+{\rm{\Delta }}{S}_{E}\simeq {\rm{\Delta }}{S}_{S}+\displaystyle \sum _{i=1}^{N}{\rm{\Delta }}{S}_{E,i},\end{eqnarray} \tag{ C2 }$$

is the mutual information developed between the system and all the environments, with ΔS_S being the change in the entropy of the system and ${\rm{\Delta }}{S}_{{E}_{i}}$ the change in the von Neumann entropy of environment E_i. The second term in equation (C1), on the other hand

$$\begin{eqnarray}&&S({\rho }_{E}(\tau )| | {\rho }_{E}^{\mathrm{th}})\simeq \displaystyle \sum _{i=1}^{N}S({\rho }_{E,i}(\tau )| | {\rho }_{E,i}^{\mathrm{th}}),\end{eqnarray} \tag{ C3 }$$

is the quantum relative entropy between the final and initial states of each environment (here $S(\rho | | \sigma )$ $=\,\mathrm{tr}(\rho {\rm{ln}}\rho -\rho {\rm{ln}}\sigma )$). The positivity of Σ then follows immediately from the positivity of the mutual information and the quantum relative entropy.

Using equation (B8) one may then readily show that equation (C1) may be written as

$$\begin{eqnarray}&&{\rm{\Sigma }}={\rm{\Delta }}{S}_{S}-\displaystyle \sum _{i=1}^{N}{\beta }_{i}\delta {Q}_{i}.\end{eqnarray} \tag{ C4 }$$

Dividing by τ and taking the limit τ → 0 we may then obtain the usual expression for the entropy production rate

$$\begin{eqnarray}&&{\rm{\Pi }}=\displaystyle \frac{{\rm{d}}{S}_{S}}{{\rm{d}}t}-\displaystyle \sum _{i=1}^{N}{\beta }_{i}\displaystyle \frac{{\rm{d}}{Q}_{i}}{{\rm{d}}t},\end{eqnarray} \tag{ C5 }$$

with the last term representing the entropy flux rates to each environment.

In the particular case where there is a single environment, or when all environments have the same temperature, we may use the first law to write this as

$$\begin{eqnarray}&&{\rm{\Pi }}=\beta \left(\displaystyle \frac{{\rm{d}}W}{{\rm{d}}t}-\displaystyle \frac{{\rm{d}}F}{{\rm{d}}t}\right),\end{eqnarray} \tag{ C6 }$$

where ${F}_{S}=\langle {H}_{S}\rangle -{{TS}}_{S}$ is the free energy of the system. Thus, in this case we recover another well known expression for the entropy production rate.