Nonequilibrium thermodynamics in the strong coupling and non-Markovian regime based on a reaction coordinate mapping

Classical thermodynamics is a weak coupling theory in the sense that boundary or surface terms of the system are negligible compared to its bulk or volume properties. This becomes particularly apparent in Maxwell's colloquial description of the zeroth law of thermodynamics: 'all heat is of the same kind'[1]. By this statement he meant that the laws governing the transformation of heat are independent of how we put two different systems into contact—a conclusion which obviously holds only if the influence of this contact can be neglected. Other implications of the weak coupling approximation are, for instance, the extensiveness of internal energy or entropy if we scale the volume of a system [2].

While the weak coupling approximation can be well justified for macroscopic systems due to simple geometric arguments (the surface to volume ratio usually decreases with increasing volume), it is harder to justify in the opposite limit when the volume of the system becomes very small. This, however, is the regime where quantum and stochastic effects dominate. Then, in order to link microscopic theory with thermodynamics, one usually starts with a Hamiltonian of the form

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

where H_S (H_E) is the Hamiltonian of the system (environment) and H_I describes their interaction (the 'contact'). Assuming that the coupling H_I is small, one then performs a perturbative expansion up to second order in H_I which (under the additional assumption that the environment is memory-less) yields a closed and Markovian evolution equation for the system density matrix ${\rho }_{{\rm{S}}}$, known as a (quantum) master equation (ME). We note that similar assumptions are needed to derive (classical) Fokker–Planck or Langevin equations. MEs derived this way can then be shown to have a transparent thermodynamic interpretation [3–5] and they provide the work horse for the field of quantum and stochastic thermodynamics, see [6–10] for recent reviews.

However, many interesting physical effects cannot be captured with such a ME approach and thus, quantum and stochastic thermodynamics is still restricted to a small regime of applicability. Consequently, many groups have started to look at thermodynamics in the strong coupling and non-Markovian regime [11–37]. Though these works present important theoretical cornerstones, they are still far away from providing a satisfactory extension of thermodynamics beyond the weak coupling limit. In particular, if one wishes to address the performance of a steadily working heat engine, the general results derived in [11–21] are not of great help because they either focus on integrated changes of thermodynamic values (e.g., the total heat exchanged in a finite time instead of the rate of heat exchange) and additionally rely on an initially decorrelated system–environment state [11–13] and/or coupling only to a single thermal reservoir [11, 13–17];or they remain very formal [18–21]. Furthermore, model-specific studies are either based on simple or exactly solvable models from the field of quantum transport [22–25] and quantum Brownian motion [26–29], or spin-boson models [30–32] often in combination with specific transformations applicable only to special Hamiltonians (polaron transformations) [31, 33–36]; or the investigations are restricted to numerical studies [37].

The goal of this paper is to close the gap between the general results, which are often hard to apply in practice, and studies restricted to overly specific models. Here, we propose a framework which allows to carry over all concepts known from the weak coupling regime of thermodynamics to the strong coupling and non-Markovian regime (we will in fact see that strong coupling is easier to treat than non-Markovianity)³ . Our framework is general in the sense that it can be applied to an arbitrary (quantum or classical⁴ and even driven) system coupled linearly to an arbitrary number of harmonic oscillator heat baths. Thus, apart from not being able to treat, e.g., fermionic reservoirs at the moment, we capture many relevant situations encountered in the study of small-scale engines.

The general idea is to give up the system–environment partition as it is dictated by the microscopic Hamiltonian (1). Instead, we define a new 'supersystem' which includes this part of the environment which is responsible for strong coupling and non-Markovian effects. By construction the resulting supersystem would be coupled weakly to Markovian residual baths and can be treated within the standard framework of quantum or stochastic thermodynamics. More specifically, to achieve this idea, we will identify a collective degree of freedom in the reservoir which is then incorporated into the description of the system. This collective degree of freedom is known as a reaction coordinate (RC) [38], which will capture non-Markovian and strong coupling effects. For this 'supersystem' (original system plus RC) we will derive the ME as usual allowing for a transparent interpretation of the laws of thermodynamics. We note that a complementary analysis of quantum Otto cycles in the strong coupling regime, also employing the RC mapping, appears in a related work by Newman et al [39]. Apart from the thermodynamic applications we propose here, it has been shown that such a RC mapping and related concepts can provide a very accurate method to investigate the behaviour of open quantum systems for a variety of problems [40–51]. Even more generally, it is possible to apply this method iteratively by including several RCs and in this way one can prove that every non-Markovian environment can be mapped to a Markovian one [52–54].

Apart from adapting this general method to treat heat engines in the strong coupling and non-Markovian regime, we also consider concrete applications. In particular, we find that non-Markovian effects can significantly enhance the steady state efficiency even in the weak coupling regime; a result which—to the best of our knowledge—has not been found before. In addition, we also investigate the relation between our method and the widely used polaron transformation, showing that the RC framework provides a means to justify particular polaron transformed ME (PME), and reduces to it under special circumstances.

Outline: We will start by introducing the general technique of the RC mapping in section 2 as far as it is needed to make the present treatment self-consistent. After having established this tool, we will present the general thermodynamic framework based on it in section 3. Particular applications to devices working out of equilibrium then follow in section 4 (efficiency study of a maser heat engine in the non-Markovian regime) and section 5 (a single electron transistor (strongly) coupled to vibrations). Final remarks about the range of validity, open problems and the thermodynamic interpretation of the method are given in section 6.

We consider an arbitrary system with Hamiltonian H_S(t) coupled linearly via some system operator s to a bath of harmonic oscillators (the coupling to several baths follows straightforwardly from this treatment, see next section). The total Hamiltonian is assumed to have the typical Brownian motion form [55, 56]

$$\begin{eqnarray}&&H={H}_{{\rm{S}}}(t)+\displaystyle \frac{1}{2}\displaystyle \sum _{k}\left[{p}_{k}^{2}+{\omega }_{k}^{2}{\left({x}_{k}-\displaystyle \frac{{c}_{k}}{{\omega }_{k}^{2}}s\right)}^{2}\right]\end{eqnarray} \tag{ 2 }$$

with mass-weighted positions x_k and momenta p_k of the bath fulfilling $[{x}_{k},{p}_{l}]={\rm{i}}{\delta }_{{kl}}$ (we set ${\rm{\hslash }}\equiv 1$ throughout). It is worth to point out that the completion of the square is important for a number of reasons, e.g., to guarantee a thermodynamically stable Hamiltonian for all coupling strengths c_k [57]. In the derivation of MEs one often neglects the quadratic system 'renormalization' term $\tfrac{1}{2}{\sum }_{k}\tfrac{{c}_{k}^{2}}{{\omega }_{k}^{2}}{s}^{2}$ from the beginning, though its contribution is, in principle, of the same order as the Lamb shift term.

An important result of the microscopic theory of Brownian motion is that the effect of the bath on the system can be captured solely by one special function known as the spectral density (SD) of the bath:

$$\begin{eqnarray}&&{J}_{0}(\omega )\equiv \displaystyle \frac{\pi }{2}\displaystyle \sum _{k}\displaystyle \frac{{c}_{k}^{2}}{{\omega }_{k}}\delta (\omega -{\omega }_{k}).\end{eqnarray} \tag{ 3 }$$

The SD is a positive function for $\omega \gt 0$ and must fulfill ${J}_{0}(\omega )\to 0$ for $\omega \to 0$ and $\omega \to \infty$.

Now, the spirit of the RC method is to define the interaction with the collective degrees of the freedom of the bath as one new coordinate: the RC X₁ (see figure 1). Hence, we seek a transformation which maps

$$\begin{eqnarray}&&\displaystyle \sum _{k}{c}_{k}{x}_{k}={\lambda }_{0}{X}_{1},\end{eqnarray} \tag{ 4 }$$

where ${\lambda }_{0}$ is an unspecified parameter so far. More formally, we perform a normal mode transformation of the form

$$\begin{eqnarray}&&{\bf{X}}={\rm{\Lambda }}{\bf{x}},\;\;\;{\bf{P}}={\rm{\Lambda }}{\bf{p}}.\end{eqnarray} \tag{ 5 }$$

Here, we used a vector notation for the collection of original and transformed bath coordinates and momenta, e.g., ${\bf{x}}={({x}_{1},\ldots ,{x}_{k},\ldots ,{x}_{N})}^{{\rm{T}}}$ describes the original bath coordinates. For definiteness we considered a finite number of N bath oscillators (having the limit $N\to \infty$ in mind). Furthermore, Λ is an orthogonal N × N matrix, i.e., ${{\rm{\Lambda }}}^{-1}={{\rm{\Lambda }}}^{{\rm{T}}}$, which guarantees that $[{X}_{k},{P}_{l}]={\rm{i}}{\delta }_{{kl}}$. Thus, Λ has $\tfrac{N(N-1)}{2}$ independent components which are fixed by the requirement that the collection of residual bath oscillators (i.e., all oscillators except the RC itself) is of normal form. This leads to

$$\begin{eqnarray}&&\displaystyle \sum _{k}{\omega }_{k}^{2}{{\rm{\Lambda }}}_{{lk}}{{\rm{\Lambda }}}_{{mk}}={\delta }_{{lm}}{{\rm{\Omega }}}_{l}^{2}\end{eqnarray} \tag{ 6 }$$

for $l\ne 1$ and $m\ne 1$ and allows us to map the Hamiltonian (2) to

$$\begin{eqnarray}&&H^{\prime} ={H}_{{\rm{S}}}(t)+\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{2}{s}^{2}-{\lambda }_{0}{{sX}}_{1}+\displaystyle \frac{1}{2}({P}_{1}^{2}+{{\rm{\Omega }}}_{1}^{2}{X}_{1}^{2})-{X}_{1}\displaystyle \sum _{k\ne 1}{C}_{k}{X}_{k}+\displaystyle \frac{1}{2}\displaystyle \sum _{k\ne 1}({P}_{k}^{2}+{{\rm{\Omega }}}_{k}^{2}{X}_{k}^{2}).\end{eqnarray} \tag{ 7 }$$

Here, we used ${{\rm{\Lambda }}}_{1k}={\lambda }_{0}^{-1}{c}_{k}$ which follows from equation (4). Furthermore, we defined ${{\rm{\Omega }}}_{1}^{2}\equiv {\sum }_{k}{\omega }_{k}^{2}{{\rm{\Lambda }}}_{1k}^{2}$ and ${C}_{k}\equiv -{\sum }_{l}{\omega }_{l}^{2}{{\rm{\Lambda }}}_{{kl}}{{\rm{\Lambda }}}_{k1}$. Finally, the system gets renormalized due to $\delta {{\rm{\Omega }}}_{0}^{2}\equiv {\sum }_{k}{\omega }_{k}^{-2}{c}_{k}^{2}$ and from $[{X}_{1},{P}_{1}]={\rm{i}}$ we can deduce that ${\lambda }_{0}^{2}={\sum }_{k}{c}_{k}^{2}$.

Zoom In Zoom Out Reset image size

At this point we can already recognize an important property of the mapping. Suppose that we scale the coupling coefficients c_k by ${c}_{k}\mapsto \alpha {c}_{k}$ for some $\alpha \in {\mathbb{R}}$. Then, the only parameters influenced by this will be ${\lambda }_{0}$ and $\delta {{\rm{\Omega }}}_{0}$, i.e., all information about the overall original system bath coupling strength is captured in the system RC coupling and a system renormalization term. The remaining terms, especially the new coupling coefficients C_k, are independent of the initial coupling strength α.

Now, the crux of the matter is that we do not have to determine Λ directly; instead, the normal mode transformation can be fully fixed by knowledge of the SD ${J}_{0}(\omega )$ only [38]. To see this we first of all note that all relevant quantities of the system and RC itself can be expressed in terms of the original SD as follows:

$$\begin{eqnarray}&&\delta {{\rm{\Omega }}}_{0}^{2}=\displaystyle \frac{2}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \frac{{J}_{0}(\omega )}{\omega },\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\lambda }_{0}^{2}=\displaystyle \frac{2}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \omega {J}_{0}(\omega ),\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{{\rm{\Omega }}}_{1}^{2}=\displaystyle \frac{2}{\pi {\lambda }_{0}^{2}}{\int }_{0}^{\infty }{\rm{d}}\omega {\omega }^{3}{J}_{0}(\omega ).\end{eqnarray} \tag{ 10 }$$

The effect of the residual environment on the system and RC is itself captured by the new SD

$$\begin{eqnarray}&&{J}_{1}(\omega )\equiv \displaystyle \frac{\pi }{2}\displaystyle \sum _{k}\displaystyle \frac{{C}_{k}^{2}}{{{\rm{\Omega }}}_{k}}\delta (\omega -{{\rm{\Omega }}}_{k})\end{eqnarray} \tag{ 11 }$$

and the remaining task is to relate ${J}_{0}(\omega )$ to ${J}_{1}(\omega )$. We will here use the prescription given by Martinazzo et al [53] who have shown the following relation under reasonable mild conditions⁵

$$\begin{eqnarray}&&{J}_{1}(\omega )=\displaystyle \frac{{\lambda }_{0}^{2}{J}_{0}(\omega )}{{| {W}_{0}^{+}(\omega )| }^{2}}.\end{eqnarray} \tag{ 12 }$$

For completeness we will present a derivation of this result in appendix A. Similar methods to link the SDs can be found in [38, 48, 52, 54]. In equation (12), ${W}_{0}(z)$ denotes the Cauchy transform of ${J}_{0}(\omega )$ given by

$$\begin{eqnarray}&&{W}_{0}(z)\equiv \displaystyle \frac{1}{\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \displaystyle \frac{{J}_{0}(\omega )}{\omega -z}\end{eqnarray} \tag{ 13 }$$

(we extended ${J}_{0}(\omega )$ to negative values of ω via ${J}_{0}(-\omega )=-{J}_{0}(\omega ))$ and we introduced the notation

$$\begin{eqnarray}&&{W}_{0}^{+}(\omega )\equiv \underset{\epsilon \searrow 0}{\mathrm{lim}}{W}_{0}(\omega +{\rm{i}}\epsilon )\;\;\;(\omega \in {\mathbb{R}}).\end{eqnarray} \tag{ 14 }$$

Furthermore, one can also show that (see [53] or appendix A again)

$$\begin{eqnarray}&&\displaystyle \frac{{\lambda }_{0}^{2}}{\delta {{\rm{\Omega }}}_{0}^{2}}={{\rm{\Omega }}}_{1}^{2}-\delta {{\rm{\Omega }}}_{1}^{2},\end{eqnarray} \tag{ 15 }$$

where we denoted the frequency renormalization term of the RC by

$$\begin{eqnarray}&&\delta {{\rm{\Omega }}}_{1}^{2}=\displaystyle \sum _{k}\displaystyle \frac{{C}_{k}^{2}}{{{\rm{\Omega }}}_{k}^{2}}=\displaystyle \frac{2}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \frac{{J}_{1}(\omega )}{\omega }.\end{eqnarray} \tag{ 16 }$$

Relation (15) allows us to rewrite the Hamiltonian (7) in a Brownian motion form [49]

$$\begin{eqnarray}&&H^{\prime} ={H}_{{\rm{S}}}(t)+\displaystyle \frac{1}{2}\left[{P}_{1}^{2}+\displaystyle \frac{{\lambda }_{0}^{2}}{\delta {{\rm{\Omega }}}_{0}^{2}}{\left({X}_{1}-\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{{\lambda }_{0}}s\right)}^{2}\right]+\displaystyle \frac{1}{2}\displaystyle \sum _{k}\left[{P}_{k}^{2}+{{\rm{\Omega }}}_{k}^{2}{\left({X}_{k}-\displaystyle \frac{{C}_{k}}{{{\rm{\Omega }}}_{k}^{2}}{X}_{1}\right)}^{2}\right],\end{eqnarray} \tag{ 17 }$$

which makes its thermodynamic stability evident. Hence, the physical frequency of the RC is not given by ${{\rm{\Omega }}}_{1}$ but by the square-root of equation (15).

Finally, we note that the power of the RC mapping also comes from the fact that it can be applied iteratively. This then yields a chain of RCs where the last one is coupled to a residual environment. Remarkably, the relation between the SDs, equation (12), still carries over to this situation (replacing the index 0 by n and the index 1 by $n+1$, where n labels the different RCs) and also all other parameters can be defined in terms of the SD as in equations (8)–(10) [52–54]. Furthermore and very importantly, the fixed point of this iteration scheme is a Markovian SD [53] and the necessary conditions for convergence to a Markovian SD were worked out in [54] and are fulfilled for the situation considered here⁶ . Thus, already at this time we can conclude that the dependence on the initial coupling strength is absorbed by including only one RC (a more critical discussion of this point is shifted to section 6) while strong non-Markovianity might require several RCs.

We now imagine the situation where our system is coupled to several reservoirs labeled by ν and where the time-dependent driving is responsible for work extraction and injection. The obvious generalization of the Hamiltonian (2) to this situation is

$$\begin{eqnarray}&&H={H}_{{\rm{S}}}(t)+\displaystyle \frac{1}{2}\displaystyle \sum _{k,\nu }\left[{p}_{k\nu }^{2}+{\omega }_{k\nu }^{2}{\left({x}_{k\nu }-\displaystyle \frac{{c}_{k\nu }}{{\omega }_{k\nu }^{2}}{s}_{\nu }\right)}^{2}\right],\end{eqnarray} \tag{ 18 }$$

where the coupling to reservoir ν is mediated by the system operator s_ν which might be all different for different ν. A sketch of a possible scenario is shown in figure 2. We can then decide to include zero, one or several RCs for each reservoir depending on the coupling strength and the shape of the SD. Again, all what we need to know for this mapping are the SDs ${J}_{0}^{(\nu )}(\omega )$ for each reservoir ν. The same transformations as introduced in section 2 will carry over in exactly the same way to the situation of multiple reservoirs. It is in particular worth pointing out that each RC mapping is a unitary transformation only on the Hilbert space of bath ν, i.e., it leaves the system part and all other baths fully untouched. This feature allows us to really present a general thermodynamic framework valid for any system, which is coupled to its environment in the prescribed way.

Zoom In Zoom Out Reset image size

After having included a sufficient number of RCs, the next step is to define a new 'supersystem' consisting out of the original system and all RCs. The idea is then to treat this supersystem within the standard Born–Markov secular (BMS) approximation, assuming that the bath of the residual oscillators is in thermal equilibrium, and to derive a Markovian ME for the supersystem. This ME then has a transparent thermodynamic interpretation (as we will review below for reasons of consistency) and this is indeed the strength of our approach: by finding this part of the environment which acts as an ideal, weakly coupled and memory-less thermal bath we are able to provide a formally clean definition of heat, which has a very precise meaning in thermodynamics and does not simply equal the energy flowing into the surroundings in the general (i.e., non-Markovian and strongly coupled) case. Besides this fact, it is worth pointing out here that already including one RC can give remarkable numerical results in agreement with the formally exact hierarchical equations of motion method, as it was recently shown by Iles-Smith et al [48, 51]. Further research in this direction was also conducted in [40–46].

The standard framework of quantum thermodynamics starts with a microscopically derived ME of the form [4–6, 8–10]

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\rho (t)=-{\rm{i}}[{H}_{{\rm{S}}}^{\prime }(t),\rho (t)]+\displaystyle \sum _{\nu }{{ \mathcal L }}^{(\nu )}(t)\rho (t).\end{eqnarray} \tag{ 19 }$$

Here, $\rho (t)$ denotes the density matrix and ${H}_{{\rm{S}}}^{\prime }(t)$ the Hamiltonian of the supersystem, i.e., the system and RCs. The time-dependence of ${H}_{{\rm{S}}}^{\prime }(t)$ might result from an initial time-dependence of H_S(t). For example, in case of a single RC we have (see equation (17))

$$\begin{eqnarray}&&{H}_{{\rm{S}}}^{\prime }(t)={H}_{{\rm{S}}}(t)+\displaystyle \frac{1}{2}\left[{P}_{1}^{2}+\displaystyle \frac{{\lambda }_{0}^{2}}{\delta {{\rm{\Omega }}}_{0}^{2}}{\left({X}_{1}-\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{{\lambda }_{0}}s\right)}^{2}\right].\end{eqnarray} \tag{ 20 }$$

Furthermore, the dissipators or thermal generators ${{ \mathcal L }}^{(\nu )}(t)$ are of Lindblad form [61] and fulfill local detailed balance, i.e.

$$\begin{eqnarray}&&{{ \mathcal L }}^{(\nu )}(t){{\rm{e}}}^{-{\beta }_{\nu }{H}_{{\rm{S}}}^{\prime }(t)}=0,\end{eqnarray} \tag{ 21 }$$

where ${\beta }_{\nu }$ is the temperature of reservoir ν and which tells us that the ratio of backward to forward transition rates is given by a Boltzmann factor for the case of a Pauli rate ME for a non-degenerate supersystem Hamiltonian⁷ [4]. We remark that, strictly speaking, a ME of the form (19) can only be derived for a slow time-dependence of H_S(t). However, using techniques from Floquet theory it is also possible to derive a ME for (strong) periodic driving [61] with a similar thermodynamic interpretation [8, 10, 62]. For an arbitrary driving H_S(t) there is no guarantee to find a simple ME for the system, but a thermodynamic analysis can be still carried out [63] (see also the general treatment [12]).

We now define the internal energy and entropy of the supersystem via

$$\begin{eqnarray}&&E(t)\equiv {\rm{tr}}\{{H}_{{\rm{S}}}^{\prime }(t)\rho (t)\},\;\;\;S(t)\equiv -{\rm{tr}}\{\rho (t)\mathrm{ln}\rho (t)\}.\end{eqnarray} \tag{ 22 }$$

The first law of thermodynamics then acquires the form

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}E(t)=\dot{W}(t)+\displaystyle \sum _{\nu }{\dot{Q}}^{(\nu )}(t),\end{eqnarray} \tag{ 23 }$$

where we identified the rate of work (power) done on the supersystem and the heat flow coming from reservoir ν as⁸

$$\begin{eqnarray}&&\dot{W}(t)\equiv {\rm{tr}}\left\{\rho (t)\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}{H}_{{\rm{S}}}(t)\right\},\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{\dot{Q}}^{(\nu )}(t)\equiv {\rm{tr}}\{{H}_{{\rm{S}}}^{\prime }(t){{ \mathcal L }}^{(\nu )}(t)\rho (t)\}.\end{eqnarray} \tag{ 25 }$$

In equation (24) we used that $\tfrac{{\rm{d}}}{{\rm{d}}t}{H}_{{\rm{S}}}^{\prime }(t)=\tfrac{{\rm{d}}}{{\rm{d}}t}{H}_{{\rm{S}}}(t)$. The second law stipulates that the rate of entropy production ${\dot{S}}_{{\bf{i}}}(t)$ is always positive

$$\begin{eqnarray}&&{\dot{S}}_{{\bf{i}}}(t)=\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}S(t)-\displaystyle \sum _{\nu }{\beta }_{\nu }{\dot{Q}}^{(\nu )}(t)\geqslant 0.\end{eqnarray} \tag{ 26 }$$

It is possible to prove equation (26) by use of Spohn's inequality stating that [3]

$$\begin{eqnarray}&&-{\rm{tr}}\{[{{ \mathcal L }}^{(\nu )}(t)\rho (t)][\mathrm{ln}\rho (t)-\mathrm{ln}{\rho }_{\mathrm{eq}}^{(\nu )}]\}\geqslant 0\end{eqnarray} \tag{ 27 }$$

for every ν and ${\rho }_{\mathrm{eq}}^{(\nu )}\equiv {{\rm{e}}}^{-{\beta }_{\nu }{H}_{{\rm{S}}}^{\prime }(t)}/{Z}_{\nu }$. By summing equation (27) over ν we obtain the second law. At this point we wish to emphasize that within our theory the reservoirs ν enter additively (or separately) in the first and second law of thermodynamics, which is reminiscent of the fact that the RC mapping can be applied to each bath separately as indicated in figure 2. Especially, the temperatures (and chemical potentials) of the residual baths are still the same and well-defined.

Finally, if the system is undriven it will eventually reach a steady state and the first and second law become

$$\begin{eqnarray}&&0=\dot{W}+\displaystyle \sum _{\nu }{\dot{Q}}^{(\nu )},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&{\dot{S}}_{{\bf{i}}}=-\displaystyle \sum _{\nu }{\beta }_{\nu }{\dot{Q}}^{(\nu )}\geqslant 0.\end{eqnarray} \tag{ 29 }$$

To indicate that we are at steady state, we dropped the time dependence on all quantities and for simplicity we will exclusively focus on the steady state regime for the rest of this paper. We also note that even for an undriven system there might be still a work source present (i.e., $\dot{W}\ne 0$) by identifying a work reservoir appropriately, see section 4, or due to the possibility of particle transport ('chemical work'), see section 5.

Before we proceed to illustrate our theory with examples of heat engines working out of equilibrium, it might be worth to stress a simple consequence of our treatment at equilibrium. If the supersystem is time independent and in contact with only one reservoir at inverse temperature β, it will relax to an equilibrium state $\rho (t\to \infty )\sim {{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}$ such that the equilibrium state of the original system is [48]

$$\begin{eqnarray}&&{\rho }_{{\rm{S}}}(t\to \infty )\sim {{\rm{tr}}}_{\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}\}.\end{eqnarray} \tag{ 30 }$$

In appendix B we will demonstrate that this state is consistent to lowest order perturbation theory in the coupling to the residual bath with the conventionally used Hamiltonian of mean force as introduced by Kirkwood [64]. In particular, this state in general does not equal the canonical equilibrium state of the system alone, i.e., ${\rho }_{{\rm{S}}}(t\to \infty )\;\not\hspace{-2pt}{\sim }\;{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}}$. Experimentally, deviations from the canonical state ${{\rm{e}}}^{-\beta {H}_{{\rm{S}}}}$ might be a clear indicator for persistent system environment correlations making it necessary to go beyond the Born approximation, e.g., by using the RC mapping. Given the SD of the reservoir, it should be then possible to test the prediction (30) in an actual experiment.

4.1. Standard treatment without RC

Possibly one of the simplest heat engine one can think of consists of three time-independent levels described by the Hamiltonian

$$\begin{eqnarray}&&{H}_{{\rm{S}}}=\displaystyle \sum _{i=0}^{2}{\epsilon }_{i}| i\rangle \langle i| \end{eqnarray} \tag{ 31 }$$

with ${\epsilon }_{0}\lt {\epsilon }_{1}\lt {\epsilon }_{2}$. The idea behind this engine is the model of a simple maser which is lasing at a particular transition, say $0\leftrightarrow 1$. This lasing corresponds to 'work' output and it is achieved due to population inversion between the levels $| 0\rangle$ and $| 1\rangle$. This in turn can be mediated via a third level $| 2\rangle$ due to the presence of two heat reservoirs at different temperatures (called the 'hot' and 'cold' reservoir respectively), also see figure 3 for a sketch. Initially, this model was investigated in 1959 [65], but it is still of interest today [66–70].

Zoom In Zoom Out Reset image size

The coupling to the reservoirs is mediated by the system operators

$$\begin{eqnarray}&&{s}_{h}=\displaystyle \frac{| 0\rangle \langle 2| +| 2\rangle \langle 0| }{\sqrt{2\epsilon }},\;\;\;{s}_{c}=\displaystyle \frac{| 1\rangle \langle 2| +| 2\rangle \langle 1| }{\sqrt{2\epsilon }},\;\;\;{s}_{w}=\displaystyle \frac{| 0\rangle \langle 1| +| 1\rangle \langle 0| }{\sqrt{2\epsilon }}.\end{eqnarray} \tag{ 32 }$$

Here, in units where ${\rm{\hslash }}=1$, the parameter $\epsilon$ has units of an energy or frequency such that s_ν has the same units as the coordinates ${x}_{k,\nu }$ of the bath oscillators. However, in all numerical calculations which follow we will simply set $\epsilon \equiv 1$. Within the standard approach (BMS approximation for the system only) the thermal generators ${{ \mathcal L }}^{(\nu )}$ in equation (19) of each bath become

$$\begin{eqnarray*}{{ \mathcal L }}^{(h)} & = & \displaystyle \frac{{J}_{0}^{(h)}({{\rm{\Delta }}}_{20})}{\epsilon }\{[1+{n}_{h}({{\rm{\Delta }}}_{20})]{{ \mathcal D }}_{02}+{n}_{h}({{\rm{\Delta }}}_{20}){{ \mathcal D }}_{20}\},\\ {{ \mathcal L }}^{(c)} & = & \displaystyle \frac{{J}_{0}^{(c)}({{\rm{\Delta }}}_{21})}{\epsilon }\{[1+{n}_{c}({{\rm{\Delta }}}_{21})]{{ \mathcal D }}_{12}+{n}_{c}({{\rm{\Delta }}}_{21}){{ \mathcal D }}_{21}\},\\ {{ \mathcal L }}^{(w)} & = & \displaystyle \frac{{J}_{0}^{(w)}({{\rm{\Delta }}}_{10})}{\epsilon }\{[1+{n}_{w}({{\rm{\Delta }}}_{10})]{{ \mathcal D }}_{01}+{n}_{w}({{\rm{\Delta }}}_{10}){{ \mathcal D }}_{10}\}.\end{eqnarray*}$$

Here, we have introduced the dissipator ${{ \mathcal D }}_{{ij}}\rho \equiv | i\rangle \langle j| \rho | j\rangle \langle i| -\tfrac{1}{2}\{| j\rangle \langle j| ,\rho \}$, the Bose distribution ${n}_{\nu }({{\rm{\Delta }}}_{{ij}})\equiv {({{\rm{e}}}^{{\beta }_{\nu }{{\rm{\Delta }}}_{{ij}}}-1)}^{-1}$ and the SDs ${J}_{0}^{(\nu )}(\omega )$ of bath ν are evaluated at the transition frequency ${{\rm{\Delta }}}_{{ij}}\equiv {\epsilon }_{i}-{\epsilon }_{j}$.

Given the prescription of section 3, it is not hard to compute the thermodynamic behaviour of our system. At steady state the first law becomes $0=\dot{W}+{\dot{Q}}^{h}+{\dot{Q}}^{c}$ with $\dot{W}\equiv {\dot{Q}}^{w}$ while the second law states that $-{\beta }_{w}\dot{W}-{\beta }_{h}{\dot{Q}}^{h}-{\beta }_{c}{\dot{Q}}^{c}\geqslant 0$. To quantify the performance of work extraction (i.e., $\dot{W}\lt 0$), we introduce the efficiency of the heat engine:

$$\begin{eqnarray}&&\eta \equiv \displaystyle \frac{-\dot{W}}{{\dot{Q}}^{h}}=\displaystyle \frac{{{\rm{\Delta }}}_{10}}{{{\rm{\Delta }}}_{20}}\leqslant \displaystyle \frac{{\beta }_{c}-{\beta }_{h}}{{\beta }_{c}-{\beta }_{w}},\end{eqnarray} \tag{ 33 }$$

where the inequality is a consequence of the second law. Furthermore and very remarkably, it is possible to show that the efficiency of the maser heat engine is always given by the ratio ${{\rm{\Delta }}}_{10}/{{\rm{\Delta }}}_{20}$ independent of all other parameters (as long as $\dot{W}\lt 0$, of course).

Note that to completely justify the notion of a 'work reservoir', we should take the limit ${\beta }_{w}\to 0$ [69, 70], which complies with Sommerfelds notion of temperature as the 'work value of heat' [71]. In this limit the entropy change $-{\beta }_{w}\dot{W}$ in the work reservoir w goes to zero, the second law acquires the form $-{\beta }_{h}{\dot{Q}}^{h}-{\beta }_{c}{\dot{Q}}^{c}\geqslant 0$ and the efficiency (33) is bounded by Carnot efficiency ${\eta }_{\mathrm{Carnot}}=1-\tfrac{{T}_{c}}{{T}_{h}}$. Microscopically, however, we recognize that the Bose distribution ${n}_{w}({{\rm{\Delta }}}_{10})$ diverges for ${\beta }_{w}\to 0$ (though $\dot{W}$ remains finite), unless we additionally require that the bath SD scales like ${J}_{0}^{(w)}({{\rm{\Delta }}}_{10})={\beta }_{w}{{\rm{\Delta }}}_{10}{{\rm{\Gamma }}}_{w}$ for small ${\beta }_{w}$. In this ideal limit we then obtain

$$\begin{eqnarray}&&\underset{{\beta }_{w}\to 0}{\mathrm{lim}}{{ \mathcal L }}^{(w)}={{\rm{\Gamma }}}_{w}({{ \mathcal D }}_{01}+{{ \mathcal D }}_{10}),\end{eqnarray} \tag{ 34 }$$

i.e., the rates of upward and downward transitions are equal. However, for all numerical results reported below we take ${\beta }_{w}$ to be small but non-zero.

Now, our approach is to consider the situation where it is actually not valid to apply this BMS ME for the system only. For instance, this could be due to a structured (non-Markovian) SD. In this case, one should actually use a non-Markovian ME for the system (e.g., the Redfield equation [61]). However, the thermodynamic interpretation of non-Markovian MEs is not clear and has not been established yet. In contrast, within our approach we know that we can include a RC into our description in order to account for non-Markovian effects on the three-level system (3LS) while the 3LS and RC evolve in a Markovian way. Any deviations from the efficiency ${{\rm{\Delta }}}_{10}/{{\rm{\Delta }}}_{20}$ then indicate non-Markovian and/or strong coupling effects which we would be unable to detect by using the naive ME approach outlined in this section.

4.2. Thermodynamics with RC

For definiteness we choose the SD of the cold bath to be parameterized as

$$\begin{eqnarray}&&{J}_{0}^{(c)}(\omega )=\displaystyle \frac{{d}_{0}^{2}\gamma \omega }{{({\omega }^{2}-{\omega }_{0}^{2})}^{2}+{\gamma }^{2}{\omega }^{2}}{\rm{\Theta }}({\omega }_{{\rm{R}}}-\omega ),\end{eqnarray} \tag{ 35 }$$

while we still assume that it is safe to apply the Markov approximation with respect to the interaction with the hot and work reservoir. Thus, by following the prescription of section 2 we obtain the modified system Hamiltonian

$$\begin{eqnarray}&&{H}_{{\rm{S}}}^{\prime }=\displaystyle \sum _{i=0}^{2}{\epsilon }_{i}| i\rangle \langle i| +\displaystyle \frac{1}{2}\left[{P}_{1}^{2}+\displaystyle \frac{{\lambda }_{0}^{2}}{\delta {{\rm{\Omega }}}_{0}^{2}}{\left({X}_{1}-\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{{\lambda }_{0}}{s}_{c}\right)}^{2}\right],\end{eqnarray} \tag{ 36 }$$

and the SD ${J}_{1}^{(c)}(\omega )$ of the residual cold bath is determined by equation (12). In figure 4 we plot ${J}_{0}^{(c)}(\omega )$ and ${J}_{1}^{(c)}(\omega )$ for comparison.

Zoom In Zoom Out Reset image size

At this point it is noteworthy that similar models have been also used in quantum biology to model light-harvesting complexes [72, 73]. Indeed, guided by this motivation, a number of researchers have started to investigate light-harvesting complexes from a heat engine perspective [74, 75] though we wish to stress that the models used there are not exactly the same as the one used here. Furthermore, [74] models the work reservoir effectively as a zero temperature reservoir, which causes divergences in the standard thermodynamic formalism. In [75] the polaron transformation was used in order to access strong coupling effects. We will introduce this transformation in section 5 for a different model to compare it with our RC method.

To proceed we first of all note that the system Hamiltonian can be alternatively written as

$$\begin{eqnarray}{H}_{{\rm{S}}}^{\prime } & = & \left({{\rm{\Delta }}}_{10}+\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{4\epsilon }\right)| 1\rangle \langle 1| +\left({{\rm{\Delta }}}_{20}+\displaystyle \frac{\delta {{\rm{\Omega }}}_{0}^{2}}{4\epsilon }\right)| 2\rangle \langle 2| \\ & & -\displaystyle \frac{{\lambda }_{0}}{\sqrt{2\epsilon }}{X}_{1}(| 1\rangle \langle 2| +| 2\rangle \langle 1| )+\displaystyle \frac{1}{2}\left({P}_{1}^{2}+\displaystyle \frac{{\lambda }_{0}^{2}}{\delta {{\rm{\Omega }}}_{0}^{2}}{X}_{1}^{2}\right),\end{eqnarray} \tag{ 37 }$$

where we made use of the fact that we can choose the ground state energy of the 3LS arbitrarily and set it to ${\epsilon }_{0}\equiv 0$. This Hamiltonian describes a Rabi model (harmonic oscillator coupled to a two-level system) plus one additional energy level $| 0\rangle$. Especially note that the energy levels of state $| 1\rangle$ and $| 2\rangle$ get both shifted by the same amount $\delta {{\rm{\Omega }}}_{0}^{2}/4\epsilon$. In the weak coupling (but non-Markovian) regime—in which we wish to compare our extended model with the one treated before—these terms become negligible small.

To investigate the thermodynamic behaviour of our system, we will now use a ME which is explicit concerning the system-RC interaction, but treats the coupling to the other reservoirs perturbatively and in a Markovian way. Following standard procedures [6, 8, 9, 61] it is then possible to derive the BMS ME of the form (19) with time-independent Hamiltonian ${H}_{{\rm{S}}}^{\prime }$ and dissipators ${{ \mathcal L }}^{(\nu )}$. Thus, the thermodynamic treatment follows straightforwardly from section 3 and is formally equivalent to the model of section 4.1. Hence, the first and second law become at steady state

$$\begin{eqnarray}&&0={\dot{\tilde{Q}}}^{h}+{\dot{\tilde{Q}}}^{c}+\dot{\tilde{W}},\end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray}&&0\leqslant -{\beta }_{h}{\dot{\tilde{Q}}}^{h}-{\beta }_{c}{\dot{\tilde{Q}}}^{c}-{\beta }_{w}\dot{\tilde{W}}.\end{eqnarray} \tag{ 39 }$$

Now, however, we used a tilde on all energy flows because they are numerically different from the corresponding quantities in section 4.1. Likewise we introduce the efficiency of our engine for $\dot{\tilde{W}}\lt 0$ as

$$\begin{eqnarray}&&\tilde{\eta }=\displaystyle \frac{-\dot{\tilde{W}}}{{\dot{\tilde{Q}}}^{h}},\end{eqnarray} \tag{ 40 }$$

which is also bounded by the Carnot efficiency in the limit ${\beta }_{w}\to 0$.

Because it is not possible to diagonalize the Hamiltonian (37) in a simple manner, one has to treat the ME numerically, and technical details of the derivation are presented in the appendix C. We also note that we decided to compare the ME based only on the Born and Markov approximation (i.e., without secular approximation) with the BMS ME. The latter is usually only well justified for systems where the level spacing is much larger than the level broadening. This becomes increasingly questionable for more complex systems with many different (and not always well separated) eigen frequencies. On the other hand, the advantage of the BMS ME is that the resulting generator is of Lindblad form and allows for a mathematically clear proof of Spohn's inequality (27) [3, 8, 61]. In the numerical results reported below, however, we never found any violation of the second law even when we used the ME without secular approximation.

Figures 5– 7 show numerical results obtained from the model without RC (dashed lines, see section 4.1) and with RC (with (solid lines) and without (dots) secular approximation) for a specific choice of parameters. This choice was done to illustrate—from our point of view—interesting features of Markovian and non-Markovian heat engines. However, a detailed investigation of the model is beyond the scope of the paper and would also be questionable because the model from section 4.1 clearly is a simplified toy model. Nevertheless, we wish to remark that the qualitative behaviour of our numerical results remains the same for a wide range of parameters. Furthermore, we have focused on plots of the efficiency and work, which are both of great interest in thermodynamics: obviously, we want to have a large power output, but on the other hand, a large power output does not mean that our engine works more efficiently. In fact, if we possess a heat engine with low power output but high efficiency, we can also build a machine with high power output and the same efficiency by running several of these machines in parallel. Hence, we think efficiency is a more universal quantity on which one should put more emphasis in the study of heat engines in the strong coupling and non-Markovian regime.

Zoom In Zoom Out Reset image size

Turning to the details, we can infer from figure 5 two important points. First, our approach predicts an efficiency enhancement of 10%–20% in comparison to the standard theoretical framework from section 4.1. We here note that it makes sense to compare the model without and with RC because all parameters of the latter are completely fixed by the initial model. We simply use two different theoretical methods to study the same engine, but only the latter allows to capture, e.g., non-Markovian effects. The second point to note concerns the secular approximation. Whereas for a large parameter regime it seems to be remarkably well justified, it can also completely fail by predicting a nearly hundred times larger power output for the case ${{\rm{\Delta }}}_{21}\approx {\omega }_{0}$ (interestingly, this effect cancels when we compute the efficiency). In fact, if this resonance condition is met⁹ , many energy levels are very close to each other (though never exactly degenerate) and a naive application of the secular approximation is not justified.

In figure 6 we show the efficiency, work and heat flow as a function of γ. As we can infer from equation (35) (also see figure 4) a smaller γ implies a stronger peaked SD ${J}_{0}^{(c)}(\omega )$. Thus, γ can be seen as a measure of the non-Markovianity of the SD and figure 6 proves that this is the cause of the efficiency enhancement. This justifies our claim that non-Markovian machines can significantly outperform their Markovian counterparts even at steady state¹⁰ . Physically, the reason for this can be traced back to the Purcell effect which predicts an enhanced spontaneous emission rate for the $1\leftrightarrow 2$ transition in presence of a (resonant) cavity [76] and thus allows for a stronger population inversion between the states $| 0\rangle$ and $| 1\rangle$. By including the RC we can indeed capture this effect which is clearly beyond the 'naive' Markovian treatment of section 4.1. Furthermore, figure 6 also shows that we recover the results from section 4.1 in the limit of large γ. In fact, for large γ the SD does not only become more Markovian, but also the SD ${J}_{1}^{(c)}(\omega )$ of the residual cold bath is directly proportional to γ, see equation (C.14), such that the RC evolves on time-scales much shorter then the 3LS and can be adiabatically eliminated. Furthermore, figure 6 also shows that the BMS ME of the supersystem completely fails in this regime. This behaviour can be traced back to the fact that the secular approximation does not commute with the adiabatic elimination of the RC, i.e., the time-scales involved in the coherent evolution of the system are of the same order as the dynamics of the relaxation due to the residual cold bath.

Zoom In Zoom Out Reset image size

Finally, figure 7 shows the thermodynamics as a function of the system-RC coupling strength d₀ (or the coupling strength between the 3LS and the cold reservoir, respectively). Again, we can observe a strong efficiency enhancement and an almost perfect agreement between the secular and non-secular ME. Furthermore, we observe that the efficiency decreases as a function of d₀ while the power output first increases up to a certain critical coupling strength and then starts to decrease again (note that the power output for the 3LS ME from section 4.1 reaches a constant value for ${d}_{0}\to \infty$ instead). In fact, for all parameters we have checked we always observed a decreasing efficiency as a function of d₀. Whether this is a general feature or model specific remains an open question, which might be eventually answerable within the RC framework.

Zoom In Zoom Out Reset image size

5.1. Model and RC mapping

As a second application we consider a single quantum dot in contact with two fermionic reservoirs (typically called a single electron transistor) and additionally coupled to a bath of phonons. Related models have been studied, e.g., in [33, 34, 77–80], in order to understand electronic transport through molecules where the phonon bath models molecular vibrations. Here, as well as outlining another model applicable to the RC method, we also wish to compare with the polaron transformation, a technique frequently used to access the regime of strong system–phonon coupling [31, 33–36, 79, 80]. We will discuss what we mean by that in more detail at the end of section 5.2.

The global Hamiltonian can be written as

$$\begin{eqnarray}&&H={H}_{\mathrm{dot}}+{H}_{\mathrm{el}}+{H}_{\mathrm{el}\mbox{--}\mathrm{dot}}+{H}_{\mathrm{ph}}.\end{eqnarray} \tag{ 41 }$$

The electronic part of the system is described by

$$\begin{eqnarray}&&{H}_{\mathrm{dot}}=\epsilon {d}^{\dagger }d,\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{el}}=\displaystyle \sum _{\nu }\displaystyle \sum _{k}{\epsilon }_{k\nu }{c}_{k\nu }^{\dagger }{c}_{k\nu },\end{eqnarray} \tag{ 43 }$$

$$\begin{eqnarray}&&{H}_{\mathrm{el}\mbox{--}\mathrm{dot}}=\displaystyle \sum _{\nu }\displaystyle \sum _{k}({t}_{k\nu }{c}_{k\nu }{d}^{\dagger }+{t}_{k\nu }^{* }{{dc}}_{k\nu }^{\dagger }),\end{eqnarray} \tag{ 44 }$$

where ${d}^{(\dagger )}$ and ${c}_{k,\nu }^{(\dagger )}$ are fermionic annihilation (creation) operators of the dot and the reservoir with associated on-site energy and lead energy ${\epsilon }_{k\nu }$. Furthermore, the dot is connected to two leads $\nu \in \{L,R\}$ with tunneling amplitudes ${t}_{k\nu }$. The interaction of the dot with the phonon bath is assumed to be

$$\begin{eqnarray}&&{H}_{\mathrm{ph}}=\displaystyle \frac{1}{2}\displaystyle \sum _{q}\left[{p}_{q}^{2}+{\omega }_{q}^{2}{\left({x}_{q}-\displaystyle \frac{{h}_{q}}{{\omega }_{q}^{2}}{d}^{\dagger }d\right)}^{2}\right].\end{eqnarray} \tag{ 45 }$$

Here, we denoted the coupling coefficients to the phonon bath by h_q because we already use ${c}_{k\nu }$ for the electrons in the fermionic reservoirs.

Next, again in order to overcome the limitations of the usual BMS ME we employ the RC mapping, this time to the phonon bath. This transforms the system and phonon part to (see equation (7))

$$\begin{eqnarray*}{H}_{\mathrm{dot}}+{H}_{\mathrm{ph}} & = & \left(\epsilon +\displaystyle \frac{1}{2}\displaystyle \sum _{q}\displaystyle \frac{{h}_{q}^{2}}{{\omega }_{q}^{2}}\right){d}^{\dagger }d-{\lambda }_{0}{d}^{\dagger }{{dX}}_{1}+\displaystyle \frac{1}{2}({P}_{1}^{2}+{{\rm{\Omega }}}_{1}^{2}{X}_{1}^{2})\\ & & -{X}_{1}\displaystyle \sum _{q}{C}_{q}{X}_{q}+\displaystyle \frac{1}{2}\displaystyle \sum _{q}({P}_{q}^{2}+{{\rm{\Omega }}}_{q}^{2}{X}_{q}^{2}).\end{eqnarray*}$$

We now identify our supersystem to be

$$\begin{eqnarray}&&{H}_{{\rm{S}}}^{\prime }=\tilde{\epsilon }{d}^{\dagger }d-{\lambda }_{0}{d}^{\dagger }{{dX}}_{1}+\displaystyle \frac{1}{2}({P}_{1}^{2}+{{\rm{\Omega }}}_{1}^{2}{X}_{1}^{2})\end{eqnarray} \tag{ 46 }$$

with $\tilde{\epsilon }\equiv \epsilon +\tfrac{1}{2}{\sum }_{q}{h}_{q}^{2}/{\omega }_{q}^{2}$ whereas the interaction with the phonon bath is described via ${H}_{{\rm{I}}}^{\mathrm{ph}}=-{X}_{1}{\sum }_{q}{C}_{q}{X}_{q}$ and the residual phonon bath Hamiltonian becomes ${H}_{{\rm{B}}}^{\mathrm{ph}}=\tfrac{1}{2}{\sum }_{q}({P}_{q}^{2}+{{\rm{\Omega }}}_{q}^{2}{X}_{q}^{2})$. Note that this does not exactly correspond to the mapping (17), but now the Hamiltonian is closer to the literature to which we wish to compare our results [33]. Furthermore, the spectrum of ${H}_{{\rm{S}}}^{\prime }$ is still bounded from below for all coupling strengths ${\lambda }_{0}$.

Below it will be convenient to work with the bosonic ladder operators, which are related to position and momentum operators via

$$\begin{eqnarray}&&{X}_{q}=\displaystyle \frac{1}{\sqrt{2{{\rm{\Omega }}}_{q}}}({a}_{q}+{a}_{q}^{\dagger }),\;\;\;{P}_{q}={\rm{i}}\sqrt{\displaystyle \frac{{{\rm{\Omega }}}_{q}}{2}}({a}_{q}-{a}_{q}^{\dagger }).\end{eqnarray} \tag{ 47 }$$

In terms of these operators we have

$$\begin{eqnarray}&&{H}_{{\rm{S}}}^{\prime }=\tilde{\epsilon }{d}^{\dagger }d+{{\rm{\Omega }}}_{1}\left({a}_{1}^{\dagger }{a}_{1}+\displaystyle \frac{1}{2}\right)-\displaystyle \frac{{\lambda }_{0}{d}^{\dagger }d}{\sqrt{2{{\rm{\Omega }}}_{1}}}({a}_{1}+{a}_{1}^{\dagger }),\end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray}&&{H}_{{\rm{I}}}^{\mathrm{ph}}=-({a}_{1}+{a}_{1}^{\dagger })\displaystyle \sum _{q}\displaystyle \frac{{C}_{q}}{2\sqrt{{{\rm{\Omega }}}_{1}{{\rm{\Omega }}}_{q}}}({a}_{q}+{a}_{q}^{\dagger }),\end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray}&&{H}_{{\rm{B}}}^{\mathrm{ph}}=\displaystyle \sum _{q}{{\rm{\Omega }}}_{q}\left({a}_{q}^{\dagger }{a}_{q}+\displaystyle \frac{1}{2}\right).\end{eqnarray} \tag{ 50 }$$

5.2. ME description

To deduce the BMS ME we need to diagonalize the system Hamiltonian (48). This can be accomplished via the unitary transformation (we introduce $\lambda \equiv {\lambda }_{0}/\sqrt{2{{\rm{\Omega }}}_{1}}$ for brevity)

$$\begin{eqnarray}&&U=\mathrm{exp}\left[\displaystyle \frac{\lambda }{{{\rm{\Omega }}}_{1}}{d}^{\dagger }d({a}_{1}-{a}_{1}^{\dagger })\right].\end{eqnarray} \tag{ 51 }$$

It transforms operators according to ${{Ua}}_{1}{U}^{\dagger }={a}_{1}+\tfrac{\lambda }{{{\rm{\Omega }}}_{1}}{d}^{\dagger }d$ and ${{UdU}}^{\dagger }=d{{\rm{e}}}^{-\tfrac{\lambda }{{{\rm{\Omega }}}_{1}}({a}_{1}-{a}_{1}^{\dagger })}$. Hence, applying (51) to ${H}_{{\rm{S}}}^{\prime }$ yields

$$\begin{eqnarray}&&{{UH}}_{{\rm{S}}}^{\prime }{U}^{\dagger }=\bar{\epsilon }{d}^{\dagger }d+{{\rm{\Omega }}}_{1}\left({a}_{1}^{\dagger }{a}_{1}+\displaystyle \frac{1}{2}\right)\end{eqnarray} \tag{ 52 }$$

with $\bar{\epsilon }\equiv \tilde{\epsilon }-{\lambda }^{2}/{{\rm{\Omega }}}_{1}$. This Hamiltonian is obviously in diagonal form with eigenstates $| n,m\rangle$, where $n\in \{0,1\}$ ($m\in \{0,1,2,...\}$) quantifies the electronic (bosonic) occupation, and eigenenergies ${E}_{{nm}}=\bar{\epsilon }n+{{\rm{\Omega }}}_{1}\left(m+\tfrac{1}{2}\right)$.

Given the spectral decomposition of the Hamiltonian it is simple to deduce the BMS ME, especially if the spectrum is non-degenerate, which will generically be the case (unless $\bar{\epsilon }$ coincides with ${{\rm{\Omega }}}_{1}$). To completely specify our model we choose to define and parametrize the SD of the fermionic leads as

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{\nu }(\omega )\equiv 2\pi \displaystyle \sum _{k}{| {t}_{k\nu }| }^{2}\delta (\omega -{\epsilon }_{k\nu })=\displaystyle \frac{{{\rm{\Gamma }}}_{\nu }{\delta }_{\nu }^{2}}{{(\omega -{\varepsilon }_{\nu })}^{2}+{\delta }_{\nu }^{2}}.\end{eqnarray} \tag{ 53 }$$

Here, ${\delta }_{\nu }$ quantifies the width and ${\varepsilon }_{\nu }$ the maximum of the Lorentzian function. Furthermore, the SD of the residual phonon bath is choosen to be ohmic with exponential cutoff

$$\begin{eqnarray}&&{J}_{1}(\omega )={J}_{\mathrm{ph}}\omega {{\rm{e}}}^{-| \omega | /{\omega }_{{\rm{R}}}},\end{eqnarray} \tag{ 54 }$$

which for large cutoff frequencies ${\omega }_{{\rm{R}}}$ corresponds to an initial SD of the form (35) (as discussed at the end of appendix C). Further technical details concerning the derivation of the BMS ME are provided in appendix D. A sketch of the resulting dynamics is provided in figure 8.

Zoom In Zoom Out Reset image size

Before we proceed to show numerical results, let us briefly explain an alternative way to treat strong system-phonon interaction and to which we refer as the polaron ME (PME). Similarly to the RC mapping, also the PME starts with a unitary transformation, which acts on the bath and system Hilbert space though. This transformation usually has the form of a generalized displacement operation such as (51). Then, within this displaced reference frame it is possible to treat the coupling as a small perturbation again because the strong part is absorbed in this new frame. Consequently, it is possible to derive a ME valid for strong coupling for the system part only (in our case here the quantum dot). In contrast, for our treatment we explicitly include the RC as a part of the (super) system and use the transformation (51) later on only to formally diagonalize the supersystem.

5.3. Results

In figure 9 we plot the matter current from left to right versus the difference in chemical potentials $V={\mu }_{{\rm{L}}}-{\mu }_{{\rm{R}}}$ of the electronic reservoirs. First, we see that at low electronic temperatures, the current displays multiple steps, which occur at ${V}_{i}=2{\rm{\Delta }}{E}_{i}$, where ${\rm{\Delta }}{E}_{i}$ are the transition frequencies of the system. Thus, the steps enable one to deduce the renormalized electron energy $\bar{\epsilon }$ and the phonon frequency ${{\rm{\Omega }}}_{1}$ from the electronic current. Furthermore, if we truncate the phonon Hilbert space at a small cutoff number ${N}_{{\rm{c}}}$, we see that only few plateaus are visible since the number of possible transitions is bounded. Most notably, however, when both ${N}_{{\rm{c}}}$ is large enough and the coupling between the RC and the residual oscillators is large, we reproduce earlier results based on a PME [33] (solid curves). In fact, for a large coupling between the RC and its environment, the RC will thermalize on much shorter time-scales as compared to the dot-lead evolution. This is exactly the regime of the PME in which it is assumed ad hoc that the environment in the polaron frame is equilibrated with respect to the dot state. We thus see that the RC method gives us a way to physically and mathematically justify the PME, but in general will be also applicable beyond that regime.

Zoom In Zoom Out Reset image size

Finally, we want to turn to the question whether it is possible to use our setup as a thermoelectric device, i.e., whether we can pump electrons against the bias due to a temperature gradient between electronic leads and phonon reservoir. From figure 9 we see that we have zero current at zero bias even for different temperatures of the electronic and phononic reservoirs (blue and turquoise). This is due to the fact that we chose ${\delta }_{{\rm{L}}}={\delta }_{{\rm{R}}}$ and ${\varepsilon }_{{\rm{L}}}={\varepsilon }_{{\rm{R}}}$ (see equation (53)), which makes our setup symmetric under exchanging the labels L and R at zero bias. This makes thermoelectric transport impossible and will be changed now.

To quantify the irreversibility of our device we take a look at the entropy production (29)

$$\begin{eqnarray}{\dot{S}}_{{\bf{i}}} & = & -{\beta }_{{\rm{L}}}{\dot{Q}}^{{\rm{L}}}-{\beta }_{{\rm{R}}}{\dot{Q}}^{{\rm{R}}}-{\beta }_{\mathrm{ph}}{\dot{Q}}^{\mathrm{ph}}\\ & = & {\beta }_{\mathrm{el}}({\mu }_{{\rm{L}}}-{\mu }_{{\rm{R}}}){I}_{{\rm{M}}}+({\beta }_{\mathrm{el}}-{\beta }_{\mathrm{ph}}){I}_{{\rm{E}}}\geqslant 0.\end{eqnarray} \tag{ 55 }$$

For the last line we assumed equal temperatures in the electronic leads ${\beta }_{{\rm{L}}}={\beta }_{{\rm{R}}}={\beta }_{{\rm{el}}}$ and we used matter conservation ${I}_{{\rm{M}}}\equiv +{I}_{{\rm{M}}}^{{\rm{L}}}=-{I}_{{\rm{M}}}^{{\rm{R}}}$ and energy conservation ${I}_{{\rm{E}}}\equiv {I}_{{\rm{E}}}^{{\rm{ph}}}=-({I}_{{\rm{E}}}^{{\rm{L}}}+{I}_{{\rm{E}}}^{{\rm{R}}})$.

When the electronic and phononic temperatures are different (e.g. ${\beta }_{{\rm{el}}}\gt {\beta }_{{\rm{ph}}}$), heat will flow between the electronic and phononic reservoirs, and the device can use a fraction of the heat to produce positive power $P=-({\mu }_{{\rm{L}}}-{\mu }_{{\rm{R}}}){I}_{{\rm{M}}}$ by transporting charges against the bias. The positivity of the entropy production leads to an upper bound for the efficiency of heat-to-power conversion

$$\begin{eqnarray}&&\eta =\displaystyle \frac{-({\mu }_{{\rm{L}}}-{\mu }_{{\rm{R}}}){I}_{{\rm{M}}}}{{I}_{{\rm{E}}}},\end{eqnarray} \tag{ 56 }$$

which is defined for $({\mu }_{{\rm{L}}}-{\mu }_{{\rm{R}}}){I}_{{\rm{M}}}\lt 0$. As expected, one can show that the upper bound is given by Carnot efficiency.

To break the inherent left–right symmetry of our model, it is necessary to use different energy dependencies of the electronic tunneling rates. We therefore consider the limit ${\delta }_{{\rm{L}}}={\delta }_{{\rm{R}}}=\delta$ but ${\varepsilon }_{{\rm{L}}}\ne {\varepsilon }_{{\rm{R}}}$. To make the effect rather strong we also have to consider δ smaller than $| {\varepsilon }_{\nu }|$, but the use of the system as a thermoelectric device does not require that ${J}_{{\rm{ph}}}\gg {\rm{\Gamma }}$. In figure 10 we plot the generated power and heat current from the hot phonon reservoir into the system versus the bias voltage. First, we see that there exists a region where the generated power becomes positive, meaning that charges are transported against the bias. This is only possible when heat flows from the phonon reservoir to the two electronic reservoirs, and within the region of positive power we can see that the efficiency (inset) in equation (56) becomes finite and reaches for our chosen parameters about half Carnot efficiency. We remark that a naive treatment of the quantum dot alone within a BMS ME would always predict I_E = 0 because $[{H}_{\mathrm{dot}},{H}_{\mathrm{ph}}]=0$ and thus, no transport of electrons against the bias would be possible in this framework. Having non-zero energy transport ${I}_{{\rm{E}}}\ne 0$ is indeed a higher order effect (i.e., beyond second order perturbation theory).

Zoom In Zoom Out Reset image size

We have presented a framework which allows to investigate thermal machines beyond the standard weak coupling and Markovian assumption. This treatment is general in the sense that it can be applied to any system where strong coupling effects and non-Markovianity are caused by a linear coupling with a bosonic bath. Our framework is by construction thermodynamically consistent because we simply apply the standard framework of quantum thermodynamics to an enlarged system. By this procedure we also automatically avoid spurious effects like efficiencies beyond Carnot. Thus, in a certain sense the RC mapping helps us to find the correct border for which it is allowed to partition the 'Universe' into a 'system' and a 'bath' part, i.e., the partition for which it is justified to apply the Born approximation.

Then, within this framework we have observed that non-Markovianity can indeed strongly enhance the efficiency of a heat engine even in the weak coupling limit. Furthermore, we also observed that strong coupling decreases the efficiency, though it is not clear at the moment whether this is a general feature. As another application we considered a single electron transistor coupled to vibrations, which is an important model to understand molecular transport. Besides studying its efficiency as a thermoelectric device we have demonstrated that the widely used PME follows from our treatment as a special case and hence, justifying the PME.

Though we believe that we have introduced a very useful and practical framework, one should also critically question it. Especially, there might be (at first sight) two weak points on which we would like to comment.

First, though it was shown in [53, 54] that the RC mapping is guaranteed to give a quasi-Ohmic (Markovian) SD after we have included a sufficent number of RCs, it is not guaranteed that the resulting final SD is also weak in the sense that it is justified to consider only second order perturbation theory in the coupling to the residual baths. Though we have seen in section 2 that the overall coupling strength of the original system to the environment can be captured by using only one RC, the resulting final SD in general depends on the shape of the initial SD (but not on its 'absolute value') and might be still large compared to parameters of the system, which has to be checked in each case separately. Still, the RC method allows us to go beyond the standard perturbative approach in a reasonable way and furthermore, by applying a conceptually similar yet technically different mapping, one can show that also the resulting SD becomes small [81].

Second, one might object that the information about what happens to the original system itself is somehow lost. Indeed, it is true that we can only give expressions for the heat flows, the power and the entropy production for the supersystem, including the RCs, but we do not know them for any particular subsystem. However, we have already stated in section 3 that the definition of heat becomes unambiguous only in case of a weakly coupled and memory-less thermal bath. In addition, we would like to defend this approach by the following remarks. First of all, the division of a thermodynamic system into subsystems (e.g., system and reservoirs) is in fact something which is only possible in the weak coupling limit. The very definition of 'system' becomes ambiguous beyond that regime (also see the discussion on the zeroth law in [8]). Second, accessing the statistics of the system alone would require to incorporate explicit measurements in our framework, a problem which was so far completely ignored also in previous work on strong coupling and non-Markovian thermodynamics. In fact, even in the weak coupling regime one already relies on abstract tools such as full counting statistics to access the statistics of energy exchanges with a heat bath [6]. However, work is usually regarded as a deterministic form of energy and as such should be easily accessible in experiments, for instance, by counting electrons as required for section 5. If the RC mapping is not directly applied to the 'work reservoir', our framework is able to predict changed work statistics which should be measurable in real experiments. Third and finally, we point out that there has been progress to understand the energetics of arbitrary multipartite systems locally [82, 83] and for special situations this is also possible for the entropic balances [84, 85] (note that the meaning of 'bipartite' in [84, 85] is more restrictive than in [82, 83]). This opens up a new interpretation of thermodynamics in the strong coupling and non-Markovian regime by recognizing the role played by the non-Markovian environment as an effective feedback controller who acts back on the system based on the information stored about it. Thus, advances in the thermodynamic understanding of multipartite systems will directly yield to new insights in the field of strong coupling and non-Markovian thermodynamics.

It is a pleasure to thank Javier Cerrillo, Mauro Cirio, Massimiliano Esposito, Jake Illes-Smith and Ahsan Nazir for valuable discussions as well as many participants of the fourth COST quantum thermodynamics conference in Erice. In addition, we explicitly wish to thank David Newman and Ahsan Nazir for sharing their related manuscript on thermodynamics and the reaction coordinate method prior to publication [39]. Furthermore, PS is especially indebted to Rocco Martinazzo for very helpful correspondence concerning the RC mapping. Financial support of the DFG (SCHA 1646/3-1, SFB 910, and GRK 1558) is gratefully acknowledged. NL is partially supported by the FY2015 Incentive Research Project.

We will here prove equation (12) following the way of Martinazzo et al [53]; for other derivations see [38, 48, 52, 54]. Because the system is completely arbitrary in our treatment, we will choose for the moment without loss of generality a particle with position q and momentum p moving in a potential V(q) and coupled via the operator s  =  q to the bath. The equations of motion according to the original Hamiltonian (2) then take on the form

$$\begin{eqnarray}&&\ddot{q}=-\displaystyle \frac{\partial V}{\partial q}+\displaystyle \sum _{k}{c}_{k}{x}_{k}-\displaystyle \sum _{k}\displaystyle \frac{{c}_{k}^{2}}{{\omega }_{k}^{2}}q,\end{eqnarray} \tag{ A.1 }$$

$$\begin{eqnarray}&&{\ddot{x}}_{k}=-{\omega }_{k}^{2}{x}_{k}+{c}_{k}q.\end{eqnarray} \tag{ A.2 }$$

After Fourier transformation according to the definition

$$\begin{eqnarray}&&\hat{f}(z)\equiv {\int }_{-\infty }^{\infty }{\rm{d}}t{{\rm{e}}}^{{\rm{i}}{zt}}f(t)\;\;\;({\mathfrak{I}}(z)\gt 0),\end{eqnarray} \tag{ A.3 }$$

we obtain

$$\begin{eqnarray}&&-{z}^{2}\hat{q}=-\displaystyle \frac{\widehat{\partial V}}{\partial q}+\displaystyle \sum _{k}{c}_{k}{\hat{x}}_{k}-\displaystyle \sum _{k}\displaystyle \frac{{c}_{k}^{2}}{{\omega }_{k}^{2}}\hat{q},\end{eqnarray} \tag{ A.4 }$$

$$\begin{eqnarray}&&-{z}^{2}{\hat{x}}_{k}=-{\omega }_{k}^{2}{\hat{x}}_{k}+{c}_{k}\hat{q}.\end{eqnarray} \tag{ A.5 }$$

Eliminating ${\hat{x}}_{k}$ we can write

$$\begin{eqnarray}&&-\displaystyle \frac{\widehat{\partial V}}{\partial q}={\hat{L}}_{0}(z)\hat{q}.\end{eqnarray} \tag{ A.6 }$$

with the Fourier space propagator

$$\begin{eqnarray}{\hat{L}}_{0}(z) & = & -{z}^{2}-\displaystyle \sum _{k}\displaystyle \frac{{c}_{k}^{2}}{{\omega }_{k}^{2}-{z}^{2}}+\displaystyle \sum _{k}\displaystyle \frac{{c}_{k}^{2}}{{\omega }_{k}^{2}},\\ & \equiv & -{z}^{2}-{W}_{0}(z)+\delta {{\rm{\Omega }}}_{0}^{2}.\end{eqnarray} \tag{ A.7 }$$

Here, we have introduced the Cauchy transform of ${J}_{0}(\omega )$:

$$\begin{eqnarray}&&{W}_{0}(z)\equiv \displaystyle \frac{2}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega \displaystyle \frac{\omega {J}_{0}(\omega )}{{\omega }^{2}-{z}^{2}}=\displaystyle \frac{1}{\pi }{\int }_{-\infty }^{\infty }{\rm{d}}\omega \displaystyle \frac{{J}_{0}(\omega )}{\omega -z},\end{eqnarray} \tag{ A.8 }$$

where we used ${J}_{0}(-\omega )=-{J}_{0}(\omega )$. As Leggett noticed [86], the SD of the bath is linked to the propagator via

$$\begin{eqnarray}&&{J}_{0}(\omega )=-\underset{\epsilon \searrow 0}{\mathrm{lim}}{\mathfrak{I}}[{\hat{L}}_{0}(\omega +{\rm{i}}\epsilon )]\equiv -{\mathfrak{I}}[{\hat{L}}_{0}^{+}(\omega )]\;\;\;(\omega \in {\mathbb{R}}).\end{eqnarray} \tag{ A.9 }$$

Especially within our notation we have ${J}_{0}(\omega )={\mathfrak{I}}[{W}_{0}^{+}(\omega )]$ what can be directly proven by use of the identity

$$\begin{eqnarray}&&\delta (\omega ^{\prime} -\omega )=\underset{\epsilon \searrow 0}{\mathrm{lim}}\displaystyle \frac{1}{\pi }\displaystyle \frac{\epsilon }{{(\omega ^{\prime} -\omega )}^{2}+{\epsilon }^{2}}.\end{eqnarray} \tag{ A.10 }$$

For the next step we look at the transformed Hamiltonian (7) to derive

$$\begin{eqnarray}&&\ddot{q}=-\displaystyle \frac{\partial V}{\partial q}+{\lambda }_{0}{X}_{1}-\delta {{\rm{\Omega }}}_{0}^{2}q,\end{eqnarray} \tag{ A.11 }$$

$$\begin{eqnarray}&&{\ddot{X}}_{1}=-{{\rm{\Omega }}}_{1}^{2}{X}_{1}+{\lambda }_{0}q+\displaystyle \sum _{k}{C}_{k}{X}_{k},\end{eqnarray} \tag{ A.12 }$$

$$\begin{eqnarray}&&{\ddot{X}}_{k}=-{{\rm{\Omega }}}_{k}^{2}{X}_{k}+{C}_{k}{X}_{1}.\end{eqnarray} \tag{ A.13 }$$

Playing the same game as above we can deduce that the Fourier space propagator for the system coordinate is

$$\begin{eqnarray}&&{\hat{L}}_{0}(z)=-{z}^{2}-\displaystyle \frac{{\lambda }_{0}^{2}}{{{\rm{\Omega }}}_{1}^{2}-{z}^{2}-{W}_{1}(z)}+\delta {{\rm{\Omega }}}_{0}^{2}\end{eqnarray} \tag{ A.14 }$$

which must be the same as equation (A.7). Furthermore, ${W}_{1}(z)$ is defined analogously to equation (A.8) with ${J}_{0}(\omega )$ replaced by ${J}_{1}(\omega )$. Then, by comparison with equation (A.7) we see that

$$\begin{eqnarray}&&{W}_{0}(z)=\displaystyle \frac{{\lambda }_{0}^{2}}{{{\rm{\Omega }}}_{1}^{2}-{z}^{2}-{W}_{1}(z)}.\end{eqnarray} \tag{ A.15 }$$

Rearranging terms we obtain

$$\begin{eqnarray}&&{W}_{1}(z)={{\rm{\Omega }}}_{1}^{2}-{z}^{2}-\displaystyle \frac{{\lambda }_{0}^{2}}{{W}_{0}(z)}\end{eqnarray} \tag{ A.16 }$$

from which we can directly deduce the relation (12). Furthermore, by noting that ${W}_{i}(0)=\delta {{\rm{\Omega }}}_{i}^{2}$ (i = 0, 1) we can also directly verify equation (15).

The Hamiltonian (or potential) of mean force is an elegant way to express the exact reduced system state of a thermal equilibrium system–bath state [64], which was also used, e.g., in [14, 15, 17]. The central idea is to introduce the Hamiltonian of mean force

$$\begin{eqnarray}&&{H}^{* }\equiv -\displaystyle \frac{1}{\beta }\mathrm{ln}\displaystyle \frac{{{\rm{tr}}}_{{\rm{E}}}\{{{\rm{e}}}^{-\beta ({H}_{{\rm{S}}}+{H}_{{\rm{I}}}+{H}_{{\rm{E}}})}\}}{{{\rm{tr}}}_{{\rm{E}}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{E}}}}\}}\end{eqnarray} \tag{ B.1 }$$

such that the reduced system state can be expressed as

$$\begin{eqnarray}&&{\rho }_{{\rm{S}}}=\displaystyle \frac{{{\rm{e}}}^{-\beta {H}^{* }}}{{Z}_{{\rm{S}}}^{* }}\end{eqnarray} \tag{ B.2 }$$

with the partition function ${Z}_{{\rm{S}}}^{* }={{\rm{tr}}}_{{\rm{S}}}\{{{\rm{e}}}^{-\beta {H}^{* }}\}$. Indeed, it is straightforward to show that equation (B.2) coincides with the reduced equilibirum state of system and environment, i.e.,

$$\begin{eqnarray}&&{\rho }_{{\rm{S}}}={{\rm{tr}}}_{{\rm{E}}}\left\{\displaystyle \frac{{{\rm{e}}}^{-\beta ({H}_{{\rm{S}}}+{H}_{{\rm{I}}}+{H}_{{\rm{E}}})}}{{{\rm{tr}}}_{\mathrm{SE}}\{{{\rm{e}}}^{-\beta ({H}_{{\rm{S}}}+{H}_{{\rm{I}}}+{H}_{{\rm{E}}})}\}}\right\}.\end{eqnarray} \tag{ B.3 }$$

To see this it suffices to note that by definition ${Z}_{{\rm{S}}}^{* }={{\rm{tr}}}_{\mathrm{SE}}\{{{\rm{e}}}^{-\beta ({H}_{{\rm{S}}}+{H}_{{\rm{I}}}+{H}_{{\rm{E}}})}\}/{{\rm{tr}}}_{{\rm{E}}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{E}}}}\}$.

To explore the connection of equation (B.2) with the equilibrium state (30) stated in the main text, we perform the RC mapping on the global system–bath Hamiltonian:

$$\begin{eqnarray}&&{H}_{{\rm{S}}}+{H}_{{\rm{I}}}+{H}_{{\rm{E}}}={H}_{{\rm{S}}}^{\prime }+{H}_{{\rm{I}}}^{\prime }+{H}_{{\rm{E}}}^{\prime }\end{eqnarray} \tag{ B.4 }$$

with ${H}_{{\rm{S}}}^{\prime }={H}_{{\rm{S}}}+{H}_{{\rm{S}}\mbox{--}\mathrm{RC}}+{H}_{\mathrm{RC}}$ describing the system, system-RC coupling and the RC respectively, and ${H}_{{\rm{I}}}^{\prime }$ describes the coupling to the residual bath described by ${H}_{{\rm{E}}}^{\prime }$. Since this mapping is exact, we can express the Hamiltonian of mean force (and consequently the partition function Z_S) in terms of the transformed Hamiltonian:

$$\begin{eqnarray}&&{{\rm{e}}}^{-\beta {H}^{* }}=\displaystyle \frac{{{\rm{tr}}}_{\mathrm{RC},E^{\prime} }\{{{\rm{e}}}^{-\beta ({H}_{{\rm{S}}}^{\prime }+{H}_{{\rm{I}}}^{\prime }+{H}_{{\rm{E}}}^{\prime })}\}}{{{\rm{tr}}}_{\mathrm{RC},E^{\prime} }\{{{\rm{e}}}^{-\beta ({H}_{\mathrm{RC}}+{H}_{{\rm{I}}}^{\prime }+{H}_{{\rm{E}}}^{\prime })}\}}.\end{eqnarray} \tag{ B.5 }$$

Now, if it is justified to regard the coupling ${H}_{{\rm{I}}}^{\prime }$ to the residual bath as weak compared to all other contributions, we obtain to lowest (i.e., zeroth) order in H_I the Hamiltonian of mean force and partition function

$$\begin{eqnarray}&&{{\rm{e}}}^{-\beta {H}^{* }}\approx \displaystyle \frac{{{\rm{tr}}}_{\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}\}}{{{\rm{tr}}}_{\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{\mathrm{RC}}}\}},\;\;\;{Z}_{{\rm{S}}}^{* }\approx \displaystyle \frac{{{\rm{tr}}}_{{\rm{S}},\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}\}}{{{\rm{tr}}}_{\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{RC}}}\}}.\end{eqnarray} \tag{ B.6 }$$

Here, we used $[{H}_{\mathrm{RC}},{H}_{{\rm{E}}}^{\prime }]=0$, which is guaranteed by construction. Thus, to lowest order

$$\begin{eqnarray}&&{\rho }_{{\rm{S}}}\approx \displaystyle \frac{{{\rm{tr}}}_{\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}\}}{{{\rm{tr}}}_{{\rm{S}},\mathrm{RC}}\{{{\rm{e}}}^{-\beta {H}_{{\rm{S}}}^{\prime }}\}}\end{eqnarray} \tag{ B.7 }$$

as described by equation (30) showing the consistency of our approach with standard results from equilibrium statistical mechanics. Especially, note that it is computationally relatively cheap to compute equation (30) as compared to the exact result (B.2).

We here provide details for the derivation of the Born–Markov ME based on the system Hamiltonian ${H}_{{\rm{S}}}^{\prime }$ (37), which treats the system-RC coupling non-perturbatively while we are aiming at a perturbative treatment of the coupling to the remaining reservoirs. The coupling Hamiltonian is

$$\begin{eqnarray}&&{H}_{{\rm{I}}}=\displaystyle \sum _{\nu }{H}_{{\rm{I}}}^{(\nu )}=\displaystyle \sum _{\nu }{s}_{\nu }\otimes {B}_{\nu }\end{eqnarray} \tag{ C.1 }$$

with s_h and s_w given in equation (32) and ${s}_{c}={X}_{1}$ due to the RC mapping. Furthermore, for $\nu \in \{h,w\}$ the coupling operators of the bath are ${B}_{\nu }=-{\sum }_{k}{c}_{k\nu }{x}_{k\nu }$ and the free bath Hamiltonian reads ${H}_{{\rm{B}}}^{(\nu )}=\tfrac{1}{2}{\sum }_{k}({p}_{k\nu }^{2}+{\omega }_{k\nu }^{2}{x}_{k\nu }^{2})$. For $\nu =c$ the form of the operators remains the same but we have to substitute ${x}_{{kc}}\to {X}_{{kc}}$, ${p}_{{kc}}\to {P}_{{kc}}$, ${c}_{{kc}}\to {C}_{{kc}}$ and ${\omega }_{{kc}}\to {{\rm{\Omega }}}_{{kc}}$.

Note that we are here neglecting renormalization terms of the form $\tfrac{1}{2}{\sum }_{k}\tfrac{{c}_{k\nu }^{2}}{{\omega }_{k\nu }^{2}}{s}_{\nu }^{2}$ which are of second order in the system–bath coupling. We will later on neglect any Lamb shift terms as well which are also of second order. This is consistent because the heat flows (25) are itself already of second order in the coupling. The contribution due to the renormalization and Lamb shift terms would then be of fourth order in total, which is beyond the validity of our perturbative approach.

After applying the Born and Markov approximation, the formal starting point of the ME is the second order equation in the interaction picture [8, 9, 61] (the interaction picture is defined by $\tilde{A}(t)\equiv {{\rm{e}}}^{+{\rm{i}}({H}_{{\rm{S}}}^{\prime }+{H}_{{\rm{B}}})t}A{{\rm{e}}}^{-{\rm{i}}({H}_{{\rm{S}}}^{\prime }+{H}_{{\rm{B}}})t}$)

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\tilde{\rho }(t) & = & -\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau {{\rm{tr}}}_{\nu }\{{\tilde{H}}_{{\rm{I}}}^{(\nu )}(t){\tilde{H}}_{{\rm{I}}}^{(\nu )}(t-\tau )\tilde{\rho }(t){R}_{0}^{(\nu )}\\ & & -{\tilde{H}}_{{\rm{I}}}^{(\nu )}(t)\tilde{\rho }(t){R}_{0}^{(\nu )}{\tilde{H}}_{{\rm{I}}}^{(\nu )}(t-\tau )+{\rm{h.c.}}\}.\end{eqnarray} \tag{ C.2 }$$

Here, $\tilde{\rho }(t)$ is the density matrix of the supersystem (in the interaction picture) and ${R}_{0}^{(\nu )}\sim {{\rm{e}}}^{-{\beta }_{\nu }{H}_{{\rm{B}}}^{(\nu )}}$ the equilibrium density matrix of the reservoir. Furthermore, the fact that the ME additively decomposes into contributions from each reservoir ν is due to the fact that ${{\rm{tr}}}_{\nu }\{{H}_{{\rm{I}}}^{(\nu )}{R}_{0}^{(\nu )}\}=0$.

After introducing the bath correlation function

$$\begin{eqnarray}&&{C}_{\nu }(\tau )\equiv {{\rm{tr}}}_{\nu }\{{\tilde{B}}_{\nu }(\tau ){B}_{\nu }{R}_{0}^{(\nu )}\},\end{eqnarray} \tag{ C.3 }$$

the ME (C.2) can be expressed as

$$\begin{eqnarray*}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\tilde{\rho }(t) & = & -\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau \{{C}_{\nu }(\tau )[{\tilde{s}}_{\nu }(t){\tilde{s}}_{\nu }(t-\tau )\tilde{\rho }(t)-{\tilde{s}}_{\nu }(t-\tau )\tilde{\rho }(t){\tilde{s}}_{\nu }(t)]\\ & & +{C}_{\nu }^{* }(\tau )[\tilde{\rho }(t){\tilde{s}}_{\nu }(t-\tau ){\tilde{s}}_{\nu }(t)-{\tilde{s}}_{\nu }(t)\tilde{\rho }(t){\tilde{s}}_{\nu }(t-\tau )]\}.\end{eqnarray*}$$

If we split the correlation function into real and imaginary part, ${C}_{\nu }(\tau )={\mathfrak{R}}{C}_{\nu }(\tau )+{\rm{i}}{\mathfrak{I}}{C}_{\nu }(\tau )$, we can write

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\tilde{\rho }(t) & = & -\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau {\mathfrak{R}}{C}_{\nu }(\tau )[{\tilde{s}}_{\nu }(t),[{\tilde{s}}_{\nu }(t-\tau ),\tilde{\rho }(t)]]\\ & & -{\rm{i}}\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau {\mathfrak{I}}{C}_{\nu }(\tau )[{\tilde{s}}_{\nu }(t),\{{\tilde{s}}_{\nu }(t-\tau ),\tilde{\rho }(t)\}].\end{eqnarray} \tag{ C.4 }$$

Because

$$\begin{eqnarray}{C}_{\nu }(\tau ) & = & \displaystyle \sum _{k}\displaystyle \frac{{c}_{k\nu }^{2}}{2{\omega }_{k\nu }}\{[1+{n}_{\nu }({\omega }_{k\nu })]{{\rm{e}}}^{-i{\omega }_{k\nu }\tau }+{n}_{\nu }({\omega }_{k\nu }){{\rm{e}}}^{{\rm{i}}{\omega }_{k\nu }\tau }\}\\ & = & \displaystyle \frac{1}{\pi }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\omega {J}^{(\nu )}(\omega )\{[1+{n}_{\nu }(\omega )]{{\rm{e}}}^{-{\rm{i}}\omega \tau }+{n}_{\nu }(\omega ){{\rm{e}}}^{{\rm{i}}\omega \tau }\},\end{eqnarray} \tag{ C.5 }$$

we obtain

$$\begin{eqnarray}&&{\mathfrak{R}}{C}_{\nu }(\tau )=\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega {J}^{(\nu )}(\omega )\mathrm{cos}(\omega \tau )\mathrm{coth}\displaystyle \frac{{\beta }_{\nu }\omega }{2},\end{eqnarray} \tag{ C.6 }$$

$$\begin{eqnarray}&&{\mathfrak{I}}{C}_{\nu }(\tau )=-\displaystyle \frac{1}{\pi }{\int }_{0}^{\infty }{\rm{d}}\omega {J}^{(\nu )}(\omega )\mathrm{sin}(\omega \tau ).\end{eqnarray} \tag{ C.7 }$$

Here, we used that $1+2{n}_{\nu }(\omega )=\mathrm{coth}\tfrac{{\beta }_{\nu }\omega }{2}$. Finally, after leaving the interaction picture the ME (C.4) becomes

$$\begin{eqnarray}\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\rho (t) & = & -{\rm{i}}[{H}_{{\rm{S}}}^{\prime },\rho (t)]\\ & & -\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau {\mathfrak{R}}{C}_{\nu }(\tau )[{s}_{\nu },[{\tilde{s}}_{\nu }(-\tau ),\rho (t)]]\\ & & -{\rm{i}}\displaystyle \sum _{\nu }{\displaystyle \int }_{0}^{\infty }{\rm{d}}\tau {\mathfrak{I}}{C}_{\nu }(\tau )[{s}_{\nu },\{{\tilde{s}}_{\nu }(-\tau ),\rho (t)\}].\end{eqnarray} \tag{ C.8 }$$

This form is still not very useful for numerical implementation. To achieve this goal we write

$$\begin{eqnarray}&&{s}_{\nu }(-\tau )=\displaystyle \sum _{{kl}}{s}_{\nu }^{{kl}}{{\rm{e}}}^{-{\rm{i}}{\omega }_{{kl}}\tau }| k\rangle \langle l| \end{eqnarray} \tag{ C.9 }$$

with ${H}_{{\rm{S}}}^{\prime }| k\rangle ={E}_{k}| k\rangle$ and ${\omega }_{{kl}}\equiv {E}_{k}-{E}_{l}$. Upon using the identity

$$\begin{eqnarray}&&{\int }_{0}^{\infty }{\rm{d}}\tau {{\rm{e}}}^{\pm {\rm{i}}\omega \tau }=\pi \delta (\omega )\pm {\rm{i}}{ \mathcal P }\displaystyle \frac{1}{\omega }\end{eqnarray} \tag{ C.10 }$$

and neglecting the principal value (Lamb shift) terms, we arrive at a ME of the form

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\rho (t)=-{\rm{i}}[{H}_{{\rm{S}}}^{\prime },\rho (t)]-\displaystyle \sum _{\nu }[{s}_{\nu },[{\chi }_{\nu },\rho (t)]+\displaystyle \sum _{\nu }[{s}_{\nu },\{{{\rm{\Theta }}}_{\nu },\rho (t)\}].\end{eqnarray} \tag{ C.11 }$$

The new operators appearing in this equation are given as

$$\begin{eqnarray}&&{\chi }_{\nu }=\displaystyle \frac{1}{2}\displaystyle \sum _{{kl}}{J}^{(\nu )}({\omega }_{{kl}})\mathrm{coth}\displaystyle \frac{{\beta }_{\nu }{\omega }_{{kl}}}{2}{s}_{\nu }^{{kl}}| k\rangle \langle l| ,\end{eqnarray} \tag{ C.12 }$$

$$\begin{eqnarray}&&{{\rm{\Theta }}}_{\nu }=\displaystyle \frac{1}{2}\displaystyle \sum _{{kl}}{J}^{(\nu )}({\omega }_{{kl}}){s}_{\nu }^{{kl}}| k\rangle \langle l| .\end{eqnarray} \tag{ C.13 }$$

Finally, the model is specified by choosing the SDs of the baths as discussed in the main text. For the work bath we again use ${J}^{(w)}(\omega )={J}_{0}^{(w)}(\omega )={\beta }_{w}{{\rm{\Delta }}}_{10}{{\rm{\Gamma }}}_{w}$ to mimic a constant SD in the infinite temperature limit. The SD of the hot bath is parameterized by ${J}^{(h)}(\omega )={J}_{0}^{(h)}(\omega )=\tfrac{{{\rm{\Gamma }}}_{h}\omega }{{{\rm{\Delta }}}_{20}}{\rm{\Theta }}({\omega }_{{\rm{R}}}-\omega )$ and for the cold bath ${J}^{(c)}(\omega )={J}_{1}^{(c)}(\omega )$ is given by the relation (12) with ${J}_{0}^{(c)}(\omega )$ given in equation (35). We are especially interested in the regime of a large cutoff frequency ${\omega }_{{\rm{R}}}\gg 1$ because this allows by virtue of the residue theorem to evaluate ${J}_{1}^{(c)}(\omega )$ exactly. Note that the residue theorem requires ${J}_{0}(\omega )$ to be analytic (except for isolated poles), which is—strictly speaking—never the case for a hard cutoff ${\rm{\Theta }}({\omega }_{{\rm{R}}}-\omega )$. However, the discrepancy with the true solution vanishs for ${\omega }_{{\rm{R}}}\to \infty$. Then, for $4{\omega }_{0}^{2}\gt {\gamma }^{2}$, it follows that

$$\begin{eqnarray}&&{J}_{1}^{(c)}(\omega )=\gamma \omega {\rm{\Theta }}({\omega }_{{\rm{R}}}-\omega ).\end{eqnarray} \tag{ C.14 }$$

Furthermore, we also have ${\lambda }_{0}={d}_{0}$ and $\delta {{\rm{\Omega }}}_{0}={d}_{0}/{\omega }_{0}$ such that the frequency of the RC is $\tfrac{{\lambda }_{0}}{\delta {{\rm{\Omega }}}_{0}}={\omega }_{0}$.

The secular ME requires to perform an additional approximation on top of the Born–Markov ME (C.11). This can be done by averaging the generator of the Born–Markov ME in the interaction picture in time (similar to a rotating wave approximation) [6, 61] or by dynamical coarse graining of the time-evolution [9, 87]. We will skip any details here because the secular ME is also reviewed in appendix D.

For a non-degenerate system Hamiltonian, it is well-known that the BMS ME yields a simple rate equation ('Pauli ME') in the eigenbasis of H_S [4, 9, 61]. This can be put into the form

$$\begin{eqnarray}&&{\dot{P}}_{k}=\displaystyle \sum _{l}{W}_{{kl}}{P}_{l},\end{eqnarray} \tag{ D.1 }$$

where P_k is the probability to find the system in state $| k\rangle$ and the transition rate from energy eigenstate l to k is given by

$$\begin{eqnarray}&&{W}_{{kl}}=\displaystyle \sum _{\alpha \beta }{\gamma }_{\alpha \beta }({E}_{l}-{E}_{k})\langle k| {A}_{\beta }| l\rangle {\langle k| {A}_{\alpha }^{\dagger }| l\rangle }^{* }.\end{eqnarray} \tag{ D.2 }$$

Here, we assumed a general interaction Hamiltonian of the form ${H}_{{\rm{I}}}={\sum }_{\alpha }{A}_{\alpha }\otimes {B}_{\alpha }$. Furthermore, the ${\gamma }_{\alpha \beta }$ denote the Fourier transforms of the reservoir correlation functions

$$\begin{eqnarray}&&{\gamma }_{\alpha \beta }(\omega )={\int }_{-\infty }^{\infty }{\rm{d}}\tau {{\rm{tr}}}_{{\rm{B}}}\{{B}_{\alpha }(\tau ){B}_{\beta }{R}_{0}\}{{\rm{e}}}^{+{\rm{i}}\omega \tau }.\end{eqnarray} \tag{ D.3 }$$

For our model from section 5, we identify the coupling operators

$$\begin{eqnarray}&&{A}_{1}=d,\qquad {B}_{1}=\displaystyle \sum _{k}{t}_{k{\rm{L}}}{c}_{k{\rm{L}}}^{\dagger },\end{eqnarray} \tag{ D.4 }$$

$$\begin{eqnarray}&&{A}_{2}={d}^{\dagger },\qquad {B}_{2}=\displaystyle \sum _{k}{t}_{k{\rm{L}}}^{* }{c}_{k{\rm{L}}},\end{eqnarray} \tag{ D.5 }$$

$$\begin{eqnarray}&&{A}_{3}=d,\qquad {B}_{3}=\displaystyle \sum _{k}{t}_{k{\rm{R}}}{c}_{k{\rm{R}}}^{\dagger },\end{eqnarray} \tag{ D.6 }$$

$$\begin{eqnarray}&&{A}_{4}={d}^{\dagger },\qquad {B}_{4}=\displaystyle \sum _{k}{t}_{k{\rm{R}}}^{* }{c}_{k{\rm{R}}},\end{eqnarray} \tag{ D.7 }$$

$$\begin{eqnarray}&&{A}_{5}=\displaystyle \frac{1}{\sqrt{2{{\rm{\Omega }}}_{1}}}({a}_{1}+{a}_{1}^{\dagger })\end{eqnarray} \tag{ D.8 }$$

$$\begin{eqnarray}&&{B}_{5}=-\displaystyle \sum _{q}\displaystyle \frac{{C}_{q}}{\sqrt{2{{\rm{\Omega }}}_{q}}}({a}_{q}+{a}_{q}^{\dagger })\end{eqnarray} \tag{ D.9 }$$

and the non-vanishing correlation functions yield

$$\begin{eqnarray}&&{\gamma }_{12}(\omega )={{\rm{\Gamma }}}_{{\rm{L}}}(-\omega ){f}_{{\rm{L}}}(-\omega ),\end{eqnarray} \tag{ D.10 }$$

$$\begin{eqnarray}&&{\gamma }_{21}(\omega )={{\rm{\Gamma }}}_{{\rm{L}}}(+\omega )[1-{f}_{{\rm{L}}}(+\omega )],\end{eqnarray} \tag{ D.11 }$$

$$\begin{eqnarray}&&{\gamma }_{34}(\omega )={{\rm{\Gamma }}}_{{\rm{R}}}(-\omega ){f}_{{\rm{R}}}(-\omega ),\end{eqnarray} \tag{ D.12 }$$

$$\begin{eqnarray}&&{\gamma }_{43}(\omega )={{\rm{\Gamma }}}_{{\rm{R}}}(+\omega )[1-{f}_{{\rm{R}}}(+\omega )],\end{eqnarray} \tag{ D.13 }$$

$$\begin{eqnarray}&&{\gamma }_{55}(\omega )={J}_{1}(+\omega )[1+{n}_{{\rm{B}}}(+\omega )].\end{eqnarray} \tag{ D.14 }$$

Here, we have used the definition of the fermionic SD (53) and the SD of the residual oscillators is as usual defined by

$$\begin{eqnarray}&&{J}_{1}(\omega )=\displaystyle \frac{\pi }{2}\displaystyle \sum _{q}\displaystyle \frac{{C}_{q}^{2}}{{{\rm{\Omega }}}_{q}}\delta (\omega -{{\rm{\Omega }}}_{q}).\end{eqnarray} \tag{ D.15 }$$

From equation (D.2) and equations (D.10)–(D.14) we see that the rates additively decompose into left, right, and phonon contributions and the total rate matrix in equation (D.1) has the structure $W={W}^{{\rm{L}}}+{W}^{{\rm{R}}}+{W}^{\mathrm{ph}}$, where for $k\ne l$ we have

$$\begin{eqnarray}&&{W}_{{kl}}^{{\rm{L}}}={\gamma }_{12}({E}_{l}-{E}_{k}){| \langle k| {d}^{\dagger }| l\rangle | }^{2}+{\gamma }_{21}({E}_{l}-{E}_{k}){| \langle k| d| l\rangle | }^{2},\end{eqnarray} \tag{ D.16 }$$

$$\begin{eqnarray}&&{W}_{{kl}}^{{\rm{R}}}={\gamma }_{34}({E}_{l}-{E}_{k}){| \langle k| {d}^{\dagger }| l\rangle | }^{2}+{\gamma }_{43}({E}_{l}-{E}_{k}){| \langle k| d| l\rangle | }^{2},\end{eqnarray} \tag{ D.17 }$$

$$\begin{eqnarray}&&{W}_{{kl}}^{\mathrm{ph}}={\gamma }_{55}({E}_{l}-{E}_{k})\displaystyle \frac{1}{2{{\rm{\Omega }}}_{1}}{| \langle k| ({a}_{1}+{a}_{1}^{\dagger })| l\rangle | }^{2}.\end{eqnarray} \tag{ D.18 }$$

To make the transition rates explicit, we have to evaluate the matrix elements of the system coupling operators, too. The electronic jumps can be separated into pure electronic transitions and bosonic excitations of the RC. Denoting by $| \widetilde{n^{\prime} ,m^{\prime} }\rangle \equiv U| n,m\rangle$ the basis in the original frame, we obtain

$$\begin{eqnarray}{| \langle n^{\prime} ,m^{\prime} | {A}_{1}| n,m\rangle | }^{2} & = & {\left|\left\langle \widetilde{n^{\prime} ,m^{\prime} }\left|d{{\rm{e}}}^{-\tfrac{\lambda }{{{\rm{\Omega }}}_{1}}({a}_{1}-{a}_{1}^{\dagger })}\right|\widetilde{n,m}\right\rangle \right|}^{2}\\ & = & {\delta }_{n^{\prime} ,0}{\delta }_{n,1}{\left|\left\langle \widetilde{m^{\prime} }\left|{{\rm{e}}}^{-\tfrac{\lambda }{{{\rm{\Omega }}}_{1}}({a}_{1}-{a}_{1}^{\dagger })}\right|\tilde{m}\right\rangle \right|}^{2},\end{eqnarray} \tag{ D.19 }$$

$$\begin{eqnarray}&&{| \langle n^{\prime} ,m^{\prime} | {A}_{2}| n,m\rangle | }^{2}={\delta }_{n^{\prime} ,1}{\delta }_{n,0}{\left|\left\langle \widetilde{m^{\prime} }\left|{{\rm{e}}}^{+\tfrac{\lambda }{{{\rm{\Omega }}}_{1}}({a}_{1}-{a}_{1}^{\dagger })}\right|\tilde{m}\right\rangle \right|}^{2},\end{eqnarray} \tag{ D.20 }$$

whereas the transitions triggered by the phonon reservoir simply yield

$$\begin{eqnarray}{| \langle n^{\prime} ,m^{\prime} | {A}_{5}| n,m\rangle | }^{2} & = & \displaystyle \frac{{\delta }_{n,n^{\prime} }}{2{{\rm{\Omega }}}_{1}}{| \langle \widetilde{m^{\prime} }| ({a}_{1}+{a}_{1}^{\dagger })| \tilde{m}\rangle | }^{2}\\ & = & \displaystyle \frac{{\delta }_{n,n^{\prime} }}{2{{\rm{\Omega }}}_{1}}[{\delta }_{m^{\prime} ,m+1}(m+1)+{\delta }_{m^{\prime} ,m-1}m].\end{eqnarray} \tag{ D.21 }$$

Having the rates at hand, we are now finally able to compute the quantities shown in section 5.3. A visualization of the resulting rate equation, which has a highly connected structure, is also provided in figure 8.