Markovian master equations for quantum thermal machines: local vs global approach

The study of quantum thermal machines, and more generally of open quantum systems, often relies on master equations. Two approaches are mainly followed. On the one hand, there is the widely used, but often criticized, local approach, where machine sub-systems locally couple to thermal baths. On the other hand, in the more established global approach, thermal baths couple to global degrees of freedom of the machine. There has been debate as to which of these two conceptually different approaches should be used in situations out of thermal equilibrium. Here we compare the local and global approaches against an exact solution for a particular class of thermal machines. We consider thermodynamically relevant observables, such as heat currents, as well as the quantum state of the machine. Our results show that the use of a local master equation is generally well justified. In particular, for weak inter-system coupling, the local approach agrees with the exact solution, whereas the global approach fails for non-equilibrium situations. For intermediate coupling, the local and the global approach both agree with the exact solution and for strong coupling, the global approach is preferable. These results are backed by detailed derivations of the regimes of validity for the respective approaches.


Introduction
Master equations are a powerful tool to study open quantum systems [1,2]. They allow for a description of the relevant degrees of freedom only, which evolve under the influence of all other degrees of freedom that are not of immediate interest. These other degrees of freedom are collectively called the environment. A particularly simple situation occurs when the system can be described by a time-local master equation with constant dissipation rates [3][4][5][6]. This results in Markovian evolution, where knowledge of the density matrix at a given time is sufficient to predict all future observables, which implies an environment that has no memory. Here we refer to this type of master equations as Markovian.
Compared to the complete problem of describing all the degrees of freedom of system and environment together, a Markovian master equation governing only the system degrees of freedom is an immense simplification. Such a drastic reduction of complexity usually comes at a price. In this case, the price comes in the form of strong approximations which are not always justified. Studying these approximations is thus of utmost importance and indeed, there is a large body of literature that addresses these issues [5, (for a recent review, see Ref. [6]). However, a large number of these studies focus on an environment that drives the system towards equilibrium. The recent rise in interest in out-of-equilibrium quantum systems, and in particular in quantum thermodynamics, calls for revisiting the question of the validity of the widely used Markovian quantum master equations [15,16,18,[20][21][22][23][24][25][26].
The standard description of these systems crucially relies on Markovian master equations to predict the relevant observables, such as heat currents and power. Two main approaches are followed in the literature. The first is a local approach, where the thermal baths couple locally to sub-systems of the machine. The second is a global approach, where thermal baths couple to the global eigenmodes of the machine. As the two approaches are conceptually different, there has been considerable debate about which one should be used in order to accurately describe thermal machines, and more generally out-of-equilibrium systems. Since the global approach describes equilibrium situations accurately (see below), while the local in some cases does not, there has been incentive to use the global approach out of equilibrium as well. Furthermore, the local approach is often believed to be more phenomenological in nature [13,14,19,28,57] and it was even argued that it is unphysical in certain regimes [27,58,59].
The goal of the present work is to discuss these questions in depth. We will consider a system for which the full unitary dynamics of the machine and the thermal baths can be solved exactly. This allows us to evaluate the performance of local and global master equations for the machine against the exact dynamics. In addition, we give detailed derivations of the local and the global approaches and discuss the involved approximations. Specifically, we consider a heat engine introduced by Kosloff [36], which can be implemented in superconducting circuits [44]. The machine consists of two subsystems (oscillators), which couple to different thermal baths, and to each other via an energy conserving interaction. In case the two oscillators have different frequencies, the machine requires an external driving field, making the Hamiltonian time-dependent. The entire system (machine plus baths) consists only of harmonic oscillators with quadratic interactions. Therefore, the system can be described exactly at the level of covariance matrices, and can be treated numerically even for relatively large baths with arbitrary precision. This exact numerical solution serves as a benchmark for evaluating the performance of both the local and global master equations. Focusing on an out-ofequilibrium steady-state regime, we discuss relevant thermodynamical observables, such as heat currents and power, as well as the quantum state (density matrix) of the machine degrees of freedom.
Our results demonstrate an overall excellent agreement between the predictions of the local approach and the exact solution. In the weak inter-system coupling regime, we see that the local approach provides an accurate description of the system, capturing both the thermodynamical observables and the quantum state. On the contrary, the global approach fails in this regime. Moving to the regime of intermediate coupling, we find that both approaches provide good descriptions of the system. Notably, the local approach still reliably captures all thermodynamical features of the machine. For very strong inter-system coupling strengths, the local approach starts to fail while the global approach still yields a faithful description of the system. We provide a detailed derivation of the regime of validity for each approach.
In the final part of the paper, we briefly discuss the case of finite dimensional machines. In particular, we consider the two-qubit entangler of Ref. [51], which is analogous to the heat engine setup considered in the first part, but with the two machine oscillators replaced by two qubits with equal level spacing (i.e. no external drive). While solving the total system (including the baths) exactly is unfortunately out of reach in this case, we can still compare the local and global approaches. We find very similar behavior to the results of the first part. In particular, the global approach still fails in the weak coupling regime, while for intermediate coupling, the two approaches agree well.
The paper is structured as follows. Section 2 gives a more detailed introduction to the local and global approaches. Sections 3-7 are devoted to the heat engine. In Sec. 3, we introduce the system. The different master equations and the respective approximations are discussed in Sec. 4, and the exact numerics are discussed in Sec. 5. The observables which are investigated are introduced in Sec. 6 and the results are given in Sec. 7. The qubit entangler is then discussed in Sec. 8 before we conclude in Sec. 9.

Local vs global
Before going into details, we provide a short introduction to the two commonly used Markovian master equations which we discuss. In the local approach, the thermal baths couple to the eigenstates of sub-systems of the machine, while the global approach is based on a secular approximation. For time-independent Hamiltonians, the global approach corresponds to baths that couple to the delocalized eigenstates of the system Hamiltonian. In an equilibrium situation (i.e. baths at equal temperatures and timeindependent Hamiltonian) the global master equation results in the desired steady state which is given by a Gibbs state with respect to the system Hamiltonian [1,60] (for a discussion on deviations from the Gibbs state due to the finite coupling between system and bath, see, e.g., Ref. [61,62]). The local approach on the other hand results in a product of Gibbs states with respect to the sub-system Hamiltonians [60]. For finite interactions, the global master equation is therefore usually considered superior to the local master equation (cf. Fig. 2). The situation drastically changes if we move away from equilibrium. Clearly, if the sub-systems do not interact, each sub-system should thermalize to its respective bath. While the local master equation yields correct results in this case, the global approach fails, resulting in a finite energy current through the system even in the absence of an interaction (cf. Fig. 3) There therefore exist limiting cases where we expect one approach to clearly be superior to the other. The goal of the present work is to connect these dots by comparing the local and the global master equations to exact numerics for a wide range of parameters in and out of equilibrium. Our main interest thus lies in the observables related to the machine operation, i.e. heat currents, powers, and efficiencies. Moreover, we also compare the quantum states of the machine.
We note that the validity of local and global approaches has been investigated before, even for out-of-equilibrium systems. Ref. [16] provides a detailed study of the different approximations, however, the authors only consider the state of the system and not the energy flows. Ref. [11] shows that the global approach neglects terms that influence the current through a two-terminal electric conductor. In Ref. [15], it was shown that the global approach can erroneously result in a vanishing heat current through a spin chain while a local approach shows good agreement with results obtained from a Redfield equation. The validity of the Redfield equation was shown in Ref. [23], where a system analogous to our heat engine in the absence of an external drive is investigated. Refs. [58,59,63] argue that a local master equation can violate the second law of thermodynamics when a non-energy preserving interaction between the subsystems is considered. These violations were however shown to be of the order of terms that are dropped when deriving a local master equation [22,64]. In Ref. [22], it was shown that in the absence of degenerate subspaces, the local approach can be understood as the zeroth order of a perturbation series in the inter-system interaction. However, most thermal machines cited above crucially rely on such degenerate subspaces. Finally, Ref. [20] investigates heat currents through a system that is equivalent to our qubit entangler, without providing a benchmark. Their results agree with the results presented in Sec. 8. All these previous works motivate our detailed investigation which compares the different master equations for a wide range of parameters. Figure 1. Sketch of the heat engine. The system consists of two harmonic oscillators with different frequencies. Each oscillator couples to a thermal bath modeled by a collection of harmonic oscillators. An external field mediates a weak coupling between the two oscillators. A thermal gradient can result in a heat flow from the hot to the cold bath, injecting some of the energy into the external field which can in principle be used to charge a battery. A thermoelectric implementation of this machine in superconducting circuits is proposed in Ref. [44]. In this work, we compare different Markovian master equations for the description of this machine, using exact numerics as a benchmark.

Heat engine
The heat engine we consider [36,44] is sketched in Fig. 1 and consists of two harmonic oscillators, with frequencies Ω c and Ω h , which each couple to a bath. Here the subscript labels both the oscillators as well as the corresponding bath (where c stands for cold and h for hot). In this setup, the oscillators constitute the weakly interacting sub-systems mentioned above. The frequencies of the oscillators differ from each other by In order to extract power from the machine, we consider a time-dependent external field with frequency E which mediates a coupling between the harmonic oscillators. The purpose of the heat engine is then to use a heat flow from the hot bath to the cold bath in order to increase the power of this external field. In Ref. [44], this external field is provided by a voltage and the power is directly related to an electrical current flowing against the voltage. The total Hamiltonian of the heat engine (including heat baths) can then be written as (throughout this paper we set = 1) whereâ α (b k,α ) denote annihilation operators of the system (baths), g denotes the interaction strength between the harmonic oscillators of the system, ω k,α are the frequencies of the bath modes, and γ k,α denote the interaction strengths between system and bath modes. The superscript l denotes the laboratory frame. The external drive accounts for the energy that is needed to convert a photon with frequency Ω c into a photon with frequency Ω h and vice versa. In contrast to Refs. [58,59,63], which lack an external field, we find that the presence of such a field ensures that the local master equation does not result in any violations of the laws of thermodynamics.
In the following, we will work in a rotating frame defined by the transformation This results in the time-independent Hamiltonian Note that the interaction termsV α are invariant under the transformation.
If not explicitly stated otherwise, all equations are given in the rotating frame.

Master equations
In this section, we consider different Markovian master equations which are used to describe the evolution of the reduced density matrix of the system. These equations allow for analytic expressions of observables such as power, heat, and efficiency which are then compared to the ones obtained from the exact dynamics. The standard way of deriving Markovian master equations is to first perform Born-Markov approximations. This procedure is discussed in detail elsewhere (see for instance Refs. [1,6,16,65]). We therefore only summarize the approximations and give the resulting expression for the system under consideration. The approximations are: • Born approximation: TreatingV α perturbatively to lowest order, • Markov approximation: Assuming invariance ofρ(t) on time-scales of the order of τ B , whereρ(t) is the reduced density matrix in the interaction picture and τ B denotes the bath-correlation time that will be introduced below. In the interaction picture (and the rotating frame), the Born-Markov approximations result in the following master equation, Here operators in the interaction picture are given bỹ whereÂ denotes the operator in the Schrödinger picture and We further introduced the bath correlation functions where we assumed the baths to be in thermal states. The Bose-Einstein distribution is given by The bath correlation functions are usually peaked around τ = 0 and decay for large times. This decay defines the bath-correlation time τ B such that C j,α (τ B ) C 2,α (0) (note that C 1,α (0) vanishes as T α → 0). For an explicit evaluation of the bath correlation functions and a discussion on the relevance of τ B (including a discussion of the zero temperature limit), we refer the reader to Ref. [16].
The Born-Markov equation represents the starting point of our analysis. Since it does not guarantee positive evolution [1], further approximations are usually made to obtain a master equation in Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form [3,4] ensuring completely positive dynamics. In the following, we will discuss two popular approximations, the local and the global approach, in some detail.

Local master equation
The Markov approximation that was made to obtain Eq. (5) is responsible for the fact that the density matrix under the integral is independent of τ . This approximation is valid as long as the characteristic time over whichρ(t) varies is much larger than τ B . In the same spirit, we can make the approximatioñ in the integral of Eq. (5). Note that we make this approximation in the rotating frame, where the fast oscillations with frequency Ω α are encoded in the bath correlation functions [cf. Eq. (8)]. In our model, we havẽ whereᾱ = α. The approximation in Eq. (10) is therefore expected to be good as long as gτ B 1. For reasonably small values of g, this approximation is therefore completely consistent with the Markov approximation.
This approximation directly results in the local master equation (in the Schrödinger picture) without the need of a secular approximation where we neglected a physically irrelevant constant and introduced as well as the Lindblad superoperators and the energy damping rate where ρ α (ω) denotes the spectral density. The renormalized Hamiltonian is given bȳ where P denotes the Cauchy principal value. We note that this renormalization should be small in order for the Born-Markov approximations to be valid (this can be understood by noting thatρ(t) will have terms that oscillate with frequency Σ α ). In the following, we will thus neglect this renormalization. As shown in Sec. 7, excellent agreement between the local master equation and exact numerics is found for a wide range of parameters without taking into account the renormalization of the Hamiltonian; see also Ref. [16]. As discussed above, the local master equation is justified as long as gτ B 1. With the help of Eqs. (8), this temporal inequality can be translated into an inequality involving the Fourier transform of the bath correlation functions. If τ B g 1, then we have exp(igτ )C j,α (τ ) C j,α (τ ) for all τ since we can approximate C j,α (τ > τ B ) 0. It is straightforward to show that this is fulfilled as long as the Fourier transform of the bath correlation functions can be approximated as constant over the energy scale of g. With the help of Eqs. (8) and (15), this results in the inequalities Here we used the fact that the main contribution for the bath correlation functions comes from the low-frequency terms. For a bosonic bath with Ohmic spectrum, the last equations are fulfilled as long as g Ω α , which is usually the case for systems where the local approach is employed [38,44,46,51]. For fermionic baths, equations analogous to Eqs. (17) can be derived. These are usually not fulfilled at low temperatures due to the step-like behavior of the Fermi-Dirac distribution [23,66]. Note that the validity of the Markov approximation, which underlies both the local as well as the global master equation, requires conditions obtained from Eq. (17) by exchanging g with κ α (Ω α ).
It is sometimes stated that the local approach is only valid for interaction strengths much smaller than the induced broadening g κ α (Ω α ) (see for instance Ref. [23]). As we show in Sec. 7, the local approach gives reliable predictions even for interactions that are several times the broadening. This is in complete agreement with Eq. (17).
We note that the local master equation is also obtained in the so-called singular coupling limit [1,5,10,12], where the bath correlation functions in Eqs. (8) tend to a delta function. This limit is often dismissed as being unrealistic. However we stress that the bath correlation functions only have to behave like delta functions on time-scales of the order of 1/g for the local master equation to be valid. In the weak coupling limit, which is often the regime of interest for thermal machines, 1/g can naturally be much bigger than τ B .
In order to solve the master equation, we make use of its bi-linearity (in creation and annihilation operators) which implies that a Gaussian state remains Gaussian at all times. Furthermore, there are no terms in the master equation which result in a displacement of the state. We can therefore restrict the analysis to states which have â α = 0. Then the state is fully described by its covariance matrix. From the local master equation in Eq. (12), one can derive the following differential equations for the covariance matrix elements where In the steady state, these differential equations are solved by the time-independent expressions

Global master equation
The second approximation we consider is the secular approximation. This approximation consists of dropping all the terms in the Born-Markov master equation [cf. Eq. (5)] which oscillate as a function of time t. The secular approximation is expected to hold over time-scales much bigger than the inverse frequencies of the oscillating terms and obviously becomes better as these frequencies increase. In order to identify the oscillating terms, we write the interaction picture operators in Eq. (11) as where we introduced the operatorŝ The secular approximation is obtained by plugging Eq. (21) into Eq. (5) and dropping all terms that oscillate with exp[2igt]. This results in the master equation with the frequencies Ω α,± = Ω α ± g.
The renormalized Hamiltonian reads and As for the local approach, we will neglect the renormalization of the Hamiltonian in the following. The secular approximation results in a master equation where the eigenmodes of the Hamiltonianâ σ are the relevant degrees of freedom. The action of the baths decouples into a bath for each eigenmode, which drives the mode towards thermal equilibrium characterized by the occupation numbers where we introduced We note that in the limit g → 0, the local master equation is no longer recovered. Because the secular approximation is no longer justified in this limit, we expect the global master equation to break down. This is also consistent with the fact that we expect the secular approximation to be valid as long as the frequency of the neglected oscillating terms is much bigger than the linewidths, i.e. κ α g. Since the Markov approximation requires κ α Ω α , and the local approach is valid for g Ω α , there is an overlap between the regimes of validity of the local and the global master equation. This implies that the local and the global master equation together are enough to describe the system for all parameters which allow for a Markovian master equation.
In the global approach, the covariance matrix is governed by the differential equations In the steady state, we find that the state is a product of thermal states with occupation numbers â † ±â± = n σ as expected.

Exact numerics
In this section we briefly describe how we obtain exact numerics which are used as a benchmark when comparing the different master equations. To this end, we simulate the unitary evolution generated by the Hamiltonian in Eq. (2) for big but finite baths. The key element here is that the Hamiltonian in Eq. (2) is quadratic, so that Gaussian states (such as thermal states) remain Gaussian throughout the whole evolution. As such, they can be fully characterized by only their first and second moments (see, e.g., [67]). This allows us to characterize the time-evolved state of the whole system, including the thermal baths, by a matrix of size ∼ N , where N is the total number of oscillators involved.
We consider baths made up of n + 1 harmonic oscillators, so that the total size of system and baths is N = 2(n + 2). The bath modes are chosen to be uniformly spread over a range (0, ω c ), where ω c is a cutoff frequency. That is, for k = 0, .., n. This definesĤ α in Eq. (2). Let us now turn our attention toV α , and hence to the couplings γ k,α . First note that the action of the baths in the Markovian master equations, Eqs. (14) and (23), is captured by the spectral density given in Eq. (15). It is then common to use an ad hoc form for the spectral density in the continuum limit, instead of specifying the coupling constants γ k,α . A common choice for the spectral density is an Ohmic spectrum, which holds for low frequencies, ω ≤ ω c . In order to relate this approach to the γ k,α 's, which are necessary to simulate the full Hamiltonian in the finite-baths scenario, let us integrate Eq. (15), obtaining, It is now convenient to discretize the integral, from where it immediately follows, which becomes increasingly accurate with increasing n. In the particular case of an Ohmic distribution, we obtain, Equation (35), together with Eq. (31), provides a simple recipe for building the discrete version of the Hamiltonian in Eq. (2) for a given spectral density. The specific choice of ω c and n for our simulations, as well as the dependence of the results on this choice, is discussed in Appendix A. The initial state for the simulations is taken to be of the form, whereτ βα are thermal states,τ with inverse temperature β α = 1/(k B T α ) and we note thatĤ l α is in the laboratory frame. Here, the initial state of the machineρ S can be an arbitrary Gaussian state. In our simulations, we takeρ S = e −β h Ω hâ † hâ h /Z h ⊗ e −βcΩcâ † câc /Z c . We note that with this choice, the state given in Eq. (37) is Gaussian.
Once the Hamiltonian and the initial state are defined, we consider the closed unitary dynamics of the full compound under the Hamiltonian in Eq. (2) (it is convenient to work in the rotating frame, where the Hamiltonian is time independent). The dynamics can be derived by considering the Heisenberg equations of motion of the operatorsâ α ,â † α , {b k,α ,b † k,α }. For details on the derivation, we refer the reader to Ref. [16].
In order to simulate the steady state using finite degrees of freedom, one needs to let the whole compound evolve for a time t that satisfies τ eq t τ rec , where τ eq is the equilibration time of the system and τ rec is the recurrence time of the bath. From Eqs. (12) and (18), we infer the equilibration time to be 1/τ eq ≈ max{κ h , κ c } (see also [68], where the equilibration time is discussed explicitly for finite systems). On the other hand, the recurrence time scales linearly with the number of oscillators in the bath [16]. Hence, by taking a sufficiently large bath (in the simulations we take ∼ 400 oscillators), we can ensure that τ eq τ rec . In our simulations, we take t ≈ 20τ eq .

Observables and reduced states
In this work, we are particularly interested in the energy flows that traverse the quantum thermal machine in a non-equilibrium situation. However, to compare to previous works [16], and to further assess the validity of the different master equations, we also consider the obtained steady states.
6.1. Heat currents, power, and efficiency 6.1.1. Local master equation As we consider a heat engine, the main quantity of interest is the power that is produced. For our system, it is defined as [33] where the superscript denotes the laboratory frame and · · · , denotes the ensemble average in the rotating frame. Note that positive power implies that energy leaves the system. In addition to the power, we consider the heat currents that enter (or leave) the system. To this end, we write where L l α is a super-operator that groups all the dissipative terms which arise from bath α in the laboratory frame. Note that under the Born-Markov approximation, the dissipators of different baths can be added [69]. The heat currents are then defined as where the dissipator in the rotating frame is related to the dissipator in the laboratory frame by With our sign convention, a positive heat current implies energy entering the system from the bath. In the steady state, the first law of thermodynamics thus reads The efficiency is defined as the ratio of the obtained power, divided by the heat that originates from the hot bath In the regime where the system operates as a heat engine (i.e. P > 0 and J h > 0), the second law of thermodynamics forces the efficiency to remain below the Carnot limit For the local master equation in Eq. (12), the heat currents in Eq. (41) can be written as From Eqs. (18), we can infer that the second term in the heat current decays exponentially in time. In the steady state, from Eqs. (20) and (39), we find for the power and the heat current resulting in the efficiency We note that this efficiency fulfills Eq. (45). When the frequencies are chosen such that the efficiency is above the Carnot efficiency, then we find P < 0 and J h < 0. As long as η < η C , our machine is thus a heat engine, η = η C denotes the point of reversibility, where all the energy currents vanish, and for η > η C , the machine acts as a refrigerator, using power to induce a heat current from the cold bath to the hot bath [44]. Finally, we note that for Ω h = Ω c , where there is no external power, heat always flows from the hot bath to the cold bath as dictated by the second law of thermodynamics. In contrast to models which consider non-energy preserving interactions, the local approach does not violate the laws of thermodynamics when including an external field that provides the energy to convert photons of frequency Ω c into photons of frequency Ω h .

Global master equation
In the global master equation, the bath couples to global states which are dressed by the external field. Therefore, the dissipative terms include the external field and the definitions introduced in the last subsection are no longer valid. To define heat and work in the global approach, we follow Refs. [12,31,35,70,71].
To this end, we first write the dissipator in the global master equation [cf. Eq. (23)] as the sum of four dissipators given by We further introduce the thermal stateŝ which constitute the steady states of the respective operators, i.e. L α,σρα,σ = 0. The heat currents in the steady state (denotedρ) are then defined as and the power is given by the first law, cf. Eq. (43). We note that for the global approach, we enforce the first law while in the local approach, it follows from the expressions for heat currents and power. An explicit calculation results in [70] Note that if Eqs. (17) are fulfilled, the last expression reduces to which is the expression obtained by the local approach in the limit g κ α . As expected, if the approximations leading to both the local and the global master equation are justified, the two approaches give the same result. Whenever Eqs. (17) are not fulfilled, the global approach implies that the heat engine makes use of two channels labeled by the subscript σ. Each of these channels has an efficiency where J α = J α,+ + J α,− and J c is obtained from J h by exchanging the labels c ↔ h

Exact numerics
When dealing with the full Hamiltonian, we define heat currents as the energy lost by the bath, where we note that bothĤ l α andˆ l (t), the unitarily time evolved state of the whole compound, are taken in the lab frame. The power can be obtained through Eq. (39).

Reduced states
In addition to the energy currents, we consider the reduced state of the system as obtained by the different solutions. Since the considered Hamiltonian is bi-linear in bosonic annihilation and creation operators, it suffices to consider the covariance matrices. It will be convenient to work in thex,p basis, withx α = 1/2Ω α (â † α +â α ) andp α = i Ω α /2(â † α −â α ). Defining the vectorr = (x h ,x c ,p h ,p c ), the covariance matrix is given by, For the master equations, the covariance matrix can be obtained straightforwardly from Sec. 4 [cf. Eqs. (20), and the discussion around Eq. (30)]. In the case of the exact numerics, one needs to consider the relevant entries of the covariance matrix of the whole compound (see, e.g., Ref [16]).
In the next section, we compare the reduced states obtained from the master equations to the reduced state obtained through exact numerics. As a measure of distinguishability between two statesρ andσ, we consider the fidelity, defined as, For two-mode Gaussian states, represented by covariance matrices C 1 and C 2 , and vanishing first order moments, the fidelity is given by [72], Here a = det(C 1 + C 2 ), b = 2 4 det(JC 1 JC 2 − I/4), c = 2 4 det(C 1 + iJ/2) det(C 2 + iJ/2), and the matrix elements of J are given by J kl = −i [r k ,r l ] .

Results
Our results for the heat engine are illustrated in Figs. 2-7. We divide our results into three regimes. First we discuss the equilibrium regime, where the global approach shows excellent agreement with numerics for all values of the inter-system interaction g. Then we discuss the presence of a thermal bias but no external field. In this case, the global approach breaks down for small values of g. Finally, we focus on the heat engine regime which requires both a thermal bias as well as an external field. Again, the global approach breaks down for small values of g as expected. In all regimes, the local approach performs well for g Ω α (for both α = c, h), which is where Eqs. (17) are fulfilled for our system. As expected, the global approach performs well for g κ α (for both α = c, h), which is where the secular approximation is well justified. In the following, all energies are given in units of Ω c and all energy currents are given in units of Ω 2 c .

Equilibrium
We first present our results for the equilibrium case, where T c = T h and Ω h = Ω c . Fig. 2 (a) shows the fidelities between the steady states obtained from Eqs. (12) and (23), and the steady state obtained from exact numerics. As expected, the global approach yields an accurate description of the steady state while the local approach gets progressively worse as the inter-system interaction g increases. The same conclusions can be drawn from Fig. 2 (b), where observables such as occupation numbers are plotted. We note that Fig. 2 goes up to g = Ω c /2 and thus covers interaction strengths which are much higher than the ones usually considered when the local master equation is employed. We note that in the limit g → 0, the local and the global master equations result in the same steady state which is given by a product of thermal states with respect to the local oscillator Hamiltonians. If we also have κ c = κ h , the two master equations coincide. For energy damping rates that differ, i.e. κ c = κ h , we expect the two approaches to result in different predictions for the transient regime.
These results confirm that in an equilibrium situation, the global approach is indeed preferable over the local approach, at least when one is interested in the steady state properties. One might be tempted to believe that this conclusion carries over to the out-of-equilibrium regime. That this is not the case is illustrated below.

Thermal bias
We now turn to the case of a thermal bias T c = T h , but still no external field, i.e. Ω h = Ω c . Our results for this regime are illustrated in Fig. 3. The heat current is plotted in Fig. 3 (a). For all models, we find J h = −J c (first law) and J h ≥ 0 (second law). As discussed above, the global approach breaks down in this case for interactions g κ α . In this limit, the secular approximation is no longer justified and the global approach gives the unphysical result of a finite heat current in the limit g → 0. The local approach on the other hand predicts the heat current extremely well up to g = Ω c /2.
The inset of Fig. 3 (a) shows the fidelities with respect to the numerical solution. As expected, the local approach reliably reproduces the steady state for small values of g while the global approach works well for large values of g. Figure 3 (b) shows other observables such as the occupation numbers, leading to the same conclusions. Note that  while the local approach shows very good agreement with numerics, the global approach breaks down when g κ α . (b) Power and efficiency. Again we find excellent agreement between the local approach and exact numerics. In particular, the numerical value for the efficiency is very close to the universal value η = 1 − Ω c /Ω h obtained from the local master equation. Note that when g becomes comparable to Ω c , the local approach starts to deviate from the exact numerics and the global approach becomes preferable. Parameters: Ω c = 1, Ω h = 2 κ h = κ c = 0.05, k B T c = 0.5, k B T h = 5. Parameters numerics: ω c = 3, n = 400, t = 20/κ. the occupation numbers obtained from the local approach differ quite a bit from the ones obtained from exact numerics for large interactions. Nevertheless, the heat current is still captured very well by the local approach.

Thermal bias and external field
Finally, we consider the regime where the considered system performs as a heat engine. This requires both a thermal bias T c = T h and an external field Ω h = Ω c . In the thermoelectric realization of Ref. [44], this situation corresponds to the presence of a thermal and a voltage bias. In Fig. 4, the energy flows through the heat engine are plotted. While the local approach agrees extremely well with exact numerics, the global approach fails for small g. Note that all models fulfill J c +J h = P (first law), and η < η C (second law). For the global approach, we note that it is crucial to use the definitions of energy currents given in Sec. 6.1.2. If one uses instead definitions similar to the local approach [see Eq. (39)], the global approach results in incorrect heat currents, leading in particular to P = 0 and J h = −J c . This is consistent with the results of Ref. [15], which discusses the absence of any currents as expressed through system observables in the global approach.
In Fig. 5 we compare the steady states against the exact numerics, using again fidelity as a figure of merit. Similarly to the case of a thermal bias without external field, we find that the two approaches give a faithful description in their respective regime of validity.
For completeness, Fig. 6 illustrates the heat currents as a function of temperature in  the presence and absence of an external field. These results strengthen the conclusions drawn above: The global approach is valid for g κ α while the local approach is valid for g Ω α . In particular, the insets show that even at reasonably strong interaction . Heat engine performance as a function of Ω h . We find good agreement between the local approach, the global approach, and the exact numerics except for efficiencies close to the Carnot efficiency (here η C = 0.5). The large differences in efficiencies result from small differences in power and heat currents. Parameters: Parameters numerics: ω c = 3, n = 400, t = 20/κ. strengths g, the local approach agrees well with numerics.
Finally, we further illustrate the heat engine performance in Fig. 7 which shows the power and efficiency as a function of Ω h , which determines the external field frequency (given by Ω h − Ω c ). We find good agreement in the power as well as the efficiency between the local approach, the global approach, and the exact numerics. Only when the machine is operated close to the Carnot point (Ω h /T h = Ω c /T c ) do the efficiencies deviate considerably. This is due to the fact that the power and the heat current become very small. Small differences in the energy flows then translate into large differences in the efficiency.

Qubit entangler
To complete our discussion, we also consider a quantum thermal machine featuring finite-dimensional systems. Specifically we consider a quantum thermal machine consisting of two interacting qubits coupled to separate bosonic thermal baths, as shown in Fig. 8, analogous to the setup of Fig. 1 for harmonic oscillators. This machine can generate entanglement between the two qubits in the steady state, as shown in Ref. [51]. In the following, however, we focus on comparing the steady states obtained from local and global master equations in the same spirit as above.
We denote the eigenstates of the free Hamiltonians of the qubits by |0 , |1 , and set the ground state energy to zero. The Hamiltonian of the system is then given by where Ω c , Ω h are the energy gaps of the qubit, and g is the interaction strength. As in Sec. 4, we compare local and global master equation models for the evolution of the system. We will focus on the degenerate case where Ω c = Ω h = Ω. Both the local and global master equations can be written in Gorini-Kossakowski-Sudarshan-Lindblad (GKSL) form with constant rates [3,4], and so lead to Markovian (specifically semigroup) evolution where D is defined in Eq. (14),L α,k are jump operators, and Γ α,k ,Γ α,k the corresponding rates. We consider bosonic baths for which Here, κ α (ε) are the bath coupling strengths, and ε is the (absolute) energy difference associated with the jump induced byL α,ε . Thus theL † α,ε correspond to jumps from lower to higher energies, absorbing energy from the bath α, whileL α,ε correspond to jumps decreasing the system energy, dissipating energy into the bath α.
Local and global master equations for the two-qubit machine can be derived using the same techniques as in Sec. 4. The system-bath coupling can be taken to have the same form as in Eq. (2), with the system annihilation and creation operators replaced byÂ c =σ − ⊗ I 2 ,Â h = I 2 ⊗σ − andÂ † c =σ + ⊗ I 2 ,Â † h = I 2 ⊗σ + respectively, wherê σ − = |0 1| andσ + = |1 0|. We take the spectral density of the baths to be Ohmic, as before. The bath coupling strengths are then linear in energy for some constants ν α . We denote κ α = κ α (Ω). For the local master equation, there are just two jump operators, both corresponding to transitions with energy Ω. They are given bŷ For the global master equation, the jump operators are found by diagonalizing the system Hamiltonian. Denoting the eigenvalues and eigenstates ofĤ S by λ and |ϕ λ respectively, one hasL In the degenerate case, Ω c = Ω h = Ω, the eigenvalues are 0, Ω ± g, and 2Ω, and the corresponding eigenstates are The possible transition energies are 2g, Ω ± g, and 2Ω. However, only transitions with energies Ω±g can be induced by the system-bath coupling considered here. The non-zero jump operators areL We note that these are indeed global in the sense that they involve transitions to and from the non-separable states |ϕ Ω±g .
We can compute the steady state solutions of both the local and globel qubit master equations. As before, we compare them varying the interaction strength. In Fig. 9 (a) we show the heat current as a function of the interaction strength between the two qubits. As for the harmonic oscillators, the global approach predicts a constant heat current, which is clearly unphysical as g → 0, while the local model predicts a vanishing heat current in this limit, as expected. The two models agree well for intermediate coupling strength (i.e. g/κ α 10 in this case, note that κ α is taken one order of magnitude smaller than in the previous sections). In Fig. 9 (b) we show the fidelity between the steady states of the two models. Again, we see that the states agree well unless the coupling is weak.

Conclusions
We investigated the accuracy of the local and global master equations for predicting thermodynamic quantities as well as system steady states describing a quantum heat engine. Exact numerics were used to benchmark the results. We found that the two approaches work very well in their respective regimes of validity which are (for bosonic baths with Ohmic spectral density): • Local approach: g Ω α , • Global approach: g κ α .
More generally, the condition under which the local approach gives a faithful description is given by Eq. (17). Since the Markov approximation, which underlies both approaches, requires κ α Ω α , the two regimes of validity overlap. As expected, we find good agreement between the local and the global approach in this region of parameter space.
We note that the local approach is by no means more phenomenological than the global approach. Indeed, the approximation leading to the local master equation is completely analogous to the Markov approximation and has a well defined regime of validity.
Finally, we also investigated a qubit entangler. For this system, no benchmark is available. However, the similarity to the results obtained for the heat engine strongly suggests that similar conclusions with respect to the applicability of the local and the global master equation are valid. We therefore conjecture that our results are qualitatively valid for a variety of baths and system Hamiltonians. As long as the bathcorrelation time is much shorter than any inverse inter-system interaction strength, the local approach is valid. The global approach is valid as long as the inter-system interaction strengths are much stronger than the system-bath interaction strengths.
We therefore conclude that the local approach provides a valid description for thermal machines that consist of weakly interacting sub-systems.
Note added -During the writing of this manuscript, we became aware of related work [63]. There the authors also compare local and global master equations, but in the absence of external fields, finding good agreement with the results presented here. In all the figures of the main text, we take ω c = 3 and n = 400. In order to see how sensitive our results are to this choice, we plot the heat current as a function of n (the inset) and ω c , and compare it with the analytic results (using the local approach) in Fig. A1. It is clearly observed that the results are independent of n. On the other hand, we see a small dependence on the values of ω c . For ω c ≈ Ω α , the results do not closely match the analytics, which is expected because energetically possible transitions are not captured by the bath. For ω c ≥ 2Ω α , very good agreement is obtained, with small differences that increase with ω c . This is due to the renormalization of the Hamiltonian [cf. Eqs. (16) and (27)], which is neglected in the analytic calculations. The choice ω c = 3 is hence large enough to capture the different energy transitions, and at the same time not too large so that the effect of the renormalization can be neglected. Similar considerations hold for the case of non-degenerate frequencies of the oscillators of the system. We hence conclude that our numerical benchmark is quite robust to the choice of n and ω c .