Reconciliation of quantum local master equations with thermodynamics

Gabriele de Chiara, Gabriel Landi, Adam Hewgill, Brendan Reid, Alessandro

Ferraro, Augusto Roncaglia, Mauro Antezza

► To cite this version:

Gabriele de Chiara, Gabriel Landi, Adam Hewgill, Brendan Reid, Alessandro Ferraro, et al. Reconciliation of quantum local master equations with thermodynamics. New Journal of Physics, 2018, 20 (11), pp.113024. 10.1088/1367-2630/aaecee. hal-01925150

HAL Id: hal-01925150

https://hal.science/hal-01925150

Submitted on 13 Oct 2020

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PAPER • OPEN ACCESS

# Reconciliation of quantum local master equations with thermodynamics

To cite this article: Gabriele De Chiara et al 2018 *New J. Phys.* **20** 113024

View the [article online](#) for updates and enhancements.

## Recent citations

- Transport in boundary-driven quantum spin systems: one-way street for the energy current
  Deborah Oliveira et al

- Probing the dynamics of an open quantum system via a single qubit
  P C López Vázquez and T Gorin

- Analysis of the conditional mutual information in ballistic and diffusive non-equilibrium steady-states
  William T B Malouf et al
**New Journal of Physics**

The open access journal at the forefront of physics

<figure><img src="image_2.png" /></figure>

IOP Institute of Physics

Published in partnership
with: Deutsche Physikalische
Gesellschaft and the Institute
of Physics

<figure><img src="image_3.png" /></figure>

CrossMark

OPEN ACCESS

RECEIVED

6 September 2018

REVISED

26 October 2018

ACCEPTED FOR PUBLICATION

31 October 2018

PUBLISHED

16 November 2018

Original content from this
work may be used under
the terms of the *Creative
Commons Attribution 3.0
licence*.

Any further distribution of
this work must maintain
attribution to the
author(s) and the title of
the work, journal citation
and DOI.

<figure><img src="image_4.png" /></figure>

PAPER

Reconciliation of quantum local master equations with thermodynamics

Gabriele De Chiara<sup>1,2,7</sup> Ⓓ, Gabriel Landi<sup>3</sup>, Adam Hewgill<sup>1</sup>, Brendan Reid<sup>1</sup>, Alessandro Ferraro<sup>1</sup>,<br>
Augusto J Roncaglia<sup>4</sup> and Mauro Antezza<sup>2,5,6</sup> Ⓓ

1 Centre for Theoretical Atomic, Molecular and Optical Physics, Queen's University Belfast, Belfast BT7 1NN, United Kingdom

2 Kavli Institute of Theoretical Physics (KITP), University of California, Santa Barbara CA 93106-4030, United States of America

3 Instituto de Física da Universidade de São Paulo, 05314-970 São Paulo, Brazil

4 Departamento de Física, Facultad de Ciencias Exactas y Naturales, Universidad de Buenos Aires and IFIBA, CONICET, Ciudad Universitaria, 1428 Buenos Aires, Argentina

5 Laboratoire Charles Coulomb (L2C), UMR 5221 CNRS-Université de Montpellier, F-34095 Montpellier, France

6 Institut Universitaire de France, 1 rue Descartes, F-75231 Paris Cedex 05, France

7 Author to whom any correspondence should be addressed.

E-mail: g.dechiara@qub.ac.uk

**Keywords:** open quantum systems, quantum thermodynamics, master equations, quantum harmonic oscillators

Abstract

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 modeling of transport in mesoscopic systems, consists in using *local* master equations (LMEs) 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 LMEs 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.

**1. Introduction**

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_1, ..., 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
<figure><img src="image_5.png" /></figure>

**Figure 1.** We consider in this paper a system S composed of several sub-systems $S_1, \dots, S_N$ (in the figure $N = 2$). Each subsystem $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$.

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.

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 $\tau$ (see figure 1(b)). In the limit $\tau \to 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.

## 2. Steady-state of two harmonic oscillators subject to a LME

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_1$ and $S_2$, each coupled to its own local heat bath via a LME (see figure 1(a)). The Hamiltonian is taken to be

$$H_S = H_1 + H_2 + H_I, \qquad (1)$$

where

$$H_I = \frac{1}{2}(p_I^2 + \omega_I^2 x_I^2) = \omega_I (a_I^\dagger a_I + 1/2), \qquad (2)$$

and

$$H_I = \epsilon (a_1^\dagger a_2 + a_1 a_2^\dagger). \qquad (3)$$

Here $a_i$ are the annihilation operators of each harmonic oscillator and $(x_i, p_i)$ the corresponding position and momentum ($\hbar = k_B = 1$). Both oscillators are assumed to have unit mass, but different frequencies $\omega_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, $\omega_1 = \omega_2$. More general interactions will be discussed in section 4.

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

$$\frac{d\rho_S}{dt} = -i[H_S, \rho_S] + D_1(\rho_S) + D_2(\rho_S), \qquad (4)$$

where

$$D_i(\rho_S) = \gamma(n_i + 1)\mathcal{L}[a_i, \rho] + \gamma n_i \mathcal{L}[a_i^\dagger, \rho]. \qquad (5)$$

Here $\gamma$ is the dissipation rate, $n_i$ is the average excitation number for the Bose-Einstein distribution

$$n_i = \frac{1}{e^{\beta_i \omega_i} - 1}, \qquad (6)$$
with inverse temperature $\beta_i = (T_i)^{-1}$ and

$$ \mathcal{L}[A, \rho] = A\rho A^{\dagger} - \frac{1}{2}\{A^{\dagger}A, \rho\}, \qquad (7) $$

is the dissipator in Lindblad form. Henceforth, we always assume that $T_1 < T_2$ (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/\omega_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:

$$ \langle a_1^\dagger a_1 \rangle = n_1 + \frac{2\epsilon^2}{\Delta^2}(n_2 - n_1), \qquad (8) $$

$$ \langle a_2^\dagger a_2 \rangle = n_2 - \frac{2\epsilon^2}{\Delta^2}(n_2 - n_1), \qquad (9) $$

$$ \langle a_2^\dagger a_1 \rangle = \frac{\epsilon}{\Delta^2}(\omega_1 - \omega_2 + i\gamma)(n_1 - n_2), \qquad (10) $$

where $\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

$$ \langle H_1 \rangle = \langle H_1 \rangle_{\text{th}} + \frac{2\omega_1 \epsilon^2}{\Delta^2}(n_2 - n_1), \qquad (11) $$

where $\langle H_1 \rangle_{\text{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

$$ \dot{Q}_{1 \to 2} = -i \langle [H_S, H_1] \rangle = -i \omega_1 \epsilon \langle a_1^\dagger a_2 - a_1 a_2^\dagger \rangle = \frac{2\gamma\epsilon^2\omega_1}{\Delta^2}(n_2 - n_1), \qquad (12) $$

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

The relevant point to highlight from these results is that the sign of the $\epsilon^2$-terms in equations (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 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_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

$$ n_1 > n_2 \quad \leftrightarrow \quad \frac{T_1}{\omega_1} > \frac{T_2}{\omega_2}, \qquad (13) $$

then there will be a current flowing from the cold to the hot bath ($\dot{Q}_{1 \to 2} < 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

$$ \frac{d\langle H_S \rangle}{dt} = \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. \qquad (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.
<table><caption>**Table 1.** 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. So $\dot{Q}_i > 0$ means energy entered S through $E_i$ and $\dot{W}_{ext} > 0$ means work was performed on the system by an external agent.</caption><thead><tr><th>Mode of operation</th><th>Q̇<sub>1</sub></th><th>Q̇<sub>2</sub></th><th>Ẇ</th><th>Operation regime</th></tr></thead><tbody><tr><td><strong>Refrigerator</strong><br>Consumes work to make heat flow from cold to hot.</td><td>Q̇<sub>1</sub> > 0</td><td>Q̇<sub>2</sub> < 0</td><td>Ẇ > 0</td><td>ω<sub>1</sub> < ω<sub>2</sub>T<sub>1</sub>/T<sub>2</sub></td></tr><tr><td><strong>Engine</strong><br>Uses hot bath to produce useful work, dumping the remainder in the cold bath.</td><td>Q̇<sub>1</sub> < 0</td><td>Q̇<sub>2</sub> > 0</td><td>Ẇ < 0</td><td>ω<sub>2</sub>T<sub>1</sub>/T<sub>2</sub> < ω<sub>1</sub> < ω<sub>2</sub></td></tr><tr><td><strong>Accelerator (oven)</strong><br>Consumes work to heat the cold reservoir more than it would with spontaneous thermal conduction.</td><td>Q̇<sub>1</sub> < 0</td><td>Q̇<sub>2</sub> > 0</td><td>Ẇ > 0</td><td>ω<sub>1</sub> > ω<sub>2</sub></td></tr></tbody></table>

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

$$ \frac{d\langle H_S \rangle}{dt} = \dot{Q}_1 + \dot{Q}_2 + \dot{W}_{ext}, \qquad (15) $$

where

$$ \dot{Q}_i = \operatorname{tr}[D_i(\rho_S)H_i] = \gamma\omega_i(n_i - \langle a_i^\dagger a_i \rangle), \qquad (16) $$

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

$$ \dot{W}_{ext} = \operatorname{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, \qquad (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

$$ \dot{Q}_1 = \frac{2\gamma\epsilon^2}{\Delta^2}\omega_1(n_1 - n_2), \qquad (18) $$

$$ \dot{Q}_2 = -\frac{2\gamma\epsilon^2}{\Delta^2}\omega_2(n_1 - n_2), \qquad (19) $$

$$ \dot{W}_{ext} = -\frac{2\gamma\epsilon^2}{\Delta^2}(\omega_1 - \omega_2)(n_1 - n_2). \qquad (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}_{ext}$ tends to zero as $\epsilon \to 0$, in which case global detailed balance is recovered. Second, the work required is proportional to $\gamma\epsilon^2$. As discussed in [52], the LME can be shown to be correct for small interactions $\epsilon$ up to the same order (see comment in section 2.4 in that paper). Third, in the steady-state $d\langle H_S \rangle/dt = 0$ so that we should have

$$ \dot{W}_{ext} = -\dot{Q}_1 - \dot{Q}_2, \qquad (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.

We begin with the case of a refrigerator, $\dot{Q}_1 > 0$. According to equation (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 $\dot{W}_{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 COP of the refrigerator is the ratio of the heat extracted from the cold environment divided by the work invested:

$$ \text{COP} = \frac{\dot{Q}_1}{\dot{W}_{ext}} = \frac{\omega_1}{\omega_2 - \omega_1} \qquad (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
<figure><img src="image_6.png" /><figcaption>Figure 2. The three regimes of operation in table 1, for the steady-state of two harmonic oscillators. (a)–(c) $\dot{Q}_1$, $\dot{Q}_2$ and $\dot{W}_{ext}$, equations (18)–(20), as a function of $\omega_1/\omega_2$. (d) $\dot{Q}_1/\dot{W}_{ext}$, which follows precisely the Otto value $\omega_1/(\omega_2 - \omega_1)$. (e) $\dot{W}_{ext}/\dot{Q}_2$ which also follows the Otto value $\omega_1/\omega_2 - 1$. All curves are plotted for $\gamma e^2/\Delta^2 = \omega_2$, assuming $T_1/T_2 = 1/2$. The vertical lines represent $\omega_1/\omega_2 = T_1/T_2$ and $\omega_1 = \omega_2$.</figcaption></figure>

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 $\omega_1$ in the interval

$$ \frac{T_1}{T_2} \omega_2 < \omega_1 < \omega_2, \qquad (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:

$$ \eta = \frac{|\dot{W}_{ext}|}{\dot{Q}_2} = 1 - \frac{\omega_1}{\omega_2} \qquad (24) $$

which is again the Otto efficiency for a thermal machine. Finally, for $\omega_2 < \omega_1$ we necessarily have $n_1 < n_2$ and therefore $\dot{Q}_1 < 0$, $\dot{Q}_2 > 0$ and $\dot{W}_{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: $|\dot{Q}_1| = \dot{W}_{ext} + \dot{Q}_2 > \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 $\omega_1/\omega_2 = T_1/T_2$, 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.

### 3. The method of 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 $\tau$ 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 $\tau \to 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_1, \dots, S_N$, where each $S_i$ interacts with its own environment $E_i$. The Hamiltonian of the system is taken to be
$$H_S = \sum_{i=1}^{N} H_{S_i} + H_I, \qquad (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

$$H_{\text{tot}} = H_S + \sum_{i=1}^{N} H_{E_i} + \frac{1}{\sqrt{\tau}} \sum_{i=1}^{N} V_i, \qquad (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 $\tau$, 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 $\rho_S(t)$ will then evolve according to the map

$$\rho_S((n + 1)\tau) = \text{tr}_E\{e^{-i\tau H_{\text{tot}}} \rho_S(n\tau) \rho_E e^{i\tau H_{\text{tot}}}\}, \qquad (27)$$

where

$$\rho_E = \bigotimes_{i=1}^{N} \frac{e^{-\beta_i H_{E_i}}}{\text{tr} e^{-\beta_i H_{E_i}}}, \qquad (28)$$

is the thermal state of the environments.

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

$$\frac{d\rho_S}{dt} = -i[H_S, \rho_S] + \sum_{i=1}^{N} D_i(\rho_S), \qquad (29)$$

where $H_S$ is given in (25) and

$$D_i(\rho_S) = -\frac{1}{2} \text{tr}_{E_i}[V_i, [V_i, \rho_S \rho_{E_i}]]. \qquad (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_i$ [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] = -\omega A$, for some $\omega > 0$ (usually referred to as a Bohr, or transition, frequency). Due to this algebra, $A$ and $A^\dagger$ function as lowering and raising operators for the $H$, causing transitions between energy levels $E$ and $E'$ separated by an energy $\omega$.

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

$$V_i = \sum_k g_{i,k} (L_{i,k}^{\dagger} A_{i,k} + L_{i,k} A_{i,k}^{\dagger}), \qquad (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:

$$[H_{S_i}, L_{i,k}] = -\omega_{i,k} L_{i,k}, \qquad (32)$$

$$[H_{E_i}, A_{i,k}] = -\omega_{i,k} A_{i,k}, \qquad (33)$$

where, for each $S_i$, $\{\omega_{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 $\{\omega_{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 $\omega$ (the natural frequency of the oscillator). Similarly, a system-environment interaction containing $a^2$ would induce transitions with Bohr frequencies $2\omega$.
In an interaction such as (31), however, there is also the additional assumption that whenever there is a transition of $+\omega_{i,k}$ in the system, the corresponding transition in $E_i$ occurs with energy $-\omega_{i,k}$. Hence, it follows that

$$[H_{S_i} + H_{E_i}, V_i] = 0, \qquad (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

$$[H_S + H_E, V] = [H_I, V] \neq 0, \qquad (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

$$D_i(\rho_S) = \sum_k \gamma_{i,k} \mathcal{L}[L_{i,k}, \rho_S] + \gamma_{i,k}^\dagger \mathcal{L}[L_{i,k}^\dagger, \rho_S], \qquad (36)$$

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

$$\begin{aligned} \gamma_{i,k}^- &= g_{i,k}^2 \langle A_{i,k} A_{i,k}^\dagger \rangle, \\ \gamma_{i,k}^\dagger &= g_{i,k}^2 \langle A_{i,k}^\dagger A_{i,k} \rangle. \end{aligned} \qquad (37)$$

Since the $A_{i,k}$ are eigenoperators of the bath which is in equilibrium at inverse temperature $\beta$, it follows that

$$\frac{\gamma_{i,k}^{\dagger}}{\gamma_{i,k}^{-}} = e^{-\beta_i \omega_{i,k}}, \qquad (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_1$ and $a_2$ treated in section 2. We assume that the environmental units $E_1$ and $E_2$ are themselves harmonic oscillators with operators $b_1$ and $b_2$<fnref n="8" />, and we take the total Hamiltonian to be

$$H_{\text{tot}} = H_S + \sum_{i=1,2} \{ \omega_i b_i^\dagger b_i + g (a_i^\dagger b_i + b_i^\dagger a_i) \}, \qquad (39)$$

where $H_S$ is given in equation (1). Since the frequency $\omega_1$ of $S_1$ and $E_1$ are the same, the operators $a_1$ and $b_1$ 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_1$ 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}^\dagger$ given by

$$\gamma_i^\dagger = 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), \qquad (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 $\delta 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 $\tau$ and taking the limit $\tau \to 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:

$$\dot{Q}_i = \sum_k \omega_{i,k} \{ \gamma_{i,k}^\dagger \langle L_{i,k} L_{i,k}^\dagger \rangle - \gamma_{i,k} \langle L_{i,k}^\dagger L_{i,k} \rangle \}. \qquad (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.

<footnote n="8">⁸ Note that it is possible to arrive at the same results by assuming the environmental units to be qubits.</footnote>
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

$$
\delta W_{\text{ext}} = \int_{n\tau}^{(n+1)\tau} \left\langle \frac{\partial H_{\text{tot}}}{\partial t} \right\rangle dt. \qquad (42)
$$

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

$$
\dot{W}_{\text{ext}} = \frac{1}{2} \sum_{i,k} \{ \gamma_{i,k} \langle L_{i,k}^{\dagger} F_{i,k} \rangle - \gamma_{i,k}^{\dagger} \langle F_{i,k} L_{i,k}^{\dagger} \rangle \} + \text{c.c.}, \quad (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$_{1}$ ] = − ϵa$_{2}$ , 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.

4. Quantum harmonic oscillators

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

$$
H_S = \sum_{i=1}^{N} \omega_i a_i^\dagger a_i + \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^\dagger) \}. \quad (44)
$$

When η = 0 we recover the excitation-conserving tight-binding model. Conversely, η = ϵ corresponds to a position–position coupling x$_{i}$x$_{i}$$_{+1}$ . The total system evolves according to the master equation

$$
\frac{d\rho_S}{dt} = -i[H_S, \rho_S] + D_1(\rho_S) + D_N(\rho_S), \qquad (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

$$
\dot{Q}_i = \gamma \omega_i (n_i - \langle a_i^\dagger a_i \rangle), \quad i = 1, N. \tag{46}
$$

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

$$
\begin{equation}
\begin{split}
\dot{W} = \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 \\
& + \frac{\gamma \eta}{2} \langle a_1 a_2 + a_1 a_2^\dagger + a_N a_{N-1} + a_N a_{N-1}^\dagger \rangle.
\end{split}
\tag{47}
\end{equation}
$$

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 η ≠ 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):
<figure><img src="image_7.png" /><figcaption>Figure 3. Similar to figure 2, but examining the influence of the counter-rotating term η in the Hamiltonian (44). The solid curves are the same as figure 2. The dashed curves correspond to equations (48)–(50) for η = 0 . 3 ω 2 (dashed) and η = 0 . 6 ω 2 (dotted-dashed). The curves were plotted for ϵ $^{2}$/γ = ω 2 , assuming T 1 /T 2 = 1 / 2. The vertical lines correspond to ω 1 /ω 2 = T 1 /T 2 and ω 1 = ω 2 .</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Similar to figure 2, but for a chain of N = 2 , 5 and 20 oscillators (solid, dashed and dotted-dashed, respectively) following a linear frequency profile (equation (51)). The parameters were γ = 2 ω$_{N}$ , ϵ = 0 . 25 ω$_{N}$ , η = 0 and T 1 /T 2 = 1 / 2.</figcaption></figure>

$$
\dot{Q}_1 = \frac{2\omega_1}{\gamma} [\epsilon^2(n_1 - n_2) - \eta^2(n_1 + n_2 + 1)], \quad (48)
$$

$$
\dot{Q}_2 = -\frac{2\omega_2}{\gamma} [\epsilon^2(n_1 - n_2) + \eta^2(n_1 + n_2 + 1)], \quad (49)
$$

$$
\begin{equation}
\begin{split}
W_{\text{ext}} = -\frac{2}{\gamma} & [\epsilon^2(n_1 - n_2)(\omega_1 - \omega_2) + \\
& \qquad - \eta^2(\omega_1 + \omega_2)(n_1 + n_2 + 1)].
\end{split}
\tag{50}
\end{equation}
$$

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.

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

$$
\omega_i = \frac{(N - i)\omega_1 + (i - 1)\omega_N}{N - 1}. \qquad (51)
$$

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 ω$_{1}$ 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]:

$$
\frac{d\langle a_1^\dagger a_1 \rangle}{dt} = -i\epsilon \langle a_1^\dagger a_2 - a_1 a_2^\dagger \rangle + \gamma (n_1 - \langle a_1^\dagger a_1 \rangle), \quad (52)
$$

$$
\frac{d\langle a_i^\dagger a_i \rangle}{dt} = -i\epsilon \langle a_i^\dagger a_{i+1} - a_i a_{i+1}^\dagger \rangle, \quad i = 2, \dots, N-1 \quad (53)
$$
$$ \frac{d\langle a_N^\dagger a_N \rangle}{dt} = 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). \quad (54) $$

Hence, at steady-state we have:

$$ \gamma(n_1 - \langle a_1^\dagger a_1 \rangle) = -\gamma(n_N - \langle a_N^\dagger a_N \rangle) \quad (55) $$

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

$$ \dot{Q}_1 = \gamma\omega_1(n_1 - \langle a_1^\dagger a_1 \rangle) = -\frac{\omega_1}{\omega_N}\dot{Q}_N. \quad (56) $$

Finally, using the first law of thermodynamics: $\dot{W}_{\text{ext}} + \dot{Q}_1 + \dot{Q}_N = 0$ we obtain:

$$ \frac{\dot{Q}_1}{\dot{W}_{\text{ext}}} = \frac{\omega_1}{\omega_N - \omega_1}; \quad \frac{\dot{W}_{\text{ext}}}{\dot{Q}_N} = -1 + \frac{\omega_1}{\omega_N}. \quad (57) $$

We remark that these ratios coincide with those of Otto thermal machines and refrigerators operating with a single harmonic oscillator driven between frequencies $\omega_1$ and $\omega_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).

## 5. Summary and Conclusions

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.

## Acknowledgments

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.

## Appendix A. The Lyapunov equation

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:

$$ Y^T = (x_1, p_1, x_2, p_2). \quad (A1) $$
We introduce the drift matrix containing the couplings and decay constants:

$$K = \begin{pmatrix} -\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{pmatrix}, \qquad (A2)$$

where:

$$\lambda = \frac{2\epsilon}{\sqrt{\omega_1\omega_2}}, \qquad \mu = 2\epsilon\sqrt{\omega_1\omega_2} \qquad (A3)$$

and the system's covariance matrix V with entries:

$$V(i, j) = \frac{1}{2} \langle \{Y(i), Y(j)\} \rangle \qquad (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:

$$\dot{V} = KV + VK^T + D, \qquad (A5)$$

where

$$D = \mathrm{diag} \left[ \frac{\gamma(2n_1 + 1)}{2\omega_1}, \frac{\gamma\omega_1(2n_1 + 1)}{2}, \right. \\ \left. \times \frac{\gamma(2n_2 + 1)}{2\omega_2}, \frac{\gamma\omega_2(2n_2 + 1)}{2} \right].$$

The steady-state is then determined by the algebraic equation

$$KV + VK^T + D = 0. \qquad (A6)$$

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

$$V = \begin{pmatrix} \frac{2n_1 + 1}{2\omega_1} & 0 & 0 & 0 \\ 0 & \frac{(2n_1 + 1)\omega_1}{2} & 0 & 0 \\ 0 & 0 & \frac{2n_2 + 1}{2\omega_2} & 0 \\ 0 & 0 & 0 & \frac{(2n_2 + 1)\omega_2}{2} \end{pmatrix}. \qquad (A7)$$

When $\epsilon \neq 0$, one finds instead the following expression for the entries of V (recall that $V_{ij} = V_{ji}$):

$$V_{11} = \frac{1}{2\omega_1} \left( \frac{4\epsilon^2(n_2 - n_1)}{\Delta^2} + 2n_1 + 1 \right), \qquad (A8)$$

$$V_{22} = \frac{1}{2}\omega_1 \left( \frac{4\epsilon^2(n_2 - n_1)}{\Delta^2} + 2n_1 + 1 \right), \qquad (A9)$$

$$V_{33} = \frac{1}{2\omega_2} \left( \frac{4\epsilon^2(n_1 - n_2)}{\Delta^2} + 2n_2 + 1 \right), \qquad (A10)$$

$$V_{44} = \frac{1}{2}\omega_2 \left( \frac{4\epsilon^2(n_1 - n_2)}{\Delta^2} + 2n_2 + 1 \right), \qquad (A11)$$

$$V_{13} = \frac{\epsilon(n_1 - n_2)(\omega_1 - \omega_2)}{\sqrt{\omega_1 \omega_2} \Delta^2}, \qquad (A12)$$

$$V_{14} = \frac{\gamma\omega_2\epsilon(n_2 - n_1)}{\sqrt{\omega_1\omega_2}\Delta^2}, \qquad (A13)$$

$$V_{23} = \frac{\gamma\omega_1\epsilon(n_1 - n_2)}{\sqrt{\omega_1\omega_2}\Delta^2}, \qquad (A14)$$

$$V_{24} = \frac{\epsilon(n_1 - n_2)(\omega_1 - \omega_2)\sqrt{\omega_1\omega_2}}{\Delta^2}, \qquad (A15)$$
$$V_{12} = V_{34} = 0, \qquad (A16)$$

which are tantamount to equations (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

$$\rho_{SE}' = e^{-i\tau H_{tot}} \rho_S \rho_E e^{i\tau H_{tot}},$$

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 $O$ (from either the system or from the environments) due to this interaction is given by

$$\langle O \rangle' = \langle e^{i\tau H_{tot}} O e^{-i\tau H_{tot}} \rangle, \qquad (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 $\tau$

$$\Delta H_S = -\frac{\tau}{2} \sum_{i=1}^{N} \langle [V_i, [V_i, H_S]] \rangle, \qquad (B2)$$

$$\Delta H_{E_i} = -\frac{\tau}{2} \langle [V_i, [V_i, H_{E_i}]] \rangle. \qquad (B3)$$

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

$$\Delta H_S = \frac{\tau}{2} \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 \}$$

$$\Delta H_{E_i} = -\tau \sum_k g_{i,k}^2 \omega_{i,k} \{ \langle A_{i,k}^\dagger A_{i,k} \rangle \langle L_{i,k}^\dagger L_{i,k} \rangle - \langle A_{i,k} A_{i,k}^\dagger \rangle \langle L_{i,k}^\dagger L_{i,k} \rangle \}.$$

Here the averages over the environment operators are always computed with respect to the same state $\rho_E$, whereas the average over system operators are computed over the instantaneous state $\rho_S(t)$.

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

$$\frac{d\langle H_S \rangle}{dt} = \frac{1}{2} \sum_{i,k} \{ \gamma_{i,k}^{\pm} \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 \}. \qquad (B4)$$

Moreover, with respect to the change in energy of the environments, we define the heat rates as $\dot{Q}_i = -\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

$$[H_S, L_{i,k}] = -\omega_{i,k} L_{i,k} + [H_i, L_{i,k}], \qquad (B5)$$

we see that $d\langle H_S \rangle/dt$ 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

$$H_{tot} = H_S + H_E + \frac{\lambda(t)}{\sqrt{\tau}} \sum_{i=1}^{N} V_i, \qquad (B6)$$

where $\lambda(t)$ has the value 1 when $t \in [n\tau, (n + 1)\tau]$ 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

$$\delta W = \sum_{i=1}^{N} \frac{\langle V_i \rangle - \langle V_i \rangle'}{\sqrt{\tau}}, \qquad (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
$$ \delta W = \frac{\tau}{2} \sum_{i=1}^{N} \langle [V_i, [V_i, H_S + H_{E_i}]] \rangle = \Delta H_S + \sum_{i=1}^{N} \Delta H_{E_i} \quad (B8) $$

Hence, taking again the limit $\tau \to 0$, one identifies the work rate as $\dot{W} = \frac{d\langle H_S \rangle}{dt} - \sum_i \dot{Q}_i$.

## Appendix C. Entropy production in LMEs

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 $\tau$, 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

$$ \Sigma := \mathcal{I}(S: E), +S(\rho_E(\tau) || \rho_E^{th}) \ge 0, \quad (C1) $$

where

$$ \mathcal{I}(S: E) = \Delta S_S + \Delta S_E \approx \Delta S_S + \sum_{i=1}^{N} \Delta S_{E,i}, \quad (C2) $$

is the mutual information developed between the system and all the environments, with $\Delta S_S$ being the change in the entropy of the system and $\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

$$ S(\rho_E(\tau) || \rho_E^{th}) \approx \sum_{i=1}^{N} S(\rho_{E,i}(\tau) || \rho_{E,i}^{th}), \quad (C3) $$

is the quantum relative entropy between the final and initial states of each environment (here $S(\rho||\sigma) = \text{tr}(\rho \ln\rho - \rho \ln\sigma)$). The positivity of $\Sigma$ 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

$$ \Sigma = \Delta S_S - \sum_{i=1}^{N} \beta_i \delta Q_i. \quad (C4) $$

Dividing by $\tau$ and taking the limit $\tau \to 0$ we may then obtain the usual expression for the entropy production rate

$$ \Pi = \frac{dS_S}{dt} - \sum_{i=1}^{N} \beta_i \frac{dQ_i}{dt}, \quad (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

$$ \Pi = \beta \left( \frac{dW}{dt} - \frac{dF}{dt} \right), \quad (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.

## ORCID iDs

Gabriele De Chiara https://orcid.org/0000-0003-3265-9021

Mauro Antezza https://orcid.org/0000-0003-4540-5864

## References

[1] Giazotto F and Martinez-Pérez M J 2012 Nature 492 401

[2] Dutta B, Peltonen J T, Antonenko D S, Meschke M, Skvortsov M A, Kubala B, König J, Winkelmann C B, Courtois H and Pekola J P 2017 Phys. Rev. Lett. 119 077701

[3] Caruso F, Chin A W, Datta A, Huelga S F and Plenio M B 2009 J. Chem. Phys. 131 105106

[4] Schlawin F, Dorfman K E, Fingerhut B P and Mukamel S 2013 Nat. Commun. 4 1782

[5] Leggio B, Messina R and Antezza M 2015a Europhys. Lett. 110 40002

[6] Doyeux P, Messina R, Leggio B and Antezza M 2017 Phys. Rev. A 95 012138

[7] Brantut J-P, Grenier C, Meineke J, Stadler D, Krinner S, Kollath C, Esslinger T and Georges A 2013 Science 342 713

[8] Brunelli M et al 2018 Phys. Rev. Lett. 121 160604

[9] Landig R, Hruby L, Dogra N, Landini M, Mottl R, Donner T and Esslinger T 2016 Nature 532 476

[10] Barreiro J T, Müller M, Schindler P, Nigg D, Monz T, Chwalla M, Hennrich M, Roos C F, Zoller P and Blatt R 2011 Nature 470 486
[11] Bermudez A, Bruderer M and Plenio M B 2013 *Phys. Rev. Lett.* 111 040601

[12] Freitas N, Martinez E A and Paz J P 2016 *Phys. Scr.* 91 013007

[13] Jarzynski C 1997 *Phys. Rev. Lett.* 78 2690

[14] Crooks G E 1998 *J. Stat. Phys.* 90 1481

[15] Tasaki H 2000 arXiv:cond-mat/0009244

[16] Mukamel S 2003 *Phys. Rev. Lett.* 90 170604

[17] Jarzynski C and Wójcik D K 2004 *Phys. Rev. Lett.* 92 230602

[18] Campisi M, Hänggi P and Talkner P 2011 *Rev. Mod. Phys.* 83 771

[19] Millen J and Xuereb A 2016 *New J. Phys.* 18 011002

[20] Vinjanampathy S and Anders J 2016 *Contemp. Phys.* 57 545

[21] Goold J, Huber M, Riera A, del Rio L and Skrzypczyk P 2016 *J. Phys. A: Math. Theor.* 49 143001

[22] Scovil H E D and Schulz-DuBois E O 1959 *Phys. Rev. Lett.* 2 262

[23] Geva E and Kosloff R 1992 *J. Chem. Phys.* 96 3054

[24] He J, Chen J and Hua B 2002 *Phys. Rev.* E 65 036145

[25] Quan H T, Liu Y-x, Sun C P and Nori F 2007 *Phys. Rev.* E 76 031105

[26] Linden N, Popescu S and Skrzypczyk P 2010 *Phys. Rev. Lett.* 105 130401

[27] Kosloff R and Levy A 2014 *Annu. Rev. Phys. Chem.* 65 365

[28] Correa L A, Palao J P, Alonso D and Adesso G 2014 *Sci. Rep.* 4 3949

[29] Zheng Y and Poletti D 2014 *Phys. Rev.* E 90 012145

[30] Uzdin R, Levy A and Kosloff R 2015 *Phys. Rev. X* 5 031044

[31] Leggio B, Bellomo B and Antezza M 2015b *Phys. Rev. A* 91 012117

[32] Leggio B, Doyeux P, Messina R and Antezza M 2015c *Europhys. Lett.* 112 40004

[33] Leggio B and Antezza M 2016 *Phys. Rev.* E 93 022122

[34] Doyeux P, Leggio B, Messina R and Antezza M 2016 *Phys. Rev.* E 93 022134

[35] Abah O and Lutz E 2016 *Europhys. Lett.* 113 60002

[36] Niedenzu W, Gelbwaser-Klimovsky D, Kofman A G and Kurizki G 2016 *New J. Phys.* 18 083012

[37] Kosloff R and Rezek Y 2017 *Entropy* 19 136

[38] Roulet A, Nimmrichter S, Arrazola J M, Seah S and Scarani V 2017 *Phys. Rev.* E 95 062131

[39] Reid B, Pigeon S, Antezza M and De Chiara G 2017 *Europhys. Lett.* 120 60006

[40] Scopa S, Landi G T and Karevski D 2018 *Phys. Rev. A* 97 062121

[41] Hewgill A, Ferraro A and De Chiara G 2018 *Phys. Rev. A* 98 042102

[42] Alicki R and Lendi K 1987 *Lecture Notes in Physics* (Berlin: Springer)

[43] Gardiner C and Zoller P 2004 *Quantum Noise* 3rd edn (Berlin: Springer)

[44] Breuer H-P and Petruccione F 2007 *The Theory of Open Quantum Systems* (Oxford: Oxford University Press)

[45] Rivas A and Huelga S F 2012 *Open Quantum Systems* (Berlin: Springer)

[46] Caldeira A O and Leggett A J 1981 *Phys. Rev. Lett.* 46 211

[47] Leggett A J, Chakravarty S, Dorsey A T, Fisher M P A, Garg A and Zwerger W 1987 *Rev. Mod. Phys.* 59 1

[48] Caldeira A O 2014 *An Introduction to Macroscopic Quantum Phenomena and Quantum Dissipation* 1st edn (Cambridge: Cambridge University Press)

[49] Rivas Á, Plato A D K, Huelga S F and Plenio M B 2010 *New J. Phys.* 12 113032

[50] Purkayastha A, Dhar A and Kulkarni M 2016 *Phys. Rev. A* 93 062114

[51] Santos J P and Landi G T 2016 *Phys. Rev. E* 94 062143

[52] González J O, Correa L A, Nocerino G, Palao J P, Alonso D and Adesso G 2017 *Open Syst. Inf. Dyn.* 24 1740010

[53] Hofer P P, Perarnau-Llobet M, Miranda L D M, Haack G, Silva R, Brask J B and Brunner N 2017 *New J. Phys.* 19 123037

[54] Wichterich H, Henrich M J, Breuer H-P, Gemmer J and Michel M 2007 *Phys. Rev. E* 76 031115

[55] Karevski D and Platini T 2009 *Phys. Rev. Lett.* 102 207207

[56] Karevski D, Popkov V and Schütz G M 2013 *Phys. Rev. Lett.* 110 047201

[57] Asadian A, Manzano D, Tiersch M and Briegel H J 2013 *Phys. Rev. E* 87 012109

[58] Landi G T, Novais E, de Oliveira M J and Karevski D 2014 *Phys. Rev. E* 90 042142

[59] Landi G T and Karevski D 2015 *Phys. Rev. B* 91 174422

[60] Schuab L, Pereira E and Landi G T 2016 *Phys. Rev. E* 94 042122

[61] Guimarães P H, Landi G T and de Oliveira M J 2016 *Phys. Rev. E* 94 032139

[62] Cui J, Cirac J I and Bañuls M C 2015 *Phys. Rev. Lett.* 114 220601

[63] Prosen T 2011a *Phys. Rev. Lett.* 107 137201

[64] Prosen T 2011b *Phys. Rev. Lett.* 106 217206

[65] Prosen T and Pižorn I 2008 *Phys. Rev. Lett.* 101 1

[66] Prosen T and Žnidarič M 2012 *Phys. Rev. B* 86 125118

[67] Prosen T and Žnidarič M 2013 *Phys. Rev. Lett.* 111 124101

[68] Nicacio F, Ferraro A, Imparato A, Paternostro M and Semião F L 2015 *Phys. Rev. E* 91 042116

[69] Campbell S, De Chiara G and Paternostro M 2016 *Sci. Rep.* 6 19730

[70] Nicacio F, Paternostro M and Ferraro A 2016 *Phys. Rev. A* 94 052129

[71] Türkpençe D, Altintas F, Paternostro M and Müstecaplıoğlu O E 2017 *Europhys. Lett.* 117 50002

[72] Chanda T, Das T, Sadhukhan D, Pal A K, Sen(De) A and Sen U 2018 *Phys. Rev. A* 97 062324

[73] Pereira E 2018 *Phys. Rev. E* 97 022115

[74] Cresser J 1992 *J. Mod. Opt.* 39 2187

[75] Levy A and Kosloff R 2014 *Europhys. Lett.* 107 20004

[76] Correa L A, Palao J P, Adesso G and Alonso D 2013 *Phys. Rev. E* 87 042131

[77] Trushechkin A S and Volovich I V 2016 *Europhys. Lett.* 113 30005

[78] Stockburger J T and Motz T 2017 *Fortschr. Phys.* 65 1600067

[79] Deçordi G and Vidiella-Barranco A 2017 *Opt. Commun.* 387 366

[80] Mitchison M T and Plenio M B 2018 *New J. Phys.* 20 033005

[81] Naseem M T, Xuereb A and Müstecaplıoğlu Ö E 2018 arXiv:1806.08175

[82] Shammah N, Ahmed S, Lambert N, De Liberato S and Nori F 2018 arXiv:1805.05129
[83] Katz G and Kosloff R 2016 *Entropy* **18** 186

[84] Seah S, Nimmrichter S and Scarani V 2018 *Phys. Rev. E* **98** 012131

[85] Scarani V, Ziman M, Štelmachovič P, Gisin N and Bužek V 2002 *Phys. Rev. Lett.* **88** 097905

[86] Ziman M, Štelmachovič P, Bužek V, Hillery M, Scarani V and Gisin N 2002 *Phys. Rev. A* **65** 042105

[87] Giovannetti V and Palma G M 2012 *Phys. Rev. Lett.* **108** 040401

[88] Barra F 2015 *Sci. Rep.* **5** 14873

[89] Strasberg P, Schaller G, Brandes T and Esposito M 2017 *Phys. Rev. X* **7** 021003

[90] Manzano G, Galve F, Zambrini R and Parrondo J M R 2016 *Phys. Rev. E* **93** 052120

[91] Ockeloen-Korppi C F, Damskägg E, Pirkkalainen J-M, Asjad M, Clerk A A, Massel F, Woolley M J and Sillanpää M A 2018 *Nature* **556** 478

[92] Martinez E A and Paz J P 2013 *Phys. Rev. Lett.* **110** 130406

[93] Buffoni L, Solfanelli A, Verrucchi P, Cuccoli A and Campisi M 2018 arXiv:1806.07814

[94] Ciccarello F, Palma G M and Giovannetti V 2013 *Phys. Rev. A* **87** 040103

[95] Lorenzo S, Farace A, Ciccarello F, Palma G M and Giovannetti V 2015a *Phys. Rev. A* **91** 022121

[96] Lorenzo S, McCloskey R, Ciccarello F, Paternostro M and Palma G M 2015b *Phys. Rev. Lett.* **115** 120403

[97] Ciccarello F 2017 *Quantum Meas. Quantum Metrol.* **4** 53

[98] Pezzutto M, Paternostro M and Omar Y 2016 *New J. Phys.* **18** 123018

[99] Pezzutto M, Paternostro M and Omar Y 2018 arXiv:1806.10075

[100] Cusumano S, Cavina V, Keck M, De Pasquale A and Giovannetti V 2018 *Phys. Rev. A* **98** 032119

[101] Campbell S, Ciccarello F, Palma G M and Vacchini B 2018 *Phys. Rev. A* **98** 012142

[102] Gross J A, Caves C M, Milburn G J and Combes J 2018 *Quantum Sci. Technol.* **3** 024005

[103] Jin J, Giovannetti V, Fazio R, Sciarrino F, Mataloni P, Crespi A and Osellame R 2015 *Phys. Rev. A* **91** 012122

[104] Campisi M and Goold J 2017 *Phys. Rev. E* **95** 1

[105] Garttner M, Bohnet J G, Safavi-Naini A, Wall M L, Bollinger J J and Rey A M 2017 *Nat. Phys.* **13** 781

[106] Katz S, Lebowitz J L and Spohn H 1983 *Phys. Rev. B* **28** 1655

[107] Wang Y D and Clerk A A 2013 *Phys. Rev. Lett.* **110** 253601

[108] Woolley M J and Clerk A A 2014 *Phys. Rev. A* **89** 063805