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

We propose a method to study the thermodynamic behaviour of small systems beyond the weak coupling and Markovian approximation, which is different in spirit from conventional approaches. The idea is to redefine the system and environment such that the effective, redefined system is again coupled weakly to Markovian residual baths and thus, allows to derive a consistent thermodynamic framework for this new system-environment partition. To achieve this goal we make use of the reaction coordinate mapping, which is a general method in the sense that it can be applied to an arbitrary (quantum or classical and even time-dependent) system coupled linearly to an arbitrary number of harmonic oscillator reservoirs. The core of the method relies on an appropriate identification of a part of the environment (the reaction coordinate), which is subsequently included as a part of the system. We demonstrate the power of this concept by showing that non-Markovian effects can significantly enhance the steady state efficiency of a three-level-maser heat engine, even in the regime of weak system-bath coupling. Furthermore, we show for a single electron transistor coupled to vibrations that our method allows one to justify master equations derived in a polaron transformed reference frame.


Introduction
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

Reaction coordinate mapping
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] ( ) ( ) with mass-weighted positions x k and momenta p k of the bath fulfilling [ ] d = x p , i k l kl (we set  º 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 å w s k c 1 2 2 k k 2 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: The SD is a positive function for w > 0 and must fulfill ( ) w  J 0 0 for w  0 and w  ¥. 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 1 (see figure 1). Hence, we seek a transformation which maps where l 0 is an unspecified parameter so far. More formally, we perform a normal mode transformation of the form . 5 Here, we used a vector notation for the collection of original and transformed bath coordinates and momenta, e.g., , , , ,  Here, we used l L =c k k 1 0 1 which follows from equation (4). Furthermore, we defined w W º å L k k k 1 At this point we can already recognize an important property of the mapping. Suppose that we scale the coupling coefficients c k by  a c c k k for some  a Î . Then, the only parameters influenced by this will be l 0 Figure 1. Sketch of the RC mapping. Before the mapping (left figure) the system can be visualized as being coupled to a large number of harmonic oscillators, see equation (2). After the mapping (right figure) the system couples to the RC only, but in turn the RC is coupled to a large number of residual oscillators as described by the Hamiltonian (17). and dW 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 ( ) w J 0 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: Furthermore, one can also show that (see [53] or appendix A again) which makes its thermodynamic stability evident. Hence, the physical frequency of the RC is not given by W 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][53][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 6 . 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.
should be continuous and strictly positive for ( ) w w Î 0, R and zero for  w w R , where w R denotes a cutoff frequency [53]. 6 In the theory of Brownian motion, Markovian behaviour is ensured by an Ohmic SD which scales linearly with ω up to a high enough frequency cutoff w R and then falls off to zero [55]. We stress, however, that the correct definition of non-Markovianity in the quantum mechanical context is non-trivial, may not be guaranteed by this condition alone, and is under much debate at the moment, see, e.g., [58][59][60].

Nonequilibrium thermodynamics within the RC framework
We now imagine the situation where our system is coupled to several reservoirs labeled by ν and where the timedependent driving is responsible for work extraction and injection. The obvious generalization of the Hamiltonian(2) to this situation is 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 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.
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][41][42][43][44][45][46].
The standard framework of quantum thermodynamics starts with a microscopically derived ME of the form  Figure 2. Sketch of a system coupled to a hot reservoir (red, blurred oscillators) and a cold reservoir (blue oscillators). After the mapping we have, as an example, included one RC for each reservoir to account for non-Markovian and strong coupling effects. Note, however, that we do not have to include a RC if the reservoir is weakly coupled and Markovian, or we might have to include two or more RCs in case of strong non-Markovianity. After the mapping we then treat the system and RCs as one new system as indicated by the shaded grey box.
where b n 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 7 [4]. We remark that, strictly speaking, a ME of the form (19) can only be derived for a slow timedependence 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 The first law of thermodynamics then acquires the form where we identified the rate of work (power) done on the supersystem and the heat flow coming from reservoir ν as 8 . The second law stipulates that the rate of entropy productioṅ ( ) S t i is always positive˙(  . 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 becomė˙( 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.,˙¹ W 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 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., ( ) r  ¥~bt e H S S . Experimentally, deviations from the canonical state b e H S might be a clear indicator for persistent system environment correlations making it necessary to go beyond the Born approximation, e.g., by 7 If we also allow for particle transport by coupling the system to a particle reservoir with chemical potential m n , we have the relation instead where ¢ N S is the particle number operator of the supersystem. 8 In presence of particle transport the heat flow˙( ) flowing into the supersystem. The first law then predicts energy conservation, , and particle number conservation, using the RC mapping. Given the SD of the reservoir, it should be then possible to test the prediction (30) in an actual experiment.

Application I: three-level-maser heat engine 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 The idea behind this engine is the model of a simple maser which is lasing at a particular transition, say « 0 1. This lasing corresponds to 'work' output and it is achieved due to population inversion between the levels | ñ 0 and | ñ 1 . This in turn can be mediated via a third level | ñ 2 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][67][68][69][70].
The coupling to the reservoirs is mediated by the system operators Here, in units where  = 1, the parameter  has units of an energy or frequency such that s ν has the same units as the coordinates n x k, of the bath oscillators. However, in all numerical calculations which follow we will simply set  º 1. Within the standard approach (BMS approximation for the system only) the thermal generators ( )  n in equation (19) of each bath become Here, we have introduced the dissipator To quantify the performance of work extraction (i.e.,˙< W 0), we introduce the efficiency of the heat engine: Figure 3. Sketch of the maser heat engine where the working medium comprises three discrete levels and each transition is coupled to a separate reservoir called the hot ('h', red), cold ('c', blue) and work ('w', green) reservoir. The black arrows indicate the direction in which we define energy flows to be positive.
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 D D 10 20 independent of all other parameters (as long as˙< W 0, of course). Note that to completely justify the notion of a 'work reservoir', we should take the limit b  0 w [69,70], which complies with Sommerfelds notion of temperature as the 'work value of heat' [71]. In this limit the entropy changeḃ -W w in the work reservoir w goes to zero, the second law acquires the forṁ˙ however, we recognize that the Bose distribution ( ) D n w 10 diverges for b  0 w (thoughẆ remains finite), unless we additionally require that the bath SD scales like ( ) for small b w . In this ideal limit we then obtain i.e., the rates of upward and downward transitions are equal. However, for all numerical results reported below we take b 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 D D 10 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.

Thermodynamics with RC
For definiteness we choose the SD of the cold bath to be parameterized as 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  (determined by equation (12), orange, dashed) over w w 0 for two different values of γ and d 0 : 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 . This Hamiltonian describes a Rabi model (harmonic oscillator coupled to a two-level system) plus one additional energy level | ñ 0 . Especially note that the energy levels of state | ñ 1 and | ñ 2 get both shifted by the same amount  dW 4 0 2 . 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 S and dissipators ( )  n . 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˜˜˜( 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˜< W 0 as˜( h which is also bounded by the Carnot efficiency in the limit b  0 w . 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.
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 w D » 21 0 (interestingly, this effect cancels when we compute the efficiency). In fact, if this resonance condition is met 9 , 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 ( ) ( ) w J c 0 . 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 10 . Physically, the reason for this can be traced back to the Purcell effect which predicts an enhanced spontaneous emission rate for the « 1 2 transition in presence of a (resonant) cavity [76] and thus allows for a stronger population inversion between the states | ñ 0 and | ñ 1 . 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 ( ) ( ) w J c 1 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.
Finally, figure 7 shows the thermodynamics as a function of the system-RC coupling strength d 0 (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 0 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 instead). In fact, for all parameters we have checked we always observed a decreasing efficiency as a function of d 0 . Whether this is a general feature or model specific remains an open question, which might be eventually answerable within the RC framework.

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 9 Note that w 0 equals the frequency of the RC for large cutoff frequency w R , see appendix C. 10 Note that [37] also demonstrates an increased ability to extract work due to non-Markovian effects. However, [37] did not study the efficiency and focused on transient effects, making it unclear whether a steadily working heat engine can be better than its Markovian counterpart. studied, e.g., in [33,34,[77][78][79][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 The electronic part of the system is described by Here, we denoted the coupling coefficients to the phonon bath by h q because we already use n c k 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)) We now identify our supersystem to be˜(  (17), but now the Hamiltonian is closer to the literature to which we wish to compare our results [33]. Furthermore, the spectrum of ¢ H S is still bounded from below for all coupling strengths l 0 .
Below it will be convenient to work with the bosonic ladder operators, which are related to position and momentum operators via In terms of these operators we have˜( 0, 1, 2 ,... ) quantifies the electronic (bosonic) occupation, and eigenenergies ( ) 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 coincides with W 1 ). To completely specify our model we choose to define and parametrize the SD of the fermionic leads as Here, d n quantifies the width and e n the maximum of the Lorentzian function. Furthermore, the SD of the residual phonon bath is choosen to be ohmic with exponential cutoff which for large cutoff frequencies w 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. Before we proceed to show numerical results, let us briefly explain an alternative way to treat strong systemphonon 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.

Results
In figure 9 we plot the matter current from left to right versus the difference in chemical potentials m m = -V L R of the electronic reservoirs. First, we see that at low electronic temperatures, the current displays multiple steps, which occur at = D V E 2 i i , where DE i are the transition frequencies of the system. Thus, the steps enable one to deduce the renormalized electron energy and the phonon frequency W 1 from the electronic current. Furthermore, if we truncate the phonon Hilbert space at a small cutoff number N c , we see that only few plateaus are visible since the number of possible transitions is bounded. Most notably, however, when both N c is , . The phonon bath (red, vertical arrows) can induce transitions between states | | ñ « + ñ n m n m , , 1whereas the electronic leads (blue, diagonal arrows) always change the occupation number n of the quantum dot. Furthermore, the jump of an electron in or out of the system can be also accompanied by multi-phonon transitions as exemplarily indicated with thin, dashed, magenta lines. 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.
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 d d = L R and e e = L 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)˙˙( . 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 d d d = = L R but e e ¹ L R . To make the effect rather strong we also have to consider δ smaller than | | e n , but the use of the system as a thermoelectric device does not require that  G J ph . 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 and thus, no transport of electrons against the bias would be possible in this framework. Having non-zero energy transport ¹ I 0 E is indeed a higher order effect (i.e., beyond second order perturbation theory).

Final remarks
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.
Eliminatingx k we can writeˆ( I what can be directly proven by use of the identity For the next step we look at the transformed Hamiltonian (7) to derive Playing the same game as above we can deduce that the Fourier space propagator for the system coordinate iŝ Rearranging terms we obtain from which we can directly deduce the relation (12). Furthermore, by noting that ( ) d = W W 0 i i 2 (i=0, 1) we can also directly verify equation (15).

Appendix B. Non-canonical equilibriums states and the potential of mean force
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  . 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.

Appendix D. BMS ME for application II
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 where P k is the probability to find the system in state | ñ k and the transition rate from energy eigenstate l to k is given by