A thermodynamically consistent Markovian master equation beyond the secular approximation

For non-interacting particles, the heat current can be written using a Landauer-like formula [42]. For bosons (with vanishing chemical potential) it takes on the form

$$\begin{equation}{J}^{\text{t}}={\int }_{0}^{\infty }\frac{\mathrm{d}\omega }{2\pi }\mathcal{T}(\omega )\omega [{n}_{\text{B}}^{\text{h}}(\omega )-{n}_{\text{B}}^{\text{c}}(\omega )],\end{equation} \tag{ 8 }$$

where we introduced the Bose–Einstein distribution

$$\begin{equation}{n}_{\text{B}}^{\alpha }(\omega )=\frac{1}{{\text{e}}^{\frac{\omega }{{k}_{\text{B}}{T}_{\alpha }}}-1},\end{equation} \tag{ 9 }$$

and the transmission function $\mathcal{T}(\omega )$ can be obtained using non-equilibrium Green's functions [43, 44]. For the present system, it reads (see for instance reference [45] for the analogous case of Fermions)

$$\begin{equation}\mathcal{T}(\omega )=\frac{{g}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}}{\vert \left(\omega -{{\Omega}}_{\text{h}}+\text{i}\frac{{\kappa }_{\text{h}}}{2}\right)\left(\omega -{{\Omega}}_{\text{c}}+\text{i}\frac{{\kappa }_{\text{c}}}{2}\right)-{g}^{2}{\vert }^{2}}.\end{equation} \tag{ 10 }$$

The Landauer-like formula has an intuitive interpretation. At each energy ω, bosons traverse the system with the rate $\mathcal{T}(\omega )\mathrm{d}\omega$. The transmission function has two peaks located at Ω₊ and Ω₋. For small g and Δ, these peaks merge as illustrated in the insets of figure 2(b).

Zoom In Zoom Out Reset image size

In this approach, the dissipators act locally on the two oscillators

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{cl}}={\kappa }_{\alpha }\left\{{n}_{\text{B}}^{\alpha }({{\Omega}}_{\alpha })\mathcal{D}[{\hat{a}}_{\alpha }^{{\dagger}}]+[{n}_{\text{B}}^{\alpha }({{\Omega}}_{\alpha })+1]\mathcal{D}[{\hat{a}}_{\alpha }]\right\},\end{equation} \tag{ 11 }$$

where cl stands for conventional local. Here we used the GKLS superoperators

$$\begin{equation}\mathcal{D}[\hat{A}]\hat{\rho }=\hat{A}\hat{\rho }{\hat{A}}^{{\dagger}}-\frac{1}{2}\left\{{\hat{A}}^{{\dagger}}\hat{A},\hat{\rho }\right\},\end{equation} \tag{ 12 }$$

with {⋅, ⋅} denoting the anti-commutator. The heat current in this approach is defined as

$$\begin{equation}{J}^{\text{cl}}=\mathrm{Tr}\left\{{\hat{H}}_{\text{S}}{\mathcal{L}}_{\text{h}}^{\text{cl}}\hat{\rho }\right\}.\end{equation} \tag{ 13 }$$

In steady state, this reduces to

$$\begin{equation}{J}^{\text{cl}}=\frac{{\kappa }_{\text{c}}{{\Omega}}_{\text{h}}+{\kappa }_{\text{h}}{{\Omega}}_{\text{c}}}{{({\kappa }_{\text{c}}+{\kappa }_{\text{h}})}^{2}}\frac{4{g}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}})-{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}})]}{4{g}^{2}+{\kappa }_{\text{c}}{\kappa }_{\text{h}}+16{{\Delta}}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}/{({\kappa }_{\text{c}}+{\kappa }_{\text{h}})}^{2}}.\end{equation} \tag{ 14 }$$

It is important to note that there is a microscopic derivation that results in this local approach, which is valid for the current system as long as ${n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{h}})\simeq {n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{c}})$ [13], which is the case for

$$\begin{equation}g,{\Delta}\ll \bar{{\Omega}}.\end{equation} \tag{ 15 }$$

When this inequality is satisfied, we may expect that the laws of thermodynamics hold, see figure 2. When this approximation is not justified, the second law of thermodynamics is not guaranteed as illustrated in figure 3. The reason for this is that the Bose–Einstein distributions are evaluated at Ω_h and Ω_c respectively in equation (14). This implies that the bosons change their energy when traversing the system. From the Landauer-like formula (cf equation (8)), it is apparent that this cannot happen. We note that while equation (15) ensures $\vert {n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{h}})-{n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{c}})\vert \ll 1$, the relative error in the heat current can still be sizable when the occupation numbers (and thus the heat current) become small. This is the case because equation (15) does not ensure $\vert {n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{h}})-{n}_{\text{B}}^{\alpha }({{\Omega}}_{\text{c}})\vert \ll {n}_{\text{B}}^{\alpha }({{\Omega}}_{{\alpha}})$.

Zoom In Zoom Out Reset image size

Here we present a thermodynamically consistent local approach that follows from our general master equation introduced below in section 3, where we motivate this approach microscopically and show how it can be generalized to more complicated systems. As in the previous approach, the dissipators act locally in this approach

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{l}}={\kappa }_{\alpha }\left\{{n}_{\text{B}}^{\alpha }(\bar{{\Omega}})\mathcal{D}[{\hat{a}}_{\alpha }^{{\dagger}}]+[{n}_{\text{B}}^{\alpha }(\bar{{\Omega}})+1]\mathcal{D}[{\hat{a}}_{\alpha }]\right\},\end{equation} \tag{ 16 }$$

the only difference being that the Bose–Einstein distributions are now all evaluated at the frequency $\bar{{\Omega}}$. This approach has the same regime of validity as the previous local approach. Indeed, equation (15) implies

$$\begin{equation}\vert {n}_{\text{B}}^{\alpha }(\bar{{\Omega}})-{n}_{\text{B}}^{\alpha }({{\Omega}}_{\alpha })\vert \ll 1,\end{equation} \tag{ 17 }$$

and we expect the dissipators in equations (11) and (16) to result in approximately the same dynamics as illustrated in figure 2. If equation (15) is not satisfied, the two local approaches may result in different results, see figure 3.

To obtain a thermodynamically consistent description, the approximation in equation (15) is exploited in the definition of the heat current which reads

$$\begin{equation}{J}^{\text{l}}=\mathrm{Tr}\left\{{\hat{H}}_{\text{TD}}{\mathcal{L}}_{\text{h}}^{\text{l}}\hat{\rho }\right\},\end{equation} \tag{ 18 }$$

where we introduced a separate Hamiltonian for the thermodynamic bookkeeping (TD)

$$\begin{equation}{\hat{H}}_{\text{TD}}=\bar{{\Omega}}({\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{c}}+{\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{h}}).\end{equation} \tag{ 19 }$$

Whenever equation (15) holds, this Hamiltonian provides the (approximately) correct energy flows. Consider the case where Δ = 0. Then, ${\hat{H}}_{\text{TD}}$ is simply obtained by dropping the coupling between the oscillators. The thermodynamic bookkeeping thus neglects both the system-bath couplings, as well as the coupling between the oscillators. Of course, all these couplings are crucial for the dynamics, where they appear either in the Hamiltonian or in the dissipators. In section 3, we illustrate how such a thermodynamic Hamiltonian can be constructed in general, providing thermodynamic consistency of the corresponding master equation at the price of a slightly reduced energy resolution of thermodynamic quantities.

In steady state, the heat current reads

$$\begin{equation}{J}^{\text{l}}=\frac{\bar{{\Omega}}}{{\kappa }_{\text{c}}+{\kappa }_{\text{h}}}\frac{4{g}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}[{n}_{\text{B}}^{\text{h}}(\bar{{\Omega}})-{n}_{\text{B}}^{\text{c}}(\bar{{\Omega}})]}{4{g}^{2}+{\kappa }_{\text{c}}{\kappa }_{\text{h}}+16{{\Delta}}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}/{({\kappa }_{\text{c}}+{\kappa }_{\text{h}})}^{2}}.\end{equation} \tag{ 20 }$$

We note that for Δ = 0, we find J^l = J^cl in the steady state. However, in the transient regime the two local approaches may result in different heat currents even for Δ = 0. The difference between the two approaches is in this case the energy associated to the coupling term (proportional to g), which is dropped in the local approach. We note that the heat current associated to this coupling term has been interpreted as a quantum contribution before [19].

Comparing the local master equation with the Landauer-like formula (cf equation (8)), one can see that transmission is approximated to happen at a single energy $\bar{{\Omega}}$. This is a good approximation when equation (15) holds. In this case, the transmission function is only non-zero around $\bar{{\Omega}}$ and the Bose–Einstein distributions may be assumed constant across all energies where transmission is non-zero. Interestingly, the local approach agrees well with the transmission approach for the present system, even when equation (15) is not satisfied, see figure 3(b). This is in stark contrast to the conventional local approach.

In the global approach, the dissipators are introducing jumps between the eigenstates of the Hamiltonian

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{g}}=\sum\limits _{\sigma =\pm }\enspace {\kappa }_{\alpha }^{\sigma }\left\{{n}_{\text{B}}^{\alpha }({{\Omega}}_{\sigma })\mathcal{D}[{\hat{a}}_{\sigma }^{{\dagger}}]+[{n}_{\text{B}}^{\alpha }({{\Omega}}_{\sigma })+1]\mathcal{D}[{\hat{a}}_{\sigma }]\right\},\end{equation} \tag{ 21 }$$

where the coupling strengths reflect the spatial distribution of the eigenstates

$$\begin{equation}{\kappa }_{\text{h}}^{+}={\kappa }_{\text{h}}\enspace {\mathrm{cos}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{h}}^{-}={\kappa }_{\text{h}}\enspace {\mathrm{sin}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{c}}^{+}={\kappa }_{\text{c}}\enspace {\mathrm{sin}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{c}}^{-}={\kappa }_{\text{c}}\enspace {\mathrm{cos}}^{2}(\theta /2).\end{equation} \tag{ 22 }$$

The heat current in the global approach reads

$$\begin{equation}{J}^{\text{g}}=\mathrm{Tr}\left\{{\hat{H}}_{\text{S}}{\mathcal{L}}_{\text{h}}^{\text{g}}\hat{\rho }\right\},\end{equation} \tag{ 23 }$$

which in steady state reduces to

$$\begin{equation}{J}^{\text{g}}=\sum\limits _{\sigma =\pm }\enspace {{\Omega}}_{\sigma }\frac{{\kappa }_{\text{c}}^{\sigma }{\kappa }_{\text{h}}^{\sigma }}{{\kappa }_{\text{c}}^{\sigma }+{\kappa }_{\text{h}}^{\sigma }}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\sigma })-{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\sigma })].\end{equation} \tag{ 24 }$$

This approach is thermodynamically consistent and we will always find heat flowing from hot to cold. As it relies on the secular approximation, it requires the condition

$$\begin{equation}{\kappa }_{\text{c}},{\kappa }_{\text{h}}\ll \sqrt{{{\Delta}}^{2}+{g}^{2}}.\end{equation} \tag{ 25 }$$

In the global master equation, the transmission is approximated to happen at the two energies Ω_±. If equation (25) holds, the transmission function given in equation (10) consists of two narrow peaks located at Ω_±, see the right inset in figure 2(b). Over the width of these peaks, the Bose–Einstein distributions can be assumed constant and the global master equation is valid. Note that if equation (25) does not hold, the peaks overlap and the secular approximation breaks down, see figure 2. We also note that equations (15) and (25) may both be valid. In this case, both the global as well as the local approach are expected to give accurate results.

In the PERLind approach, the dissipators read

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{p}}=\sum\limits _{k=-1,1}\mathcal{D}[{\hat{J}}_{\alpha ,k}],\end{equation} \tag{ 26 }$$

with the jump operators

$$\begin{equation}\begin{aligned}\hfill {\hat{J}}_{\text{c,}-1}& =\sqrt{{\kappa }_{\text{c}}^{+}[{n}_{\text{B}}^{\text{c}}({{\Omega}}_{+})+1]}{\hat{a}}_{+}+\sqrt{{\kappa }_{\text{c}}^{-}[{n}_{\text{B}}^{\text{c}}({{\Omega}}_{-})+1]}{\hat{a}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{c,}\;1}& =\sqrt{{\kappa }_{\text{c}}^{+}{n}_{\text{B}}^{\text{c}}({{\Omega}}_{+})}{\hat{a}}_{+}^{{\dagger}}+\sqrt{{\kappa }_{\text{c}}^{-}{n}_{\text{B}}^{\text{c}}({{\Omega}}_{-})}{\hat{a}}_{-}^{{\dagger}},\hfill \\ \hfill {\hat{J}}_{\text{h,}-1}& =\sqrt{{\kappa }_{\text{h}}^{+}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{+})+1]}{\hat{a}}_{+}-\sqrt{{\kappa }_{\text{h}}^{-}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{-})+1]}{\hat{a}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{h,}\;1}& =\sqrt{{\kappa }_{\text{h}}^{+}{n}_{\text{B}}^{\text{h}}({{\Omega}}_{+})}{\hat{a}}_{+}^{{\dagger}}-\sqrt{{\kappa }_{\text{h}}^{-}{n}_{\text{B}}^{\text{h}}({{\Omega}}_{-})}{\hat{a}}_{-}^{{\dagger}}.\hfill \end{aligned}\end{equation} \tag{ 27 }$$

In this approach, the heat current is defined as

$$\begin{equation}{J}^{\text{p}}=\mathrm{Tr}\left\{{\hat{H}}_{\text{S}}{\mathcal{L}}_{\text{h}}^{\text{p}}\hat{\rho }\right\}.\end{equation} \tag{ 28 }$$

For the present system, the heat current in the PERLind approach agrees perfectly with the result obtained from the transmission approach for all considered parameter values. This is expected as we consider scenarios where equation (6) is fulfilled and the Born–Markov approximations are justified. Clearly, the PERLind approach has strong advantages. However, it also has its disadvantages. First, the analytical expressions quickly become unwieldy, which is the case already for the simple system considered here. Second, the PERLind approach may violate the second law of thermodynamics [29]. When heat is defined by equation (28), a necessary and sufficient condition for the second law to hold is given by [46, 47]

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{p}}{\text{e}}^{-{\beta }_{\alpha }{\hat{H}}_{\text{S}}}=0.\end{equation} \tag{ 29 }$$

This condition is not fulfilled as shown in reference [30]. However, as long as the Born–Markov approximations are justified, any second law violations should become vanishingly small. It is an open question if a different definition for heat, that is consistent with the Born–Markov approximation, could salvage the second law.

As we have shown above, all approaches have their advantages and disadvantages. The PERLind approach provides the most accurate description, while the conventional local approach is the most simple, not requiring diagonalization of the Hamiltonian. Both approaches may however result in violations of the second law. The global and the local approach both respect the laws of thermodynamics. Furthermore, together they provide an accurate description for all parameter values as they are justified in complementary, but overlapping, parameter regimes (cf equations (6), (15) and (25), assuming the Born–Markov approximations are justified). In the following, we illustrate how these thermodynamically consistent approaches follow from a unified framework that can be extended to scenarios where neither a global nor a local description is justified.

We consider the general scenario described by the Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{tot}}(t)={\hat{H}}_{\text{S}}(t)+\sum\limits _{\alpha }\left({\hat{H}}_{\alpha }+{\hat{V}}_{\alpha }\right)={\hat{H}}_{\text{S}}(t)+{\hat{H}}_{\text{B}}+\hat{V},\end{equation} \tag{ 30 }$$

where the first term describes the Hamiltonian of the system (which may be time-dependent) and the second and third term describe the thermal reservoirs (labeled by α) and their coupling to the system respectively. The system exchanges energy and particles with the reservoirs, such that energy changes can be divided into heat and work. The average heat that leaves bath α during the time interval [0, t] is given by

$$\begin{equation}\langle {Q}_{\alpha }\rangle =-\mathrm{Tr}\left\{({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha }){\hat{\rho }}_{\text{tot}}(t)\right\}+\mathrm{Tr}\left\{({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha }){\hat{\rho }}_{\text{tot}}(0)\right\},\end{equation} \tag{ 31 }$$

where ${\hat{N}}_{\alpha }$ denotes the particle number operator for reservoir α and μ_α its chemical potential. The average work provided by reservoir α is

$$\begin{equation}\langle {W}_{\alpha }\rangle =-{\mu }_{\alpha }\left(\mathrm{Tr}\left\{{\hat{N}}_{\alpha }{\hat{\rho }}_{\text{tot}}(t)\right\}-\mathrm{Tr}\left\{{\hat{N}}_{\alpha }{\hat{\rho }}_{\text{tot}}(0)\right\}\right).\end{equation} \tag{ 32 }$$

In addition, the time-dependence of the system Hamiltonian results in the external average power

$$\begin{equation}{P}_{\text{ext}}(t)=\mathrm{Tr}\left\{[{\partial }_{t}{\hat{H}}_{\text{S}}(t)]{\hat{\rho }}_{\text{tot}}(t)\right\}.\end{equation} \tag{ 33 }$$

Before deriving a Markovian description, we discuss the laws of thermodynamics as they hold for the general scenario.

3.2.1. The 0th law

For a large environment in thermal equilibrium (i.e. described by a single inverse temperature β and chemical potential μ) and a time-independent system Hamiltonian, the reduced state of the system tends to [48, 49]

$$\begin{equation}{\hat{\rho }}_{\text{S}}={\mathrm{Tr}}_{\text{B}}\left\{{\text{e}}^{-\beta ({\hat{H}}_{\text{tot}}-\mu {\hat{N}}_{\mathrm{t}\mathrm{o}\mathrm{t}})}/Z\right\},\quad Z=\mathrm{Tr}\left\{{\text{e}}^{-\beta ({\hat{H}}_{\text{tot}}-\mu {\hat{N}}_{\mathrm{t}\mathrm{o}\mathrm{t}})}\right\},\end{equation} \tag{ 34 }$$

where Tr_B denotes the trace over the reservoir degrees of freedom. In the weak system-bath coupling limit, equation (34) reduces to the Gibbs state

$$\begin{equation}{\hat{\rho }}_{\text{S}}=\frac{{\text{e}}^{-\beta ({\hat{H}}_{\text{S}}-\mu {\hat{N}}_{\text{S}})}}{\mathrm{Tr}\left\{{\text{e}}^{-\beta ({\hat{H}}_{\text{S}}-\mu {\hat{N}}_{\text{S}})}\right\}}.\end{equation} \tag{ 35 }$$

3.2.2. The 1st law

For future reference, we first introduce the heat current and power provided by bath α as

$$\begin{equation}{J}_{\alpha }(t)={\partial }_{t}\langle {Q}_{\alpha }\rangle ,\hspace{25.0pt}{P}_{\alpha }(t)={\partial }_{t}\langle {W}_{\alpha }\rangle .\end{equation} \tag{ 36 }$$

The first law of thermodynamics is then given by

$$\begin{equation}{\partial }_{t}U(t)={P}_{\text{ext}}(t)+\sum\limits _{\alpha }[{J}_{\alpha }(t)+{P}_{\alpha }(t)],\quad U=\mathrm{Tr}\left\{[{\hat{H}}_{\text{S}}(t)+\hat{V}]{\hat{\rho }}_{\text{tot}}(t)\right\}.\end{equation} \tag{ 37 }$$

In the weak system-bath coupling limit, the coupling energy can be neglected and $U=\mathrm{Tr}\left\{{\hat{H}}_{\text{S}}(t){\hat{\rho }}_{\text{tot}}(t)\right\}$ reduces to the usual internal energy of the system.

3.2.3. The 2nd law

To express the second law of thermodynamics, we impose the initial condition

$$\begin{equation}{\hat{\rho }}_{\text{tot}}(0)={\hat{\rho }}_{\text{S}}(0)\underset{\alpha }{\bigotimes}{\hat{\tau }}_{\alpha },\quad {\hat{\tau }}_{\alpha }=\frac{{\text{e}}^{-\beta ({\hat{H}}_{\alpha }-\mu {\hat{N}}_{\alpha })}}{Tr\left\{{\text{e}}^{-{\beta }_{\alpha }({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha })}\right\}},\end{equation} \tag{ 38 }$$

i.e. all reservoirs are in local thermal equilibrium and uncorrelated with the system and the other reservoirs. With this initial condition, the second law of thermodynamics can be written as [50]

$$\begin{equation}{\Sigma}(t)\equiv {\Delta}S(t)-\sum\limits _{\alpha }\frac{\langle {Q}_{\alpha }\rangle }{{T}_{\alpha }}=S\left[{\hat{\rho }}_{\text{tot}}(t){\Vert}{\hat{\rho }}_{\text{S}}(t)\underset{\alpha }{\bigotimes}{\hat{\tau }}_{\alpha }\right]\geqslant 0,\end{equation} \tag{ 39 }$$

where ${\hat{\rho }}_{\text{S}}={\mathrm{Tr}}_{\text{B}}\left\{{\hat{\rho }}_{\text{tot}}(t)\right\}$ denotes the reduced state of the system, $S[{\hat{\rho }}_{1}{\Vert}\hat{{\rho }_{2}}]=\mathrm{Tr}\left\{{\hat{\rho }}_{1}(\mathrm{ln}\enspace {\hat{\rho }}_{1}-\mathrm{ln}\enspace {\hat{\rho }}_{2})\right\}$ is the quantum relative entropy (which is by definition positive for positive definite density matrices ${\hat{\rho }}_{1},{\hat{\rho }}_{2}$ with unity trace), and ΔS denotes the change in the system's von Neumann entropy

$$\begin{equation}{\Delta}S(t)=-{k}_{\text{B}}\enspace \mathrm{Tr}\left\{{\hat{\rho }}_{\text{S}}(t)\mathrm{ln}\enspace {\hat{\rho }}_{\text{S}}(t)\right\}+{k}_{\text{B}}\enspace \mathrm{Tr}\left\{{\hat{\rho }}_{\text{S}}(0)\mathrm{ln}\enspace {\hat{\rho }}_{\text{S}}(0)\right\}.\end{equation} \tag{ 40 }$$

Equation (39) has an intuitive interpretation: the entropy production Σ denotes the information that is lost when describing system and reservoirs by the state ${\hat{\rho }}_{\text{S}}(t){\bigotimes}_{\alpha }{\hat{\tau }}_{\alpha }$. In this description, correlations between the system and the reservoirs, as well as any displacement from equilibrium of the reservoirs are neglected [51].

While equation (39) is always non-negative, the same is not necessarily true for the entropy production rate

$$\begin{equation}\dot{{\Sigma}}(t)\equiv {\partial }_{t}{\Delta}S(t)-\sum\limits _{\alpha }\frac{{J}_{\alpha }(t)}{{T}_{\alpha }}.\end{equation} \tag{ 41 }$$

A negative entropy production rate can be understood as information backflow and is a hallmark of non-Markovian behavior [52]. For systems amenable to a Markovian description, we expect $\dot{{\Sigma}}\geqslant 0$, as the initial condition in equation (38) can be translated in time without altering the dynamics.

Throughout this work, we are mainly interested in average thermodynamic quantities. In order to obtain a Markovian description of these, it is nevertheless instructive to start with the full probability distribution for the heat and work exchanged with the reservoirs. With the initial condition given in equation (38), this distribution can be written as

$$\begin{align}\hfill P(\boldsymbol{Q},\boldsymbol{W})& =\sum\limits _{\boldsymbol{E},{\boldsymbol{E}}^{\prime },\boldsymbol{N},{\boldsymbol{N}}^{\prime }}{P}_{t}({\boldsymbol{E}}^{\prime },{\boldsymbol{N}}^{\prime }\vert \boldsymbol{E},\boldsymbol{N}){P}_{0}(\boldsymbol{E},\boldsymbol{N})\hfill \\ \hfill & \quad \times \prod\limits _{\alpha }\delta ({W}_{\alpha }-{\mu }_{\alpha }{N}_{\alpha }^{\prime }+{\mu }_{\alpha }{N}_{\alpha })\delta ({Q}_{\alpha }+{W}_{\alpha }-{E}_{\alpha }^{\prime }+{E}_{\alpha }),\hfill \end{align} \tag{ 42 }$$

where we grouped Q_α , W_α , E_α , and N_α into vectors Q , W , E , and N and similarly for ${E}_{\alpha }^{\prime }$ and ${N}_{\alpha }^{\prime }$. The joint probability for each bath α having energy E_α and particle number N_α at time t = 0 is given by

$$\begin{equation}{P}_{0}(\boldsymbol{E},\boldsymbol{N})=\prod\limits _{\alpha }\frac{{\text{e}}^{-{\beta }_{\alpha }({E}_{\alpha }-{\mu }_{\alpha }{N}_{\alpha })}}{\mathrm{Tr}\left\{{\text{e}}^{-{\beta }_{\alpha }({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha })}\right\}}.\end{equation} \tag{ 43 }$$

The conditional probability that the reservoirs have energies ${E}_{\alpha }^{\prime }$ and particle number ${N}_{\alpha }^{\prime }$ at time t, given that their energies and particle numbers where E_α and N_α initially reads

$$\begin{equation}{P}_{t}({\boldsymbol{E}}^{\prime },{\boldsymbol{N}}^{\prime }\vert \boldsymbol{E},\boldsymbol{N})=\mathrm{Tr}\left\{\hat{M}{\hat{\rho }}_{\text{S}}(0){\hat{M}}^{{\dagger}}\right\},\quad \hat{M}={\mathrm{Tr}}_{\text{B}}\left\{\vert \boldsymbol{E},\boldsymbol{N}\rangle \langle {\boldsymbol{E}}^{\prime },{\boldsymbol{N}}^{\prime }\vert \hat{U}(t)\right\},\end{equation} \tag{ 44 }$$

where ${\hat{H}}_{\alpha }\vert \boldsymbol{E},\boldsymbol{N}\rangle ={E}_{\alpha }\vert \boldsymbol{E},\boldsymbol{N}\rangle$, ${\hat{N}}_{\alpha }\vert \boldsymbol{E},\boldsymbol{N}\rangle ={N}_{\alpha }\vert \boldsymbol{E},\boldsymbol{N}\rangle$, and we introduced the time-evolution operator

$$\begin{equation}\hat{U}(t)=\mathcal{T}{\text{e}}^{-\text{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }{\hat{H}}_{\text{tot}}({t}^{\prime })},\end{equation} \tag{ 45 }$$

with $\mathcal{T}$ denoting the time-ordering operator. We note that the probability distribution given in (42) can in principle be measured by applying projective measurements on the reservoirs at the initial and final time. Importantly, while such measurements may not be experimentally feasible, they do not influence the dynamics of the system due to the chosen initial condition.

The moment generating function is provided by the Fourier transform of the probability distribution in (42)

$$\begin{equation}{\Lambda}(\boldsymbol{\lambda },\boldsymbol{\chi })\equiv \int \mathrm{d}\boldsymbol{Q}\enspace \mathrm{d}\boldsymbol{W}P(\boldsymbol{Q},\boldsymbol{W}){\text{e}}^{-\text{i}\boldsymbol{\lambda }\cdot \boldsymbol{Q}-i\boldsymbol{\chi }\cdot \boldsymbol{W}}=\mathrm{Tr}\left\{{\hat{\rho }}_{\text{tot}}(\boldsymbol{\lambda },\boldsymbol{\chi };t)\right\},\end{equation} \tag{ 46 }$$

where we introduced

$$\begin{equation}{\hat{\rho }}_{\text{tot}}(\boldsymbol{\lambda },\boldsymbol{\chi };t)=\hat{U}(\boldsymbol{\lambda },\boldsymbol{\chi };t){\hat{\rho }}_{\text{tot}}(0){\hat{U}}^{{\dagger}}(-\boldsymbol{\lambda },-\boldsymbol{\chi };t)\end{equation} \tag{ 47 }$$

with the modified time-evolution operator

$$\begin{equation}\hat{U}(\boldsymbol{\lambda },\boldsymbol{\chi };t)={\text{e}}^{-\frac{\text{i}}{2}\sum\limits _{\alpha }[{\lambda }_{\alpha }({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha })+{\chi }_{\alpha }{\mu }_{\alpha }{\hat{N}}_{\alpha }]}\hat{U}(t){\text{e}}^{\frac{\text{i}}{2}\sum\limits _{\alpha }[{\lambda }_{\alpha }({\hat{H}}_{\alpha }-{\mu }_{\alpha }{\hat{N}}_{\alpha })+{\chi }_{\alpha }{\mu }_{\alpha }{\hat{N}}_{\alpha }]}.\end{equation} \tag{ 48 }$$

The quantities λ_α and χ_α are known as counting fields and allow us to keep track of work and heat exchanged with the reservoirs (for the use of counting fields in master equations, see references [4, 53]). From the moment generating function, we recover the average values for heat and work given in section 3.1 as

$$\begin{equation}\langle {Q}_{\alpha }\rangle =\text{i}{\partial }_{{\lambda }_{\alpha }}{\Lambda}(\boldsymbol{\lambda },0){\vert }_{\boldsymbol{\lambda }=0},\hspace{25.0pt}\langle {W}_{\alpha }\rangle =\text{i}{\partial }_{{\chi }_{\alpha }}{\Lambda}(0,\boldsymbol{\chi }){\vert }_{\boldsymbol{\chi }=0}.\end{equation} \tag{ 49 }$$

We note that we do not include the external power in the probability distribution for heat and work as this goes beyond the scope of this paper. Indeed, fluctuations of the external power may not necessarily be described by a positive probability distribution without introducing a measurement scheme that potentially alters the dynamics of the system [54–56].

We are now in a position to derive a Markovian master equation, keeping track of the thermodynamic bookkeeping associated to the heat and work exchanged with the bath. Together with a secular approximation, resulting in a global master equation, such a procedure has been applied before [57–59] (for a similar approach based on path-integrals, see [60]). Here we will make an approximation that is different from the secular approximation, allowing for treating (near) degeneracies while still ensuring thermodynamic consistency. We follow standard procedure [1, 4] and write the system-bath coupling as

$$\begin{equation}\hat{V}=\sum\limits _{\alpha ,k}\enspace {\hat{S}}_{\alpha ,k}{\hat{B}}_{\alpha ,k},\end{equation} \tag{ 50 }$$

where the operators ${\hat{S}}_{\alpha ,k}\enspace ({\hat{B}}_{\alpha ,k})$ act only on the system (reservoir). We note that we do not assume that these operators are Hermitian. To determine the jump operators that will enter the Markovian master equation, we need the Fourier coefficients of the operators ${\hat{S}}_{\alpha ,k}$ in the interaction picture, i.e.

$$\begin{equation}{\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{S}}_{\alpha ,k}{\hat{U}}_{\text{S}}(t)=\sum\limits _{j}\enspace {\text{e}}^{-\text{i}{\omega }_{j}t}{\hat{S}}_{\alpha ,k;j},\quad {\hat{U}}_{\text{S}}(t)=\mathcal{T}{\text{e}}^{-\text{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }{\hat{H}}_{\text{S}}({t}^{\prime })}.\end{equation} \tag{ 51 }$$

For a time-independent Hamiltonian, the frequencies ω_j denote the energy gaps in the Hamiltonian and the operators ${\hat{S}}_{\alpha ,k;j}$ are ladder operators that induce transitions between the corresponding states. For a periodic Hamiltonian with period t_p = 2π/ϖ, the ${\hat{S}}_{\alpha ,k;j}$ are ladder operators of an averaged Hamiltonian defined by

$$\begin{equation}{\hat{U}}_{\text{S}}({t}_{\text{p}})={\text{e}}^{-\text{i}{\hat{H}}_{\text{av}}{t}_{\text{p}}}.\end{equation} \tag{ 52 }$$

The frequencies then fulfill ω_j = ν_j + l_j ϖ, where ν_j denotes an energy gap of ${\hat{H}}_{\text{av}}$ and l_j denotes an integer, corresponding to the exchange of photons with the driving field [57, 61, 62].

Standard application of Born–Markov approximations then results in the Redfield equation including counting fields (in the interaction picture)

$$\begin{align}\hfill {\partial }_{t}{\hat{\rho }}_{\text{S}}(t)& =-\sum\limits _{\alpha ,k,{k}^{\prime },j,{j}^{\prime }}{\text{e}}^{\text{i}({\omega }_{j}-{\omega }_{{j}^{\prime }})t}{\int }_{0}^{\infty }\mathrm{d}s\enspace \mathcal{I}(s,t),\hfill \\ \hfill \mathcal{I}(s,t)& ={\text{e}}^{\text{i}{\omega }_{{j}^{\prime }}s}{C}_{k,{k}^{\prime }}^{\alpha }(s){\hat{S}}_{\alpha ,k;j}^{{\dagger}}{\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }}{\hat{\rho }}_{\text{S}}(t)+{\text{e}}^{-\text{i}{\omega }_{j}s}{C}_{k,{k}^{\prime }}^{\alpha }(-s){\hat{\rho }}_{\text{S}}(t){\hat{S}}_{\alpha ,k;j}^{{\dagger}}{\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }}\hfill \\ \hfill & \quad -{\text{e}}^{\text{i}{\mu }_{\alpha }{n}_{\alpha ,k}({\lambda }_{\alpha }-{\chi }_{\alpha })}\left[{\text{e}}^{\text{i}{\omega }_{{j}^{\prime }}s}{C}_{k,{k}^{\prime }}^{\alpha }(s+{\lambda }_{\alpha })+{\text{e}}^{-\text{i}{\omega }_{j}s}{C}_{k,{k}^{\prime }}^{\alpha }(-s+{\lambda }_{\alpha })\right]{\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }}{\hat{\rho }}_{\text{S}}(t){\hat{S}}_{\alpha ,k;j}^{{\dagger}},\hfill \end{align} \tag{ 53 }$$

where we suppressed the counting field dependence of the density matrix for ease of notation. We note that the term Redfield equation sometimes refers to the non-Markovian equation that is obtained before taking the time-integration to infinity [1] while we include this limit here. Here we used global particle conservation, which ensures

$$\begin{equation}[{\hat{S}}_{\alpha ,k;j},{\hat{N}}_{\text{S}}]={n}_{\alpha ,k}{\hat{S}}_{\alpha ,k;j},\hspace{25.0pt}[{\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }},{\hat{N}}_{\text{S}}]={n}_{\alpha ,k}{\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }},\end{equation} \tag{ 54 }$$

for all pairs k and k' that appear together in equation (53), i.e. ${\hat{S}}_{\alpha ,k;j}$ and ${\hat{S}}_{\alpha ,{k}^{\prime };{j}^{\prime }}$ change the particle number by the same amount, such that no superpositions of particle numbers are created. We further introduced the bath correlation functions

$$\begin{equation}{C}_{k,{k}^{\prime }}^{\alpha }(s)=\mathrm{Tr}\left\{{\text{e}}^{\text{i}s{\hat{H}}_{\alpha }}{\hat{B}}_{\alpha ,k}^{{\dagger}}\enspace {\text{e}}^{-\text{i}s{\hat{H}}_{\alpha }}{\hat{B}}_{\alpha ,{k}^{\prime }}{\hat{\tau }}_{\alpha }\right\}.\end{equation} \tag{ 55 }$$

These correlation functions define a bath correlation time τ_B by their characteristic decay time. The Markov approximation is generally justified when the bath correlation time is much shorter than the characteristic time-scale over which ${\hat{\rho }}_{\text{S}}$ changes in the interaction picture, τ_S. The time τ_S describes the relaxation time of the system and is determined by the inverse of the system-bath coupling. Note however that the counting fields λ_α enter the argument of the bath correlation times. For the counting-field dependent density matrix, the Markov approximation is only justified if ${C}_{k,{k}^{\prime }}^{\alpha }(\pm \tau +{\lambda }_{\alpha })\simeq 0$ for τ ≳ τ_S. This implies the regime of validity

$$\begin{equation}{\tau }_{\text{B}}\ll {\tau }_{\text{S}},\quad \vert {\lambda }_{\alpha }\vert \ll {\tau }_{\text{S}}.\end{equation} \tag{ 56 }$$

The last inequality has a very important consequence. It implies that only the low frequency components of the heat distribution are to be trusted. Due to the uncertainty principle between Fourier conjugate variables [63], this implies that energy-differences in heat of the order of 1/τ_S cannot be resolved. Thus, whenever a Markovian description is employed, the heat exchanged with the reservoirs suffers from a limited energy-resolution. With this in mind, it is no surprise that Markovian descriptions may result in thermodynamic inconsistencies. There is however a straightforward solution to the problem: as our resolution of heat is finite, we may change the definition of heat such that it fulfills the laws of thermodynamics while remaining the same within our limited resolution.

It is well known that the Redfield equation does not preserve positivity of the density matrix and can thus result in negative probabilities. Multiple schemes have been put forward to achieve positivity. Here we introduce a novel scheme that has the same regime of validity as the Born–Markov assumptions and thus goes beyond the regime of validity of the secular approximation. As we show in the next subsection, our approach allows for a thermodynamically consistent formulation.

Equation (56) ensures that for any two transition frequencies we either have |ω_j − ω_j'| ≪ 1/τ_B or |ω_j − ω_j'| ≫ 1/τ_S. Indeed, both of these conditions may be fulfilled simultaneously as τ_B ≪ τ_S. We may thus group the transition frequencies into sets x_q , such that the first (second) inequality holds if the transition frequencies are in the same (different) set, that is

$$\begin{equation}\begin{aligned}\hfill & \vert {\omega }_{j}-{\omega }_{{j}^{\prime }}\vert \ll 1/{\tau }_{\text{B}}\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace {\omega }_{j}\in {x}_{q},\enspace {\omega }_{{j}^{\prime }}\in {x}_{{q}^{\prime }}\enspace \enspace \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\enspace \enspace q={q}^{\prime },\hfill \\ \hfill & \vert {\omega }_{j}-{\omega }_{{j}^{\prime }}\vert \gg 1/{\tau }_{\text{S}}\quad \mathrm{f}\mathrm{o}\mathrm{r}\enspace {\omega }_{j}\in {x}_{q},\enspace {\omega }_{{j}^{\prime }}\in {x}_{{q}^{\prime }}\enspace \enspace \mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}\enspace \enspace q\ne {q}^{\prime }.\hfill \end{aligned}\end{equation} \tag{ 57 }$$

It is important to note that this procedure may not always work, for instance if the ω_j form a continuum. However, in small quantum systems, where the number of transition frequencies is finite, this procedure is expected to work. It is particularly well suited for thermal machines that consist of weakly coupled sub-units. These naturally exhibit sets of near-degenerate transition frequencies (see the examples below).

In the spirit of the secular approximation, we then drop terms in equation (53) where ω_j and ω_j' belong to different sets. This is justified as the corresponding terms exhibit fast oscillations that average to zero. For transition frequencies in the same set x_q , we follow the spirit of the Markov approximation and replace

$$\begin{equation}{\text{e}}^{\text{i}{\omega }_{j}s},{\text{e}}^{\text{i}{\omega }_{{j}^{\prime }}s}\to {\text{e}}^{\text{i}{\omega }_{q}s},\hspace{25.0pt}\vert {\omega }_{q}-{\omega }_{j}\vert \ll 1/{\tau }_{\text{B}}\quad \forall \enspace {\omega }_{j}\in {x}_{q},\end{equation} \tag{ 58 }$$

within $\mathcal{I}(s)$ in equation (53), while keeping the frequency differences in the prefactor ${\text{e}}^{\text{i}({\omega }_{j}-{\omega }_{{j}^{\prime }})t}$ governing the coherent dynamics of the system. For each set, we thus choose a set frequency ω_q . All transition frequencies ω_j in the set x_q are then replaced by the set frequency, because they are virtually indistinguishable over the time-scale over which the integrand in equation (53) is finite.

This approximation results in the Lindblad master equation

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{\text{S}}(\boldsymbol{\lambda },\boldsymbol{\chi };t)=-\text{i}[{\hat{H}}_{\text{LS}},{\hat{\rho }}_{\text{S}}(\boldsymbol{\lambda },\boldsymbol{\chi };t)]+\sum\limits _{\alpha }\enspace {\tilde{\mathcal{L}}}_{\alpha }^{{\chi }_{\alpha },{\lambda }_{\alpha }}{\hat{\rho }}_{\text{S}}(\boldsymbol{\lambda },\boldsymbol{\chi };t),\end{equation} \tag{ 59 }$$

with

$$\begin{equation}{\tilde{\mathcal{L}}}_{\alpha }^{{\chi }_{\alpha },{\lambda }_{\alpha }}\hat{\rho }=\sum\limits _{k;q}\enspace {{\Gamma}}_{k}^{\alpha }({\omega }_{q})\left[{\text{e}}^{-\text{i}{\lambda }_{\alpha }{\omega }_{q}-\text{i}({\chi }_{\alpha }-{\lambda }_{\alpha }){\mu }_{\alpha }{n}_{\alpha ,k}}{\hat{S}}_{\alpha ,k;q}(t)\hat{\rho }{\hat{S}}_{\alpha ,k;q}^{{\dagger}}(t)-\frac{1}{2}\left\{{\hat{S}}_{\alpha ,k;q}^{{\dagger}}(t){\hat{S}}_{\alpha ,k;q}(t),\hat{\rho }\right\}\right],\end{equation} \tag{ 60 }$$

where the tilde denotes the interaction picture and we introduced the jump operators

$$\begin{equation}{\hat{S}}_{\alpha ,k;q}(t)=\sum\limits _{\left\{j\vert {\omega }_{j}\in {x}_{q}\right\}}{\text{e}}^{-\text{i}{\omega }_{j}t}{\hat{S}}_{\alpha ,k;j},\end{equation} \tag{ 61 }$$

and the Lamb-shift Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{LS}}=\sum\limits _{\alpha ,k;q}{{\Delta}}_{k}^{\alpha }({\omega }_{q}){\hat{S}}_{\alpha ,k;q}^{{\dagger}}(t){\hat{S}}_{\alpha ,k;q}(t),\end{equation} \tag{ 62 }$$

as well as the quantities

$$\begin{equation}{{\Gamma}}_{k}^{\alpha }(\omega )={\int }_{-\infty }^{\infty }\mathrm{d}s\enspace {\text{e}}^{\text{i}\omega s}{C}_{k,k}^{\alpha }(s)\hspace{25.0pt}{{\Delta}}_{k}^{\alpha }(\omega )=-\frac{\text{i}}{2}{\int }_{-\infty }^{\infty }\mathrm{d}s\enspace {\text{e}}^{\text{i}\omega s}\enspace \text{sign}(s){C}_{k,k}^{\alpha }(s).\end{equation} \tag{ 63 }$$

For simplicity, we assumed ${C}_{k,{k}^{\prime }}^{\alpha }\propto {\delta }_{k,{k}^{\prime }}$ to derive equation (59). We note that relaxing this assumption is straightforward.

In the absence of counting fields, we may use the Kubo–Martin–Schwinger condition [1] (which is slightly complicated by our definition of the bath-correlation functions) to write the master equation as

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{\text{S}}=-\text{i}[{\hat{H}}_{\text{LS}}(t),{\hat{\rho }}_{\text{S}}]+\sum\limits _{\alpha }\enspace {\tilde{\mathcal{L}}}_{\alpha }{\hat{\rho }}_{\text{S}},\end{equation} \tag{ 64 }$$

with

$$\begin{equation}{\tilde{\mathcal{L}}}_{\alpha }=\sum\limits _{\left\{q\vert {\omega }_{q} > 0\right\}}\sum\limits _{k}\enspace {{\Gamma}}_{k}^{\alpha }({\omega }_{q})\left\{\mathcal{D}[{\hat{S}}_{\alpha ,k;q}(t)]+{\text{e}}^{-{\beta }_{\alpha }({\omega }_{q}-{\mu }_{\alpha }{n}_{\alpha ,k})}\mathcal{D}[{\hat{S}}_{\alpha ,k;q}^{{\dagger}}(t)]\right\}.\end{equation} \tag{ 65 }$$

For a time-independent Hamiltonian, the master equation in the Schrödinger picture is given by

$$\begin{equation}{\partial }_{t}{\hat{\rho }}_{\text{S}}=-\text{i}[{\hat{H}}_{\text{S}}+{\hat{H}}_{\text{LS}},{\hat{\rho }}_{\text{S}}]+\sum\limits _{\alpha }\enspace {\mathcal{L}}_{\alpha }{\hat{\rho }}_{\text{S}},\end{equation} \tag{ 66 }$$

with

$$\begin{equation}{\mathcal{L}}_{\alpha }=\sum\limits _{\left\{q\vert {\omega }_{q} > 0\right\}}\sum\limits _{k}\enspace {{\Gamma}}_{k}^{\alpha }({\omega }_{q})\left\{\mathcal{D}[{\hat{S}}_{\alpha ,k;q}]+{\text{e}}^{-{\beta }_{\alpha }({\omega }_{q}-{\mu }_{\alpha }{n}_{\alpha ,k})}\mathcal{D}[{\hat{S}}_{\alpha ,k;q}^{{\dagger}}]\right\}.\end{equation} \tag{ 67 }$$

where ${\hat{S}}_{\alpha ,k;q}\equiv {\hat{S}}_{\alpha ,k;q}(0)$ and ${\hat{H}}_{\text{LS}}\equiv {\hat{H}}_{\text{LS}}(0)$. For time-dependent Hamiltonians, one has to be more careful because in general ${\hat{S}}_{\alpha ,k;q}(t)\ne {\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{S}}_{\alpha ,k;q}{\hat{U}}_{\text{S}}(t)$, see also appendix A.

There are two simple limits for the master equation in equation (66). If all pairs of transition frequencies fulfill |ω_j − ω_j'| ≫ 1/τ_S, then each transition frequency may be associated to a separate set x_j . We then recover the secular approximation where $\mathcal{D}[{\hat{S}}_{\alpha ,k;q}(t)]=\mathcal{D}[{\hat{S}}_{\alpha ,k;j}]$ (the time-dependence of the jump operators drops out in the interaction picture). If all pairs of transition frequencies fulfill |ω_j − ω_j'| ≪ 1/τ_B, then all transition frequencies may be grouped into a single set x₀ such that ${\hat{S}}_{\alpha ,k;q}(t)={\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{S}}_{\alpha ,k}{\hat{U}}_{\text{S}}(t)$ (the time-dependence of the jump operators drops out upon returning to the Schrödinger picture). We then recover a local master equation, which has the appeal that the Hamiltonian does not need to be diagonalized in order to identify the jump operators. This local approach differs from the conventional local approach as only a single frequency enters the bath distribution function.

Importantly, the counting fields λ_α have an influence on the regime of validity of equation (59), just as for the Markov approximation. Indeed, for the approximation to be valid, |ω_q − ω_j | does not only have to be much smaller than 1/τ_B but also has to be much smaller than 1/|λ_α | (in complete analogy to equation (56)). Just like the Markov approximation, this frequency-grouping thus limits the energy-resolution for heat exchanged with the reservoirs. As a consequence, energy changes of the order of |ω_j − ω_j'|, where both frequencies are within the same set, can no longer be resolved. This finite resolution in energy may result in the false conclusion that particles change their energy when traversing the system. In the conventional local approach, this results in thermodynamic inconsistencies when using standard definitions for heat, see section 2.

We note that dropping terms that involve frequencies from different sets is reminiscent of the partial secular approximation developed in references [15, 34, 35, 64, 65]. Here, as well as in reference [41], an additional coarse-graining in energy (replacing frequencies with set-frequencies) results in a master equation in GKLS form.

We now turn to the question of how to utilize the finite energy-resolution to ensure thermodynamic consistency. To do this, we introduce a second Hamiltonian, ${\hat{H}}_{\text{TD}}$, that provides the correct thermodynamic bookkeeping under the approximations that resulted in our master equation. For simplicity, we consider here a time-independent system Hamiltonian. The time-dependent case is slightly more complicated and treated in appendix A. Due to the frequency grouping outlined in the last subsection, all frequencies in the set x_q are effectively replaced by the frequency ω_q from the point of view of the reservoirs. To ensure a consistent thermodynamic bookkeeping, the Hamiltonian ${\hat{H}}_{\text{TD}}$ needs to fulfill

$$\begin{equation}[{\hat{S}}_{\alpha ,k;j},{\hat{H}}_{\text{TD}}]={\omega }_{q}{\hat{S}}_{\alpha ,k;j},\end{equation} \tag{ 68 }$$

for all frequencies ω_j ∈ x_q . To fulfill this, the thermodynamic Hamiltonian can be obtained from ${\hat{H}}_{\text{S}}$ by changing its eigenvalues such that all frequencies ω_j → ω_q for ω_j ∈ x_q . Such a rescaling is expected to always be possible if the frequency grouping is possible as discussed above.

With the thermodynamic Hamiltonian at hand, we define the internal energy of the system as

$$\begin{equation}U=\mathrm{Tr}\left\{{\hat{H}}_{\text{TD}}\hat{\rho }\right\}.\end{equation} \tag{ 69 }$$

Furthermore, from equations (59), (36), and (49) we find that the heat current and power provided by bath α can be cast into

$$\begin{equation}{J}_{\alpha }=\mathrm{Tr}\left\{({\hat{H}}_{\text{TD}}-{\mu }_{\alpha }{\hat{N}}_{\text{S}}){\mathcal{L}}_{\alpha }{\hat{\rho }}_{\text{S}}\right\},\hspace{25.0pt}{P}_{\alpha }={\mu }_{\alpha }\enspace \mathrm{Tr}\left\{{\hat{N}}_{\text{S}}{\mathcal{L}}_{\alpha }{\hat{\rho }}_{\text{S}}\right\},\end{equation} \tag{ 70 }$$

where ${\mathcal{L}}_{\alpha }$ is defined in equation (67). We note that while ${\hat{H}}_{\text{TD}}$ determines the thermodynamic bookkeeping, it is still ${\hat{H}}_{\text{S}}$ that determines the kinetics, i.e. enters the master equation in equation (59). For future reference, we also introduce the total power produced by the system as

$$\begin{equation}{P}_{\text{S}}(t)=-\sum\limits _{\alpha }\enspace {P}_{\alpha }(t)-{P}_{\text{ext}}(t).\end{equation} \tag{ 71 }$$

This will be the quantity of interest when we consider heat engines.

3.6.1. The 0th law

Using equation (68), it is straightforward to show that

$$\begin{equation}{\mathcal{L}}_{\alpha }{\text{e}}^{-{\beta }_{\alpha }({\hat{H}}_{\text{TD}}-{\mu }_{\alpha }{\hat{N}}_{\text{S}})}=0.\end{equation} \tag{ 72 }$$

as well as $[{\hat{H}}_{\text{TD}},{\hat{H}}_{\text{LS}}]=0$. Furthermore, as ${\hat{H}}_{\text{TD}}$ is obtained by changing only the eigenvalues of ${\hat{H}}_{\text{S}}$, these two Hamiltonians also commute. When all reservoirs have the same inverse temperature β and chemical potential μ, it then follows that the Gibbs state with respect to the thermodynamic Hamiltonian,

$$\begin{equation}{\hat{\rho }}_{\text{G}}(\beta ,\mu )=\frac{{\text{e}}^{-\beta ({\hat{H}}_{\text{TD}}-\mu {\hat{N}}_{\text{S}})}}{\mathrm{Tr}\left\{{\text{e}}^{-\beta ({\hat{H}}_{\text{TD}}-\mu {\hat{N}}_{\text{S}})}\right\}},\end{equation} \tag{ 73 }$$

is the steady state of equation (66). Compared to equation (35), this state neglects the system-bath coupling, as well as the differences between ${\hat{H}}_{\text{S}}$ and ${\hat{H}}_{\text{TD}}$. This is consistent with our approximations, which imply that these differences cannot be resolved within our Markovian treatment.

3.6.2. The first law

Using $[{\hat{H}}_{\text{TD}},{\hat{H}}_{\text{S}}+{\hat{H}}_{\text{LS}}]=0$, the first law of thermodynamics follows directly from equations (69) and (70) and reads

$$\begin{equation}{\partial }_{t}U=\sum\limits _{\alpha }\enspace ({P}_{\alpha }+{J}_{\alpha }).\end{equation} \tag{ 74 }$$

3.6.3. The second law

The entropy production rate can be written as

$$\begin{align}\hfill \dot{{\Sigma}}/{k}_{\text{B}}& =-{\partial }_{t}\enspace \mathrm{Tr}\left\{{\hat{\rho }}_{\text{S}}\enspace \mathrm{ln}\enspace {\hat{\rho }}_{\text{S}}\right\}-\sum\limits _{\alpha }\enspace {\beta }_{\alpha }{J}_{\alpha }=-\sum\limits _{\alpha }\mathrm{Tr}\left\{({\mathcal{L}}_{\alpha }{\hat{\rho }}_{\text{S}})[\mathrm{ln}\enspace {\hat{\rho }}_{\text{S}}-\mathrm{ln}\enspace {\hat{\rho }}_{\text{G}}({\beta }_{\alpha },{\mu }_{\alpha })]\right\}\geqslant 0.\hfill \end{align} \tag{ 75 }$$

Here we used equations (66), (70), and the last inequality is known as Spohns inequality [46, 66] and relies on ${\hat{\rho }}_{\text{G}}({\beta }_{\alpha },{\mu }_{\alpha })$ being a fixed point of ${\mathcal{L}}_{\alpha }$, cf equation (72).

The master equation in equation (66), together with the thermodynamic bookkeeping introduced in equations (69) and (70), thus provide a thermodynamically consistent description. We stress that the only assumption that went into equation (66) is the one that justifies the standard Born–Markov approximations, i.e. τ_B ≪ τ_S. However, when considering the full probability distribution for heat and work exchanged with the reservoirs, cf (59), then we further have restrictions on the counting fields λ_α which are the Fourier transform variables of heat. These restrictions imply that we lose energy-resolution on the scale of 1/τ_S as well as $\vert {\omega }_{j}-{\omega }_{j}^{\prime }\vert$ for pairs of frequencies that are close to each other (and grouped into one set). Thermodynamic consistency is obtained by using an appropriate definition of heat that is consistent with our limited energy-resolution (i.e. we may freely add terms of order $\vert {\omega }_{j}-{\omega }_{j}^{\prime }\vert$ to heat, as our results are not to be trusted to this order in the first place).

In this section, we illustrate the derivation of our master equation for a simple but non-trivial example of non-interacting electrons. Here, as well as for all time-independent examples, we provide all equations in the Schrödinger picture. The system under consideration is sketched in figure 4 and consists of two coupled single-level quantum dots

$$\begin{equation}{\hat{H}}_{\text{S}}={{\Omega}}_{\text{L}}{\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{{\Omega}}_{\text{R}}{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}+g({\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{R}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{L}}),\end{equation} \tag{ 76 }$$

with the standard anti-commutation relations

$$\begin{equation}\left\{{\hat{d}}_{\alpha },{\hat{d}}_{\beta }^{{\dagger}}\right\}={\delta }_{\alpha ,\beta },\hspace{25.0pt}\left\{{\hat{d}}_{\alpha },{\hat{d}}_{\beta }\right\}=0.\end{equation} \tag{ 77 }$$

In equation (76), Ω_α denote the on-site energies and g the coupling strength. We note that this system is a fermionic version of the system considered in section 2 and thus shares many of its properties. Due to the finite chemical potential of the reservoirs, it can however feature as a heat engine, leveraging a heat current to drive a particle current against a chemical potential bias. The reservoirs and their coupling to the system are described by

$$\begin{equation}{\hat{H}}_{\alpha }=\sum\limits _{l}\enspace {\varepsilon }_{\alpha ,l}{\hat{c}}_{\alpha ,l}^{{\dagger}}{\hat{c}}_{\alpha ,l},\hspace{25.0pt}{\hat{V}}_{\alpha }=\sum\limits _{l}\enspace {t}_{\alpha ,l}({\hat{d}}_{\alpha }^{{\dagger}}{\hat{c}}_{\alpha ,l}+{\hat{c}}_{\alpha ,l}^{{\dagger}}{\hat{d}}_{\alpha }),\end{equation} \tag{ 78 }$$

with α = L, R. From equations (50) and (78), we identify the following system and bath operators

$$\begin{equation}{\hat{S}}_{\alpha ,1}={\hat{d}}_{\alpha }^{{\dagger}},\hspace{25.0pt}{\hat{S}}_{\alpha ,-1}={\hat{d}}_{\alpha },\end{equation} \tag{ 79 }$$

and

$$\begin{equation}{\hat{B}}_{\alpha ,1}=\sum\limits _{l}\enspace {t}_{\alpha ,l}{\hat{c}}_{\alpha ,l},\hspace{25.0pt}{\hat{B}}_{\alpha ,-1}=\sum\limits _{l}\enspace {t}_{\alpha ,l}{\hat{c}}_{\alpha ,l}^{{\dagger}}.\end{equation} \tag{ 80 }$$

To derive the master equation and judge its validity, the bath correlation functions need to be inspected. We defer this discussion to appendix B.1, and simply state the conclusion that the Born Markov approximations are valid whenever

$$\begin{equation}{\kappa }_{\alpha }\ll \mathrm{max}\left\{{k}_{\text{B}}{T}_{\beta },\enspace \vert {\omega }_{j}-{\mu }_{\beta }\vert \right\},\end{equation} \tag{ 81 }$$

where the inequality should hold for all choices of α, β, and j and we assumed an energy-independent bath spectral density. Physically, equation (81) ensures that the rates

$$\begin{equation}{{\Gamma}}_{1}^{\alpha }(-\omega )={\kappa }_{\alpha }{n}_{\text{F}}^{\alpha }(\omega ),\hspace{25.0pt}{{\Gamma}}_{-1}^{\alpha }(\omega )={\kappa }_{\alpha }[1-{n}_{\text{F}}^{\alpha }(\omega )],\end{equation} \tag{ 82 }$$

with the Fermi–Dirac distribution

$$\begin{equation}{n}_{\text{F}}^{\alpha }(\omega )=\frac{1}{{\text{e}}^{\frac{\omega -{\mu }_{\alpha }}{{k}_{\text{B}}{T}_{\alpha }}}+1},\end{equation} \tag{ 83 }$$

is flat over the energy scale κ_α . At high temperatures, this is ensured because the Fermi–Dirac distribution becomes flat. For large |ω_j − μ_β |, the eigenenergies of the double quantum dot lie far away from the chemical potential, where the Fermi–Dirac distribution is flat and takes on the value zero or one.

Zoom In Zoom Out Reset image size

To identify the transition frequencies and the jump operators, we require the Fourier coefficients of the system operators given in equation (79)

$$\begin{equation}\begin{aligned}\hfill & {\text{e}}^{\text{i}{\hat{H}}_{\text{S}}t}{\hat{d}}_{\text{R}}{\text{e}}^{-\text{i}{\hat{H}}_{\text{S}}t}={\text{e}}^{-\text{i}{{\Omega}}_{-}t}\enspace \text{cos}(\theta /2){\hat{d}}_{-}+{\text{e}}^{-\text{i}{{\Omega}}_{+}t}\enspace \text{sin}(\theta /2){\hat{d}}_{+},\hfill \\ \hfill & {\text{e}}^{\text{i}{\hat{H}}_{\text{S}}t}{\hat{d}}_{\text{L}}{\text{e}}^{-\text{i}{\hat{H}}_{\text{S}}t}=-{\text{e}}^{-\text{i}{{\Omega}}_{-}t}\enspace \text{sin}(\theta /2){\hat{d}}_{-}+{\text{e}}^{-\text{i}{{\Omega}}_{+}t}\enspace \text{cos}(\theta /2){\hat{d}}_{+},\hfill \end{aligned}\end{equation} \tag{ 84 }$$

where we used a similar notation to section 2, i.e. the diagonalized Hamiltonian reads

$$\begin{equation}{\hat{H}}_{\text{S}}={{\Omega}}_{+}{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}+{{\Omega}}_{-}{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-},\end{equation} \tag{ 85 }$$

where ${{\Omega}}_{\pm }=\bar{{\Omega}}\pm \sqrt{{{\Delta}}^{2}+{g}^{2}}$, with $\bar{{\Omega}}=({{\Omega}}_{\text{R}}+{{\Omega}}_{\text{L}})/2$, Δ = (Ω_R − Ω_L)/2, and $\mathrm{cos}(\theta )={\Delta}/\sqrt{{{\Delta}}^{2}+{g}^{2}}$. Comparing equation (84), and their Hermitian conjugates, to equation (51), we find the transition frequencies

$$\begin{equation}\left\{{\omega }_{j}\right\}=\left\{{{\Omega}}_{+},\enspace {{\Omega}}_{-},\enspace -{{\Omega}}_{+},\enspace -{{\Omega}}_{-}\right\},\end{equation} \tag{ 86 }$$

and the corresponding jump operators

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{R},-1;j}\right\}=\left\{\mathrm{sin}(\theta /2){\hat{d}}_{+},\enspace \mathrm{cos}(\theta /2){\hat{d}}_{-},\enspace 0,\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},-1;j}\right\}=\left\{\mathrm{cos}(\theta /2){\hat{d}}_{+},\enspace -\mathrm{sin}(\theta /2){\hat{d}}_{-},\enspace 0,\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{R},1;j}\right\}=\left\{0,\enspace 0,\enspace \mathrm{sin}(\theta /2){\hat{d}}_{+}^{{\dagger}},\enspace \mathrm{cos}(\theta /2){\hat{d}}_{-}^{{\dagger}}\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},1;j}\right\}=\left\{0,\enspace 0,\enspace \mathrm{cos}(\theta /2){\hat{d}}_{+}^{{\dagger}},\enspace -\mathrm{sin}(\theta /2){\hat{d}}_{-}^{{\dagger}}\right\}.\hfill \end{aligned}\end{equation} \tag{ 87 }$$

We note that the frequencies Ω_± only feature in the jump operators with k = −1 (corresponding to electrons leaving the system), while the frequencies −Ω_± only feature in the jump operators with k = 1 (corresponding to electrons entering the system). This implies that we may consider only the frequencies Ω_± when deciding how to group the transition frequencies into sets. We may group Ω₊ and Ω₋ in the same set, or we may assign them to different sets. For the frequencies −Ω_±, we then choose an equivalent grouping.

Grouping Ω_± into different sets is justified as long as their difference $2\sqrt{{{\Delta}}^{2}+{g}^{2}}$ is much larger than max{κ_L, κ_R}, which corresponds to 1/τ_S. This grouping results in a different set for each frequency, such that {ω_q } = {ω_j } and $\left\{{\hat{S}}_{\alpha ,k;q}\right\}=\left\{{\hat{S}}_{\alpha ,k;j}\right\}$. Having identified the set frequencies and jump operators, we may then use equation (67) to obtain the well-known dissipator in the secular approximation

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{g}}=\sum\limits _{\sigma =\pm }{\kappa }_{\alpha }^{\sigma }\left\{{n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma })\mathcal{D}[{\hat{d}}_{\sigma }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma })]\mathcal{D}[{\hat{d}}_{\sigma }]\right\},\end{equation} \tag{ 88 }$$

with the coupling strengths

$$\begin{equation}{\kappa }_{\text{L}}^{+}={\kappa }_{\text{L}}\enspace {\mathrm{cos}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{L}}^{-}={\kappa }_{\text{L}}\enspace {\mathrm{sin}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{R}}^{+}={\kappa }_{\text{R}}\enspace {\mathrm{sin}}^{2}(\theta /2),\hspace{7.11317pt}{\kappa }_{\text{R}}^{-}={\kappa }_{\text{R}}\enspace {\mathrm{cos}}^{2}(\theta /2).\end{equation} \tag{ 89 }$$

As the jump operators in equation (87) are ladder operators of ${\hat{H}}_{\text{S}}$, no rescaling is required for obtaining the thermodynamic Hamiltonian and we find ${\hat{H}}_{\text{TD}}={\hat{H}}_{\text{S}}$, resulting in the standard definition for heat currents (cf equation (70)).

Alternatively, we may group Ω₊ and Ω₋ in the same set. In this case, we need to choose a set frequency. The average value $\bar{{\Omega}}$ is a natural choice. Using the same arguments that result in equation (81), this grouping is justified as long as $2\sqrt{{{\Delta}}^{2}+{g}^{2}}\ll \mathrm{max}\left\{{k}_{\text{B}}{T}_{\beta },\enspace \vert {\omega }_{j}-{\mu }_{\beta }\vert \right\}$. In this case, we end up with two sets with the set frequencies

$$\begin{equation}\left\{{\omega }_{q}\right\}=\left\{\bar{{\Omega}},\enspace -\bar{{\Omega}}\right\},\end{equation} \tag{ 90 }$$

and the jump operators

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{R},-1;q}\right\}=\left\{{\hat{d}}_{\text{R}},\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},-1;q}\right\}=\left\{{\hat{d}}_{\text{L}},\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{R},1;q}\right\}=\left\{0,\enspace {\hat{d}}_{\text{R}}^{{\dagger}}\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},1;q}\right\}=\left\{0,\enspace {\hat{d}}_{\text{L}}^{{\dagger}}\right\},\hfill \end{aligned}\end{equation} \tag{ 91 }$$

which are simply obtained by adding the jump operators in equation (87) which belong to the same set. Inserting these quantities into equation (60) results in the local approach with the dissipator

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{l}}={\kappa }_{\alpha }\left\{{n}_{\text{F}}^{\alpha }(\bar{{\Omega}})\mathcal{D}[{\hat{d}}_{\alpha }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }(\bar{{\Omega}})]\mathcal{D}[{\hat{d}}_{\alpha }]\right\},\end{equation} \tag{ 92 }$$

for α = L, R.

The thermodynamic Hamiltonian is obtained by rescaling all transition frequencies to their set frequency. In the present case this rescaling is obtained by ${{\Omega}}_{\pm }\to \bar{{\Omega}}$ resulting in

$$\begin{equation}{\hat{H}}_{\text{TD}}=\bar{{\Omega}}({\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}+{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-})=\bar{{\Omega}}({\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}).\end{equation} \tag{ 93 }$$

Depending on how the transition frequencies are grouped, we obtain either the global or the local approach. As long as the Born–Markov approximations are justified, i.e. as long as equation (81) holds, at least one of the two procedures is justified. As the parameters of the system are changed, we may thus need to change from the local to the global approach to maintain an adequate description of the system. This is a general drawback of our master equation: as parameters are varied, the grouping of transition frequencies into sets may need adaptation. This results in discontinuities which are expected to be small as long as the Born–Markov approximations are justified.

In figure 5, we illustrate the global, the local, as well as the benchmark provided by the transmission approach for this system when operated as a heat engine (see appendix B.2 for analytical expressions). We find that our master equation reproduces the transmission approach well if the frequency grouping is chosen such that the local (global) approach is obtained for g below (above) 5κ. Interestingly, we find that the efficiency is generally better reproduced by the global approach (see insets in figure 5). While the local approach predicts Carnot efficiency when ${n}_{\text{F}}^{\text{L}}(\bar{{\Omega}})={n}_{\text{F}}^{\text{R}}(\bar{{\Omega}})$, the transmission and the global approach predict a drop in efficiency. The reason for this drop is that the heat current remains finite at vanishing power when transmission occurs at more than one energy [21].

Zoom In Zoom Out Reset image size

The upper inset in figure 5(b) shows that at temperatures k_B T_α ≲ κ_α , the Born–Markov approximation breaks down if one of the transition frequencies Ω_± is close to the chemical potential μ_α , cf equation (81). In this case, the step-like behavior of the Fermi–Dirac distribution renders the assumption of a Markovian bath unjustified. By approximating transmission to occur only at the transition frequencies of the system, the Markovian master equations predict a step-like behavior of heat currents and power. The transmission approach shows that the step-like feature is smeared out by the transmission function, which is a continuous function of energy.

In this section, we provide an example of a time-dependent system. To this end, we consider the heat engine introduced by Kosloff in 1984 [67] with the system Hamiltonian

$$\begin{equation}{\hat{H}}_{\text{S}}(t)={{\Omega}}_{\text{c}}{\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{c}}+{{\Omega}}_{\text{h}}{\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{h}}+g({\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{c}}{\text{e}}^{-\text{i}\varpi t}+{\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{h}}{\text{e}}^{\text{i}\varpi t}),\end{equation} \tag{ 94 }$$

where the bosonic annihilation and creation operators fulfill the standard commutation relations given in equation (2), and the frequency of the external drive reads

$$\begin{equation}\varpi ={{\Omega}}_{\text{h}}-{{\Omega}}_{\text{c}-2{\Delta}}.\end{equation} \tag{ 95 }$$

This Hamiltonian describes two bosonic modes with different frequencies that are coupled via a time-dependent term with a detuning quantified by Δ. For a physical implementation of this Hamiltonian based on a superconducting circuit, see reference [68]. The reservoirs and their coupling to the system are described by

$$\begin{equation}{\hat{H}}_{\alpha }=\sum\limits _{l}\enspace {\varepsilon }_{\alpha ,l}{\hat{b}}_{\alpha ,l}^{{\dagger}}{\hat{b}}_{\alpha ,l},\hspace{25.0pt}{\hat{V}}_{\alpha }=\sum\limits _{l}\enspace {t}_{\alpha ,l}({\hat{a}}_{\alpha }^{{\dagger}}{\hat{b}}_{\alpha ,l}+{\hat{b}}_{\alpha ,l}^{{\dagger}}{\hat{a}}_{\alpha }),\end{equation} \tag{ 96 }$$

with α = c, h and we consider a vanishing chemical potential for the reservoirs, μ_α = 0.

From equations (50) and (96), we identify the following system and bath operators

$$\begin{equation}{\hat{S}}_{\alpha ,1}={\hat{a}}_{\alpha }^{{\dagger}},\hspace{25.0pt}{\hat{S}}_{\alpha ,-1}={\hat{a}}_{\alpha },\end{equation} \tag{ 97 }$$

and

$$\begin{equation}{\hat{B}}_{\alpha ,1}=\sum\limits _{l}\enspace {t}_{\alpha ,l}{\hat{b}}_{\alpha ,l},\hspace{25.0pt}{\hat{B}}_{\alpha ,-1}=\sum\limits _{l}\enspace {t}_{\alpha ,l}{\hat{b}}_{\alpha ,l}^{{\dagger}}.\end{equation} \tag{ 98 }$$

We find that the Born–Markov approximations are justified when (see appendix C.1 and references [12, 13])

$$\begin{equation}{\kappa }_{\alpha }\ll {\omega }_{j},\end{equation} \tag{ 99 }$$

where the inequality should hold for all values of α and j. This ensures that the Bose–Einstein distribution is flat around the transition energies over the energy scale κ_α . Assuming a flat bath spectral density, the same holds for the rates

$$\begin{equation}{{\Gamma}}_{1}^{\alpha }(-\omega )={\kappa }_{\alpha }{n}_{\text{B}}^{\alpha }(\omega ),\hspace{25.0pt}{{\Gamma}}_{-1}^{\alpha }(\omega )={\kappa }_{\alpha }[1-{n}_{\text{B}}^{\alpha }(\omega )],\end{equation} \tag{ 100 }$$

where ω > 0 as we are dealing with Bosons with zero chemical potential.

The time-dependence in equation (94) can be removed by a suitable unitary transformation, see reference [13], where this model was used to compare local and global master equations. Here, we use the model to illustrate how a time-dependent Hamiltonian affects our master equation and we thus remain in the lab frame. As discussed in section 3, the average Hamiltonian defined in equation (52) is an important quantity. For our system, we find (see appendix C.2)

$$\begin{equation}{\hat{H}}_{\text{av}}={{\Omega}}_{\text{c}}^{+}{\hat{a}}_{+}^{{\dagger}}{\hat{a}}_{+}+{{\Omega}}_{\text{c}}^{-}{\hat{a}}_{-}^{{\dagger}}{\hat{a}}_{-},\end{equation} \tag{ 101 }$$

where the eigenmodes are given in equations (4) and (5) and we introduced the frequencies

$$\begin{equation}{{\Omega}}_{\text{c}}^{\pm }={{\Omega}}_{\text{c}}+{\Delta}\pm \sqrt{{{\Delta}}^{2}+{g}^{2}},\hspace{25.0pt}{{\Omega}}_{\text{h}}^{\pm }={{\Omega}}_{\text{h}}-{\Delta}\pm \sqrt{{{\Delta}}^{2}+{g}^{2}}.\end{equation} \tag{ 102 }$$

We note that this choice for the average Hamiltonian is not unique. Indeed, ${\hat{H}}_{\text{av}}^{\prime }={{\Omega}}_{\text{h}}^{+}{\hat{a}}_{+}^{{\dagger}}{\hat{a}}_{+}+{{\Omega}}_{\text{h}}^{-}{\hat{a}}_{-}^{{\dagger}}{\hat{a}}_{-}$ does also fulfill the defining relation given in equation (52).

The Fourier coefficients of the bath operators given in equation (98) are determined by

$$\begin{equation}\begin{aligned}\hfill & {\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{a}}_{\text{c}}{\hat{U}}_{\text{S}}(t)={\text{e}}^{-\text{i}t{{\Omega}}_{\text{c}}^{+}}\enspace \mathrm{sin}(\theta /2){\hat{a}}_{+}+{\text{e}}^{-\text{i}t{{\Omega}}_{\text{c}}^{-}}\enspace \mathrm{cos}(\theta /2){\hat{a}}_{-},\hfill \\ \hfill & {\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{a}}_{\text{h}}{\hat{U}}_{\text{S}}(t)={\text{e}}^{-\text{i}t{{\Omega}}_{\text{h}}^{+}}\enspace \mathrm{cos}(\theta /2){\hat{a}}_{+}-{\text{e}}^{-\text{i}t{{\Omega}}_{\text{h}}^{-}}\enspace \mathrm{sin}(\theta /2){\hat{a}}_{-},\hfill \end{aligned}\end{equation} \tag{ 103 }$$

where θ is defined in equation (5) (where equation (95) determines Δ). From these equations (and their Hermitian conjugates) we can identify the transition rates and jump operators. As we have seen for the fermionic heat engine, it is sufficient to consider only the frequencies and operators related to particles leaving the system (the others can be treated separately and equivalently). These transition frequencies read

$$\begin{equation}\left\{{\omega }_{j}\right\}=\left\{{{\Omega}}_{\text{c}}^{+},\enspace {{\Omega}}_{\text{c}}^{-},\enspace {{\Omega}}_{\text{h}}^{+},\enspace {{\Omega}}_{\text{h}}^{-}\right\}.\end{equation} \tag{ 104 }$$

As discussed above, these can be written as ν_j + lϖ, where ν_j denotes a transition frequency of the average Hamiltonian (i.e. ${{\Omega}}_{\text{c}}^{\pm }$ for equation (101)) and l = 0, 1. The corresponding jump operators read

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{c},-1;j}\right\}=\left\{\mathrm{sin}(\theta /2){\hat{a}}_{+},\enspace \mathrm{cos}(\theta /2){\hat{a}}_{-},\enspace 0,\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{h},-1;j}\right\}=\left\{0,\enspace 0,\enspace \mathrm{cos}(\theta /2){\hat{a}}_{+},\enspace -\mathrm{sin}(\theta /2){\hat{a}}_{-}\right\}.\hfill \end{aligned}\end{equation} \tag{ 105 }$$

Again, we have two choices for grouping the frequencies which result in the local and global approach respectively. We may group ${{\Omega}}_{\alpha }^{\pm }$ into different sets or assign them to the same set. Grouping them into different sets results in a single frequency per set such that {ω_q } = {ω_j } and $\left\{{\hat{S}}_{\alpha ,k;q}(t)\right\}=\left\{{\hat{S}}_{\alpha ,k;j}\enspace \mathrm{exp}(-\text{i}{\omega }_{j}t)\right\}$ (cf equation (61)). Inserting these quantities into equation (65), we find the dissipator (in the interaction picture)

$$\begin{equation}{\tilde{\mathcal{L}}}_{\alpha }^{\text{g}}=\sum\limits _{\sigma =\pm }\frac{{\kappa }_{\alpha }^{\sigma }}{2}\left\{{n}_{\text{B}}^{\alpha }({{\Omega}}_{\alpha }^{\sigma })\mathcal{D}[{\hat{a}}_{\sigma }^{{\dagger}}]+[{n}_{\text{B}}^{\alpha }({{\Omega}}_{\alpha }^{\sigma })+1]\mathcal{D}[{\hat{a}}_{\sigma }]\right\},\end{equation} \tag{ 106 }$$

with ${\kappa }_{\alpha }^{\sigma }$ being defined in equation (22). As outlined in appendix A, a single-thermodynamic Hamiltonian is usually not sufficient for the time-dependent scenario. Furthermore, even when a secular approximation is performed, a rescaling needs to be performed [62]. In the present case, the rescaling is obtained starting from ${\hat{H}}_{\text{av}}$ in equation (101) and rescaling the relevant transition frequencies to ω_q . It turns out that it is sufficient to consider two thermodynamic Hamiltonians, one for each reservoir

$$\begin{equation}{\hat{H}}_{\text{TD}}^{\alpha }={{\Omega}}_{\alpha }^{+}{\hat{a}}_{+}^{{\dagger}}{\hat{a}}_{+}+{{\Omega}}_{\alpha }^{-}{\hat{a}}_{-}^{{\dagger}}{\hat{a}}_{-}.\end{equation} \tag{ 107 }$$

The expression for the heat current given in equation (A.5) then reduces to

$$\begin{equation}{J}_{\alpha }=\mathrm{Tr}\left\{{\hat{H}}_{\text{TD}}^{\alpha }{\tilde{\mathcal{L}}}_{\alpha }^{\text{g}}{\hat{\rho }}_{\text{S}}\right\}.\end{equation} \tag{ 108 }$$

Note that the two thermodynamic Hamiltonians correspond to two possible choices for ${\hat{H}}_{\text{av}}$.

In this scenario, the external power can no longer be accessed by the standard expression given in equation (33). As a consequence of the secular approximation, this quantity evaluates to zero [13]. The time-averaged power can however be recovered by relying on the first law as in equation (A.6).

Alternatively, we may group ${{\Omega}}_{\alpha }^{\pm }$ into the same set. This grouping is justified as long as g, Δ ≪ Ω_c, Ω_h (in analogy to equation (99)). A natural choice for the set frequencies then reads

$$\begin{equation}\left\{{\omega }_{q}\right\}=\left\{{{\Omega}}_{\text{c}}+{\Delta},\enspace {{\Omega}}_{\text{h}}-{\Delta}\right\},\end{equation} \tag{ 109 }$$

with the jump operators

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{c},-1;q}(t)\right\}=\left\{{\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{a}}_{\text{c}}{\hat{U}}_{\text{S}}(t),\enspace 0\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{h},-1;q}(t)\right\}=\left\{0,\enspace {\hat{U}}_{\text{S}}^{{\dagger}}(t){\hat{a}}_{\text{h}}{\hat{U}}_{\text{S}}(t)\right\}.\hfill \end{aligned}\end{equation} \tag{ 110 }$$

Upon returning to the Schrödinger picture, these quantities result in the local dissipator

$$\begin{equation}\begin{aligned}\hfill & {\mathcal{L}}_{\text{c}}^{\text{l}}={\kappa }_{\text{c}}\left\{{n}_{\text{B}}({{\Omega}}_{\text{c}}+{\Delta})\mathcal{D}[{\hat{a}}_{\text{c}}^{{\dagger}}]+[{n}_{\text{B}}({{\Omega}}_{\text{c}}+{\Delta})+1]\mathcal{D}[{\hat{a}}_{\text{c}}]\right\},\hfill \\ \hfill & {\mathcal{L}}_{\text{h}}^{\text{l}}={\kappa }_{\text{h}}\left\{{n}_{\text{B}}({{\Omega}}_{\text{h}}-{\Delta})\mathcal{D}[{\hat{a}}_{\text{h}}^{{\dagger}}]+[{n}_{\text{B}}({{\Omega}}_{\text{h}}-{\Delta})+1]\mathcal{D}[{\hat{a}}_{\text{h}}]\right\}.\hfill \end{aligned}\end{equation} \tag{ 111 }$$

Two thermodynamic Hamiltonians can again be obtained by re-scaling the transition frequencies of the average Hamiltonian in equation (101). However, due to the local structure of the jump operators, we can find a single thermodynamic Hamiltonian in the Schrödinger picture (cf appendix A)

$$\begin{equation}{\hat{H}}_{\text{TD}}=({{\Omega}}_{\text{c}}+{\Delta}){\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{c}}+({{\Omega}}_{\text{h}}-{\Delta}){\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{h}}.\end{equation} \tag{ 112 }$$

From equation (A.9), we recover the usual expression for the power produced by the system

$$\begin{equation}{P}_{\text{S}}=-{P}_{\text{ext}}=\text{i}\enspace \mathrm{Tr}\left\{[{\hat{H}}_{\text{TD}},{\hat{H}}_{\text{S}}(t)]{\hat{\rho }}_{\text{S}}\right\}=-\mathrm{Tr}\left\{[{\partial }_{t}{\hat{H}}_{\text{S}}(t)]{\hat{\rho }}_{\text{S}}\right\}.\end{equation} \tag{ 113 }$$

For Δ = 0, the thermodynamic Hamiltonian in equation (112) was used before in reference [56] and has an intuitive explanation: for the thermodynamic bookkeeping, we neglect the coupling energy between the bosonic modes, just as we neglect the coupling energy between system and bath.

In figure 6, power, heat current, and efficiency of the bosonic heat engine for a resonant drive are illustrated. The local and global approaches for this scenario were already compared in reference [13], where exact numerics for finite reservoirs was used as a benchmark. In agreement with these results, we find that the local and global approaches together reproduce exact solutions over the full range of parameters where the Born–Markov approximation is justified. This implies that our master equation, which reduces either to the local or to the global approach, describes this scenario very well. We note that the PERLind approach reproduces the transmission approach extremely well. The disadvantage of the PERLind approach is that it is not thermodynamically consistent. The disadvantage of our master equation is that it exhibits a discontinuity when the frequency grouping is adapted (i.e. when changing from the local to the global approach).

Zoom In Zoom Out Reset image size

In figure 7, the performance of the heat engine at finite detuning Δ is investigated. For small detuning, we find that the global approach correctly reproduces the g → 0 limit, as an exact degeneracy is no longer reached in this limit (panel (a)). As Δ increases, the low-g behavior is increasingly well captured by the global approach, because the secular approximation becomes better. For a fixed coupling strength g, we find a decrease in the heat engine performance as Δ is increased (panel (b)). This behavior is expected, because energy transfer between the drive and the system works best on resonance.

Zoom In Zoom Out Reset image size

In the previous examples, the derived master equation reduces to the well known global or local approaches. In this last example, we consider a system that has more transition frequencies, such that the master equation may differ from both the local as well as the global approach. To this end, we consider a spinless double quantum dot with Coulomb interactions described by the Hamiltonian

$$\begin{align}\hfill {\hat{H}}_{\text{S}}& ={\Omega}{\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{\Omega}{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}+g({\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{R}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{L}})+U{\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}\hfill \\ \hfill & ={{\Omega}}_{-}{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}+{{\Omega}}_{+}{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}+U{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}.\hfill \end{align} \tag{ 114 }$$

Apart from the interaction term, the system as well as the bath properties are identical to section 4 (with Ω ≡ Ω_L = Ω_R and Ω_± = Ω ± g) where all notation used in this section is defined. From equation (114), we see that the Hamiltonian is diagonal in the occupation basis of the ± modes. The transition frequencies and jump operators are obtained from the equations

$$\begin{align}\hfill {\text{e}}^{\text{i}{\hat{H}}_{\text{S}}t}{\hat{d}}_{\text{R}}{\text{e}}^{-\text{i}{\hat{H}}_{\text{S}}t}& =\frac{{\text{e}}^{-\text{i}{{\Omega}}_{-}t}}{\sqrt{2}}(1-{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}){\hat{d}}_{-}+\frac{{\text{e}}^{-\text{i}({{\Omega}}_{-}+U)t}}{\sqrt{2}}{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}{\hat{d}}_{-}\hfill \\ \hfill & \quad +\frac{{\text{e}}^{-\text{i}{{\Omega}}_{+}t}}{\sqrt{2}}(1-{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}){\hat{d}}_{+}+\frac{{\text{e}}^{-\text{i}({{\Omega}}_{+}+U)t}}{\sqrt{2}}{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}{\hat{d}}_{+},\hfill \\ \hfill {\text{e}}^{\text{i}{\hat{H}}_{\text{S}}t}{\hat{d}}_{\text{L}}{\text{e}}^{-\text{i}{\hat{H}}_{\text{S}}t}& =-\frac{{\text{e}}^{-\text{i}{{\Omega}}_{-}t}}{\sqrt{2}}(1-{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}){\hat{d}}_{-}-\frac{{\text{e}}^{-\text{i}({{\Omega}}_{-}+U)t}}{\sqrt{2}}{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}{\hat{d}}_{-}\hfill \\ \hfill & \quad +\frac{{\text{e}}^{-\text{i}{{\Omega}}_{+}t}}{\sqrt{2}}(1-{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}){\hat{d}}_{+}+\frac{{\text{e}}^{-\text{i}({{\Omega}}_{+}+U)t}}{\sqrt{2}}{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}{\hat{d}}_{+}.\hfill \end{align} \tag{ 115 }$$

Considering only transition frequencies that correspond to particles leaving the system, we find

$$\begin{equation}\left\{{\omega }_{j}\right\}=\left\{{{\Omega}}_{-},\enspace {{\Omega}}_{+},\enspace {{\Omega}}_{-}+U,\enspace {{\Omega}}_{+}+U\right\},\end{equation} \tag{ 116 }$$

and the jump operators

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{R},-1;j}\right\}=\left\{\frac{1-{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}}{\sqrt{2}}{\hat{d}}_{-},\enspace \frac{1-{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}}{\sqrt{2}}{\hat{d}}_{+},\enspace \frac{{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}}{\sqrt{2}}{\hat{d}}_{-},\enspace \frac{{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}}{\sqrt{2}}{\hat{d}}_{+}\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},-1;j}\right\}=\left\{-\frac{1-{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}}{\sqrt{2}}{\hat{d}}_{-},\enspace \frac{1-{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}}{\sqrt{2}}{\hat{d}}_{+},\enspace -\frac{{\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}}{\sqrt{2}}{\hat{d}}_{-},\enspace \frac{{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}}{\sqrt{2}}{\hat{d}}_{+}\right\}.\hfill \end{aligned}\end{equation} \tag{ 117 }$$

In the global approach, we choose a different set for each transition frequency. This results in the dissipator

$$\begin{align}\hfill {\mathcal{L}}_{\alpha }^{\text{g}}& =\sum\limits _{\sigma =\pm }\frac{{\kappa }_{\alpha }}{2}\left\{{n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma })\mathcal{D}[(1-{\hat{d}}_{\bar{\sigma }}^{{\dagger}}{\hat{d}}_{\bar{\sigma }}){\hat{d}}_{\sigma }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma })]\mathcal{D}[(1-{\hat{d}}_{\bar{\sigma }}^{{\dagger}}{\hat{d}}_{\bar{\sigma }}){\hat{d}}_{\sigma }]\right.\hfill \\ \hfill & \quad \left.\times {n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma }+U)\mathcal{D}[{\hat{d}}_{\bar{\sigma }}^{{\dagger}}{\hat{d}}_{\bar{\sigma }}{\hat{d}}_{\sigma }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({{\Omega}}_{\sigma }+U)]\mathcal{D}[{\hat{d}}_{\bar{\sigma }}^{{\dagger}}{\hat{d}}_{\bar{\sigma }}{\hat{d}}_{\sigma }]\right\},\hfill \end{align} \tag{ 118 }$$

where $\bar{\sigma }\ne \sigma$. As usual, the thermodynamic Hamiltonian in the global approach is equal to ${\hat{H}}_{\text{S}}$. This frequency grouping is justified when

$$\begin{equation}\vert {\omega }_{j}-{\omega }_{{j}^{\prime }}\vert \gg {\kappa }_{\alpha },\end{equation} \tag{ 119 }$$

for all j, j', and α.

In the local approach, we group all transition frequencies into a single set. A natural set frequency is then the average frequency ω_q = Ω + U/2. The jump operators are given by summing the jump operators in equation (117) over j, which results in the local dissipator

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{l}}={\kappa }_{\alpha }\left\{{n}_{\text{F}}^{\alpha }({\Omega}+U/2)\mathcal{D}[{\hat{d}}_{\alpha }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({\Omega}+U/2)]\mathcal{D}[{\hat{d}}_{\alpha }]\right\}.\end{equation} \tag{ 120 }$$

The thermodynamic Hamiltonian, obtained by rescaling all transitions to Ω + U/2, is given by

$$\begin{equation}{\hat{H}}_{\text{TD}}=({\Omega}+U/2)({\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}})=({\Omega}+U/2)({\hat{d}}_{+}^{{\dagger}}{\hat{d}}_{+}+{\hat{d}}_{-}^{{\dagger}}{\hat{d}}_{-}).\end{equation} \tag{ 121 }$$

We note that apart from the choice of the transition frequency, the dissipator and thermodynamic Hamiltonian are identical to the one obtained for U = 0, cf equations (92) and (93). The frequency grouping for the local approach is justified as long as

$$\begin{equation}\vert {n}_{\text{F}}^{\alpha }({\omega }_{j})-{n}_{\text{F}}^{\alpha }({\omega }_{{j}^{\prime }})\vert \ll 1,\end{equation} \tag{ 122 }$$

for all j, j', and α.

In addition to the local and global approaches, we consider the low-g frequency grouping, where Ω_± can be replaced by Ω. This results in the set frequencies

$$\begin{equation}\left\{{\omega }_{q}\right\}=\left\{{\Omega},\enspace {\Omega}+U\right\},\end{equation} \tag{ 123 }$$

with the corresponding jump operators

$$\begin{equation}\begin{aligned}\hfill & \left\{{\hat{S}}_{\text{R},-1;q}\right\}=\left\{(1-{\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}){\hat{d}}_{\text{R}},\enspace {\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}{\hat{d}}_{\text{R}}\right\},\hfill \\ \hfill & \left\{{\hat{S}}_{\text{L},-1;q}\right\}=\left\{(1-{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}){\hat{d}}_{\text{L}},\enspace {\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}{\hat{d}}_{\text{L}}\right\}.\hfill \end{aligned}\end{equation} \tag{ 124 }$$

Note that the jump operators still locally change the number of electrons. However, the jumps are now dependent on the occupancy of the other dot due to the Coulomb interaction. Inserting these quantities into equation (60) results in the dissipator

$$\begin{align}\hfill {\mathcal{L}}_{\alpha }^{\text{sl}}& ={\kappa }_{\alpha }\left\{{n}_{\text{F}}^{\alpha }({\Omega})\mathcal{D}[(1-{\hat{d}}_{\bar{\alpha }}^{{\dagger}}{\hat{d}}_{\bar{\alpha }}){\hat{d}}_{\alpha }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({\Omega})]\mathcal{D}[(1-{\hat{d}}_{\bar{\alpha }}^{{\dagger}}{\hat{d}}_{\bar{\alpha }}){\hat{d}}_{\alpha }]\right.\hfill \\ \hfill & \quad \left.+{n}_{\text{F}}^{\alpha }({\Omega}+U)\mathcal{D}[{\hat{d}}_{\bar{\alpha }}^{{\dagger}}{\hat{d}}_{\bar{\alpha }}{\hat{d}}_{\alpha }^{{\dagger}}]+[1-{n}_{\text{F}}^{\alpha }({\Omega}+U)]\mathcal{D}[{\hat{d}}_{\bar{\alpha }}^{{\dagger}}{\hat{d}}_{\bar{\alpha }}{\hat{d}}_{\alpha }]\right\},\hfill \end{align} \tag{ 125 }$$

where $\bar{\alpha }\ne \alpha$. The thermodynamic Hamiltonian corresponding to this frequency grouping reads

$$\begin{equation}{\hat{H}}_{\text{TD}}={\Omega}({\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}})+U{\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}.\end{equation} \tag{ 126 }$$

This frequency grouping is justified for

$$\begin{equation}\vert {n}_{\text{F}}^{\alpha }({{\Omega}}_{+})-{n}_{\text{F}}^{\alpha }({{\Omega}}_{-})\vert \ll 1,\quad \text{and}\quad \vert {n}_{\text{F}}^{\alpha }({{\Omega}}_{+}+U)-{n}_{\text{F}}^{\alpha }({{\Omega}}_{-}+U)\vert \ll 1.\end{equation} \tag{ 127 }$$

To ensure that the frequencies in different sets obey |ω_j − ω_j'| ≫ κ_α , one may expect the additional condition U ≫ g, κ_α . However, because no coherences can build up between states with a different total number of electrons in the system, this additional condition is not required. The semi-local approach thus enjoys a strictly larger regime of validity than the local approach (cf equation (122)), remaining valid even for a vanishing interaction strength U = 0.

Due to the interaction term, the transmission approach no longer applies. As a benchmark, we instead use the PERLind approach, see appendix D. As mentioned above, this approach is expected to provide accurate results whenever the Born–Markov approximations are justified, an expectation that was confirmed in the examples we considered above.

Figure 8 illustrates heat current and power in the interacting double quantum dot. We find that at finite interaction strengths, the semi-local approach should be employed instead of the local approach for small coupling g. In particular, comparing figure 8(a) with figure 5(a) we find a very similar interplay between the semi-local and global approach at U ≠ 0 as we observed between the local and global approaches for U = 0. We furthermore find that the semi-local approach agrees with the PERLind approach for any value of the interaction strength, as long as the coupling g is sufficiently small to respect equation (127), see figure 8(b). These results highlight the importance to go beyond the local and global approaches for systems that have more than two competing transition frequencies.

Zoom In Zoom Out Reset image size

The bath operators given in equation (80) result in the bath correlation functions (cf equation (55))

$$\begin{equation}\begin{aligned}\hfill & {C}_{1}^{\alpha }=\frac{{\kappa }_{\alpha }}{2}{\text{e}}^{\text{i}s{\mu }_{\alpha }}\left[\delta (s)-\text{i}\frac{{k}_{\text{B}}{T}_{\alpha }}{\mathrm{sinh}(\pi s{k}_{\text{B}}{T}_{\alpha })}\right],\hfill \\ \hfill & {C}_{-1}^{\alpha }=\frac{{\kappa }_{\alpha }}{2}{\text{e}}^{-\text{i}s{\mu }_{\alpha }}\left[\delta (s)-\text{i}\frac{{k}_{\text{B}}{T}_{\alpha }}{\mathrm{sinh}(\pi s{k}_{\text{B}}{T}_{\alpha })}\right],\hfill \end{aligned}\end{equation} \tag{ B.1 }$$

where we introduced the bath spectral density

$$\begin{equation}\frac{{\kappa }_{\alpha }}{2\pi }=\sum\limits _{l}\enspace {t}_{\alpha ,l}^{2}\delta (\omega -{\varepsilon }_{\alpha ,l}),\end{equation} \tag{ B.2 }$$

quantifying the system-bath coupling and which we assume to be constant as a function of ω. We note that ${C}_{k}^{\alpha }\equiv {C}_{k,k}^{\alpha }$ and ${C}_{k,{k}^{\prime }}^{\alpha }\propto {\delta }_{k,{k}^{\prime }}$. The bath correlation functions decay on the time-scale β_α . The bath correlation time can thus be identified by τ_B = max{β_L, β_R}. The system (in the interaction picture) changes on the time-scale τ_S = min{1/κ_L, 1/κ_R}. Thus, the Born–Markov approximations are expected to hold for temperatures that are much higher than the coupling between system and bath. Furthermore, from the Redfield equation in equation (53), together with equation (B.1), we find that the integral over s involves factors exp[is(ω_j ± μ_α )]. Whenever |ω_j ± μ_α | ≫ k_B T_α , the rapid oscillations of these factors imply that only values of s ≲ 1/(ω_j ± μ_α ) contribute to the integral. For larger values of s, the term k_B T_α /sinh(πsk_B T_α ) varies on a time-scale much slower than the oscillations. This implies that the Born–Markov approximations are also justified in the regime |ω_j ± μ_α | ≫ k_B T_α and min{|ω_j ± μ_L|, |ω_j ± μ_R|} ≫ 1/τ_S. Now in case 1/τ_S ≪ |ω_j ± μ| (dropping the bath index for ease of notation), either we have k_B T ≪ |ω_j ± μ|, such that the Born–Markov approximations are justified, or we have k_B T ≳|ω_j ± μ| which implies 1/τ_S ≪ k_B T, also ensuring the Born–Markov approximation. In the first case, the bath-correlation functions oscillate to zero over times much smaller than τ_S, in the latter case the bath correlation functions decay much faster than τ_S. As a consequence, the Born–Markov approximations are fulfilled whenever κ is much smaller than either k_B T, or |ω_j ± μ| as expressed in equation (81).

The Fourier transforms of the bath correlation functions are given by equation (82). We reach the same conclusions about the validity of the Born–Markov approximation by demanding that these functions are flat over the energy scale κ (the inverse of τ_S). For κ ≪ k_B T, the Fermi–Dirac distribution is everywhere slowly varying compared to κ. For κ ≪ |ω − μ| and k_B T ≪ |ω − μ|, the Fermi function is either very close to zero or one, depending on the sign of ω − μ, and remains so over the energy scale of κ.

For the double dot system with non-interacting fermions the Landauer-like formula provides the following steady state heat current (out of the left lead)

$$\begin{equation}{J}_{\text{L}}^{\text{t}}=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega \mathcal{T}(\omega )(\omega -{\mu }_{\text{L}})[{n}_{\text{F}}^{\text{L}}(\omega )-{n}_{\text{F}}^{\text{R}}(\omega )],\end{equation} \tag{ B.3 }$$

and power

$$\begin{equation}{P}_{\text{S}}^{\text{t}}=\frac{({\mu }_{\text{R}}-{\mu }_{\text{L}})}{2\pi }{\int }_{-\infty }^{\infty }\mathrm{d}\omega \mathcal{T}(\omega )[{n}_{\text{F}}^{\text{L}}(\omega )-{n}_{\text{F}}^{\text{R}}(\omega )],\end{equation} \tag{ B.4 }$$

with the transmission function

$$\begin{equation}\mathcal{T}(\omega )=\frac{{g}^{2}{\kappa }_{\text{L}}{\kappa }_{\text{R}}}{\vert \left(\omega -{{\Omega}}_{\text{L}}+\text{i}\frac{{\kappa }_{\text{L}}}{2}\right)\left(\omega -{{\Omega}}_{\text{R}}+\text{i}\frac{{\kappa }_{\text{R}}}{2}\right)-{g}^{2}{\vert }^{2}}.\end{equation} \tag{ B.5 }$$

Note that the transmission function is the same as for the bosonic system, cf equation (10).

For the global and local approaches, the definitions for the heat current and the power are given in equations (70) and (71), with ${\hat{N}}_{\text{S}}={\hat{d}}_{\text{L}}^{{\dagger}}{\hat{d}}_{\text{L}}+{\hat{d}}_{\text{R}}^{{\dagger}}{\hat{d}}_{\text{R}}$. For the local approach, we find

$$\begin{equation}{J}_{\text{L}}^{\text{l}}=\frac{(\bar{{\Omega}}-{\mu }_{\text{L}})}{{\kappa }_{\text{L}}+{\kappa }_{\text{R}}}\frac{4{g}^{2}{\kappa }_{\text{L}}{\kappa }_{\text{R}}[{n}_{\text{F}}^{\text{L}}(\bar{{\Omega}})-{n}_{\text{F}}^{\text{R}}(\bar{{\Omega}})]}{4{g}^{2}+{\kappa }_{\text{L}}{\kappa }_{\text{R}}+16{{\Delta}}^{2}{\kappa }_{\text{L}}{\kappa }_{\text{R}}/{({\kappa }_{\text{L}}+{\kappa }_{\text{R}})}^{2}},\end{equation} \tag{ B.6 }$$

and

$$\begin{equation}{P}_{\text{S}}^{\text{l}}\equiv -{P}_{\text{L}}^{\text{l}}-{P}_{\text{R}}^{\text{l}}=\frac{({\mu }_{\text{R}}-{\mu }_{\text{L}})}{(\bar{{\Omega}}-{\mu }_{\text{L}})}{J}^{\text{l}}.\end{equation} \tag{ B.7 }$$

This tight-coupling condition between heat current and power (which can be inferred from the form of ${\hat{H}}_{\text{TD}}$, cf equation (93)) is a direct result of approximating transmission to happen at a single energy.

For the global approach, we find

$$\begin{equation}{J}_{\text{L}}^{\text{g}}=\sum\limits _{\sigma =\pm }({{\Omega}}_{\sigma }-{\mu }_{\text{L}})\frac{{\kappa }_{\text{L}}^{\sigma }{\kappa }_{\text{R}}^{\sigma }}{{\kappa }_{\text{L}}^{\sigma }+{\kappa }_{\text{R}}^{\sigma }}[{n}_{\text{F}}^{\text{L}}({{\Omega}}_{\sigma })-{n}_{\text{F}}^{\text{R}}({{\Omega}}_{\sigma })],\end{equation} \tag{ B.8 }$$

and

$$\begin{equation}{P}_{\text{S}}^{\text{g}}\equiv -{P}_{\text{L}}^{\text{g}}-{P}_{\text{R}}^{\text{g}}=({\mu }_{\text{R}}-{\mu }_{\text{L}})\sum\limits _{\sigma =\pm }\frac{{\kappa }_{\text{L}}^{\sigma }{\kappa }_{\text{R}}^{\sigma }}{{\kappa }_{\text{L}}^{\sigma }+{\kappa }_{\text{R}}^{\sigma }}[{n}_{\text{F}}^{\text{L}}({{\Omega}}_{\sigma })-{n}_{\text{F}}^{\text{R}}({{\Omega}}_{\sigma })].\end{equation} \tag{ B.9 }$$

Finally, in the PERLind approach the dissipators read

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{p}}=\sum\limits _{k=-1,1}\mathcal{D}[{\hat{J}}_{\alpha ,k}],\end{equation} \tag{ B.10 }$$

with the jump operators

$$\begin{equation}\begin{aligned}\hfill {\hat{J}}_{\text{R,}-1}& =\sqrt{{\kappa }_{\text{R}}^{+}[1-{n}_{\text{F}}^{\text{R}}({{\Omega}}_{+})]}{\hat{d}}_{+}+\sqrt{{\kappa }_{\text{R}}^{-}[1-{n}_{\text{F}}^{\text{R}}({{\Omega}}_{-})]}{\hat{d}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{R,}1}& =\sqrt{{\kappa }_{\text{R}}^{+}{n}_{\text{F}}^{\text{R}}({{\Omega}}_{+})}{\hat{d}}_{+}^{{\dagger}}+\sqrt{{\kappa }_{\text{R}}^{-}{n}_{\text{F}}^{\text{R}}({{\Omega}}_{-})}{\hat{d}}_{-}^{{\dagger}},\hfill \\ \hfill {\hat{J}}_{\text{L,}-1}& =\sqrt{{\kappa }_{\text{L}}^{+}[1-{n}_{\text{F}}^{\text{L}}({{\Omega}}_{+})]}{\hat{d}}_{+}-\sqrt{{\kappa }_{\text{L}}^{-}[1-{n}_{\text{F}}^{\text{L}}({{\Omega}}_{-})]}{\hat{d}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{L,}1}& =\sqrt{{\kappa }_{\text{L}}^{+}{n}_{\text{F}}^{\text{L}}({{\Omega}}_{+})}{\hat{d}}_{+}^{{\dagger}}-\sqrt{{\kappa }_{\text{L}}^{-}{n}_{\text{F}}^{\text{L}}({{\Omega}}_{-})}{\hat{d}}_{-}^{{\dagger}}.\hfill \end{aligned}\end{equation} \tag{ B.11 }$$

In this approach, heat current and power are defined as

$$\begin{equation}{J}_{\text{L}}^{\text{p}}=\mathrm{Tr}\left\{\left({\hat{H}}_{\text{S}}-{\mu }_{\text{L}}{\hat{N}}_{\text{S}}\right){\mathcal{L}}_{\text{L}}^{\text{p}}\hat{\rho }\right\},\end{equation} \tag{ B.12 }$$

and

$$\begin{equation}{P}_{\text{S}}^{\text{p}}=-{\mu }_{\text{L}}\enspace \mathrm{Tr}\left\{{\hat{N}}_{\text{S}}{\mathcal{L}}_{\text{L}}^{\text{p}}\hat{\rho }\right\}-{\mu }_{\text{R}}\enspace \mathrm{Tr}\left\{{\hat{N}}_{\text{S}}{\mathcal{L}}_{\text{R}}^{\text{p}}\hat{\rho }\right\}.\end{equation} \tag{ B.13 }$$

The bath operators in equation (98) result in the bath correlation functions (cf equation (55))

$$\begin{equation}\begin{aligned}\hfill & {C}_{1}^{\alpha }={\int }_{0}^{\infty }\frac{\mathrm{d}\omega }{2\pi }{\kappa }_{\alpha }{\text{e}}^{\text{i}s\omega }{n}_{\text{B}}^{\alpha }(\omega ),\hfill \\ \hfill & {C}_{-1}^{\alpha }={\int }_{0}^{\infty }\frac{\mathrm{d}\omega }{2\pi }{\kappa }_{\alpha }{\text{e}}^{-\text{i}s\omega }[{n}_{\text{B}}^{\alpha }(\omega )+1],\hfill \end{aligned}\end{equation} \tag{ C.1 }$$

where κ_α is given by equation (B.2) and we have ${C}_{k}^{\alpha }\equiv {C}_{k,k}^{\alpha }$ and ${C}_{k,{k}^{\prime }}^{\alpha }\propto {\delta }_{k,{k}^{\prime }}$. We note that these integrals only converge if κ_α → 0 for ω → 0. Here we do not give a specific form of κ_α and directly proceed to investigating the Fourier transforms of the bath correlation functions. For a detailed discussion on the bath correlation functions for an ohmic spectral density, we refer to reference [12]. By Fourier transforming equation (C.1), we find the rates given in equation (100). We note that even if κ_α needs to vanish as ω → 0, we can still assume it to be constant over all energies that the system probes (i.e. all energies where the transmission function is non-zero).

Similarly to the fermionic case, the Born–Markov approximations are justified when the quantities in equation (100) remain constant over the energy interval κ_α around the transition frequencies of the system. This is the case if κ_α ≪ ω_j . We note that for very small temperatures (k_B T ≪ ω_j ), n_B(ω_j ) becomes exponentially small and the dynamics is dominated by the temperature-independent term in equation (100) which corresponds to spontaneous emission. We further note that κ_α ≪ k_B T_α is not a sufficient condition for the Born–Markov approximation to be justified. When κ_α → ω_j , then n_B(ω_j − κ_α ) diverges, no matter the temperature. The conclusions drawn here from the Fourier transform of the bath correlation functions are in complete agreement with the conclusions obtained in reference [12] by considering the bath correlation functions for an ohmic spectral density.

To derive the average Hamiltonian for the system given in equation (94), we first consider the rotating frame determined by the unitary transformation

$$\begin{equation}{\hat{U}}_{\text{r}}={\text{e}}^{\text{i}t[({{\Omega}}_{\text{c}}+{\Delta}){\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{c}}+({{\Omega}}_{\text{h}}-{\Delta}){\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{h}}]}.\end{equation} \tag{ C.2 }$$

In the rotating frame, we find the time-independent Hamiltonian

$$\begin{equation}{\tilde{H}}_{\text{S}}={\hat{U}}_{\text{r}}(t){\hat{H}}_{\text{S}}(t){\hat{U}}_{\text{r}}^{{\dagger}}(t)-\text{i}{\hat{U}}_{\text{r}}(t){\partial }_{t}{\hat{U}}_{\text{r}}^{{\dagger}}(t)={\Delta}({\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{h}}-{\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{c}})+g({\hat{a}}_{\text{h}}^{{\dagger}}{\hat{a}}_{\text{c}}+{\hat{a}}_{\text{c}}^{{\dagger}}{\hat{a}}_{\text{h}}).\end{equation} \tag{ C.3 }$$

The time-evolution operator in the lab frame fulfills the relation

$$\begin{equation}{\hat{U}}_{\text{S}}(t)=\mathcal{T}{\text{e}}^{-\text{i}{\int }_{0}^{t}\mathrm{d}{t}^{\prime }{\hat{H}}_{\text{S}}({t}^{\prime })}={\hat{U}}_{\text{r}}^{{\dagger}}(t){\tilde{U}}_{\text{S}}(t),\end{equation} \tag{ C.4 }$$

where ${\tilde{U}}_{\text{S}}(t)$ denotes the time-evolution operator for the system in the rotating frame.

To obtain the average Hamiltonian, we use the relations

$$\begin{equation}{\text{e}}^{\text{i}{t}_{\text{p}}\varpi {\hat{n}}_{j}}=1,\end{equation} \tag{ C.5 }$$

for t_p = 2π/ϖ, where ${\hat{n}}_{j}$ is any operator that only has integer eigenvalues (such as a photon number operator). With the help of this relation, one can show

$$\begin{equation}{\hat{U}}_{\text{S}}({t}_{\text{p}})={\text{e}}^{-\text{i}{t}_{\text{p}}[{{\Omega}}_{\alpha }^{+}{\hat{a}}_{+}^{{\dagger}}{\hat{a}}_{+}+{{\Omega}}_{\alpha }^{-}{\hat{a}}_{-}^{{\dagger}}{\hat{a}}_{-}]},\end{equation} \tag{ C.6 }$$

which holds both for α = c, h (the frequencies are given in equation (102)), resulting in the two choices for the average Hamiltonian discussed in the main text.

In the rotating frame discussed in the last subsection, the bosonic heat engine looks just like the system discussed in section 2, upon setting $\bar{{\Omega}}=0$ and shifting the energies of the reservoir modes by Ω_c + Δ and Ω_h − Δ respectively (see also reference [13]). We then obtain the Landauer-like formula for the heat current

$$\begin{equation}{J}_{\text{h}}^{\text{t}}=\frac{1}{2\pi }{\int }_{\varpi }^{\infty }\mathrm{d}\omega \mathcal{T}(\omega )\omega [{n}_{\text{B}}^{\text{h}}(\omega )-{n}_{\text{B}}^{\text{c}}(\omega -\varpi )],\end{equation} \tag{ C.7 }$$

and for power we find

$$\begin{equation}{P}_{\text{S}}^{\text{t}}=\frac{1}{2\pi }{\int }_{-{{\Omega}}_{\text{c}}}^{\infty }\mathrm{d}\omega \mathcal{T}(\omega )\varpi [{n}_{\text{B}}^{\text{h}}(\omega )-{n}_{\text{B}}^{\text{c}}(\omega -\varpi )],\end{equation} \tag{ C.8 }$$

with the transmission function

$$\begin{equation}\mathcal{T}(\omega )=\frac{{g}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}}{\vert \left(\omega -{{\Omega}}_{\text{h}}+\text{i}\frac{{\kappa }_{\text{c}}}{2}\right)\left(\omega -{{\Omega}}_{\text{c}}-\varpi +\text{i}\frac{{\kappa }_{\text{h}}}{2}\right)-{g}^{2}{\vert }^{2}}.\end{equation} \tag{ C.9 }$$

For the local approach, we find the heat current

$$\begin{equation}{J}_{\text{h}}^{\text{l}}=\frac{{{\Omega}}_{\text{h}}-{\Delta}}{{\kappa }_{\text{c}}+{\kappa }_{\text{h}}}\frac{4{g}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}-{\Delta})-{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}+{\Delta})]}{4{g}^{2}+{\kappa }_{\text{c}}{\kappa }_{\text{h}}+16{{\Delta}}^{2}{\kappa }_{\text{c}}{\kappa }_{\text{h}}/{({\kappa }_{\text{c}}+{\kappa }_{\text{h}})}^{2}},\end{equation} \tag{ C.10 }$$

and power

$$\begin{equation}{P}^{\text{l}}=\left(1-\frac{{{\Omega}}_{\text{c}}+{\Delta}}{{{\Omega}}_{\text{h}}-{\Delta}}\right){J}_{\text{h}}^{\text{l}}.\end{equation} \tag{ C.11 }$$

The global approach results in the heat current

$$\begin{equation}{J}_{\text{h}}^{\text{g}}=\sum\limits _{\sigma =\pm }{{\Omega}}_{\text{h}}^{\sigma }\frac{{\kappa }_{\text{c}}^{\sigma }{\kappa }_{\text{h}}^{\sigma }}{{\kappa }_{\text{c}}^{\sigma }+{\kappa }_{\text{h}}^{\sigma }}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}^{\sigma })-{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{\sigma })],\end{equation} \tag{ C.12 }$$

and power

$$\begin{equation}{P}_{\text{S}}^{\text{g}}\equiv {J}_{\text{h}}+{J}_{\text{c}}=\sum\limits _{\sigma =\pm }\varpi \frac{{\kappa }_{\text{c}}^{\sigma }{\kappa }_{\text{h}}^{\sigma }}{{\kappa }_{\text{c}}^{\sigma }+{\kappa }_{\text{h}}^{\sigma }}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}^{\sigma })-{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{\sigma })].\end{equation} \tag{ C.13 }$$

Finally, for the PERLind approach, the dissipators read (in the rotating frame)

$$\begin{equation}{\mathcal{L}}_{\alpha }^{\text{p}}=\sum\limits _{k=-1,1}\mathcal{D}[{\hat{J}}_{\alpha ,k}],\end{equation} \tag{ C.14 }$$

with the jump operators

$$\begin{equation}\begin{aligned}\hfill {\hat{J}}_{\text{c,}-1}& =\sqrt{{\kappa }_{\text{c}}^{+}[{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{+})+1]}{\hat{a}}_{+}+\sqrt{{\kappa }_{\text{c}}^{-}[{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{-})+1]}{\hat{a}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{c,}\;1}& =\sqrt{{\kappa }_{\text{c}}^{+}{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{+})}{\hat{a}}_{+}^{{\dagger}}+\sqrt{{\kappa }_{\text{c}}^{-}{n}_{\text{B}}^{\text{c}}({{\Omega}}_{\text{c}}^{-})}{\hat{a}}_{-}^{{\dagger}},\hfill \\ \hfill {\hat{J}}_{\text{h,}-1}& =\sqrt{{\kappa }_{\text{h}}^{+}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}^{+})+1]}{\hat{a}}_{+}-\sqrt{{\kappa }_{\text{h}}^{-}[{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}^{-})+1]}{\hat{a}}_{-},\hfill \\ \hfill {\hat{J}}_{\text{h,}\;1}& =\sqrt{{\kappa }_{\text{h}}^{+}{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{h}}^{+})}{\hat{a}}_{+}^{{\dagger}}-\sqrt{{\kappa }_{\text{h}}^{-}{n}_{\text{B}}^{\text{h}}({{\Omega}}_{\text{g}}^{-})}{\hat{a}}_{-}^{{\dagger}}.\hfill \end{aligned}\end{equation} \tag{ C.15 }$$

In this approach, the heat currents can be obtained from equation (53) by applying the approximations discussed in the supplemental material of reference [30]. For the present system, we obtain

$$\begin{equation}{J}_{\alpha }^{\text{p}}=\sum\limits _{\sigma ,{\sigma }^{\prime }=\pm }\frac{{{\Omega}}_{\alpha }^{\sigma }+{{\Omega}}_{\alpha }^{{\sigma }^{\prime }}}{2}\mathrm{Tr}\left\{\left({\hat{J}}_{\alpha ,1;\sigma }^{{\dagger}}{\hat{J}}_{\alpha ,1;{\sigma }^{\prime }}-{\hat{J}}_{\alpha ,-1;\sigma }^{{\dagger}}{\hat{J}}_{\alpha ,-1;{\sigma }^{\prime }}\right)\tilde{\rho }\right\},\end{equation} \tag{ C.16 }$$

where

$$\begin{equation}{\hat{J}}_{\alpha ,k}={\hat{J}}_{\alpha ,k;+}+{\hat{J}}_{\alpha ,k;-},\end{equation} \tag{ C.17 }$$

and ${\hat{J}}_{\alpha ,k;\sigma }$ can be inferred from equation (C.15) as the term that is proportional to ${\hat{a}}_{\sigma }^{({\dagger})}$. The power is given by

$$\begin{equation}{P}_{\text{S}}^{\text{p}}={J}_{\text{h}}^{\text{p}}+{J}_{\text{c}}^{\text{p}}.\end{equation} \tag{ C.18 }$$