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 (RC) 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 RC), 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.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][7][8][9][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 . 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][12][13][14][15][16][17][18][19][20][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][12][13] and/or coupling only to a single thermal reservoir [11,[13][14][15][16][17];or they remain very formal [18][19][20][21]. Furthermore, model-specific studies are either based on simple or exactly solvable models from the field of quantum transport [22][23][24][25] and quantum Brownian motion [26][27][28][29], or spin-boson models [30][31][32] often in combination with specific transformations applicable only to special Hamiltonians (polaron transformations) [31,[33][34][35][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 studi...