Markovian master equations for quantum thermal machines: local versus global approach

Master equations are a powerful tool to study open quantum systems [1, 2]. They allow for a description of the relevant degrees of freedom only, which evolve under the influence of all other degrees of freedom that are not of immediate interest. These other degrees of freedom are collectively called the environment. A particularly simple situation occurs when the system can be described by a time-local master equation with constant dissipation rates [3–6]. This results in Markovian evolution, where knowledge of the density matrix at a given time is sufficient to predict all future observables, which implies an environment that has no memory. Here we refer to this type of master equations as Markovian.

Compared to the complete problem of describing all the degrees of freedom of system and environment together, a Markovian master equation governing only the system degrees of freedom is an immense simplification. Such a drastic reduction of complexity usually comes at a price. In this case, the price comes in the form of strong approximations which are not always justified. Studying these approximations is thus of utmost importance and indeed, there is a large body of literature that addresses these issues [5, 7–29] (for a recent review, see [6]). However, a large number of these studies focus on an environment that drives the system towards equilibrium. The recent rise in interest in out-of-equilibrium quantum systems, and in particular in quantum thermodynamics, calls for revisiting the question of the validity of the widely used Markovian quantum master equations [15, 16, 18, 20–26].

Quantum thermal machines are devices that perform useful tasks by exploiting thermal gradients in the environment; for recent reviews, see e.g. [30–35]. This task can for instance be the production of work [36–41], or more concretely of an electrical current [34, 42–44], the refrigeration of a quantum degree of freedom [45–50], the creation of entanglement [51–53], the determination of low temperatures [54], or the design of thermal transistors [55] and autonomous quantum clocks [56].

The standard description of these systems crucially relies on Markovian master equations to predict the relevant observables, such as heat currents and power. Two main approaches are followed in the literature. The first is a local approach, where the thermal baths couple locally to sub-systems of the machine. The second is a global approach, where thermal baths couple to the global eigenmodes of the machine. As the two approaches are conceptually different, there has been considerable debate about which one should be used in order to accurately describe thermal machines, and more generally out-of-equilibrium systems. Since the global approach describes equilibrium situations accurately (see below), while the local in some cases does not, there has been incentive to use the global approach out of equilibrium as well. Furthermore, the local approach is often believed to be more phenomenological in nature [13, 14, 19, 28, 57] and it was even argued that it is unphysical in certain regimes [27, 58, 59].

Zoom In Zoom Out Reset image size

The goal of the present work is to discuss these questions in depth. We will consider a system for which the full unitary dynamics of the machine and the thermal baths can be solved exactly. This allows us to evaluate the performance of local and global master equations for the machine against the exact dynamics. In addition, we give detailed derivations of the local and the global approaches and discuss the involved approximations. Specifically, we consider a heat engine introduced by Kosloff [36], which can be implemented in superconducting circuits [44]. The machine consists of two sub-systems (oscillators), which couple to different thermal baths, and to each other via an energy conserving interaction. In case the two oscillators have different frequencies, the machine requires an external driving field, making the Hamiltonian time-dependent. The entire system (machine plus baths) consists only of harmonic oscillators with quadratic interactions. Therefore, the system can be described exactly at the level of covariance matrices, and can be treated numerically even for relatively large baths with arbitrary precision. This exact numerical solution serves as a benchmark for evaluating the performance of both the local and global master equations. Focusing on an out-of-equilibrium steady-state regime, we discuss relevant thermodynamical observables, such as heat currents and power, as well as the quantum state (density matrix) of the machine degrees of freedom.

Our results demonstrate an overall excellent agreement between the predictions of the local approach and the exact solution. In the weak inter-system coupling regime, we see that the local approach provides an accurate description of the system, capturing both the thermodynamical observables and the quantum state. On the contrary, the global approach fails in this regime. Moving to the regime of intermediate coupling, we find that both approaches provide good descriptions of the system. Notably, the local approach still reliably captures all thermodynamical features of the machine. For very strong inter-system coupling strengths, the local approach starts to fail while the global approach still yields a faithful description of the system. We provide a detailed derivation of the regime of validity for each approach.

In the final part of the paper, we briefly discuss the case of finite dimensional machines. In particular, we consider the two-qubit entangler of [51], which is analogous to the heat engine setup considered in the first part, but with the two machine oscillators replaced by two qubits with equal level spacing (i.e. no external drive). While solving the total system (including the baths) exactly is unfortunately out of reach in this case, we can still compare the local and global approaches. We find very similar behavior to the results of the first part. In particular, the global approach still fails in the weak coupling regime, while for intermediate coupling, the two approaches agree well.

The paper is structured as follows. Section 2 gives a more detailed introduction to the local and global approaches. Sections 3–7 are devoted to the heat engine. In section 3, we introduce the system. The different master equations and the respective approximations are discussed in section 4, and the exact numerics are discussed in section 5. The observables which are investigated are introduced in section 6 and the results are given in section 7. The qubit entangler is then discussed in section 8 before we conclude in section 9.

Before going into details, we provide a short introduction to the two commonly used Markovian master equations which we discuss. In the local approach, the thermal baths couple to the eigenstates of sub-systems of the machine, while the global approach is based on a secular approximation. For time-independent Hamiltonians, the global approach corresponds to baths that couple to the delocalized eigenstates of the system Hamiltonian. In an equilibrium situation (i.e. baths at equal temperatures and time-independent Hamiltonian) the global master equation results in the desired steady state which is given by a Gibbs state with respect to the system Hamiltonian [1, 60] (for a discussion on deviations from the Gibbs state due to the finite coupling between system and bath, see, e.g., [61, 62]). The local approach on the other hand results in a product of Gibbs states with respect to the sub-system Hamiltonians [60]. For finite interactions, the global master equation is therefore usually considered superior to the local master equation. The situation drastically changes if we move away from equilibrium. Clearly, if the sub-systems do not interact, each sub-system should thermalize to its respective bath. While the local master equation yields correct results in this case, the global approach fails, resulting in a finite energy current through the system even in the absence of an interaction. There therefore exist limiting cases where we expect one approach to clearly be superior to the other. The goal of the present work is to connect these dots by comparing the local and the global master equations to exact numerics for a wide range of parameters in and out of equilibrium. Our main interest thus lies in the observables related to the machine operation, i.e. heat currents, powers, and efficiencies. Moreover, we also compare the quantum states of the machine.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We note that the validity of local and global approaches has been investigated before, even for out-of-equilibrium systems. Reference [16] provides a detailed study of the different approximations, however, the authors only consider the state of the system and not the energy flows. Reference [11] shows that the global approach neglects terms that influence the current through a two-terminal electric conductor. In [15], it was shown that the global approach can erroneously result in a vanishing heat current through a spin chain while a local approach shows good agreement with results obtained from a Redfield equation. The validity of the Redfield equation was shown in [23], where a system analogous to our heat engine in the absence of an external drive is investigated. References [58, 59, 63] argue that a local master equation can violate the second law of thermodynamics when a non-energy preserving interaction between the sub-systems is considered. These violations were however shown to be of the order of terms that are dropped when deriving a local master equation [22, 64]. In [22], it was shown that in the absence of degenerate subspaces, the local approach can be understood as the zeroth order of a perturbation series in the inter-system interaction. However, most thermal machines cited above crucially rely on such degenerate subspaces. Finally, [20] investigates heat currents through a system that is equivalent to our qubit entangler, without providing a benchmark. Their results agree with the results presented in section 8. All these previous works motivate our detailed investigation which compares the different master equations for a wide range of parameters.

The heat engine we consider [36, 44] is sketched in figure 1 and consists of two harmonic oscillators, with frequencies ${{\rm{\Omega }}}_{c}$ and ${{\rm{\Omega }}}_{h}$, which each couple to a bath. Here the subscript labels both the oscillators as well as the corresponding bath (where c stands for cold and h for hot). In this setup, the oscillators constitute the weakly interacting sub-systems mentioned above. The frequencies of the oscillators differ from each other by

$$\begin{eqnarray}&&{ \mathcal E }={{\rm{\Omega }}}_{h}-{{\rm{\Omega }}}_{c}\geqslant 0.\end{eqnarray} \tag{ 1 }$$

In order to extract power from the machine, we consider a time-dependent external field with frequency ${ \mathcal E }$ which mediates a coupling between the harmonic oscillators. The purpose of the heat engine is then to use a heat flow from the hot bath to the cold bath in order to increase the power of this external field. In [44], this external field is provided by a voltage and the power is directly related to an electrical current flowing against the voltage.

The total Hamiltonian of the heat engine (including heat baths) can then be written as (throughout this paper we set ${\hslash }=1$)

$$\begin{eqnarray}{\hat{H}}^{l}(t) & = & {\hat{H}}_{S}^{l}(t)+\displaystyle \sum _{\alpha =c,h}[{\hat{H}}_{\alpha }^{l}+{\hat{V}}_{\alpha }],\\ {\hat{H}}_{S}^{l}(t) & = & \displaystyle \sum _{\alpha =h,c}{{\rm{\Omega }}}_{\alpha }{\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }+g({\hat{a}}_{c}^{\dagger }{\hat{a}}_{h}{{\rm{e}}}^{{\rm{i}}{ \mathcal E }t}+{\rm{h}}{\rm{.}}{\rm{c}}.),\\ {\hat{H}}_{\alpha }^{l} & = & \displaystyle \sum _{k}{\omega }_{k,\alpha }{\hat{b}}_{k,\alpha }^{\dagger }{\hat{b}}_{k,\alpha },\hspace{1cm}{\hat{V}}_{\alpha }=\displaystyle \sum _{k}{\gamma }_{k,\alpha }({\hat{b}}_{k,\alpha }^{\dagger }{\hat{a}}_{\alpha }+{\hat{a}}_{\alpha }^{\dagger }{\hat{b}}_{k,\alpha }),\end{eqnarray} \tag{ 2 }$$

where ${\hat{a}}_{\alpha }$ (${\hat{b}}_{k,\alpha }$) denote annihilation operators of the system (baths), g denotes the interaction strength between the harmonic oscillators of the system, ${\omega }_{k,\alpha }$ are the frequencies of the bath modes, and ${\gamma }_{k,\alpha }$ denote the interaction strengths between system and bath modes. The superscript l denotes the laboratory frame. The external drive accounts for the energy that is needed to convert a photon with frequency ${{\rm{\Omega }}}_{c}$ into a photon with frequency ${{\rm{\Omega }}}_{h}$ and vice versa. In contrast to [58, 59, 63], which lack an external field, we find that the presence of such a field ensures that the local master equation does not result in any violations of the laws of thermodynamics.

In the following, we will work in a rotating frame defined by the transformation

$$\begin{eqnarray}&&{\hat{U}}_{r}(t)=\exp \left[{\rm{i}}{t}\displaystyle \sum _{\alpha =h,c}{{\rm{\Omega }}}_{\alpha }\left({\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }+\displaystyle \sum _{k}{\hat{b}}_{k,\alpha }^{\dagger }{\hat{b}}_{k,\alpha }\right)\right].\end{eqnarray} \tag{ 3 }$$

This results in the time-independent Hamiltonian

$$\begin{eqnarray}\hat{H} & = & {\hat{U}}_{r}(t){\hat{H}}^{l}(t){\hat{U}}_{r}^{\dagger }(t)-{\rm{i}}{\hat{U}}_{r}(t){\partial }_{t}{\hat{U}}_{r}^{\dagger }(t)={\hat{H}}_{S}+\displaystyle \sum _{\alpha =c,h}[{\hat{H}}_{\alpha }+{\hat{V}}_{\alpha }],\\ {\hat{H}}_{S} & = & g({\hat{a}}_{c}^{\dagger }{\hat{a}}_{h}+{\rm{h}}{\rm{.}}{\rm{c}}.),\hspace{1cm}{\hat{H}}_{\alpha }=\displaystyle \sum _{k}({\omega }_{k,\alpha }-{{\rm{\Omega }}}_{\alpha }){\hat{b}}_{k,\alpha }^{\dagger }{\hat{b}}_{k,\alpha }.\end{eqnarray} \tag{ 4 }$$

Note that the interaction terms ${\hat{V}}_{\alpha }$ are invariant under the transformation.

If not explicitly stated otherwise, all equations are given in the rotating frame.

In this section, we consider different Markovian master equations which are used to describe the evolution of the reduced density matrix of the system. These equations allow for analytic expressions of observables such as power, heat, and efficiency which are then compared to the ones obtained from the exact dynamics. The standard way of deriving Markovian master equations is to first perform Born–Markov approximations. This procedure is discussed in detail elsewhere (see for instance [1, 6, 16, 65]). We therefore only summarize the approximations and give the resulting expression for the system under consideration. The approximations are:

• Born approximation: treating ${\hat{V}}_{\alpha }$ perturbatively to lowest order,
• Markov approximation: assuming invariance of $\tilde{\rho }(t)$ on time-scales of the order of ${\tau }_{B}$,

where $\tilde{\rho }(t)$ is the reduced density matrix in the interaction picture and ${\tau }_{B}$ denotes the bath-correlation time that will be introduced below. In the interaction picture (and the rotating frame), the Born–Markov approximations result in the following master equation,

$$\begin{eqnarray}{\partial }_{t}\tilde{\rho }(t) & = & \displaystyle \sum _{\alpha =h,c}{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau \{{C}_{1,\alpha }(\tau )[{\tilde{a}}_{\alpha }^{\dagger }(t-\tau )\tilde{\rho }(t){\tilde{a}}_{\alpha }(t)-{\tilde{a}}_{\alpha }(t){\tilde{a}}_{\alpha }^{\dagger }(t-\tau )\tilde{\rho }(t)]\\ & & +{C}_{2,\alpha }(\tau )[{\tilde{a}}_{\alpha }(t-\tau )\tilde{\rho }(t){\tilde{a}}_{\alpha }^{\dagger }(t)-{\tilde{a}}_{\alpha }^{\dagger }(t){\tilde{a}}_{\alpha }(t-\tau )\tilde{\rho }(t)]\}+{\rm{h}}{\rm{.}}{\rm{c}}.\end{eqnarray} \tag{ 5 }$$

Here operators in the interaction picture are given by

$$\begin{eqnarray}&&\tilde{A}(t)={\hat{U}}_{S}^{\dagger }(t)\hat{A}{\hat{U}}_{S}(t),\end{eqnarray} \tag{ 6 }$$

where $\hat{A}$ denotes the operator in the Schrödinger picture and

$$\begin{eqnarray}&&{\hat{U}}_{S}(t)={{\rm{e}}}^{-{\rm{i}}{\hat{H}}_{S}t}.\end{eqnarray} \tag{ 7 }$$

We further introduced the bath correlation functions

$$\begin{eqnarray}{C}_{1,\alpha }(\tau ) & = & \displaystyle \sum _{k}{\gamma }_{k,\alpha }^{2}{n}_{B}^{\alpha }({\omega }_{k,\alpha }){{\rm{e}}}^{{\rm{i}}({\omega }_{k,\alpha }-{{\rm{\Omega }}}_{\alpha })\tau },\\ {C}_{2,\alpha }(\tau ) & = & \displaystyle \sum _{k}{\gamma }_{k,\alpha }^{2}[1+{n}_{B}^{\alpha }({\omega }_{k,\alpha })]{{\rm{e}}}^{-{\rm{i}}({\omega }_{k,\alpha }-{{\rm{\Omega }}}_{\alpha })\tau },\end{eqnarray} \tag{ 8 }$$

where we assumed the baths to be in thermal states. The Bose–Einstein distribution is given by

$$\begin{eqnarray}&&{n}_{B}^{\alpha }(\omega )=\displaystyle \frac{1}{{{\rm{e}}}^{\omega /({k}_{B}{T}_{\alpha })}-1}.\end{eqnarray} \tag{ 9 }$$

The bath correlation functions are usually peaked around $\tau =0$ and decay for large times. This decay defines the bath-correlation time ${\tau }_{B}$ such that ${C}_{j,\alpha }({\tau }_{B})\ll {C}_{2,\alpha }(0)$ (note that ${C}_{1,\alpha }(0)$ vanishes as ${T}_{\alpha }\to 0$). For an explicit evaluation of the bath correlation functions and a discussion on the relevance of ${\tau }_{B}$ (including a discussion of the zero temperature limit), we refer the reader to [16].

The Born–Markov equation represents the starting point of our analysis. Since it does not guarantee positive evolution [1], further approximations are usually made to obtain a master equation in Gorini–Kossakowski–Sudarshan–Lindblad (GKSL) form [3, 4] ensuring completely positive dynamics. In the following, we will discuss two popular approximations, the local and the global approach, in some detail.

4.1. Local master equation

The Markov approximation that was made to obtain equation (5) is responsible for the fact that the density matrix under the integral is independent of τ. This approximation is valid as long as the characteristic time over which $\tilde{\rho }(t)$ varies is much larger than ${\tau }_{B}$. In the same spirit, we can make the approximation

$$\begin{eqnarray}&&{\tilde{a}}_{\alpha }(t-\tau )\simeq {\tilde{a}}_{\alpha }(t),\end{eqnarray} \tag{ 10 }$$

in the integral of equation (5). Note that we make this approximation in the rotating frame, where the fast oscillations with frequency ${{\rm{\Omega }}}_{\alpha }$ are encoded in the bath correlation functions (see equation (8)). In our model, we have

$$\begin{eqnarray}&&{\tilde{a}}_{\alpha }(t)={\hat{a}}_{\alpha }\cos ({gt})-{\rm{i}}{\hat{a}}_{\bar{\alpha }}\sin ({gt}),\end{eqnarray} \tag{ 11 }$$

where $\bar{\alpha }\ne \alpha$. The approximation in equation (10) is therefore expected to be good as long as $g{\tau }_{B}\ll 1$. For reasonably small values of g, this approximation is therefore completely consistent with the Markov approximation.

This approximation directly results in the local master equation (in the Schrödinger picture) without the need of a secular approximation

$$\begin{eqnarray}&&{\partial }_{t}\hat{\rho }(t)=-{\rm{i}}[{\bar{H}}_{S},\hat{\rho }(t)]+\displaystyle \sum _{\alpha =h,c}\{{\bar{{\rm{\Gamma }}}}_{\alpha }{ \mathcal D }[{\hat{a}}_{\alpha }]\hat{\rho }(t)+{{\rm{\Gamma }}}_{\alpha }{ \mathcal D }[{\hat{a}}_{\alpha }^{\dagger }]\hat{\rho }(t)\},\end{eqnarray} \tag{ 12 }$$

where we neglected a physically irrelevant constant and introduced

$$\begin{eqnarray}&&{\bar{{\rm{\Gamma }}}}_{\alpha }={\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha })[{n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha })+1],\hspace{1cm}{{\rm{\Gamma }}}_{\alpha }={\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha }){n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha }).\end{eqnarray} \tag{ 13 }$$

as well as the Lindblad superoperators

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

and the energy damping rate

$$\begin{eqnarray}&&{\kappa }_{\alpha }(\omega )=2\pi \displaystyle \sum _{k}{\gamma }_{k,\alpha }^{2}\delta (\omega -{\omega }_{k,\alpha })=2\pi {\rho }_{\alpha }(\omega ),\end{eqnarray} \tag{ 15 }$$

where ${\rho }_{\alpha }(\omega )$ denotes the spectral density. The renormalized Hamiltonian is given by ${\bar{H}}_{S}={\hat{H}}_{S}+{\sum }_{\alpha =h,c}{{\rm{\Sigma }}}_{\alpha }{\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }$ with

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\alpha }=P{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \frac{{\rho }_{\alpha }(\omega )}{{{\rm{\Omega }}}_{\alpha }-\omega },\end{eqnarray} \tag{ 16 }$$

where P denotes the Cauchy principal value. We note that this renormalization should be small in order for the Born–Markov approximations to be valid (this can be understood by noting that $\tilde{\rho }(t)$ will have terms that oscillate with frequency ${{\rm{\Sigma }}}_{\alpha }$). In the following, we will thus neglect this renormalization. As shown in section 7, excellent agreement between the local master equation and exact numerics is found for a wide range of parameters without taking into account the renormalization of the Hamiltonian; see also [16].

As discussed above, the local master equation is justified as long as $g{\tau }_{B}\ll 1$. With the help of equations (8), this temporal inequality can be translated into an inequality involving the Fourier transform of the bath correlation functions. If ${\tau }_{B}g\ll 1$, then we have $\exp ({\rm{i}}{g}\tau ){C}_{j,\alpha }(\tau )\simeq {C}_{j,\alpha }(\tau )$ for all τ since we can approximate ${C}_{j,\alpha }(\tau \gt {\tau }_{B})\simeq 0$. It is straightforward to show that this is fulfilled as long as the Fourier transform of the bath correlation functions can be approximated as constant over the energy scale of g. With the help of equations (8) and (15), this results in the inequalities

$$\begin{eqnarray}&&| {n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha }\pm g)-{n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha })| \ll 1\,\mathrm{or}\,{n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha }),\\ &&\quad | {\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha }\pm g)-{\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha })| \ll {\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha }).\end{eqnarray} \tag{ 17 }$$

Here we used the fact that the main contribution for the bath correlation functions comes from the low-frequency terms For a bosonic bath with Ohmic spectrum, the last equations are fulfilled as long as $g\ll {{\rm{\Omega }}}_{\alpha }$, which is usually the case for systems where the local approach is employed [38, 44, 46, 51]. For fermionic baths, equations analogous to equations (17) can be derived. These are usually not fulfilled at low temperatures due to the step-like behavior of the Fermi–Dirac distribution [23, 66]. Note that the validity of the Markov approximation, which underlies both the local as well as the global master equation, requires conditions obtained from equation (17) by exchanging g with ${\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha })$.

It is sometimes stated that the local approach is only valid for interaction strengths much smaller than the induced broadening $g\ll {\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha })$ (see for instance [23]). As we show in section 7, the local approach gives reliable predictions even for interactions that are several times the broadening. This is in complete agreement with equation (17).

We note that the local master equation is also obtained in the so-called singular coupling limit [1, 5, 10, 12], where the bath correlation functions in equations (8) tend to a delta function. This limit is often dismissed as being unrealistic. However we stress that the bath correlation functions only have to behave like delta functions on time-scales of the order of $1/g$ for the local master equation to be valid. In the weak coupling limit, which is often the regime of interest for thermal machines, $1/g$ can naturally be much bigger than ${\tau }_{B}$.

In order to solve the master equation, we make use of its bi-linearity (in creation and annihilation operators) which implies that a Gaussian state remains Gaussian at all times. Furthermore, there are no terms in the master equation which result in a displacement of the state. We can therefore restrict the analysis to states which have $\langle {\hat{a}}_{\alpha }\rangle =0$. Then the state is fully described by its covariance matrix. From the local master equation in equation (12), one can derive the following differential equations for the covariance matrix elements

$$\begin{eqnarray}{\partial }_{t}\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{h}\rangle & = & 2g\mathrm{Im}\{\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle \}+{\kappa }_{h}({n}_{B}^{h}-\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{h}\rangle ),\\ {\partial }_{t}\langle {\hat{a}}_{c}^{\dagger }{\hat{a}}_{c}\rangle & = & -2g\mathrm{Im}\{\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle \}+{\kappa }_{c}({n}_{B}^{c}-\langle {\hat{a}}_{c}^{\dagger }{\hat{a}}_{c}\rangle ),\\ {\partial }_{t}\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle & = & -\displaystyle \frac{({\kappa }_{c}+{\kappa }_{h})}{2}\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle -{\rm{i}}{g}(\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{h}\rangle -\langle {\hat{a}}_{c}^{\dagger }{\hat{a}}_{c}\rangle ),\end{eqnarray} \tag{ 18 }$$

where

$$\begin{eqnarray}&&{\kappa }_{\alpha }={\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha }),\hspace{1cm}{n}_{B}^{\alpha }={n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha }).\end{eqnarray} \tag{ 19 }$$

In the steady state, these differential equations are solved by the time-independent expressions

$$\begin{eqnarray}\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{h}\rangle & = & {n}_{B}^{h}-\displaystyle \frac{4{g}^{2}{\kappa }_{c}({n}_{B}^{h}-{n}_{B}^{c})}{({\kappa }_{h}+{\kappa }_{c})({\kappa }_{c}{\kappa }_{h}+4{g}^{2})},\\ \langle {\hat{a}}_{c}^{\dagger }{\hat{a}}_{c}\rangle & = & {n}_{B}^{c}+\displaystyle \frac{4{g}^{2}{\kappa }_{h}({n}_{B}^{h}-{n}_{B}^{c})}{({\kappa }_{h}+{\kappa }_{c})({\kappa }_{c}{\kappa }_{h}+4{g}^{2})},\\ \langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle & = & \displaystyle \frac{-{\rm{i}}2g{\kappa }_{c}{\kappa }_{h}({n}_{B}^{h}-{n}_{B}^{c})}{({\kappa }_{h}+{\kappa }_{c})({\kappa }_{c}{\kappa }_{h}+4{g}^{2})}.\end{eqnarray} \tag{ 20 }$$

4.2. Global master equation

The second approximation we consider is the secular approximation. This approximation consists of dropping all the terms in the Born–Markov master equation (see equation (5)) which oscillate as a function of time t. The secular approximation is expected to hold over time-scales much bigger than the inverse frequencies of the oscillating terms and obviously becomes better as these frequencies increase. In order to identify the oscillating terms, we write the interaction picture operators in equation (11) as

$$\begin{eqnarray}&&{\tilde{a}}_{h}(t)=\displaystyle \frac{1}{\sqrt{2}}({\hat{a}}_{+}{{\rm{e}}}^{-{\rm{i}}{g}{t}}+{\hat{a}}_{-}{{\rm{e}}}^{{\rm{i}}{g}{t}}),\hspace{.75cm}{\tilde{a}}_{c}(t)=\displaystyle \frac{1}{\sqrt{2}}({\hat{a}}_{+}{{\rm{e}}}^{-{\rm{i}}{g}{t}}-{\hat{a}}_{-}{{\rm{e}}}^{{\rm{i}}{g}{t}}),\end{eqnarray} \tag{ 21 }$$

where we introduced the operators

$$\begin{eqnarray}&&{\hat{a}}_{\pm }=\displaystyle \frac{1}{\sqrt{2}}({\hat{a}}_{h}\pm {\hat{a}}_{c}).\end{eqnarray} \tag{ 22 }$$

The secular approximation is obtained by plugging equation (21) into equation (5) and dropping all terms that oscillate with $\exp [2{\rm{i}}{g}{t}]$. This results in the master equation

$$\begin{eqnarray}&&{\partial }_{t}\hat{\rho }(t)=-{\rm{i}}[{\bar{H}}_{S},\hat{\rho }(t)]+\displaystyle \frac{1}{2}\displaystyle \sum _{\sigma =\pm }\{{\bar{{\rm{\Gamma }}}}_{\sigma }{ \mathcal D }[{\hat{a}}_{\sigma }]\hat{\rho }(t)+{{\rm{\Gamma }}}_{\sigma }{ \mathcal D }[{\hat{a}}_{\sigma }^{\dagger }]\hat{\rho }(t)\},\end{eqnarray} \tag{ 23 }$$

where

$$\begin{eqnarray}&&{\bar{{\rm{\Gamma }}}}_{\sigma }=\displaystyle \sum _{\alpha =h,c}{\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma })[{n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma })+1],\hspace{.75cm}{{\rm{\Gamma }}}_{\sigma }=\displaystyle \sum _{\alpha =h,c}{\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma }){n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma }),\end{eqnarray} \tag{ 24 }$$

with the frequencies

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{\alpha ,\pm }={{\rm{\Omega }}}_{\alpha }\pm g.\end{eqnarray} \tag{ 25 }$$

The renormalized Hamiltonian reads

$$\begin{eqnarray}&&{\bar{H}}_{S}=\displaystyle \sum _{\sigma =\pm }(\sigma g+{{\rm{\Sigma }}}_{h,\sigma }+{{\rm{\Sigma }}}_{c,\sigma }){\hat{a}}_{\sigma }^{\dagger }{\hat{a}}_{\sigma },\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\alpha ,\sigma }=\displaystyle \frac{1}{2}P{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \frac{{\rho }_{\alpha }(\omega )}{{{\rm{\Omega }}}_{\alpha ,\sigma }-\omega }.\end{eqnarray} \tag{ 27 }$$

As for the local approach, we will neglect the renormalization of the Hamiltonian in the following.

The secular approximation results in a master equation where the eigenmodes of the Hamiltonian ${\hat{a}}_{\sigma }$ are the relevant degrees of freedom. The action of the baths decouples into a bath for each eigenmode, which drives the mode towards thermal equilibrium characterized by the occupation numbers

$$\begin{eqnarray}&&{n}_{\sigma }=\displaystyle \frac{{\kappa }_{h,\sigma }{n}_{B}^{h,\sigma }+{\kappa }_{c,\sigma }{n}_{B}^{c,\sigma }}{{\kappa }_{h,\sigma }+{\kappa }_{c,\sigma }},\end{eqnarray} \tag{ 28 }$$

where we introduced

$$\begin{eqnarray}&&{\kappa }_{\alpha ,\sigma }={\kappa }_{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma }),\hspace{1cm}{n}_{B}^{\alpha ,\sigma }={n}_{B}^{\alpha }({{\rm{\Omega }}}_{\alpha ,\sigma }).\end{eqnarray} \tag{ 29 }$$

We note that in the limit $g\to 0$, the local master equation is no longer recovered. Because the secular approximation is no longer justified in this limit, we expect the global master equation to break down. This is also consistent with the fact that we expect the secular approximation to be valid as long as the frequency of the neglected oscillating terms is much bigger than the linewidths, i.e. ${\kappa }_{\alpha }\ll g$. Since the Markov approximation requires ${\kappa }_{\alpha }\ll {{\rm{\Omega }}}_{\alpha }$, and the local approach is valid for $g\ll {{\rm{\Omega }}}_{\alpha }$, there is an overlap between the regimes of validity of the local and the global master equation. This implies that the local and the global master equation together are enough to describe the system for all parameters which allow for a Markovian master equation.

In the global approach, the covariance matrix is governed by the differential equations

$$\begin{eqnarray}{\partial }_{t}\langle {\hat{a}}_{\sigma }^{\dagger }{\hat{a}}_{\sigma }\rangle & = & \displaystyle \frac{1}{2}\displaystyle \sum _{\alpha =h,c}{\kappa }_{\alpha ,\sigma }[{n}_{B}^{\alpha ,\sigma }-\langle {\hat{a}}_{\sigma }^{\dagger }{\hat{a}}_{\sigma }\rangle ],\\ {\partial }_{t}\langle {\hat{a}}_{+}^{\dagger }{\hat{a}}_{-}\rangle & = & \left[2{\rm{i}}{g}-\displaystyle \frac{1}{4}\displaystyle \sum _{\alpha ,\sigma }{\kappa }_{\alpha ,\sigma }\right]\langle {\hat{a}}_{+}^{\dagger }{\hat{a}}_{-}\rangle .\end{eqnarray} \tag{ 30 }$$

In the steady state, we find that the state is a product of thermal states with occupation numbers $\langle {\hat{a}}_{\pm }^{\dagger }{\hat{a}}_{\pm }\rangle ={n}_{\sigma }$ as expected.

In this section we briefly describe how we obtain exact numerics which are used as a benchmark when comparing the different master equations. To this end, we simulate the unitary evolution generated by the Hamiltonian in equation (2) for big but finite baths. The key element here is that the Hamiltonian in equation (2) is quadratic, so that Gaussian states (such as thermal states) remain Gaussian throughout the whole evolution. As such, they can be fully characterized by only their first and second moments (see e.g. [67]). This allows us to characterize the time-evolved state of the whole system, including the thermal baths, by a matrix of size $\sim N$, where N is the total number of oscillators involved.

We consider baths made up of $n+1$ harmonic oscillators, so that the total size of system and baths is $N=2(n+2)$. The bath modes are chosen to be uniformly spread over a range $(0,{\omega }_{{c}})$, where ${\omega }_{{c}}$ is a cutoff frequency. That is,

$$\begin{eqnarray}&&{\omega }_{k,\alpha }=\displaystyle \frac{k}{n}{\omega }_{{c}},\end{eqnarray} \tag{ 31 }$$

for $k=0,\mathrm{..},\,n$. This defines ${\hat{H}}_{\alpha }$ in equation (2). Let us now turn our attention to ${\hat{V}}_{\alpha }$, and hence to the couplings ${\gamma }_{k,\alpha }$. First note that the action of the baths in the Markovian master equations, equations (14) and (23), is captured by the spectral density given in equation (15). It is then common to use an ad hoc form for the spectral density in the continuum limit, instead of specifying the coupling constants ${\gamma }_{k,\alpha }$. A common choice for the spectral density is an Ohmic spectrum,

$$\begin{eqnarray}&&{\rho }_{\alpha }(\omega )\propto \omega ,\end{eqnarray} \tag{ 32 }$$

which holds for low frequencies, $\omega \leqslant {\omega }_{{c}}$. In order to relate this approach to the ${\gamma }_{k,\alpha }$'s, which are necessary to simulate the full Hamiltonian in the finite-baths scenario, let us integrate equation (15), obtaining,

$$\begin{eqnarray}&&{\int }_{0}^{{\omega }_{{c}}}{\rho }_{\alpha }(\omega ){\rm{d}}\omega =\displaystyle \sum _{k=0}^{n}{\gamma }_{k,\alpha }^{2}.\end{eqnarray} \tag{ 33 }$$

It is now convenient to discretize the integral,

$$\begin{eqnarray}&&{\int }_{0}^{{\omega }_{{c}}}{\rho }_{\alpha }(\omega ){\rm{d}}\omega \approx \displaystyle \sum _{k=0}^{n}{\rho }_{\alpha }({\omega }_{k,\alpha })\displaystyle \frac{{\omega }_{{c}}}{n},\end{eqnarray} \tag{ 34 }$$

from where it immediately follows,

$$\begin{eqnarray}&&{\gamma }_{k,\alpha }^{2}\approx {\rho }_{\alpha }({\omega }_{k,\alpha })\displaystyle \frac{{\omega }_{c}}{n},\end{eqnarray} \tag{ 35 }$$

which becomes increasingly accurate with increasing n. In the particular case of an Ohmic distribution, we obtain,

$$\begin{eqnarray}&&{\gamma }_{k,\alpha }\propto \sqrt{k}\displaystyle \frac{{\omega }_{c}}{n}.\end{eqnarray} \tag{ 36 }$$

Equation (35), together with equation (31), provides a simple recipe for building the discrete version of the Hamiltonian in equation (2) for a given spectral density. The specific choice of ${\omega }_{c}$ and n for our simulations, as well as the dependence of the results on this choice, is discussed in the appendix section.

The initial state for the simulations is taken to be of the form,

$$\begin{eqnarray}&&{\hat{\rho }}_{0}={\hat{\tau }}_{{\beta }_{h}}\otimes {\hat{\rho }}_{S}\otimes {\hat{\tau }}_{{\beta }_{c}}\end{eqnarray} \tag{ 37 }$$

where ${\hat{\tau }}_{{\beta }_{\alpha }}$ are thermal states,

$$\begin{eqnarray}&&{\hat{\tau }}_{{\beta }_{\alpha }}=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }{\hat{H}}_{\alpha }^{l}}}{{{ \mathcal Z }}_{B,\alpha }},\end{eqnarray} \tag{ 38 }$$

with inverse temperature ${\beta }_{\alpha }=1/({k}_{B}{T}_{\alpha })$ and we note that ${\hat{H}}_{\alpha }^{l}$ is in the laboratory frame. Here, the initial state of the machine ${\hat{\rho }}_{S}$ can be an arbitrary Gaussian state. In our simulations, we take ${\hat{\rho }}_{S}={{\rm{e}}}^{-{\beta }_{h}{{\rm{\Omega }}}_{h}{\hat{a}}_{h}^{\dagger }{\hat{a}}_{h}}/{{ \mathcal Z }}_{h}$ $\otimes \,{{\rm{e}}}^{-{\beta }_{c}{{\rm{\Omega }}}_{c}{\hat{a}}_{c}^{\dagger }{\hat{a}}_{c}}/{{ \mathcal Z }}_{c}$. We note that with this choice, the state given in equation (37) is Gaussian.

Once the Hamiltonian and the initial state are defined, we consider the closed unitary dynamics of the full compound under the Hamiltonian in equation (2) (it is convenient to work in the rotating frame, where the Hamiltonian is time independent). The dynamics can be derived by considering the Heisenberg equations of motion of the operators ${\hat{a}}_{\alpha },{\hat{a}}_{\alpha }^{\dagger },\{{\hat{b}}_{k,\alpha },{\hat{b}}_{k,\alpha }^{\dagger }\}$. For details on the derivation, we refer the reader to [16].

In order to simulate the steady state using finite degrees of freedom, one needs to let the whole compound evolve for a time t that satisfies ${\tau }_{\mathrm{eq}}\ll t\ll {\tau }_{\mathrm{rec}}$, where ${\tau }_{\mathrm{eq}}$ is the equilibration time of the system and ${\tau }_{\mathrm{rec}}$ is the recurrence time of the bath. From equations (12) and (18), we infer the equilibration time to be $1/{\tau }_{\mathrm{eq}}\approx \max \{{\kappa }_{h},{\kappa }_{c}\}$ (see also [68], where the equilibration time is discussed explicitly for finite systems). On the other hand, the recurrence time scales linearly with the number of oscillators in the bath [16]. Hence, by taking a sufficiently large bath (in the simulations we take ∼400 oscillators), we can ensure that ${\tau }_{\mathrm{eq}}\ll {\tau }_{\mathrm{rec}}$. In our simulations, we take $t\approx 20{\tau }_{\mathrm{eq}}$.

In this work, we are particularly interested in the energy flows that traverse the quantum thermal machine in a non-equilibrium situation. However, to compare to previous works [16], and to further assess the validity of the different master equations, we also consider the obtained steady states.

6.1. Heat currents, power, and efficiency

6.1.1. Local master equation

As we consider a heat engine, the main quantity of interest is the power that is produced. For our system, it is defined as [33]

$$\begin{eqnarray}&&P=-\mathrm{Tr}\{[{\partial }_{t}{\hat{H}}_{S}^{l}(t)]{\hat{\rho }}^{l}(t)\}=-2g{ \mathcal E }\mathrm{Im}\{\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}\rangle \},\end{eqnarray} \tag{ 39 }$$

where the superscript denotes the laboratory frame and $\langle \cdots \rangle$, denotes the ensemble average in the rotating frame. Note that positive power implies that energy leaves the system. In addition to the power, we consider the heat currents that enter (or leave) the system. To this end, we write

$$\begin{eqnarray}&&{\partial }_{t}{\hat{\rho }}^{l}=-{\rm{i}}[{\hat{H}}_{S}^{l}(t),{\hat{\rho }}^{l}(t)]+\displaystyle \sum _{\alpha =h,c}{{ \mathcal L }}_{\alpha }^{l}{\hat{\rho }}^{l}(t),\end{eqnarray} \tag{ 40 }$$

where ${{ \mathcal L }}_{\alpha }^{l}$ is a super-operator that groups all the dissipative terms which arise from bath α in the laboratory frame. Note that under the Born–Markov approximation, the dissipators of different baths can be added [69]. The heat currents are then defined as

$$\begin{eqnarray}&&{J}_{\alpha }=\mathrm{Tr}\{{\hat{H}}_{S}^{l}(t){{ \mathcal L }}_{\alpha }^{l}{\hat{\rho }}^{l}\}=\mathrm{Tr}\left\{\left[{\hat{H}}_{S}+\displaystyle \sum _{\alpha =h,c}{{\rm{\Omega }}}_{\alpha }{\hat{a}}_{\alpha }^{\dagger }\hat{a}\right]{{ \mathcal L }}_{\alpha }\hat{\rho }\right\},\end{eqnarray} \tag{ 41 }$$

where the dissipator in the rotating frame is related to the dissipator in the laboratory frame by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\alpha }^{l}{\hat{\rho }}^{l}(t)={\hat{U}}_{r}^{\dagger }(t)[{{ \mathcal L }}_{\alpha }\hat{\rho }(t)]{\hat{U}}_{r}(t).\end{eqnarray} \tag{ 42 }$$

With our sign convention, a positive heat current implies energy entering the system from the bath. In the steady state, the first law of thermodynamics thus reads

$$\begin{eqnarray}&&P={J}_{c}+{J}_{h}.\end{eqnarray} \tag{ 43 }$$

The efficiency is defined as the ratio of the obtained power, divided by the heat that originates from the hot bath

$$\begin{eqnarray}&&\eta =\displaystyle \frac{P}{{J}_{h}}.\end{eqnarray} \tag{ 44 }$$

In the regime where the system operates as a heat engine (i.e. $P\gt 0$ and ${J}_{h}\gt 0$), the second law of thermodynamics forces the efficiency to remain below the Carnot limit

$$\begin{eqnarray}&&\eta \lt 1-\displaystyle \frac{{T}_{c}}{{T}_{h}}={\eta }_{C}.\end{eqnarray} \tag{ 45 }$$

For the local master equation in equation (12), the heat currents in equation (41) can be written as

$$\begin{eqnarray}&&{J}_{\alpha }={\kappa }_{\alpha }\left[{{\rm{\Omega }}}_{\alpha }({n}_{B}^{\alpha }-\langle {\hat{a}}_{\alpha }^{\dagger }{\hat{a}}_{\alpha }\rangle )-\displaystyle \frac{g}{2}\langle {\hat{a}}_{h}^{\dagger }{\hat{a}}_{c}+{\hat{a}}_{c}^{\dagger }{\hat{a}}_{h}\rangle \right].\end{eqnarray} \tag{ 46 }$$

From equations (18), we can infer that the second term in the heat current decays exponentially in time. In the steady state, from equations (20) and (39), we find for the power

$$\begin{eqnarray}&&P=\displaystyle \frac{({{\rm{\Omega }}}_{h}-{{\rm{\Omega }}}_{c})4{g}^{2}{\kappa }_{c}{\kappa }_{h}({n}_{B}^{h}-{n}_{B}^{c})}{({\kappa }_{h}+{\kappa }_{c})({\kappa }_{c}{\kappa }_{h}+4{g}^{2})},\end{eqnarray} \tag{ 47 }$$

and the heat current

$$\begin{eqnarray}&&{J}_{h}={{\rm{\Omega }}}_{h}\displaystyle \frac{4{g}^{2}{\kappa }_{c}{\kappa }_{h}({n}_{B}^{h}-{n}_{B}^{c})}{({\kappa }_{h}+{\kappa }_{c})({\kappa }_{c}{\kappa }_{h}+4{g}^{2})},\end{eqnarray} \tag{ 48 }$$

resulting in the efficiency

$$\begin{eqnarray}&&\eta =1-\displaystyle \frac{{{\rm{\Omega }}}_{c}}{{{\rm{\Omega }}}_{h}}.\end{eqnarray} \tag{ 49 }$$

We note that this efficiency fulfills equation (45). When the frequencies are chosen such that the efficiency is above the Carnot efficiency, then we find $P\lt 0$ and ${J}_{h}\lt 0$. As long as $\eta \lt {\eta }_{C}$, our machine is thus a heat engine, $\eta ={\eta }_{C}$ denotes the point of reversibility, where all the energy currents vanish, and for $\eta \gt {\eta }_{C}$, the machine acts as a refrigerator, using power to induce a heat current from the cold bath to the hot bath [44]. Finally, we note that for ${{\rm{\Omega }}}_{h}={{\rm{\Omega }}}_{c}$, where there is no external power, heat always flows from the hot bath to the cold bath as dictated by the second law of thermodynamics. In contrast to models which consider non-energy preserving interactions, the local approach does not violate the laws of thermodynamics when including an external field that provides the energy to convert photons of frequency ${{\rm{\Omega }}}_{c}$ into photons of frequency ${{\rm{\Omega }}}_{h}$.

6.1.2. Global master equation

In the global master equation, the bath couples to global states which are dressed by the external field. Therefore, the dissipative terms include the external field and the definitions introduced in the last subsection are no longer valid. To define heat and work in the global approach, we follow [12, 31, 35, 70, 71]. To this end, we first write the dissipator in the global master equation (see equation (23)) as the sum of four dissipators given by

$$\begin{eqnarray}&&{{ \mathcal L }}_{\alpha ,\sigma }\hat{\rho }(t)=\displaystyle \frac{1}{2}{\kappa }_{\alpha ,\sigma }({n}_{B}^{\alpha ,\sigma }+1){ \mathcal D }[{\hat{a}}_{\sigma }]\hat{\rho }(t)+\displaystyle \frac{1}{2}{\kappa }_{\alpha ,\sigma }{n}_{B}^{\alpha ,\sigma }{ \mathcal D }[{\hat{a}}_{\sigma }^{\dagger }]\hat{\rho }(t).\end{eqnarray} \tag{ 50 }$$

We further introduce the thermal states

$$\begin{eqnarray}&&{\hat{\rho }}_{\alpha ,\sigma }=\displaystyle \frac{{{\rm{e}}}^{-{\beta }_{\alpha }{{\rm{\Omega }}}_{\alpha ,\sigma }{\hat{a}}_{\sigma }^{\dagger }{\hat{a}}_{\sigma }}}{{{ \mathcal Z }}_{\alpha ,\sigma }},\end{eqnarray} \tag{ 51 }$$

which constitute the steady states of the respective operators, i.e. ${{ \mathcal L }}_{\alpha ,\sigma }{\hat{\rho }}_{\alpha ,\sigma }=0$.

The heat currents in the steady state (denoted $\hat{\rho }$) are then defined as

$$\begin{eqnarray}&&{J}_{\alpha }=-{k}_{B}{T}_{\alpha }\displaystyle \sum _{\sigma =\pm }\mathrm{Tr}\{[{{ \mathcal L }}_{\alpha ,\sigma }\hat{\rho }]\mathrm{ln}{\hat{\rho }}_{\alpha ,\sigma }\},\end{eqnarray} \tag{ 52 }$$

and the power is given by the first law, see equation (43). We note that for the global approach, we enforce the first law while in the local approach, it follows from the expressions for heat currents and power.

An explicit calculation results in [70]

$$\begin{eqnarray}&&{J}_{h}=\displaystyle \frac{1}{2}\displaystyle \sum _{\sigma =\pm }{{\rm{\Omega }}}_{h,\sigma }\displaystyle \frac{{\kappa }_{h,\sigma }{\kappa }_{c,\sigma }}{{\kappa }_{h,\sigma }+{\kappa }_{c,\sigma }}({n}_{B}^{h,\sigma }-{n}_{B}^{c,\sigma }).\end{eqnarray} \tag{ 53 }$$

Note that if equations (17) are fulfilled, the last expression reduces to

$$\begin{eqnarray}&&{J}_{h}={{\rm{\Omega }}}_{h}\displaystyle \frac{{\kappa }_{h}{\kappa }_{c}}{{\kappa }_{h}+{\kappa }_{c}}({n}_{B}^{h}-{n}_{B}^{c}),\end{eqnarray} \tag{ 54 }$$

which is the expression obtained by the local approach in the limit $g\gg {\kappa }_{\alpha }$. As expected, if the approximations leading to both the local and the global master equation are justified, the two approaches give the same result. Whenever equations (17) are not fulfilled, the global approach implies that the heat engine makes use of two channels labeled by the subscript σ. Each of these channels has an efficiency

$$\begin{eqnarray}&&{\eta }_{\sigma }=\displaystyle \frac{{P}_{\sigma }}{{J}_{h,\sigma }}=\displaystyle \frac{{J}_{h,\sigma }+{J}_{c,\sigma }}{{J}_{h,\sigma }}=1-\displaystyle \frac{{{\rm{\Omega }}}_{c,\sigma }}{{{\rm{\Omega }}}_{h,\sigma }}\leqslant {\eta }_{C},\end{eqnarray} \tag{ 55 }$$

where ${J}_{\alpha }={J}_{\alpha ,+}+{J}_{\alpha ,-}$ and J_c is obtained from J_h by exchanging the labels $c\leftrightarrow h$ (see equation (53)). We focus on the heat engine regime where ${J}_{h,\sigma },{P}_{\sigma }\,\gt 0$. The total efficiency is then of the form $\eta =c{\eta }_{+}+(1-c){\eta }_{-}$, with $c={J}_{h,+}/({J}_{h,+}+{J}_{h,-})\leqslant 1$. The efficiency can only reach the Carnot value if one of the channels carries no energy (i.e. c = 1 or c = 0) or if both channels can simultaneously reach the Carnot point (which is the case if equations (17) hold). The fact that heat engines can only reach Carnot efficiency if they are effectively reduced to a single channel is also discussed in [38].

6.1.3. Exact numerics

When dealing with the full Hamiltonian, we define heat currents as the energy lost by the bath,

$$\begin{eqnarray}&&{J}_{\alpha }=-\displaystyle \frac{{\rm{d}}\mathrm{Tr}\{{\hat{H}}_{\alpha }^{l}{\hat{\varrho }}^{l}(t)\}}{{\rm{d}}{t}},\end{eqnarray} \tag{ 56 }$$

where we note that both ${\hat{H}}_{\alpha }^{l}$ and ${\hat{\varrho }}^{l}(t)$, the unitarily time evolved state of the whole compound, are taken in the lab frame. The power can be obtained through equation (39).

6.2. Reduced states

In addition to the energy currents, we consider the reduced state of the system as obtained by the different solutions. Since the considered Hamiltonian is bi-linear in bosonic annihilation and creation operators, it suffices to consider the covariance matrices. It will be convenient to work in the $\hat{x}$, $\hat{p}$ basis, with ${\hat{x}}_{\alpha }=\sqrt{1/2{{\rm{\Omega }}}_{\alpha }}({\hat{a}}_{\alpha }^{\dagger }+{\hat{a}}_{\alpha })$ and ${\hat{p}}_{\alpha }={\rm{i}}\sqrt{{{\rm{\Omega }}}_{\alpha }/2}({\hat{a}}_{\alpha }^{\dagger }-{\hat{a}}_{\alpha })$. Defining the vector $\hat{{\bf{r}}}=({\hat{x}}_{h},{\hat{x}}_{c},{\hat{p}}_{h},{\hat{p}}_{c})$, the covariance matrix is given by,

$$\begin{eqnarray}&&{{ \mathcal C }}_{{ij}}=\displaystyle \frac{1}{2}\mathrm{Tr}(\hat{\rho }({\hat{r}}_{i}{\hat{r}}_{j}+{\hat{r}}_{j}{\hat{r}}_{i})).\end{eqnarray} \tag{ 57 }$$

For the master equations, the covariance matrix can be obtained straightforwardly from section 4 (see equation (20), and the discussion around equation (30)). In the case of the exact numerics, one needs to consider the relevant entries of the covariance matrix of the whole compound (see e.g. [16]).

In the next section, we compare the reduced states obtained from the master equations to the reduced state obtained through exact numerics. As a measure of distinguishability between two states $\hat{\rho }$ and $\hat{\sigma }$, we consider the fidelity, defined as,

$$\begin{eqnarray}&&{ \mathcal F }(\hat{\rho },\hat{\sigma })=\mathrm{Tr}(\sqrt{\sqrt{\hat{\rho }}\hat{\sigma }\sqrt{\hat{\rho }}}).\end{eqnarray} \tag{ 58 }$$

For two-mode Gaussian states, represented by covariance matrices ${{ \mathcal C }}_{1}$ and ${{ \mathcal C }}_{2}$, and vanishing first order moments, the fidelity is given by [72],

$$\begin{eqnarray}&&{ \mathcal F }({{ \mathcal C }}_{1},{{ \mathcal C }}_{2})={[\sqrt{b}+\sqrt{c}-\sqrt{{(\sqrt{b}+\sqrt{c})}^{2}-a}]}^{-1}.\end{eqnarray} \tag{ 59 }$$

Here $a=\det ({{ \mathcal C }}_{1}+{{ \mathcal C }}_{2})$, $b={2}^{4}\det (J{{ \mathcal C }}_{1}J{{ \mathcal C }}_{2}-{\mathbb{I}}/4)$, $c={2}^{4}\det ({{ \mathcal C }}_{1}+{\rm{i}}{J}/2)\det ({{ \mathcal C }}_{2}+{\rm{i}}{J}/2)$, and the matrix elements of J are given by ${J}_{{kl}}=-{\rm{i}}\langle [{\hat{r}}_{k},{\hat{r}}_{l}]\rangle$.

Our results for the heat engine are illustrated in figures 2–7. We divide our results into three regimes. First we discuss the equilibrium regime, where the global approach shows excellent agreement with numerics for all values of the inter-system interaction g. Then we discuss the presence of a thermal bias but no external field. In this case, the global approach breaks down for small values of g. Finally, we focus on the heat engine regime which requires both a thermal bias as well as an external field. Again, the global approach breaks down for small values of g as expected. In all regimes, the local approach performs well for $g\ll {{\rm{\Omega }}}_{\alpha }$ (for both $\alpha =c,h$), which is where equations (17) are fulfilled for our system. As expected, the global approach performs well for $g\gg {\kappa }_{\alpha }$ (for both $\alpha =c,h$), which is where the secular approximation is well justified. In the following, all energies are given in units of ${{\rm{\Omega }}}_{c}$ and all energy currents are given in units of ${{\rm{\Omega }}}_{c}^{2}$.

7.1. Equilibrium

7.2. Thermal bias

7.3. Thermal bias and external field

Finally, we consider the regime where the considered system performs as a heat engine. This requires both a thermal bias ${T}_{c}\ne {T}_{h}$ and an external field ${{\rm{\Omega }}}_{h}\ne {{\rm{\Omega }}}_{c}$. In the thermoelectric realization of [44], this situation corresponds to the presence of a thermal and a voltage bias. In figure 4, the energy flows through the heat engine are plotted. While the local approach agrees extremely well with exact numerics, the global approach fails for small g. Note that all models fulfill ${J}_{c}+{J}_{h}=P$ (first law), and $\eta \lt {\eta }_{C}$ (second law). For the global approach, we note that it is crucial to use the definitions of energy currents given in section 6.1.2. If one uses instead definitions similar to the local approach (see equation (39)), the global approach results in incorrect heat currents, leading in particular to P = 0 and ${J}_{h}=-{J}_{c}$. This is consistent with the results of [15], which discusses the absence of any currents as expressed through system observables in the global approach.

Zoom In Zoom Out Reset image size

In figure 5 we compare the steady states against the exact numerics, using again fidelity as a figure of merit. Similarly to the case of a thermal bias without external field, we find that the two approaches give a faithful description in their respective regime of validity.

Zoom In Zoom Out Reset image size

For completeness, figure 6 illustrates the heat currents as a function of temperature in the presence and absence of an external field. These results strengthen the conclusions drawn above: the global approach is valid for $g\gg {\kappa }_{\alpha }$ while the local approach is valid for $g\ll {{\rm{\Omega }}}_{\alpha }$. In particular, the insets show that even at reasonably strong interaction strengths g, the local approach agrees well with numerics.

Zoom In Zoom Out Reset image size

Finally, we further illustrate the heat engine performance in figure 7 which shows the power and efficiency as a function of ${{\rm{\Omega }}}_{h}$, which determines the external field frequency (given by ${{\rm{\Omega }}}_{h}-{{\rm{\Omega }}}_{c}$). We find good agreement in the power as well as the efficiency between the local approach, the global approach, and the exact numerics. Only when the machine is operated close to the Carnot point (${{\rm{\Omega }}}_{h}/{T}_{h}={{\rm{\Omega }}}_{c}/{T}_{c}$) do the efficiencies deviate considerably. This is due to the fact that the power and the heat current become very small. Small differences in the energy flows then translate into large differences in the efficiency.

Zoom In Zoom Out Reset image size

To complete our discussion, we also consider a quantum thermal machine featuring finite-dimensional systems. Specifically we consider a quantum thermal machine consisting of two interacting qubits coupled to separate bosonic thermal baths, as shown in figure 8, analogous to the setup of figure 1 for harmonic oscillators. This machine can generate entanglement between the two qubits in the steady state, as shown in [51]. In the following, however, we focus on comparing the steady states obtained from local and global master equations in the same spirit as above.

Zoom In Zoom Out Reset image size

We denote the eigenstates of the free Hamiltonians of the qubits by $| 0\rangle$, $| 1\rangle$, and set the ground state energy to zero. The Hamiltonian of the system is then given by ${\hat{H}}_{S}={\hat{H}}_{0}+{\hat{H}}_{{\rm{int}}}$ with

$$\begin{eqnarray}{\hat{H}}_{0} & = & {{\rm{\Omega }}}_{c}| 1\rangle \langle 1| \otimes {{\mathbb{I}}}_{2}+{{\rm{\Omega }}}_{h}{{\mathbb{I}}}_{2}\otimes | 1\rangle \langle 1| ,\\ {\hat{H}}_{{\rm{int}}} & = & g(| 0,1\rangle \langle 1,0| +| 1,0\rangle \langle 0,1| ),\end{eqnarray} \tag{ 60 }$$

where ${{\rm{\Omega }}}_{c}$, ${{\rm{\Omega }}}_{h}$ are the energy gaps of the qubit, and g is the interaction strength. As in section 4, we compare local and global master equation models for the evolution of the system. We will focus on the degenerate case where ${{\rm{\Omega }}}_{c}={{\rm{\Omega }}}_{h}={\rm{\Omega }}$.

Both the local and global master equations can be written in GKSL form with constant rates [3, 4], and so lead to Markovian (specifically semigroup) evolution

$$\begin{eqnarray}&&{\partial }_{t}\hat{\rho }(t)=-{\rm{i}}[{\hat{H}}_{S},\hat{\rho }(t)]+\displaystyle \sum _{\alpha =c,h}\displaystyle \sum _{k}{{\rm{\Gamma }}}_{\alpha ,k}{ \mathcal D }[{\hat{L}}_{\alpha ,k}^{\dagger }]\rho (t)+{\bar{{\rm{\Gamma }}}}_{\alpha ,k}{ \mathcal D }[{\hat{L}}_{\alpha ,k}]\hat{\rho }(t),\end{eqnarray} \tag{ 61 }$$

where ${ \mathcal D }$ is defined in equation (14), ${\hat{L}}_{\alpha ,k}$ are jump operators, and ${{\rm{\Gamma }}}_{\alpha ,k}$, ${\bar{{\rm{\Gamma }}}}_{\alpha ,k}$ the corresponding rates. We consider bosonic baths for which

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\alpha ,\varepsilon }={\kappa }_{\alpha }(\varepsilon ){n}_{B}^{\alpha }(\varepsilon ),\hspace{1cm}{\bar{{\rm{\Gamma }}}}_{\alpha ,\varepsilon }={\kappa }_{\alpha }(\varepsilon )[{n}_{B}^{\alpha }(\varepsilon )+1].\end{eqnarray} \tag{ 62 }$$

Here, ${\kappa }_{\alpha }(\varepsilon )$ are the bath coupling strengths, and ε is the (absolute) energy difference associated with the jump induced by ${\hat{L}}_{\alpha ,\varepsilon }$. Thus the ${\hat{L}}_{\alpha ,\varepsilon }^{\dagger }$ correspond to jumps from lower to higher energies, absorbing energy from the bath α, while ${\hat{L}}_{\alpha ,\varepsilon }$ correspond to jumps decreasing the system energy, dissipating energy into the bath α.

Local and global master equations for the two-qubit machine can be derived using the same techniques as in section 4. The system-bath coupling can be taken to have the same form as in equation (2), with the system annihilation and creation operators replaced by ${\hat{A}}_{c}={\hat{\sigma }}_{-}\otimes {{\mathbb{I}}}_{2}$, ${\hat{A}}_{h}={{\mathbb{I}}}_{2}\otimes {\hat{\sigma }}_{-}$ and ${\hat{A}}_{c}^{\dagger }={\hat{\sigma }}_{+}\otimes {{\mathbb{I}}}_{2}$, ${\hat{A}}_{h}^{\dagger }={{\mathbb{I}}}_{2}\otimes {\hat{\sigma }}_{+}$ respectively, where ${\hat{\sigma }}_{-}=| 0\rangle \langle 1|$ and ${\hat{\sigma }}_{+}=| 1\rangle \langle 0|$. We take the spectral density of the baths to be Ohmic, as before. The bath coupling strengths are then linear in energy

$$\begin{eqnarray}&&{\kappa }_{\alpha }(\varepsilon )={\nu }_{\alpha }\varepsilon ,\end{eqnarray} \tag{ 63 }$$

for some constants ${\nu }_{\alpha }$. We denote ${\kappa }_{\alpha }={\kappa }_{\alpha }({\rm{\Omega }})$.

For the local master equation, there are just two jump operators, both corresponding to transitions with energy Ω. They are given by

$$\begin{eqnarray}&&{\hat{L}}_{c,{\rm{\Omega }}}={A}_{c}={\hat{\sigma }}_{-}\otimes {{\mathbb{I}}}_{2},\hspace{1cm}{\hat{L}}_{h,{\rm{\Omega }}}={A}_{h}={{\mathbb{I}}}_{2}\otimes {\hat{\sigma }}_{-}.\end{eqnarray} \tag{ 64 }$$

For the global master equation, the jump operators are found by diagonalizing the system Hamiltonian. Denoting the eigenvalues and eigenstates of ${\hat{H}}_{S}$ by λ and $| {\varphi }_{\lambda }\rangle$ respectively, one has

$$\begin{eqnarray}&&{\hat{L}}_{\alpha ,\varepsilon }=\displaystyle \sum _{\lambda -\lambda ^{\prime} =\varepsilon }| {\varphi }_{\lambda ^{\prime} }\rangle \langle {\varphi }_{\lambda ^{\prime} }| {\hat{A}}_{\alpha }| {\varphi }_{\lambda }\rangle \langle {\varphi }_{\lambda }| .\end{eqnarray} \tag{ 65 }$$

In the degenerate case, ${{\rm{\Omega }}}_{c}={{\rm{\Omega }}}_{h}={\rm{\Omega }}$, the eigenvalues are 0, ${\rm{\Omega }}\pm g$, and $2{\rm{\Omega }}$, and the corresponding eigenstates are

$$\begin{eqnarray}&&| {\varphi }_{0}\rangle =| 0,0\rangle ,\hspace{.75cm}| {\varphi }_{{\rm{\Omega }}\pm g}\rangle =(| 0,1\rangle \pm | 1,0\rangle )/\sqrt{2},\hspace{.75cm}| {\varphi }_{2{\rm{\Omega }}}\rangle =| 1,1\rangle .\end{eqnarray} \tag{ 66 }$$

The possible transition energies are $2g$, ${\rm{\Omega }}\pm g$, and $2{\rm{\Omega }}$. However, only transitions with energies ${\rm{\Omega }}\pm g$ can be induced by the system-bath coupling considered here. The non-zero jump operators are

$$\begin{eqnarray}{\hat{L}}_{c,{\rm{\Omega }}-g} & = & \displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{{\rm{\Omega }}+g}\rangle \langle {\varphi }_{2{\rm{\Omega }}}| -\displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{0}\rangle \langle {\varphi }_{{\rm{\Omega }}-g}| ,\\ {\hat{L}}_{c,{\rm{\Omega }}+g} & = & \displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{{\rm{\Omega }}-g}\rangle \langle {\varphi }_{2{\rm{\Omega }}}| +\displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{0}\rangle \langle {\varphi }_{{\rm{\Omega }}+g}| ,\\ {\hat{L}}_{h,{\rm{\Omega }}-g} & = & \displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{{\rm{\Omega }}+g}\rangle \langle {\varphi }_{2{\rm{\Omega }}}| +\displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{0}\rangle \langle {\varphi }_{{\rm{\Omega }}-g}| ,\\ {\hat{L}}_{h,{\rm{\Omega }}+g} & = & -\displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{{\rm{\Omega }}-g}\rangle \langle {\varphi }_{2{\rm{\Omega }}}| +\displaystyle \frac{1}{\sqrt{2}}| {\varphi }_{0}\rangle \langle {\varphi }_{{\rm{\Omega }}+g}| .\end{eqnarray} \tag{ 67 }$$

We note that these are indeed global in the sense that they involve transitions to and from the non-separable states $| {\varphi }_{{\rm{\Omega }}\pm g}\rangle$.

We can compute the steady state solutions of both the local and globel qubit master equations. As before, we compare them varying the interaction strength. In figure 9(a) we show the heat current as a function of the interaction strength between the two qubits. As for the harmonic oscillators, the global approach predicts a constant heat current, which is clearly unphysical as $g\to 0$, while the local model predicts a vanishing heat current in this limit, as expected. The two models agree well for intermediate coupling strength (i.e. $g/{\kappa }_{\alpha }\gtrsim 10$ in this case, note that ${\kappa }_{\alpha }$ is taken one order of magnitude smaller than in the previous sections). In figure 9(b) we show the fidelity between the steady states of the two models. Again, we see that the states agree well unless the coupling is weak.

Zoom In Zoom Out Reset image size

We investigated the accuracy of the local and global master equations for predicting thermodynamic quantities as well as system steady states describing a quantum heat engine. Exact numerics were used to benchmark the results. We found that the two approaches work very well in their respective regimes of validity which are (for bosonic baths with Ohmic spectral density):

• Local approach: $g\ll {{\rm{\Omega }}}_{\alpha }$.
• Global approach: $g\gg {\kappa }_{\alpha }$.

More generally, the condition under which the local approach gives a faithful description is given by equation (17). Since the Markov approximation, which underlies both approaches, requires ${\kappa }_{\alpha }\ll {{\rm{\Omega }}}_{\alpha }$, the two regimes of validity overlap. As expected, we find good agreement between the local and the global approach in this region of parameter space.

We note that the local approach is by no means more phenomenological than the global approach. Indeed, the approximation leading to the local master equation is completely analogous to the Markov approximation and has a well defined regime of validity.

Finally, we also investigated a qubit entangler. For this system, no benchmark is available. However, the similarity to the results obtained for the heat engine strongly suggests that similar conclusions with respect to the applicability of the local and the global master equation are valid. We therefore conjecture that our results are qualitatively valid for a variety of baths and system Hamiltonians. As long as the bath-correlation time is much shorter than any inverse inter-system interaction strength, the local approach is valid. The global approach is valid as long as the inter-system interaction strengths are much stronger than the system-bath interaction strengths.

We therefore conclude that the local approach provides a valid description for thermal machines that consist of weakly interacting sub-systems.

We acknowledge fruitful discussions with Mark Mitchison and valuable feedback from Ronnie Kosloff, Amikam Levy, Gediminas Kiršanskas, Andreas Wacker, Felipe Barra, Aashish Clerk, Kay Brandner, and Massimiliano Esposito. MP-L acknowledges support from the Alexander von Humboldt Foundation. All other authors acknowledge the Swiss National Science Foundation (Starting grant DIAQ, grant 200021_169002, Marie-Heim Vögtlin grant 164466, and QSIT). All authors are grateful for support from the EU COST Action MP1209 on Thermodynamics in the quantum regime.

In this appendix we discuss the robustness of the numerics, that we take as a benchmark, to the specific choices of n (the number of oscillators) and ${\omega }_{c}$ (the cutoff). First of all, recall that we model the bath as a collection of oscillators with frequencies,

$$\begin{eqnarray}&&{\omega }_{k,\alpha }=\displaystyle \frac{k}{n}{\omega }_{{c}},\end{eqnarray} \tag{ A.1 }$$

and coupling constants,

$$\begin{eqnarray}&&{\gamma }_{k,\alpha }={\eta }_{\alpha }\sqrt{k}\displaystyle \frac{{\omega }_{c}}{n}\end{eqnarray} \tag{ A.2 }$$

which models an Ohmic spectral density in the continuous limit. The ${\eta }_{\alpha }$ in equation (A.2) are given by ${\kappa }_{\alpha }=2\pi {\eta }_{\alpha }{{\rm{\Omega }}}_{\alpha }$.

In all the figures of the main text, we take ${\omega }_{c}=3$ and n = 400. In order to see how sensitive our results are to this choice, we plot the heat current as a function of n (the inset) and ${\omega }_{c}$, and compare it with the analytic results (using the local approach) in figure A1. It is clearly observed that the results are independent of n. On the other hand, we see a small dependence on the values of ${\omega }_{c}$. For ${\omega }_{c}\approx {{\rm{\Omega }}}_{\alpha }$, the results do not closely match the analytics, which is expected because energetically possible transitions are not captured by the bath. For ${\omega }_{c}\geqslant 2{{\rm{\Omega }}}_{\alpha }$, very good agreement is obtained, with small differences that increase with ${\omega }_{c}$. This is due to the renormalization of the Hamiltonian (see equations (16) and (27)), which is neglected in the analytic calculations. The choice ${\omega }_{c}=3$ is hence large enough to capture the different energy transitions, and at the same time not too large so that the effect of the renormalization can be neglected. Similar considerations hold for the case of non-degenerate frequencies of the oscillators of the system. We hence conclude that our numerical benchmark is quite robust to the choice of n and ${\omega }_{c}$.

Zoom In Zoom Out Reset image size