A thermodynamically consistent Markovian master equation beyond the secular approximation

Markovian master equations provide a versatile tool for describing open quantum systems when memory effects of the environment may be neglected. As these equations are of an approximate nature, they often do not respect the laws of thermodynamics when no secular approximation is performed in their derivation. Here we introduce a Markovian master equation that is thermodynamically consistent and provides an accurate description whenever memory effects can be neglected. The thermodynamic consistency is obtained through a rescaled Hamiltonian for the thermodynamic bookkeeping, exploiting the fact that a Markovian description implies a limited resolution for heat. Our results enable a thermodynamically consistent description of a variety of systems where the secular approximation breaks down.


Introduction
Theoretical descriptions of open quantum systems are crucial for understanding various scenarios, as a complete shielding from the environment is usually not feasible or not even desirable. The latter is the case in the field of quantum thermodynamics, where heat flows in out-of-equilibrium situations are the key issue of concern. As a microscopic description of the environmental degrees of freedom becomes quickly intractable, numerous approaches exist to approximate the behavior of the degrees of freedom of the system alone [1][2][3][4]. Of particular interest are Markovian master equations in so-called GKLS form, after Gorini, Kosakowski, Sudarshan [5], and Linblad [6]. By neglecting any memory effects induced by the environment, they provide a particularly tractable description that gives access to the full density matrix of the system alone (i.e., the reduced system). These equations usually rely on Born-Markov approximations which do not ensure GKLS form. For this reason additional approximations are usually employed.
The most prominent approximation is the so-called secular approximation [1], where oscillating terms are dropped from the master equation (for a recent generalization to time-dependent systems, see Ref. [7]). This approximation has the desirable feature that it ensures consistency with the laws of thermodynamics. The secular approximation has also been termed the global approach, because it uses the (delocalized) eigenstates of the Hamiltonian. This is in contrast to the local approach which is widely used in, e.g., quantum optics [8,9] (for a comparison between the global and the local approach, see Refs. [10][11][12][13][14][15][16][17][18][19][20]). The difference between the local and the global approach becomes particularly apparent when considering a system of weakly coupled components (e.g., two or three qubits or harmonic oscillators weakly coupled to each other). In this case, the local approach can be obtained by deriving a master equation for the individual, uncoupled components and simply adding the coupling term. We stress however that the local approach is by no means phenomenological, as a microscopic derivation exists [13,20]. This approach is appealing mainly for two reasons: First, it does not require diagonalization of the Hamiltonian and can thus be applied to problems where a secular approximation is difficult to implement analytically. Second, it holds for systems that consist of weakly coupled degenerate units, where the secular approximation breaks down due to the near-degeneracies. Such systems are widely studied in the field of quantum thermodynamics as they are promising for manipulating and exploiting heat flows [21,22]. However, the local approach was criticized for not being thermodynamically consistent as it may result in violations of the second law of thermodynamics [23]. These violations where shown to be small as long as the local approach is justified [24,25]. One may thus argue that sizable violations of the second law provide a useful red flag, indicating that the approach is applied outside its regime of validity. Nevertheless, thermodynamic consistency is desirable in any Markovian master equation as it allows for falling back on well established laws [26]. It was shown that the local approach can be rendered thermodynamically consistent by re-defining heat. This was motivated from a collisional model, where work is required to maintain the collisions [27], as well as directly from the master equation itself [28]. Below, we show how thermodynamic consistency of the local approach can be obtained starting from the standard microscopic system-bath picture by exploiting a crucial insight: the approximations that result in a Markovian master equation impose limitations on the energy-resolution for heat. Within this limited resolution, we may re-define heat without compromising the accuracy of the approach.
Recently, a novel approach termed PERLind (Position and Energy Resolved LINDblad equation) was introduced [29] in an attempt to interpolate between the local and global approach. Since then, multiple microscopic derivations were given [30][31][32][33], showing that the PERLind approach does not require any additional assumptions going beyond the ones already present when performing Born-Markov approximations. Compared to the global and local approaches, the PERLind approach thus has an increased regime of validity, making it a very promising approach. However, it also comes with disadvantages: it requires diagonalization of the Hamiltonian and it is not thermodynamically consistent. Furthermore, an analytical treatment is often complicated by tedious expressions as shown below. Another approach that goes beyond the secular approximation is based on introducing a coarse-graining time [15,[34][35][36][37][38]. While taking this time to infinity recovers the global master equation, a GKLS master equation can be obtained by choosing a finite coarse-graining time. In general, coarsegrained master equations are not thermodynamically consistent. However, just as for the local approach, thermodynamic consistency can be recovered by re-defining heat as motivated by a collisional model [39]. Yet an alternative approach for obtaining a GKLS master equation beyond the secular approximation is provided by truncating the Redfield equation [40].
Here we introduce a novel Markovian master equation in GKLS form that goes beyond the secular approximation. Its main merit compared to previous approaches is that it is thermodynamically consistent. Our approach is based on the same principle that is at the heart of all Markovian master equations: the environment properties are slowly changing in energy. This allows us to not only neglect the broadening of energy levels but to also treat transition energies that are close as having the same value. Employing the same approximation in the definition of heat then naturally results in a thermodynamically consistent treatment. Our approach may reduce to the global approach (when no transition energies are close) or to a thermodynamically consistent version of the local approach (when all transition energies are close). We note that for a time-independent Hamiltonian, our results are in agreement with the unified GKLS master equation which was very recently derived in an independent work [41].
The rest of this paper is structured as follows. In Sec. 2, we compare different Markovian master equations and illustrate the main results without going into technical details. In particular, we introduce a thermodynamically consistent local master equation as an example of our general master equation, which is derived in Sec. 3. Sections 4 and 5 illustrate our master equation with the examples of a fermionic and a bosonic heat engine, covering both the time-independent as well as the time-dependent case. In Sec. 6, we apply our master equation to an interacting double quantum dot, a problem where we do not have an exact solution and where our master equation may provide a description that differs from both the global and the local approach. We conclude in Sec. 7.

Comparing different master equations
In this section, we compare different master equations, in order to illustrate some of our main results without going into mathematical details. To shed light on the discussion about the thermodynamic consistency, we focus on the system that was considered in Ref. [23], where it was shown that a local approach (here referred to as conventional local approach) may violate the second law of thermodynamics. The system (sketched in Fig. 1) consists of two coupled harmonic oscillators and is described by the Hamiltonian with the standard commutation relations The two oscillators are labeled by c and h, as they are coupled to a cold and a hot bath respectively. Diagonalizing the Hamiltonian results in the eigenfrequencies and the corresponding ladder operatorŝ where the angle θ is defined through where we use the branch 0 ≤ θ ≤ π, which corresponds to g > 0, which we tacitly assume througout this work. We will now compare five different approaches to describe the heat transport through this system: (i) The transmission function approach (valid for arbitrary system-bath coupling).
(ii) The conventional local approach (may violate the second law of thermodynamics).
(iii) Our local approach (ensures the laws of thermodynamics).
(iv) The global approach (relies on the secular approximation).
(v) The PERLind approach (interpolates between local and global). As we show below, the master equation we introduce in this work reduces either to the local, or to the global approach depending on the parameters. We stress that the approach we refer to as local is thermodynamically consistent, in contrast to the conventional local approach investigated, e.g., in Ref. [23]. The transmission function approach is used as a benchmark, because it is not perturbative in the system-bath coupling. Furthermore, by resolving the energy dependence of heat transport, it indicates specific shortcomings in the other approaches. We include the local approach used in Ref. [23] to illustrate when and why the second law of thermodynamics may be violated.
In all five approaches, we take the so-called wide-band limit, assuming an energyindependent coupling between system and bath and neglecting any Lamb shift of the Hamiltonian. While this is the only assumption for the transmission approach, all master equation approaches further rely on the Born-Markov approximations. For the current scenario, these approximations are valid whenever where κ α denotes the transition rate between system and reservoir α = c, h. All master equation approaches are of the form with different dissipators L j α , where j labels the approach.

The transmission function
For non-interacting particles, the heat current can be written using a Landauer-like formula [42]. For bosons (with vanishing chemical potential) it takes on the form where we introduced the Bose-Einstein distribution and the transmission function T (ω) can be obtained using non-equilibrium Green's functions [43,44]. For the present system, it reads (see for instance Ref. [45] for the analogous case of Fermions) The Landauer-like formula has an intuitive interpretation. At each energy ω, bosons traverse the system with the rate T (ω)dω. The transmission function has two peaks located at Ω + and Ω − . For small g and ∆, these peaks merge as illustrated in the insets of Fig. 2 (b).

The conventional local approach
In this approach, the dissipators act locally on the two oscillators where cl stands for conventional local. Here we used the GKLS superoperators with {·, ·} denoting the anti-commutator. The heat current in this approach is defined as J cl = Tr Ĥ S L cl hρ .
In steady state, this reduces to It is important to note that there is a microscopic derivation that results in this local approach, which is valid for the current system as long as n α B (Ω h ) n α B (Ω c ) [13], which is the case for g, ∆ Ω .
When this inequality is satisfied, we may expect that the laws of thermodynamics hold, see Fig. 2. When this approximation is not justified, the second law of thermodynamics is not guaranteed as illustrated in Fig. 3. The reason for this is that the Bose-Einstein distributions are evaluated at Ω h and Ω c respectively in Eq. (14). This implies that the bosons change their energy when traversing the system. From the Landauer-like formula [cf. Eq. (8)], it is apparent that this cannot happen. We note that while Eq. (15) ensures |n α B (Ω h ) − n α B (Ω c )| 1, the relative error in the heat current can still be sizable when the occupation numbers (and thus the heat current) become small. This is the case because Eq. (15) does not ensure |n α . Steady state heat current for degenerate harmonic oscillators (∆ = 0). Grey: Transmission approach (J t ) which serves as a benchmark. The PERLind approach (J p ) agrees perfectly with the transmission approach for these figures. Blue: local approach (J l ). In steady state, this approach reduces to the conventional local approach (J cl ) for ∆ = 0. Red: Global approach (J g ) which relies on the secular approximation. (a) As expected, the transmission approach interpolates between the local and global approaches as a function of the coupling strength g. In agreement with Eq. (15), the local approach may be justified even for g > κ. The insets illustrate the transmission function as a function of energy for two different values of g, where the local or the global approach is well justified. (b) The master equation approaches agree well with the transmission function approach for a wide range of temperatures. The inset (where T c = 0) illustrates the regime where both temperatures become small. We note that in this regime, the relative error associated to the master equation approaches may become large but the absolute error stays small. Parameters: κ ≡ κ c = κ h = 0.02Ω, ∆ = 0, k B T h =Ω, (a) k B T c = 0.8Ω, (b) g = 5 κ, inset: T c = 0.

Our local approach
Here we present a thermodynamically consistent local approach that follows from our general master equation introduced below in Sec. 3, where we motivate this approach microscopically and show how it can be generalized to more complicated systems. As in the previous approach, the dissipators act locally in this approach the only difference being that the Bose-Einstein distributions are now all evaluated at the frequencyΩ. This approach has the same regime of validity as the previous local approach. Indeed, Eq. (15) implies and we expect the dissipators in Eqs. (11) and (16) to result in approximately the same dynamics as illustrated in Fig. 2. If Eq. (15) is not satisfied, the two local approaches may result in different results, see Fig. 3.
To obtain a thermodynamically consistent description, the approximation in Eq. (15) is exploited in the definition of the heat current which reads where we introduced a separate Hamiltonian for the thermodynamic bookkeeping (TD) Whenever Eq. (15) holds, this Hamiltonian provides the (approximately) correct energy flows. Consider the case where ∆ = 0. Then,Ĥ TD is simply obtained by dropping the coupling between the oscillators. The thermodynamic bookkeeping thus neglects both the system-bath couplings, as well as the coupling between the oscillators. Of course, all these couplings are crucial for the dynamics, where they appear either in the Hamiltonian or in the dissipators. In Sec. 3, we illustrate how such a thermodynamic Hamiltonian can be constructed in general, providing thermodynamic consistency of the corresponding master equation at the price of a slightly reduced energy resolution of thermodynamic quantities. In steady state, the heat current reads We note that for ∆ = 0, we find J l = J cl in the steady state. However, in the transient regime the two local approaches may result in different heat currents even for ∆ = 0. The difference between the two approaches is in this case the energy associated to the coupling term (proportional to g), which is dropped in the local approach. We note that the heat current associated to this coupling term has been interpreted as a quantum contribution before [19].
Comparing the local master equation with the Landauer-like formula [cf. Eq. (8)], one can see that transmission is approximated to happen at a single energyΩ. This is a good approximation when Eq. (15) holds. In this case, the transmission function is only non-zero aroundΩ and the Bose-Einstein distributions may be assumed constant across all energies where transmission is non-zero. Interestingly, the local approach agrees well with the transmission approach for the present system, even when Eq. (15) is not satisfied, see Fig. 3 (b). This is in stark contrast to the conventional local approach.

The global approach
In the global approach, the dissipators are introducing jumps between the eigenstates of the Hamiltonian where the coupling strengths reflect the spatial distribution of the eigenstates κ + h = κ h cos 2 (θ/2), κ − h = κ h sin 2 (θ/2), κ + c = κ c sin 2 (θ/2), κ − c = κ c cos 2 (θ/2). (22) The heat current in the global approach reads  . Steady state heat current as a function of the detuning between the oscillator frequencies ∆ = Ω h − Ω c . Grey: Transmission approach (J t ) which serves as a benchmark. The PERLind approach (J p ) agrees perfectly with the transmission approach for these figures. Blue: local approach (J l ). Green: Conventional local approach (J cl ). Red: Global approach (J g ). (a) For small coupling g, only a very small violation of the second law (negative heat current) is observed in the conventional local approach. (b) For larger g, the violation of the second law becomes larger. Interestingly, the local approach provides an accurate heat current, even though it is no longer microscopically justified for these parameters. Parameters: which in steady state reduces to This approach is thermodynamically consistent and we will always find heat flowing from hot to cold. As it relies on the secular approximation, it requires the condition In the global master equation, the transmission is approximated to happen at the two energies Ω ± . If Eq. (25) holds, the transmission function given in Eq. (10) consists of two narrow peaks located at Ω ± , see the right inset in Fig. 2 (b). Over the width of these peaks, the Bose-Einstein distributions can be assumed constant and the global master equation is valid. Note that if Eq. (25) does not hold, the peaks overlap and the secular approximation breaks down, see Fig. 2. We also note that Eqs. (15) and (25) may both be valid. In this case, both the global as well as the local approach are expected to give accurate results.

The PERLind approach
In the PERLind approach, the dissipators read with the jump operatorŝ In this approach, the heat current is defined as For the present system, the heat current in the PERLind approach agrees perfectly with the result obtained from the transmission approach for all considered parameter values. This is expected as we consider scenarios where Eq. (6) is fulfilled and the Born-Markov approximations are justified. Clearly, the PERLind approach has strong advantages. However, it also has its disadvantages. First, the analytical expressions quickly become unwieldy, which is the case already for the simple system considered here. Second, the PERLind approach may violate the second law of thermodynamics [29]. When heat is defined by Eq. (28), a necessary and sufficient condition for the second law to hold is given by [46,47] L p α e −βαĤ S = 0.
This condition is not fulfilled as shown in Ref. [30]. However, as long as the Born-Markov approximations are justified, any second law violations should become vanishingly small.
It is an open question if a different definition for heat, that is consistent with the Born-Markov approximation, could salvage the second law.
2.6. Which approach to use?
As we have shown above, all approaches have their advantages and disadvantages.
The PERLind approach provides the most accurate description, while the conventional local approach is the most simple, not requiring diagonalization of the Hamiltonian. Both approaches may however result in violations of the second law. The global and the local approach both respect the laws of thermodynamics. Furthermore, together they provide an accurate description for all parameter values as they are justified in complementary, but overlapping, parameter regimes [cf. Eqs. (15), (25), and (6), assuming the Born-Markov approximations are justified]. In the following, we illustrate how these thermodynamically consistent approaches follow from a unified framework that can be extended to scenarios where neither a global nor a local description is justified.

A thermodynamically consistent master equation
We first revisit the general scenario under consideration and discuss the laws of thermodynamics. We then provide a detailed derivation and discuss how thermodynamic consistency is obtained by a consistent application of the approximations on the thermodynamic bookkeeping.

The general scenario
We consider the general scenario described by the Hamiltonian where the first term describes the Hamiltonian of the system (which may be timedependent) and the second and third term describe the thermal reservoirs (labeled by α) and their coupling to the system respectively. The system exchanges energy and particles with the reservoirs, such that energy changes can be divided into heat and work. The average heat that leaves bath α during the time interval [0, t] is given by whereN α denotes the particle number operator for reservoir α and µ α its chemical potential. The average work provided by reservoir α is In addition, the time-dependence of the system Hamiltonian results in the external average power

The laws of thermodynamics
Before deriving a Markovian description, we discuss the laws of thermodynamics as they hold for the general scenario.
3.2.1. The 0th law For a large environment in thermal equilibrium (i.e., described by a single inverse temperature β and chemical potential µ) and a time-independent system Hamiltonian, the reduced state of the system tends to [48,49] where Tr B denotes the trace over the reservoir degrees of freedom. In the weak systembath coupling limit, Eq. (34) reduces to the Gibbs statê

The 1st law
For future reference, we first introduce the heat current and power provided by bath α as The first law of thermodynamics is then given by In the weak system-bath coupling limit, the coupling energy can be neglected and U = Tr{Ĥ S (t)ρ tot (t)} reduces to the usual internal energy of the system.

The 2nd law
To express the second law of thermodynamics, we impose the initial conditionρ i.e., all reservoirs are in local thermal equilibrium and uncorrelated with the system and the other reservoirs. With this initial condition, the second law of thermodynamics can be written as [50] whereρ S = Tr B {ρ tot (t)} denotes the reduced state of the system, S[ρ 1 ||ρ 2 ] = Tr{ρ 1 (lnρ 1 − lnρ 2 )} is the quantum relative entropy (which is by definition positive for positive definite density matricesρ 1 ,ρ 2 with unity trace), and ∆S denotes the change in the system's von Neumann entropy Equation (39) has an intuitive interpretation: The entropy production Σ denotes the information that is lost when describing system and reservoirs by the stateρ S (t) ατ α . In this description, correlations between the system and the reservoirs, as well as any displacement from equilibrium of the reservoirs are neglected [51]. While Eq. (39) is always non-negative, the same is not necessarily true for the entropy production rateΣ A negative entropy production rate can be understood as information backflow and is a hallmark of non-Markovian behavior [52]. For systems amenable to a Markovian description, we expectΣ ≥ 0, as the initial condition in Eq. (38) can be translated in time without altering the dynamics.

Distribution for heat and work exchanged with the reservoirs
Throughout this work, we are mainly interested in average thermodynamic quantities. In order to obtain a Markovian description of these, it is nevertheless instructive to start with the full probability distribution for the heat and work exchanged with the reservoirs. With the initial condition given in Eq. (38), this distribution can be written as where we grouped Q α , W α , E α , and N α into vectors Q, W , E, and N and similarly for E α and N α . The joint probability for each bath α having Energy E α and particle number N α at time t = 0 is given by Tr e −βα(Ĥα−µαNα) .
The conditional probability that the reservoirs have energies E α and particle number N α at time t, given that their energies and particle numbers where E α and N α initially reads with T denoting the time-ordering operator. We note that the probability distribution given in (42) can in principle be measured by applying projective measurements on the reservoirs at the initial and final time. Importantly, while such measurements may not be experimentally feasible, they do not influence the dynamics of the system due to the chosen initial condition. The moment generating function is provided by the Fourier transform of the probability distribution in (42) where we introducedρ with the modified time-evolution operator The quantities λ α and χ α are known as counting fields and allow us to keep track of work and heat exchanged with the reservoirs (for the use of counting fields in master equations, see Refs. [4,53]). From the moment generating function, we recover the average values for heat and work given in Sec. 3.1 as We note that we do not include the external power in the probability distribution for heat and work as this goes beyond the scope of this paper. Indeed, fluctuations of the external power may not necessarily be described by a positive probability distribution without introducing a measurement scheme that potentially alters the dynamics of the system [54][55][56].

Born-Markov approximations
We are now in a position to derive a Markovian master equation, keeping track of the thermodynamic bookkeeping associated to the heat and work exchanged with the bath. Together with a secular approximation, resulting in a global master equation, such a procedure has been applied before [57][58][59] (for a similar approach based on path-integrals, see [60]). Here we will make an approximation that is different from the secular approximation, allowing for treating (near) degeneracies while still ensuring thermodynamic consistency. We follow standard procedure [1,4] and write the systembath coupling asV = α,kŜ where the operatorsŜ α,k (B α,k ) act only on the system (reservoir). We note that we do not assume that these operators are Hermitian. To determine the jump operators that will enter the Markovian master equation, we need the Fourier coefficients of the operatorsŜ α,k in the interaction picture, i.e., For a time-independent Hamiltonian, the frequencies ω j denote the energy gaps in the Hamiltonian and the operatorsŜ α,k;j are ladder operators that induce transitions between the corresponding states. For a periodic Hamiltonian with period t p = 2π/ , theŜ α,k;j are ladder operators of an averaged Hamiltonian defined bŷ The frequencies then fulfill ω j = ν j + l j , where ν j denotes an energy gap ofĤ av and l j denotes an integer, corresponding to the exchange of photons with the driving field [57,61,62].
Standard application of Born-Markov approximations then results in the Redfield equation including counting fields (in the interaction picture) where we suppressed the counting field dependence of the density matrix for ease of notation. We note that the term Redfield equation sometimes refers to the non-Markovian equation that is obtained before taking the time-integration to infinity [1] while we include this limit here. Here we used global particle conservation, which ensures for all pairs k and k that appear together in Eq. (53), i.e.,Ŝ α,k;j andŜ α,k ;j change the particle number by the same amount, such that no superpositions of particle numbers are created. We further introduced the bath correlation functions These correlation functions define a bath correlation time τ B by their characteristic decay time. The Markov approximation is generally justified when the bath correlation time is much shorter than the characteristic time-scale over whichρ S changes in the interaction picture, τ S . The time τ S describes the relaxation time of the system and is determined by the inverse of the system-bath coupling. Note however that the counting fields λ α enter the argument of the bath correlation times. For the counting-field dependent density matrix, the Markov approximation is only justified if C α k,k (±τ + λ α ) 0 for τ τ S . This implies the regime of validity The last inequality has a very important consequence. It implies that only the low frequency components of the heat distribution are to be trusted. Due to the uncertainty principle between Fourier conjugate variables [63], this implies that energy-differences in heat of the order of 1/τ S cannot be resolved. Thus, whenever a Markovian description is employed, the heat exchanged with the reservoirs suffers from a limited energyresolution. With this in mind, it is no surprise that Markovian descriptions may result in thermodynamic inconsistencies. There is however a straightforward solution to the problem: as our resolution of heat is finite, we may change the definition of heat such that it fulfills the laws of thermodynamics while remaining the same within our limited resolution.

Frequency grouping for positivity
It is well known that the Redfield equation does not preserve positivity of the density matrix and can thus result in negative probabilities. Multiple schemes have been put forward to achieve positivity. Here we introduce a novel scheme that has the same regime of validity as the Born-Markov assumptions and thus goes beyond the regime of validity of the secular approximation. As we show in the next subsection, our approach allows for a thermodynamically consistent formulation. Equation (56) ensures that for any two transition frequencies we either have Indeed, both of these conditions may be fulfilled simultaneously as τ B τ S . We may thus group the transition frequencies into sets x q , such that the first (second) inequality holds if the transition frequencies are in the same (different) set, that is It is important to note that this procedure may not always work, for instance if the ω j form a continuum. However, in small quantum systems, where the number of transition frequencies is finite, this procedure is expected to work. It is particularly well suited for thermal machines that consist of weakly coupled sub-units. These naturally exhibit sets of near-degenerate transition frequencies (see the examples below).
In the spirit of the secular approximation, we then drop terms in Eq. (53) where ω j and ω j belong to different sets. This is justified as the corresponding terms exhibit fast oscillations that average to zero. For transition frequencies in the same set x q , we follow the spirit of the Markov approximation and replace within I(s) in Eq. (53), while keeping the frequency differences in the prefactor e i(ω j −ω j )t governing the coherent dynamics of the system. For each set, we thus choose a set frequency ω q . All transition frequencies ω j in the set x q are then replaced by the set frequency, because they are virtually indistinguishable over the time-scale over which the integrand in Eq. (53) is finite. This approximation results in the Lindblad master equation with L χα,λα αρ = k;q Γ α k (ω q ) e −iλαωq−i(χα−λα)µαn α,kŜ α,k;q (t)ρŜ † α,k;q (t)− 1 2 where the tilde denotes the interaction picture and we introduced the jump operatorŝ and the Lamb-shift Hamiltonian as well as the quantities For simplicity, we assumed C α k,k ∝ δ k,k to derive Eq. (59). We note that relaxing this assumption is straightforward.
In the absence of counting fields, we may use the Kubo-Martin-Schwinger condition [1] (which is slightly complicated by our definition of the bath-correlation functions) to write the master equation as For a time-independent Hamiltonian, the master equation in the Schrödinger picture is given by with whereŜ α,k;q ≡Ŝ α,k;q (0) andĤ LS ≡Ĥ LS (0). For time-dependent Hamiltonians, one has to be more careful because in generalŜ α,k;q (t) =Û † S (t)Ŝ α,k;qÛS (t), see also Appendix A. There are two simple limits for the master equation in Eq. (66). If all pairs of transition frequencies fulfill |ω j − ω j | 1/τ S , then each transition frequency may be associated to a separate set x j . We then recover the secular approximation where D[Ŝ α,k;q (t)] = D[Ŝ α,k;j ] (the time-dependence of the jump operators drops out in the interaction picture). If all pairs of transition frequencies fulfill |ω j − ω j | 1/τ B , then all transition frequencies may be grouped into a single set x 0 such thatŜ α,k;q (t) = U † S (t)Ŝ α,kÛS (t) (the time-dependence of the jump operators drops out upon returning to the Schrödinger picture). We then recover a local master equation, which has the appeal that the Hamiltonian does not need to be diagonalized in order to identify the jump operators. This local approach differs from the conventional local approach as only a single frequency enters the bath distribution function.
Importantly, the counting fields λ α have an influence on the regime of validity of Eq. (59), just as for the Markov approximation. Indeed, for the approximation to be valid, |ω q − ω j | does not only have to be much smaller than 1/τ B but also has to be much smaller than 1/|λ α | [in complete analogy to Eq. (56)]. Just like the Markov approximation, this frequency-grouping thus limits the energy-resolution for heat exchanged with the reservoirs. As a consequence, energy changes of the order of |ω j − ω j |, where both frequencies are within the same set, can no longer be resolved. This finite resolution in energy may result in the false conclusion that particles change their energy when traversing the system. In the conventional local approach, this results in thermodynamic inconsistencies when using standard definitions for heat, see Sec. 2.
We note that dropping terms that involve frequencies from different sets is reminiscent of the partial secular approximation developed in Refs. [15,34,35,64,65].
Here, as well as in Ref. [41], an additional coarse-graining in energy (replacing frequencies with set-frequencies) results in a master equation in GKLS form.

Two Hamiltonians for thermodynamic consistency
We now turn to the question of how to utilize the finite energy-resolution to ensure thermodynamic consistency. To do this, we introduce a second Hamiltonian,Ĥ TD , that provides the correct thermodynamic bookkeeping under the approximations that resulted in our master equation. For simplicity, we consider here a time-independent system Hamiltonian. The time-dependent case is slightly more complicated and treated in Appendix A. Due to the frequency grouping outlined in the last subsection, all frequencies in the set x q are effectively replaced by the frequency ω q from the point of view of the reservoirs. To ensure a consistent thermodynamic bookkeeping, the HamiltonianĤ TD needs to fulfill [Ŝ α,k;j ,Ĥ TD ] = ω qŜα,k;j , for all frequencies ω j ∈ x q . To fulfill this, the thermodynamic Hamiltonian can be obtained fromĤ S by changing its eigenvalues such that all frequencies ω j → ω q for ω j ∈ x q . Such a rescaling is expected to always be possible if the frequency grouping is possible as discussed above.
With the thermodynamic Hamiltonian at hand, we define the internal energy of the system as U = Tr{Ĥ TDρ }.
Furthermore, from Eqs. (59), (36), and (49) we find that the heat current and power provided by bath α can be cast into where L α is defined in Eq. (67). We note that whileĤ TD determines the thermodynamic bookkeeping, it is stillĤ S that determines the kinetics, i.e., enters the master equation in Eq. (59). For future reference, we also introduce the total power produced by the system as This will be the quantity of interest when we consider heat engines.
3.6.1. The 0th law Using Eq. (68), it is straightforward to show that as well as [Ĥ TD ,Ĥ LS ] = 0. Furthermore, asĤ TD is obtained by changing only the eigenvalues ofĤ S , these two Hamiltonians also commute. When all reservoirs have the same inverse temperature β and chemical potential µ, it then follows that the Gibbs state with respect to the thermodynamic Hamiltonian, is the steady state of Eq. (66). Compared to Eq. (35), this state neglects the systembath coupling, as well as the differences betweenĤ S andĤ TD . This is consistent with our approximations, which imply that these differences cannot be resolved within our Markovian treatment.
3.6.2. The first law Using [Ĥ TD ,Ĥ S +Ĥ LS ] = 0, the first law of thermodynamics follows directly from Eqs. (69) and (70) and reads 3.6.3. The second law The entropy production rate can be written aṡ Here we used Eqs. (66), (70), and the last inequality is known as Spohns inequality [46,66] and relies onρ G (β α , µ α ) being a fixed point of L α , cf. Eq. (72). The master equation in Eq. (66), together with the thermodynamic bookkeeping introduced in Eq. (69) and (70), thus provide a thermodynamically consistent description. We stress that the only assumption that went into Eq. (66) is the one that justifies the standard Born-Markov approximations, i.e., τ B τ S . However, when considering the full probability distribution for heat and work exchanged with the reservoirs, cf. (59), then we further have restrictions on the counting fields λ α which are the Fourier transform variables of heat. These restrictions imply that we lose energyresolution on the scale of 1/τ S as well as |ω j −ω j | for pairs of frequencies that are close to each other (and grouped into one set). Thermodynamic consistency is obtained by using an appropriate definition of heat that is consistent with our limited energy-resolution (i.e., we may freely add terms of order |ω j − ω j | to heat, as our results are not to be trusted to this order in the first place).

System
In this section, we illustrate the derivation of our master equation for a simple but nontrivial example of non-interacting electrons. Here, as well as for all time-independent examples, we provide all equations in the Schrödinger picture. The system under consideration is sketched in Fig. 4 and consists of two coupled single-level quantum dotsĤ with the standard anti-commutation relations In Eq. (76), Ω α denote the on-site energies and g the coupling strength. We note that this system is a fermionic version of the system considered in Sec. 2 and thus shares many of its properties. Due to the finite chemical potential of the reservoirs, it can however feature as a heat engine, leveraging a heat current to drive a particle current against a chemical potential bias. The reservoirs and their coupling to the system are described byĤ with α = L, R. From Eqs. (50) and (78), we identify the following system and bath operatorsŜ To derive the master equation and judge its validity, the bath correlation functions need to be inspected. We defer this discussion to Appendix B.1, and simply state the conclusion that the Born Markov approximations are valid whenever where the inequality should hold for all choices of α, β, and j and we assumed an energy-independent bath spectral density. Physically, Eq. (81) ensures that the rates with the Fermi-Dirac distribution is flat over the energy scale κ α . At high temperatures, this is ensured because the Fermi-Dirac distribution becomes flat. For large |ω j − µ β |, the eigenenergies of the double quantum dot lie far away from the chemical potential, where the Fermi-Dirac distribution is flat and takes on the value zero or one.
To identify the transition frequencies and the jump operators, we require the Fourier coefficients of the system operators given in Eq. (79) where we used a similar notation to Sec. 2, i.e., the diagonlized Hamiltonian readŝ where Ω ± =Ω ± ∆ 2 + g 2 , withΩ = (Ω R + Ω L )/2, ∆ = (Ω R − Ω L )/2, and cos(θ) = ∆/ ∆ 2 + g 2 . Comparing Eqs. (84), and their Hermitian conjugates, to Eq. (51), we find the transition frequencies and the corresponding jump operators We note that the frequencies Ω ± only feature in the jump operators with k = −1 (corresponding to electrons leaving the system), while the frequencies −Ω ± only feature in the jump operators with k = 1 (corresponding to electrons entering the system). This implies that we may consider only the frequencies Ω ± when deciding how to group the transition frequencies into sets. We may group Ω + and Ω − in the same set, or we may assign them to different sets. For the frequencies −Ω ± , we then choose an equivalent grouping.

Global approach
Grouping Ω ± into different sets is justified as long as their difference 2 ∆ 2 + g 2 is much larger than max{κ L , κ R }, which corresponds to 1/τ S . This grouping results in a different set for each frequency, such that {ω q } = {ω j } and {Ŝ α,k;q } = {Ŝ α,k;j }. Having identified the set frequencies and jump operators, we may then use Eq. (67) to obtain the well-known dissipator in the secular approximation with the coupling strengths As the jump operators in Eq. (87) are ladder operators ofĤ S , no rescaling is required for obtaining the thermodynamic Hamiltonian and we findĤ TD =Ĥ S , resulting in the standard definition for heat currents [cf. Eq. (70)].

Local approach
Alternatively, we may group Ω + and Ω − in the same set. In this case, we need to choose a set frequency. The average valueΩ is a natural choice. Using the same arguments that result in Eq. (81), this grouping is justified as long as 2 ∆ 2 + g 2 max{k B T β , |ω j − µ β |}. In this case, we end up with two sets with the set frequencies and the jump operators which are simply obtained by adding the jump operators in Eq. (87) which belong to the same set. Inserting these quantities into Eq. (60) results in the local approach with the dissipator for α = L, R. The thermodynamic Hamiltonian is obtained by rescaling all transition frequencies to their set frequency. In the present case this rescaling is obtained by Ω ± →Ω resulting inĤ  1 κ), the Born-Markov approximation breaks down for on-site energies close to the chemical potential µ R . In this regime, the PERLind approach differs from the transmission approach. Parameters:

Results
Depending on how the transition frequencies are grouped, we obtain either the global or the local approach. As long as the Born-Markov approximations are justified, i.e., as long as Eq. (81) holds, at least one of the two procedures is justified. As the parameters of the system are changed, we may thus need to change from the local to the global approach to maintain an adequate description of the system. This is a general drawback of our master equation: as parameters are varied, the grouping of transition frequencies into sets may need adaptation. This results in discontinuities which are expected to be small as long as the Born-Markov approximations are justified.
In Fig. 5, we illustrate the global, the local, as well as the benchmark provided by the transmission approach for this system when operated as a heat engine (see Appendix B.2 for analytical expressions). We find that our master equation reproduces the transmission approach well if the frequency grouping is chosen such that the local (global) approach is obtained for g below (above) 5κ. Interestingly, we find that the efficiency is generally better reproduced by the global approach (see insets in Fig. 5). While the local approach predicts Carnot efficiency when n L F (Ω) = n R F (Ω), the transmission and the global approach predict a drop in efficiency. The reason for this drop is that the heat current remains finite at vanishing power when transmission occurs at more than one energy [21].
The upper inset in Fig. 5 (b) shows that at temperatures k B T α κ α , the Born-Markov approximation breaks down if one of the transition frequencies Ω ± is close to the chemical potential µ α , cf. Eq. (81). In this case, the step-like behavior of the Fermi-Dirac distribution renders the assumption of a Markovian bath unjustified. By approximating transmission to occur only at the transition frequencies of the system, the Markovian master equations predict a step-like behavior of heat currents and power. The transmission approach shows that the step-like feature is smeared out by the transmission function, which is a continuous function of energy.

System
In this section, we provide an example of a time-dependent system. To this end, we consider the heat engine introduced by Kosloff in 1984 [67] with the system Hamiltonian where the bosonic annihilation and creation operators fulfill the standard commutation relations given in Eq. (2), and the frequency of the external drive reads This Hamiltonian describes two bosonic modes with different frequencies that are coupled via a time-dependent term with a detuning quantified by ∆. For a physical implementation of this Hamiltonian based on a superconducting circuit, see Ref. [68]. The reservoirs and their coupling to the system are described bŷ with α = c, h and we consider a vanishing chemical potential for the reservoirs, µ α = 0. From Eqs. (50) and (96), we identify the following system and bath operatorŝ We find that the Born-Markov approximations are justified when (see Appendix C.1 and Refs. [12,13]) where the inequality should hold for all values of α and j. This ensures that the Bose-Einstein distribution is flat around the transition energies over the energy scale κ α . Assuming a flat bath spectral density, the same holds for the rates where ω > 0 as we are dealing with Bosons with zero chemical potential. The time-dependence in Eq. (94) can be removed by a suitable unitary transformation, see Ref. [13], where this model was used to compare local and global master equations. Here, we use the model to illustrate how a time-dependent Hamiltonian affects our master equation and we thus remain in the lab frame. As discussed in Sec. 3, the average Hamiltonian defined in Eq. (52) is an important quantity. For our system, we find (see Appendix C.2) where the eigenmodes are given in Eqs. (4,5) and we introduced the frequencies where θ is defined in Eq. (5) [where Eq. (95) determines ∆]. From these equations (and their Hermitian conjugates) we can identify the transition rates and jump operators. As we have seen for the fermionic heat engine, it is sufficient to consider only the frequencies and operators related to particles leaving the system (the others can be treated separately and equivalently). These transition frequencies read As discussed above, these can be written as ν j + l , where ν j denotes a transition frequency of the average Hamiltonian [i.e. Ω ± c for Eq. (101)] and l = 0, 1. The corresponding jump operators read (105)

Global approach
Again, we have two choices for grouping the frequencies which result in the local and global approach respectively. We may group Ω ± α into different sets or assign them to the same set. Grouping them into different sets results in a single frequency per set such that {ω q } = {ω j } and {Ŝ α,k;q (t)} = {Ŝ α,k;j exp(−iω j t)} [cf. Eq. (61)]. Inserting these quantities into Eq. (65), we find the dissipator (in the interaction picture) whith κ σ α being defined in Eq. (22). As outlined in Appendix A, a single-thermodynamic Hamiltonian is usually not sufficient for the time-dependent scenario. Furthermore, even when a secular approximation is performed, a rescaling needs to be performed [62]. In the present case, the rescaling is obtained starting fromĤ av in Eq. (101) and rescaling the relevant transition frequencies to ω q . It turns out that it is sufficient to consider two thermodynamic Hamiltonians, one for each reservoir The expression for the heat current given in Eq. (A.5) then reduces to Note that the two thermodynamic Hamiltonians correspond to two possible choices for H av . In this scenario, the external power can no longer be accessed by the standard expression given in Eq. (33). As a consequence of the secular approximation, this quantity evaluates to zero [13]. The time-averaged power can however be recovered by relying on the first law as in Eq. (A.6).

Local approach
Alternatively, we may group Ω ± α into the same set. This grouping is justified as long as g, ∆ Ω c , Ω h [in analogy to Eq. (99)]. A natural choice for the set frequencies then reads with the jump operators Upon returning to the Schrödinger picture, these quantities result in the local dissipator Two thermodynamic Hamiltonians can again be obtained by re-scaling the transition frequencies of the average Hamiltonian in Eq. (101). However, due to the local structure of the jump operators, we can find a single thermodynamic Hamiltonian in the Schrödinger picture [cf. Appendix A] From Eq. (A.9), we recover the usual expression for the power produced by the system  [68]). Along the arrow, Ω h is increased. The vertical line denotes the Carnot efficiency: 1 − T c /T h . For these plots, the PERLind approach is visibly indistinguishable from the transmission approach. Parameters: For ∆ = 0, the thermodynamic Hamiltonian in Eq. (112) was used before in Ref. [56] and has an intuitive explanation: For the thermodynamic bookkeeping, we neglect the coupling energy between the bosonic modes, just as we neglect the coupling energy between system and bath.

Results
In Fig. 6, power, heat current, and efficiency of the bosonic heat engine for a resonant drive are illustrated. The local and global approaches for this scenario were already compared in Ref. [13], where exact numerics for finite reservoirs was used as a benchmark. In agreement with these results, we find that the local and global approaches together reproduce exact solutions over the full range of parameters where the Born-Markov approximation is justified. This implies that our master equation, which reduces either to the local or to the global approach, describes this scenario very well. We note that the PERLind approach reproduces the transmission approach extremely well. The disadvantage of the PERLind approach is that it is not thermodynamically consistent. The disadvantage of our master equation is that it exhibits a discontinuity when the frequency grouping is adapted (i.e., when changing from the local to the global approach).
In Fig. 7, the performance of the heat engine at finite detuning ∆ is investigated. For small detuning, we find that the global approach correctly reproduces the g → 0 limit, as an exact degeneracy is no longer reached in this limit [panel (a)]. As ∆ increases, the low-g behavior is increasingly well captured by the global approach, because the secular approximation becomes better. For a fixed coupling strength g, we find a decrease in Bosonic heat engine for a a finite detuning ∆. Steady state heat current from the hot bath J j h and power P j S as a function of (a) the coupling strength g, and (b) the detuning ∆. The superscripts l, g, t, refer to the local, global, and transmission approach respectively. For these plots, the PERLind approach is visibly indistinguishable from the transmission approach. Parameters: the heat engine performance as ∆ is increased [panel (b)]. This behavior is expected, because energy transfer between the drive and the system works best on resonance.
6. Interacting double quantum dot -beyond local and global approaches

System
In the previous examples, the derived master equation reduces to the well known global or local approaches. In this last example, we consider a system that has more transition frequencies, such that the master equation may differ from both the local as well as the global approach. To this end, we consider a spinless double quantum dot with Coulomb interactions described by the Hamiltonian Apart from the interaction term, the system as well as the bath properties are identical to Sec. 4 (with Ω ≡ Ω L = Ω R and Ω ± = Ω ± g) where all notation used in this section is defined. From Eq. (114), we see that the Hamiltonian is diagonal in the occupation basis of the ± modes. The transition frequencies and jump operators are obtained from the equations Considering only transition frequencies that correspond to particles leaving the system, we find and the jump operators (117)

Global approach
In the global approach, we choose a different set for each transition frequency. This results in the dissipator whereσ = σ. As usual, the thermodynamic Hamiltonian in the global approach is equal toĤ S . This frequency grouping is justified when for all j, j , and α.

Local approach
In the local approach, we group all transition frequencies into a single set. A natural set frequency is then the average frequency ω q = Ω + U/2. The jump operators are given by summing the jump operators in Eq. (117) over j, which results in the local dissipator The thermodynamic Hamiltonian, obtained by rescaling all transitions to Ω + U/2, is given byĤ We note that apart from the choice of the transition frequency, the dissipator and thermodynamic Hamiltonian are identical to the one obtained for U = 0, c.f., Eqs. (92) and (93). The frequency grouping for the local approach is justified as long as for all j, j , and α.

Semi-local approach
In addition to the local and global approaches, we consider the low-g frequency grouping, where Ω ± can be replaced by Ω. This results in the set frequencies with the corresponding jump operators Note that the jump operators still locally change the number of electrons. However, the jumps are now dependent on the occupancy of the other dot due to the Coulomb interaction. Inserting these quantities into Eq. (60) results in the dissipator whereᾱ = α. The thermodynamic Hamiltonian corresponding to this frequency grouping readsĤ This frequency grouping is justified for To ensure that the frequencies in different sets obey |ω j − ω j | κ α , one may expect the additional condition U g, κ α . However, because no coherences can build up between states with a different total number of electrons in the system, this additional condition is not required. The semi-local approach thus enjoys a strictly larger regime of validity than the local approach [c.f. Eq. (122)], remaining valid even for a vanishing interaction strength U = 0.

Results
Due to the interaction term, the transmission approach no longer applies. As a benchmark, we instead use the PERLind approach, see Appendix D. As mentioned above, this approach is expected to provide accurate results whenever the Born-Markov approximations are justified, an expectation that was confirmed in the examples we considered above. Figure 8 illustrates heat current and power in the interacting double quantum dot. We find that at finite interaction strengths, the semi-local approach should be employed instead of the local approach for small coupling g. In particular, comparing Fig. 8 (a) with Fig. 5 (a) we find a very similar interplay between the semi-local and global approach at U = 0 as we observed between the local and global approaches for U = 0. We furthermore find that the semi-local approach agrees with the PERLind approach for any value of the interaction strength, as long as the coupling g is sufficiently small to respect Eq. (127), see Fig. 8 (b). These results highlight the importance to go beyond the local and global approaches for systems that have more than two competing transition frequencies.

Conclusions
Markovian master equations in GKLS form are approximate descriptions for the reduced system state. The approximations involved in deriving these equations may not preserve the laws of thermodynamics (as is the case, e.g., in the PERLind approach). A comparison to the transmission approach for non-interacting particles illustrates that the source of any thermodynamic inconsistency is an inconsistent assignment of heat to the different jump operators appearing in the master equation. To shed light on this issue, we performed a microscopic derivation of the probability distribution for heat, employing the same approximations that are used for deriving master equations in GKLS form. We found that, when employing master equations, the resolution in heat is limited by the approximations that are performed. Exploiting this limited resolution, we derived a thermodynamically consistent master equation. To this end, we adapted the thermodynamic bookkeeping to enforce a consistent assignment of heat to the jump operators. The changes in heat induced by this procedure remain smaller than the resolution that the master equation allows for and are thus consistent with the performed approximations. This refined thermodynamic bookkeeping is captured by a thermodynamic Hamiltonian, which may differ from the Hamiltonian that determines the dynamics of the quantum state.
We illustrated our master equation with three different examples, including a timedependent model and a model that includes interactions. Our master equation may reduce to the well-known global approach, or to a thermodynamically consistent version of the local approach in their respective regime of validity. For systems where neither of these approaches is adequate, it provides a novel description.
Our master equation provides a thermodynamically consistent description for quantum systems that are amenable to a Markovian description. This allows future investigations on the thermodynamics of quantum devices to fully rely on the conclusions that follow from basic thermodynamics. larger than the bath correlation time τ B . In this case, a master equation may be derived by considering the Hamiltonian to be frozen, treating its time argument as any other parameter [69]. We may still use Eq. (66) but the transition frequencies ω j and ω q , as well as the operatorsŜ α,k;j become time-dependent. The thermodynamic Hamiltonian, obtained as in the time-independent case outlined in Sec. 3.6, then also becomes timedependent and the first law reads where the internal energy as well as the heat-current and power from reservoir α are still defined by Eqs. (69) and (70) respectively. The first term on the right-hand side of Eq. (A.1) denotes the power provided by the external drive. The second law still holds as discussed in Sec. 3.6. Next, we consider a time-periodic system Hamiltonian with period t p = 2π/ . In this case, the frequencies ω j are not directly transition frequencies but rather fulfill where l is an integer and ν j are the transition frequencies of the averaged Hamiltonian H av defined in Eq. (52). For clarity, we extended the index of ω j,l to explicitly include l, such that j uniquely defines the transition frequency of the averaged Hamiltonian. We stress that in this scenario, Eqs. (59) and (64) hold. Throughout this appendix, we remain in the interaction picture, because going to the Schrödinger picture is nontrivial for time-dependent Hamiltonians [cf. the comment below Eq. (67)]. Since any time-dependence may be thought of as a single-period of a periodic process, our master equation is applicable for any time-dependence. However, for a thermodynamically consistent description, we require the additional assumption which is complementary to the slow driving regime discussed above. As in the timeindependent case, we may group the frequencies into sets x q , such that |ω j,l − ω j ,l | 1/τ B for ω j,l ∈ x q , ω j ,l ∈ x q with q = q , |ω j,l − ω j ,l | 1/τ S for ω j,l ∈ x q , ω j ,l ∈ x q with q = q .
We note that due to Eq. (A.3), we may choose the sets x q such that ω j,l and ω j,l are not in the same set for l = l . In that case, all frequencies in a given set have different values for j, i.e., correspond to different transition frequencies ν j of the averaged Hamiltonian. For a given set x q , a thermodynamic HamiltonianĤ q TD may be obtained by rescaling ν j → ω q for all j that occur in the set x q . In contrast to the time-independent case, the same transition frequency ν j may feature in different sets (as it may feature in different ω j,l ). The rescaling therefore has to be done individually for each set. We note that this rescaling is even required in the secular approximation since ν j and ω j,l differ by l [62]. The heat current from bath α, obtained from Eqs. (59) and (36) can then be cast into Tr{(Ĥ q TD − µ αNS )L α;qρS }, (A.5) withL α;q given in Eq. (65).
To arrive at the master equation in Eq. (66), we dropped terms that oscillate with frequencies much larger than 1/τ S . Unfortunately, these terms may be important when evaluating the standard expression for external power given in Eq. (33). Loosely speaking, the oscillations of the density matrix cancel oscillations of ∂ tĤS (t), thereby contributing to the average power. In general, one may no longer disentangle the internal energy from the external power [70]. We may then still consider the first law upon taking a long-time average, ensuring that changes in the internal energy become negligible. Similarly, one could average over a single period when the system has reached a limit cycle. The first law then reduces to where the bar denotes any kind of average that ensures a vanishing change in internal energy. Interestingly, one may still access the internal energy (and thus the time-dependent power) when the master equation takes on a local structure. To illustrate this, let us consider a coupling Hamiltonian of the form (A.7) A local master equation is obtained when we group the transition frequencies such that each set frequency corresponds to a single reservoir, i.e., {ω q } = {ω α }. In this case, the jump operators (in the interaction picture and dropping the index q) are of the formŜ α (t) =Û † S (t)Ŝ αÛS (t) [cf. Eq. (61)]. For each jump operator, we then introduce a corresponding number operatorn α , such that [Ŝ α ,n α ] =Ŝ α . In this case, we may define a single, albeit time-dependent, thermodynamic Hamiltonian (in the interaction picture)Ĥ which fulfills the defining property [Ŝ α (t),Ĥ TD ] = ω αŜα (t). One may then still use Eq. (69) to define the internal energy and the first law reads which is valid both in the interaction, as well as in the Schrödinger picture where both the jump operators as well as the thermodynamic Hamiltonian become time-independent.
The first term on the right-hand side of Eq. (A.9) can be identified with the power provided by the external drive.
To verify the second law, we write the entropy production rate aṡ Spohns inequality ensures the positivity of the entropy production rate sinceL α;qρ q G = 0.
where we introduced the bath spectral density κ α 2π = l t 2 α,l δ(ω − ε α,l ), (B.2) quantifying the system-bath coupling and which we assume to be constant as a function of ω. We note that C α k ≡ C α k,k and C α k,k ∝ δ k,k . The bath correlation functions decay on the time-scale β α . The bath correlation time can thus be identified by τ B = max{β L , β R }. The system (in the interaction picture) changes on the time-scale τ S = min{1/κ L , 1/κ R }. Thus, the Born-Markov approximations are expected to hold for temperatures that are much higher then the coupling between system and bath. Furthermore, from the Redfield equation in Eq. (53), together with Eq. (B.1), we find that the integral over s involves factors exp[is(ω j ±µ α )]. Whenever |ω j ±µ α | k B T α , the rapid oscillations of these factors imply that only values of s 1/(ω j ± µ α ) contribute to the integral. For larger values of s, the term k B T α / sinh(πsk B T α ) varies on a time-scale much slower than the oscillations. This implies that the Born-Markov approximations are also justified in the regime |ω j ± µ α | k B T α and min{|ω j ± µ L |, |ω j ± µ R |} 1/τ S . Now in case 1/τ S |ω j ± µ| (dropping the bath index for ease of notation), either we have k B T |ω j ± µ|, such that the Born-Markov approximations are justified, or we have k B T |ω j ± µ| which implies 1/τ S k B T , also ensuring the Born-Markov approximation. In the first case, the bath-correlation functions oscillate to zero over times much smaller than τ S , in the latter case the bath correlation functions decay much faster than τ S . As a consequence, the Born-Markov apprixmations are fulfilled whenever κ is much smaller than either k B T , or |ω j ± µ| as expressed in Eq. (81).
The Fourier transforms of the bath correlation functions are given by Eq. (82). We reach the same conclusions about the validity of the Born-Markov approximation by demanding that these functions are flat over the energy scale κ (the inverse of τ S ). For κ k B T , the Fermi-Dirac distribution is everywhere slowly varying compared to κ. For κ |ω − µ| and k B T |ω − µ|, the Fermi function is either very close to zero or one, depending on the sign of ω − µ, and remains so over the energy scale of κ.

Appendix B.2. Power and heat currents
For the double dot system with non-interacting fermions the Landauer-like formula provides the following steady state heat current (out of the left lead) and power with the transmission function T (ω) = g 2 κ L κ R |(ω − Ω L + i κ L 2 )(ω − Ω R + i κ R 2 ) − g 2 | 2 .

(B.5)
Note that the transmission function is the same as for the bosonic system, cf. Eq. (10).
For the global and local approaches, the definitions for the heat current and the power are given in Eq. (70) and (71), withN S =d † Ld L +d † Rd R . For the local approach, we find and This tight-coupling condition between heat current and power [which can be inferred from the form ofĤ TD , cf. Eq. (93)] is a direct result of approximating transmission to happen at a single energy.
For the global approach, we find and with the jump operatorŝ (B.11) In this approach, heat current and power are defined as where κ α is given by Eq. (B.2) and we have C α k ≡ C α k,k and C α k,k ∝ δ k,k . We note that these integrals only converge if κ α → 0 for ω → 0. Here we do not give a specific form of κ α and directly proceed to investigating the Fourier transforms of the bath correlation functions. For a detailed discussion on the bath correlation functions for an Ohmic spectral density, we refer to Ref. [12]. By Fourier transforming Eq. (C.1), we find the rates given in Eq. (100). We note that even if κ α needs to vanish as ω → 0, we can still assume it to be constant over all energies that the system probes (i.e., all energies where the transmission function is non-zero).
Similarly to the fermionic case, the Born-Markov approximations are justified when the quantities in Eq. (100) remain constant over the energy interval κ α around the transition frequencies of the system. This is the case if κ α ω j . We note that for very small temperatures (k B T ω j ), n B (ω j ) becomes exponentially small and the dynamics is dominated by the temperature-independent term in Eq. (100) which corresponds to spontaneous emission. We further note that κ α k B T α is not a sufficient (C.9) For the local approach, we find the heat current (C.10) and power The global approach results in the heat current ], (C.12) and power Finally, for the PERLind approach, the dissipators read (in the rotating frame) (C. 15) In this approach, the heat currents can be obtained from Eq. (53) by applying the approximations discussed in the supplemental material of Ref. [30]. For the present system, we obtain with the jump operatorŝ