Reconciliation of quantum local master equations with thermodynamics

The study of open quantum systems often relies on approximate master equations derived under the assumptions of weak coupling to the environment. However when the system is made of several interacting subsystems such a derivation is in many cases very hard. An alternative method, employed especially in the modelling of transport in mesoscopic systems, consists in using {\it local} master equations containing Lindblad operators acting locally only on the corresponding subsystem. It has been shown that this approach however generates inconsistencies with the laws of thermodynamics. In this paper we demonstrate that using a microscopic model of local master equations based on repeated collisions all thermodynamic inconsistencies can be resolved by correctly taking into account the breaking of global detailed balance related to the work cost of maintaining the collisions. We provide examples based on a chain of quantum harmonic oscillators whose ends are connected to thermal reservoirs at different temperatures. We prove that this system behaves precisely as a quantum heat engine or refrigerator, with properties that are fully consistent with basic thermodynamics.

The study of open quantum systems often relies on approximate master equations derived under the assumptions of weak coupling to the environment. However when the system is made of several interacting subsystems such a derivation is in many cases very hard. An alternative method, employed especially in the modelling of transport in mesoscopic systems, consists in using local master equations containing Lindblad operators acting locally only on the corresponding subsystem. It has been shown that this approach however generates inconsistencies with the laws of thermodynamics. In this paper we demonstrate that using a microscopic model of local master equations based on repeated collisions all thermodynamic inconsistencies can be resolved by correctly taking into account the breaking of global detailed balance related to the work cost of maintaining the collisions.
We provide examples based on a chain of quantum harmonic oscillators whose ends are connected to thermal reservoirs at different temperatures. We prove that this system behaves precisely as a quantum heat engine or refrigerator, with properties that are fully consistent with basic thermodynamics.
The starting point for addressing these phenomena is the theory of open quantum systems [42][43][44][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][47][48].
The situation becomes more involved when the system S is made up of several interacting subsystems S 1 , . . . , S n and each subsystem interacts with a local environment E i , possibly held at different temperatures T i (see Fig. 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][50][51][52][53]. However, these derivations are in general quite involved since they require knowledge of the full set of eigen-values 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.
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 local master equations (LME) (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][56][57][58][59][60][61][62][63][64][65][66][67][68][69][70].
It turns out, however, that the nonlocal terms neglected in the LME may still lead to non-thermal steady states [71] 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 [72]. 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 S1 (a) FIG. 1. We consider in this paper a system S composed of several sub-systems S 1 , . . . , S N (in the figure N = 2). Each sub-system S i is connected to a local environment E i prepared in a different temperature T i . (a) The standard bosonic heat bath model: the environment is assumed to consist of an ensemble of independent quantum harmonic oscillators with different frequencies in thermal equilibrium and coupled permanently to the system. (b) In this paper we focus instead on the framework of the repeated interactions method: the environment E i is divided into a series of ancillas (in this case represented by individual bosonic modes with identical frequencies) which interact with S i sequentially. This type of method leads to local master equations (LMEs), irrespective of the internal interactions between S i and S j . of these effects, as shown in [72], 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 local master equations [49,52,53,61,[73][74][75][76][77][78][79] as well as clever possible alternatives [80,81].
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,[82][83][84][85][86][87]. 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 Fig. 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 [85]. Hence, it is possible to put forth a consistent thermodynamic framework for the repeated interactions scheme, including definitions of heat, work and entropy production [85,86].
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 nonequilibrium 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 quantum harmonic oscillators (QHO), which has been the subject of considerable interest due to potential applications in trapped ions [10,11], optomechanical systems [8,88] 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 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 [89], which states that no au-tonomous 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 Sec. II 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 Sec. III 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 Sec. IV 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 Sec. V we summarise 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 Sec. III. In Appendix C we compute the entropy production showing that it is always positive.

II. STEADY-STATE OF TWO HARMONIC OSCILLATORS SUBJECT TO A LOCAL MASTER EQUATION
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 [72], consisting of two coupled harmonic oscillators, S 1 and S 2 , each coupled to its own local heat bath via a LME (see Fig. 1(a)). The Hamiltonian is taken to be where and 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, ω 1 = ω 2 . More general interactions will be discussed in Sec.

IV
The system is also subject to two baths E 1 and E 2 , which we assume can be modeled by LMEs acting only on S 1 and S 2 . Hence, we take the system to evolve according to the Lindblad master equation where Here γ is the dissipation rate, n i is the average excitation number for the Bose-Einstein distribution, with inverse temperature β i = (T i ) −1 and is the dissipator in Lindblad form. Henceforth, we always assume that T 1 < T 2 (see Fig. 1(a)). Below, it will also be useful to keep in mind that n i is a monotonically increasing function of T i /ω i .

A. 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: where (2) we then find that the energy of oscillator 1 will be where H 1 th = ω 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 ↔ 2. Similarly, the heat flux from the first oscillator to the second will bė which is consistent with the expression obtained in Ref. [72]. The relevant point to highlight from these results is that the sign of the 2 -terms in Eqs. (11) and (12) will depend on the relative magnitudes of n 1 and n 2 and not of T 1 and T 2 . This simple fact, when combined with the definition of n i in Eq. (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 1 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 1 = T 2 . And third, if T 1 < T 2 but we happen to satisfy then there will be a current flowing from the cold to the hot bath (Q 1→2 < 0). As pointed out in Ref. [72], these results seem to violate the second law of thermodynamics, a fact which the authors attribute to the inadequacy of local master equations 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 [89].

B. Resolving the thermodynamic inconsistencies
Using Eq. (4) one finds that the total rate of change of the system Hamiltonian (1) is given by However, this change in energy cannot be associated only with heat flowing to the reservoirs [85]. 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 Sec. III, we now discuss how one can remove all the inconsistencies by including this work cost.
In particular, we find that the correct splitting of Eq. (14) is is the heat rate flowing from the system to each reservoir (see Eq. (41) for the general expression) anḋ is the required work rate (see Eq. (43) for the general expression). From Eqs. (8)-(10) we then find that in the steady-statė These are very important results, which we shall now discuss in depth. But before doing so, some consistency checks are in order. First,Ẇ ext tends to zero as → 0, in which case global detailed balance is recovered. Second, the work required is proportional to γ 2 . As discussed in Ref. [52], the local master equation can be shown to be correct for small interactions up to the same order (see comment in Sec. II D in that paper). Third, in the steady-state d H S / dt = 0 so that we should haveẆ 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 Eqs. (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 ofQ 1 ,Q 2 andẆ is presented in Table I. We also show in Fig. 2 a plot of these 3 quantities illustrating the typical operation regimes.
We begin with the case of a refrigerator,Q 1 > 0. According to Eq. (18), this requires n 1 > n 2 which is tantamount to the condition (13). But we now see that when this condition is satisfied, one also hasẆ ext > 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 coefficient of performance (COP) of the refrigerator is the ratio of the heat extracted from the cold environment divided by the work invested: which coincides with the COP of the Otto cycle. As will become clear in Sec. III, 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 Sec. IV, the introduction of counter-rotating terms only reduces the COP below the Otto value.
Increasing ω 1 in the interval I. The three modes of operation of the steady-state of the two oscillators, always assuming E 1 is the cold reservoir (T 1 < T 2 ). In this paper positive heat or work always means energy entering the system. SoQ i > 0 means energy entered S through E i andẆ ext > 0 means work was performed on the system by an external agent.
Uses hot bath to produce useful work, dumping the remainder in the cold bath.

FIG. 2. The three regimes of operation in
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: which is again the Otto efficiency for a thermal machine. Finally, for ω 2 < ω 1 we necessarily have n 1 < n 2 and thereforė Q 1 < 0,Q 2 > 0 andẆ ext > 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: Hence, it functions as an accelerator (i.e. an oven) [90], heating the cold reservoir faster than it would with spontaneous thermal conduction. Notice that at the special value ω 1 /ω 2 = T 1 /T 2 , the system operates as a thermal machine with the Carnot efficiency and zero power.
We conclude this section by emphasising that within the repeated interactions framework one is able to resolve all of the seeming thermodynamic inconsistencies of local master equations. 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, we do not violate the theorem by Martinez and Paz [89], 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.

III. THE METHOD OF REPEATED INTERACTIONS
We now turn to the microscopic derivation of LMEs using the method of repeated interactions [55,58,[82][83][84][85][86][87][91][92][93][94][95][96][97][98][99], 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 modelled by a Lindblad master equation that is, in general, local. We emphasise that a larger class of non-Markovian master equations including the global master equation can be engineered in a similar way [84] and realised experimentally as proposed in Ref. [100].
We consider here this method implemented in the general scenario of Fig. 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 1 , . . . , S N , where each S i interacts with its own environment E i . The Hamiltonian of the system is taken to be where H S i is the Hamiltonian of subsystem S i and H I summarises all interactions between subsystems. Moreover, the total Hamiltonian comprising the system, the environments and all interactions, is taken as 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,86]. Following the usual repeated interactions approach, the reduced density matrix of the system ρ S (t) will then evolve according to the map where is the thermal state of the environments. Expanding Eq. (27) in a power series in τ one then finds, up to first order, the LME, where H S is given in (25) and 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 [101,102].

A. Local vs. 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.
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 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: [ 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 2 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 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, where 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 Eq. (31) in Eq. (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 A i,k A i,q = 0 and A † i,k A i,q ∝ δ k,q . Hence, Eq. (30) reduces to where L[A, ρ] is the Lindblad dissipator defined in Eq. (7), whereas Since the A i,k are eigenoperators of the bath which is in equilibrium at inverse temperature β i , it follows that which is another manifestation of local detailed balance.

B. Example: two harmonic oscillators
To show how these results may be applied, we consider again the two harmonic oscillators a 1 and a 2 treated in Sec. II. We assume that the environmental units E 1 and E 2 are themselves harmonic oscillators with operators b 1 and b 2 [103], and we take the total Hamiltonian to be where H S is given in Eq (1). Since the frequency ω 1 of S 1 and E 1 are the same, the operators a 1 and b 1 are eigenopera- Hence, local detailed balance, Eq. (34), is satisfied. However, a 1 is not an eigenoperator of H S due to the interaction term H I = (a † 1 a 2 + a 1 a † 2 ) [Eq. (3)] so that global detailed balance is in general broken.
Applying Eq. (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 γ ± i,k given by (40) which are the coefficients in Eq. (5).

C. 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 H E i 0 − H E i τ . Dividing by τ and taking the limit τ → 0 will give the heat currentQ i to environment E i . The explicit calculation is postponed to Appendix B. As a result, we find the surprisingly simple formula: For example, in the case of the two harmonic oscillators, Eq. (41) gives precisely the heat-rate formula (16) that was used in Sec. II.
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 Again, dividing by τ and taking the limit τ → 0 yields the work rateẆ ext . As shown in Appendix B, carrying out this computation yieldṡ 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 systemsystem interaction H I . If global detailed balance, Eq. (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 1 ] = − a 2 , so that a direct application of Eq. (43) leads to expression (17). From a practical point of view, Eqs. (41) and (43) are our main results, as they offer general expressions that may be applied over a broad range of situations. 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.

IV. QUANTUM HARMONIC OSCILLATORS
We consider once again the problem studied in Sec. II, 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 When η = 0 we recover the excitation-conserving tightbinding model. Conversely, η = corresponds to a positionposition coupling x i x i+1 . The total system evolves according to the master equation where D i (ρ S ) is given in Eq (5).
Using Eq. (41) we find that the heat rate to the first and last environments will bė The work rate, on the other hand, is found using Eq. (43) and readsẆ If N = 2 and η = 0 this reduces to Eq. (17).
A. 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 Sec. II to include the counterrotating terms η. We then show that, interpreting our steadystate as a heat engine, as in Fig. 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 η 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): + η 2 (ω 1 + ω 2 )(n 1 + n 2 + 1) .
The influence of η in these thermodynamic quantities is shown in Fig. 3, which compares the results with those of Fig. 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 Eqs. (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 [88,104,105], which is precisely what is observed in Fig. 3.

B. 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 Fig. 4 for the case of a linear function profile As can be seen, in general the thermodynamic quantities depend on the system size. However, surprisingly, their ratios (Figs. 4(d) and (e)) do not, behaving precisely as in Fig. 2.

V. SUMMARY AND CONCLUSIONS
In this paper we have systematically analysed the external work required for local master equations based on the repeated interaction model. Specifically, we have provided explicit expressions, Eq. (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 local master equations, 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 realising 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 discuss here how to compute the steady-state of Gaussian preserving master equations. For concreteness, we focus on the problem discussed in Sec. II, more specifically Eq. (4). The extension to multiple oscillators (Sec. IV) 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 vec-tor: We introduce the drift matrix containing the couplings and decay constants: where: and the system's covariance matrix V with entries: 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: The steady-state is then determined by the algebraic equation In absence of coupling between the first and second oscillator, = 0, the stationary state of each oscillator is thermal at its respective temperature: When 0, one finds instead the following expression for the entries of V (recall that V i j = V ji ): which are tantamount to Eqs. (8)-(9).

Appendix B: Details of the repeated interactions derivation
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 where H tot is given in Eq. (26) (with V i → V i / √ τ, as discussed in the main text). Hence, the evolution of any observable O (from either the system or from the environments) due to this interaction is given by 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 Carrying out the computations, using Eq. (31), we find 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 γ ± i,k in Eq. (37), we may then divide both sides by τ and take the limit τ → 0, which yields Moreover, with respect to the change in energy of the environments, we define the heat rates asQ i = −∆H E i /τ, which then yields precisely Eq. (41). The heat rate may now be found in two ways. First, using Eq. (25) and noting that we see that d H S / dt can be split into two terms, one of which is precisely iQi . Hence, the remainder must be attributed to the work that has to be performed, which yields precisely Eq. (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 Eq. (26), the Hamiltonian for this interaction is more appropriately defined as 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 Eq. (42). Carrying out the integration one then finds where, once again, the primed expectation value refers to the state after the interaction. Using once again Eq. (B1) for evaluating this expectation value, we then find Hence, taking again the limit τ → 0, one identifies the work rate asẆ = d H S dt − iQi .
Appendix C: Entropy production in local master equations 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 [86], we can define the entropy production in a single stroke as Σ := I(S : E), +S (ρ E (τ)||ρ th E ) ≥ 0. where 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 ∆S E i the change in the entropy of environment E i . The second term in Eq. (C1), on the other hand, is the quantum relative entropy between the final and initial states of each environment [here S (ρ||σ) = tr(ρ ln ρ − ρ ln σ)].
The positivity of Σ then follows immediately from the positivity of the mutual information and the quantum relative entropy.
Using Eq. (B8) one may then readily show that Eq. (C1) may be written as Dividing by τ and taking the limit τ → 0 we may then obtain the usual expression for the entropy production rate 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 where F S = H S − T S S is the free energy of the system. Thus, in this case we recover another well known expression for the entropy production rate.