Local master equations bypass the secular approximation

Master equations are a vital tool to model heat flow through nanoscale thermodynamic systems. Most practical devices are made up of interacting sub-system, and are often modelled using either local master equations (LMEs) or global master equations (GMEs). While the limiting cases in which either the LME or the GME breaks down are well understood, there exists a 'grey area' in which both equations capture steady-state heat currents reliably, but predict very different transient heat flows. In such cases, which one should we trust? Here, we show that, when it comes to dynamics, the local approach can be more reliable than the global one for weakly interacting open quantum systems. This is due to the fact that the secular approximation, which underpins the GME, can destroy key dynamical features. To illustrate this, we consider a minimal transport setup and show that its LME displays exceptional points (EPs). These singularities have been observed in a superconducting-circuit realisation of the model [1]. However, in stark contrast to experimental evidence, no EPs appear within the global approach. We then show that the EPs are a feature built into the Redfield equation, which is more accurate than the LME and the GME. Finally, we show that the local approach emerges as the weak-interaction limit of the Redfield equation, and that it entirely avoids the secular approximation.


Introduction
Master equations and quantum thermodynamics go hand in hand. The former have become essential tools to make sense of the 'thermodynamics' of quantum systems. But, conversely, the early works on quantum thermodynamics [2-4] focused on the study of the mathematical properties of master equations. Nowadays the field is evolving very rapidly [5], and quantum heat devices are making the transition from theory to experiments on a wide range of platforms, including trapped ions, solid-state systems, atomic gases, single-electron systems, nanoscale thermoelectrics, and superconducting circuits [6][7][8][9][10][11][12][13][14].
The Gorini-Kossakovski-Lindblad-Sudarshan (GKLS) quantum master equation [15,16] makes it easy to draw parallels between the dissipative dynamics of a single open quantum system and the thermodynamics of macroscopic devices [17][18][19]. Namely, these equations can be derived from first principles in the limit of weak systemenvironment coupling, and may lead to thermal equilibrium [20]. We shall refer to such 'thermalising' GKLS equations as global master equations (GMEs). Furthermore, heat currents can be formally defined such that they obey the second law of (classical) thermodynamics [2,17]. However, the underlying assumptions of the global master equation require a clean timescale separation [21], which may break down for, e.g., small multipartite quantum-thermodynamic devices that interact weakly among them (see, e.g., [22,23]) and large many-body open quantum systems [24].
Alternatively, in multipartite open quantum systems, master equations have been often built heuristically by 'adding up' GKLS terms (cf. Fig. 1). These are referred-to as local master equations (LMEs). While such equations do comply, by construction, with the minimum expectation of generating a completely positive dynamics, they have been criticised for their thermodynamic deficiencies [22,[25][26][27][28][29][30][31][32][33][34][35][36]. Namely, unlike global master equations, LMEs fail to bring systems to thermal equilibrium, even in the limit of weak system-environment interactions [29,32]. When applied to quantum heat devices, they entirely miss crucial physics, such as heat leaks and internal dissipation [26,37]. They may even predict flagrant violations of the Second Law of thermodynamics, in the form of cold-to-hot stationary heat flows [27]. Surprisingly, however, the LME does prove very accurate in some caseseven more accurate than the GME [22,23]. The aim of this paper is to understand when and why is this the case. Using the most suitable master equation in each situation can make a crucial difference when studying the thermodynamics of any nanoscale heat device.
Local master equations may be understood as a rough approximation to the true dissipative dynamics, valid in the limit of weak interactions between the sub-systems. Starting from a microscopic model this may be shown in two closely related ways-either by carefully introducing a coarse-graining in the time-evolution of the open system [38][39][40][41], or by truncating a perturbative expansion of the master equation in the internal coupling strength [42,43].
The disagreement between the steady-state thermodynamic predictions of LME and GME had been illustrated before [23,28,44]. Here, we put the spotlight on situations in which they agree in the steady state, but differ during the transient dynamics. In such cases, which master equation is correct? One that is thermodynamically sound, or a truncated series expansion with limited validity and serious thermodynamic deficiencies? Strikingly, our answer is that we should always trust the latter within its error bars, which assume weak internal couplings. This holds for any multipartite quantum heat device which uses frequency filters to couple to the environment [45], e.g., a qubit or a harmonic oscillator.
For illustration, we focus on the specific model of two coupled resonators that connect two thermal reservoirs. We find that the LME exhibits a family of exceptional points [46] (EPs) in its dynamics, while the corresponding GME does not. Yet, EPs have indeed been found experimentally in a superconducting circuit described by the same model [1]. In addition, at weak internal couplings, the transient heat currents obtained from the LME agree with the much more accurate Redfield equation [47] while not with the global approach. The reason for the failures of the GME is that it is underpinned by the (somewhat crude) secular approximation, which misses relevant physics. In contrast, the LME can be obtained directly from the Redfield equation, bypassing the secular approximation. Our results thus add much needed clarity to the longstanding 'local-versus-global' debate, and explain various previously reported features of both the global and the local approach. This paper is structured as follows: We begin by describing the details of the example model in Sec. 2 and then, give an overview of open-system dynamics within the global and local approach in Sec. 3. In Sec. 4 we introduce the concept of exceptional points and discuss how to search for them, given the equations of motion of a linear open quantum system. We then illustrate the different exceptional-point structure in the parameter space of the model according to the local, global, and Redfield approaches (Sec. 5). Finally, in Sec. 5.4, we show that, in resonance, the local approach succeeds at capturing the correct heatflow dynamics, while the global master equation fails. In Sec. 6 we summarise and conclude.

The system
The system considered in what follows is sketched in Fig. 1. It is comprised of two coupled resonators with frequencies ω c and ω h , which we model as where x x x α and p p p α are the corresponding quadratures, and k denotes the inter-resonator coupling strength. Also m α = 1. From now on, we work in units of = k B = 1. This model can describe two capacitively coupled superconducting resonators studied experimentally in Ref.
[1], in search for exceptional points (see footnote [48]). Each node is weakly connected to a local heat bath. Their equilibrium temperatures are T c < T h ('c' for cold and 'h' for hot), and we denote the respective dissipation strengths by λ α . When it comes to the baths themselves, they are bosonic reservoirs in thermal equilibrium. The resonatorbath couplings are where the bath operators B B B α are and q (α) ν stands for the coordinate of the environmental mode of bath α at frequency ω ν . The couplings g (α) ν can be collected into the spectral densities where δ(ω − ω ν ) stand for Dirac deltas. In the following we will work with the standard Ohmicalgebraic spectral densities Here, the phenomenological parameter Λ sets an 'ultraviolet' cut-off on the spectrum of the bath. We note that the precise analytical form of the spectral density does not play an active role in our problem, as long as Λ is large when compared with all other relevant energy scales.
where H H H (loc) i stands for the local Hamiltonian of each sub-system and V V V , for the interactions between them, and k controls the magnitude of the latter. For the sake of our discussion, let us consider that every sub-system couples to its own independent bath, i.e., where S S S α is an arbitrary system operator so that [H H H S , S S S α ] = 0, thus allowing for energy dissipation as well as decoherence.
The effective equation of motion for any arbitrary observable O O O of such system can be cast in the standard GKLS form [15,16]. Although its microscopic derivation is textbook material [49], we provide it in Appendix A, making as few assumptions as possible. Essentially, (i) that the dissipation strengths λ α are small, (ii) that system and bath start uncorrelated, (iii) that the bath correlation functions are short-lived, (iv) and that there is a clear-cut timescale separation between (fast) coherent and (slow) dissipative processes.
Assumptions (i)-(iii) yield the Redfield equation [47] (cf. Appendix A). Constraint (iv) justifies the secular approximation, and allows to reduce it to the GKLS form where the super-operators G † α (O O O) are given by The notation · · · indicates averaging over a thermal-equilibrium state of bath α. Hence, Eq. (10) together with (3) implies that G † α is O(λ 2 α ). Finally, the non-Hermitian operators A A A As already mentioned, (8) is a global master equation. The name tag highlights the fact that the operators A A A (α) ω enable 'jumps' between eigenstates of the full multipartite H H H S , rather than between states of the sub-system H H H (loc) α coupled to bath α (see Fig. 1).
Since full diagonalisation of H H H S is needed to construct these jump operators, setting up Eq. (8) may become computationally unworkable; especially, when one wishes to scale up a many-body open quantum system. Furthermore, the GME suffers from another important issue, especially when applied to systems with a dense energy spectrum. In such cases, assumption (iv) from the list above is likely to break down [24], which could invalidate the GME's predictions [22,23,50]. On the 'plus side', constructing jump operators that fulfil Eqs. (11) is guaranteed to bring the system into a state of thermodynamic equilibrium whenever all temperatures coincide, i.e., T α = T ∀α. This is analogous to the Zeroth Law of thermodynamics [51]. Since, in addition, the dynamical map resulting from (8) is completely positive [3], it can be shown that where the steady-state heat currentsQ α account for the stationary rate of energy influx from bath α (i.e., · · · ∞ denotes here stationary average). Eq. (12) is interpreted as the Second Law of thermodynamics, since it is formally identical to the statement of Clausius' theorem [2,17]. Therefore, the global GKLS equation is particularly well suited to study the thermodynamics of (weakly dissipative) open quantum systems.

The local master equation
Another type of quantum master equations may be built by simply adding up dissipators that act locally on a specific part of the system and S S S α = ω L L L (α) ω but, now, satisfying in place of Eq. (11a). Here, the sum does not run over the Bohr frequencies ω of H H H S , but over those of α H H H (loc) α ; hence the different notation ω . Such local equation is easily scalable, as the nodes must be diagonalised individually when searching for the operators L L L (α) ω . It also shares with the global master equation its GKLS form, which means that the resulting dynamics is, again, completely positive. Importantly, 'completely positive' is often equated to 'physical', but this is not always the case [32]-it is known that LMEs violate the Zeroth Law by construction [29]. Indeed, due to Eq. (15), the dissipators L † α 'try' to pull the system towards the local thermal state ∝ exp (− α H H H (loc) α /T α ), which does not commute with the full H H H S appearing in the 'coherent-evolution' term of Eq. (13). Hence, according to the LME, the system would never thermalise, even if all temperatures T α are identical. As a result, also the Second Law as stated in Eq. (12), may be violated [27]. In fact, the local approach invariably predicts unphysical cold-tohot heat currents whenever ω c /T c < ω h /T h in our example model of Eq. (1). But the LME can also be physical provided that one applies it to open systems with specific active environments [52], as opposed to problems with passive dissipation into equilibrium reservoirs.

Comparing local and global approach
A local master equation can be motivated physically beyond merely "adding up dissipators". For instance, we may obtain the LME directly from a formal collisional model, in the limit of instantaneous collisions [53,54]. However, in this case, the thermodynamics of the LME needs to be reassessed to account for the 'cost' of those collisions with bath 'units' [52,55]. On the other hand, for microscopic Hamiltonian models, such as (2), the general Redfield equation can be simplified by coarse-graining over a relevant timescale [38,39,56]. Depending on how this averaging is done, both the LME and the GME can emerge naturally [39,40].
The LME (13) can be alternatively viewed as the lowest order in an expansion of the GME (8) [42,43]. Namely, applying perturbation theory to the eigenstates and eigenvalues of H H H S for small k and following the procedure to obtain a GME [21,42], gives Whenever G † (0) α = L † α , one can claim that the LME becomes equivalent to truncating the expansion at zeroth order in k. Since all terms G † (n) are O(λ 2 ), the LME would thus be accurate within kλ 2 -sized 'error bars'. If k is weak enough, this may be acceptable when compared with the intrinsic error of the GME, set at O(λ 3 ) (cf. Eq. (8)). In light of this interpretation, the unphysical cold-to-hot heat flows predicted by the local approach are a mere artifact of the truncated expansion [42].
However, if the degeneracy of H H H S changes depending on whether k = 0 or k = 0, the zeroth order of the global dissipator does not coincide with the local one (i.e., . This is precisely what happens to the example model (1) when ω c = ω h = ω [22,23,39]. Instead, we find that it is the zeroth order of the Redfield dissipator that coincides with the local one (i.e.,

Exceptional points
We now introduce the concept of exceptional point formally, before drawing the link to opensystem dynamics. EPs are branch-point singularities which appear under variation of parameters of non-Hermitian matrices, such as those describing the dynamics of quantum dissipative systems. EPs have been used as resource for applications such as sensitivity amplification in microresonators [57,58], laser-mode selectivity [59], and topological chirality [60,61]. More recently, a gain in signal-to-noise-ratio has been shown in EP sensors [62] highlighting the practical relevance of exceptional points.

Formal definition and witnesses
Let M(k) ∈ C N ×N be an N -by-N matrix dependent on some parameter (or set of parameters) k.
We denote the right eigenvectors of this matrix The Notice that (v j | † is here a column vector due to the Hermitian conjugation. These two indexed families of vectors form a bi-orthogonal set [63]; that is, Hermitian. Note that, in general, (v i |v i can be negative. We say that the matrix M(k) has an exceptional point for those parameter choices k resulting in (v i |v i = 0 for two or more of the indices i ∈ {1, · · · , N }. This phenomenon is called self-orthogonality and it is the hallmark of the coalescence of two or more right eigenvectors [64]. In order to locate the exceptional points of M(k) we could search for zeros of the phase rigidities φ i (k) := |(v i |v i | for i ∈ {1, · · · , N } as a function of k. However, analysing the behaviour of every single eigenvector can be time-consuming for large N . Luckily, we shall only be interested in finding where in parameter space an EP of M(k) is located, and not in which or how many eigenvectors coalesce. We can thus exploit the fact that, at an EP, the set {|v i } i=1,··· ,N does not form a complete basis. Therefore, the matrix V k = (|v 1 , · · · , |v N ) features a singularity and the norm of its inverse diverges. Consequently, the condition number of V k , denoted κ(V k ), could be a suitable witness for an EP. Namely, where the operator norm · is defined as x is an arbitrary non-zero vector, and · p stands for the p-norm Here, the parameter p could take on any real value p ≥ 1.

EPs in open quantum systems
In many cases of practical interest, an open system may be fully described by choosing a set of observables θ = (ϑ 1 , · · · , ϑ m ) T and applying the corresponding adjoint quantum master equation to each of them; i.e., where the dissipators D † α can take, e.g., the global (14), and (A.15)).
The aim is to pick observables ϑ j so that (23) becomes a closed set of equations [65]. This can then be expressed in compact form as were the resulting matrix of coefficients M D ∈ C m×m is generally non-Hermitian. At an EP, M D stops being diagonalizable, which has a detectable impact on the dynamics and thermodynamics of the open system [65,66].
Specifically, in continuous-variable settings with linear H H H S -as the model (1) studied hereit is always possible to write a set of equations like (24). For instance, we can build θ by grouping the four position and momentum operators together with the 10 distinct combinations e.g., That is, for the Hamiltonian H H H S , one finds that M D is the direct sum of sub-matrices M D,1 ∈ C 4×4 and M D,2 ∈ C 10×10 . Looking back at Eq. (24), we thus see that the dynamics of the first-order moments q i decouples from that of the second-order moments C ij . Note that the same block-diagonal structure is found for the local (L ), the global (G ), and Redfield (R) equations. In Sec. 5 below, we focus on the appearance of EPs in the sub-matrices M D,1 . The discussion about M D,2 has been deferred to Appendix B.
5 Results and discussion

Local master equation
In order to obtain the equations of motion for the first-order moments q q q of the system according to the LME, we need to know the Bohr frequencies involved, and the corresponding jump operators. From now on, we work with resonant oscillators ω h = ω c = ω. Therefore, the Bohr frequencies in Eq. (14) are {ω } = {±ω}, and the local jump operators of (15) are where a a a α is an annihilation operator, so that Replacing these into the LME (14) gives with local coefficient matrix M L ,1 where we have introduced the notation Since M L ,1 is simple enough, one can obtain an analytic expression for the exceptional points in parameter space. Looking at its eigenvalues, we see that degeneracy appears, provided that Further replacing the expressions of the decay rates γ (α) ±ω into the coefficients ∆ α for our choice of spectral density (cf. Eqs. (5) and (10)) results in the remarkably simple expression Resorting now to the condition number of the eigenvector matrix of M L ,1 , we confirm that whole family of exceptional points does lie along (33) (see leftmost panel Fig. 2). That is, adjusting resonance frequency and internal couplings, it is possible to tune the system into an EP. Interestingly, for exceptional points to appear in this system, dissipation must be asymmetric and the oscillators resonant. Note as well that, at resonance, the LME cannot give rise to unphysical cold-to-hot heat currents [27]. Moreover, in a recent experiment with coupled superconducting resonators signatures of these EPs were indeed detected [1].

Global master equation
In order to find the jump operators A A A (α) ω within the global dissipators G † α (cf. (9) and Eqs. (11)), we must rotate our system into its normal-mode quadratures Q Q Q = (η η η 1 , Π Π Π 1 , η η η 2 , Π Π Π 2 ) T , so that In resonance, the orthogonal transformation (x x x h , x x x c ) T = P (η η η 1 , η η η 2 ) T between local and global modes has the form and the normal-mode frequencies are We must decompose the system's coupling operators x x x α (cf. Eq. (2)) in eigenoperators of H H H S . Taking, for instance, x x x h , one can see that where b b b j is the annihilation operator associated with η η η j . Hence, Also note that, unlike in the local approach, now there are two open decay channels into each bath, at frequencies Ω 1 and Ω 2 , respectively. We can obtain the equations of motion for the normal-mode variables Q Q Q from (9); namely, where we have introduced the new coefficients For completeness, we can rotate Eq. (40) back to the original coordinates q q q, which gives The coefficient matrices M G ,1 and M G ,1 have the same condition number κ, since they are connected via an orthogonal transformation. As a result, they also have the same EPs, since the norms involved in the calculation of κ would remain unaffected (see Eq. (20)). In the central panel of Fig. 2, we see that the 'exceptional lines' of diverging condition number in the frequencycoupling space disappear completely, according to the global master equation. One may question the validity of the GME for such parameters. Namely, in Fig. 2 we set k ∼ λ 2 α whereas, to be on the safe side when it comes to the secular approximation, we should ensure instead that k max α λ 2 α [22]. However, as we show in Sec. 5.4, the GME does lead to the correct steady-state properties at all plotted points save for the fringe |k| 0.1 λ 2 h . Hence the disappearance of the EPs cannot be simply attributed to the global approach breaking down.
Thinking of the local coefficient matrix M L ,1 as the lowest order M it would even be tempting to disregard the EP singularities predicted by the LME as mathematical artifacts, and trust instead in the a priori more physical GME. However, as advanced in Sec. 3.3, this interpretation is not valid here. To see why, we only need to calculate M That is, even in the limit k = 0, each heat bath continues to act globally on both resonators, rather than on the node directly coupled to it. This introduces a 'heat leak' channel which explains why the GME predicts steady state heat flows across resonant nodes even at vanishing coupling, as noted in Refs. [22,23]. Eq. (43) confronts us with the fact that the LME is not, in general, a limiting case of the GME. But then, what is it? What we are after is a microscopic justification of the local equation (29) capable of explaining why it succeeds in capturing the EPs detected experimentally in [1], while the global equation (42a) fails.

Redfield equation
We now resort to the Redfield equation (given in (A.15) in Appendix A) to shed light on the nature of the LME. As outlined in Sec. 3.1, this equation is the last step in the derivation of the GME, before forcing the GKLS form by means of the secular approximation. Consequently, the Redfield equation is always more accurate than the GME. In the context of quantum thermodynamics, however, the main shortcoming of the Redfield equation is its lack of GKLS structure. Without it, it may fail to generate a completely positive dynamics [67,68]. Nonetheless, if used with caution, the Redfield equation can still yield thermodynamically sound results [69][70][71].
To be more precise, (A.15) is a partial Redfield equation [22,39,69]. That is, a simplified version of the equation in which its most rapidly oscillating terms-but not all oscillating terms-are averaged out under a coarse graining of time. Below, we also discard the Lamb shifts, defined in Appendix A. The resulting system is where δ α := 1 2 (∆ 2 ), and In spite of the cumbersome expressions, one can see that the zeroth order term M a) b) c) Figure 3: Heat currents at an exceptional point. a) Steady-state heat currents from the hot bath (red) and the cold bath (blue) as a function of k, according to the GME (solid) and the LME (dashed). b) Phase rigidity calculated according to the LME (dashed red) for the two coalescing eigenvectors, and according to the GME (solid black), for which no eigenvectors coalesce (cf. note [72]). The position of the EP according to (33) is indicated by the dotted grey line across a) and b). c) Transient of the hot heat currentQ h (t) according to the LME (solid red), the GME (dashed blue), and the Redfield equation (dotted green). The steady-state value is indicated by the dot-dashed grey line. The global approach thus deviates from the true heat-flow dynamics, set by the Redfield equation. On the contrary, the local approach remains accurate throughout. Parameters are as in Fig. 2, except for ω = 1. As initial state, we take the tensor product of the thermal state of each oscillator at its local bath's temperature.
is identical to M L ,1 in Eq. (30) (or, equivalently, R † (0) α = L † α ). One only needs to set k = 0 in all dissipative contributions of Eqs. (A.15), which is equivalent to setting k = 0 in the∆ α and δ α terms, as well as in the four matrix elements written out in Eqs. (44c).
We have thus shown that, in resonance, the LME is the low-k limit of the Redfield equation; not of the GME. In fact, it is easy to see that R † (0) α = L † α holds as well out of resonance for any multipartite model, as long as the coupling to the heat baths is mediated by 'frequency filters', such as harmonic oscillators (see Appendix C). In such settings, the LME emerges directly from the Redfield equation. Crucially, unlike the Redfield approach, the LME is guaranteed to generate a completely positive dissipative dynamics. This is why M R,1 and M L ,1 share the same pattern of EPs in parameter space (cf. Fig. 2). In addition, this explains the unlikely success of the LME over the GME at low k, when the secular approximation breaks down [22,23]. In this new light, we see that the LME simply bypasses the secular approximation. This is one of the main results of this paper. Furthermore, note from Eq. (35) that the eigenstates of H H H S do not depend on k. It then becomes clear why the local approach remains accurate even at larger couplings in resonance [23]: A small-k approximation of the dissipators R † α (or G † α ) is valid over a wider range of couplings if the expansion affects only the Bohr frequencies ω, but not the jump operators A A A (α) ω . Conversely, out of resonance, the LME loses validity at larger couplings, since the eigenstates of H H H S are then explicitly dependent on k (see Appendix C).
Next, we will show that, even when the LME and GME do agree in their steady-state predictions, the local approach can generate more accurate heat-flow dynamics. The failure of the GME at capturing the correct dynamics is precisely due to the fact that the secular approximationrequired by the GME-'washes away' relevant dynamical features. We have just illustrated this with the disappearance of the EPs.

Heat currents
The instantaneous rate of heat flow into the system from each of the baths can be obtained by generalising Eq. (12b) tȯ where the transient heat currentsQ α (t) generally do not sum up to zero, nor obey the Clausiuslike inequality (12a); they only do so at very long times, t → ∞. Directly applying this definition using the local, global, and Redfield dissipators giveṡ where Σ (α) ω ω ), and η η η i η η η j and Π Π Π i Π Π Π j are second-order moments. Therefore, in order to evaluate Eqs. (47) we must set up and solve the corresponding linear system of 10 equations with coefficient matrices M 2,D from Eq. (27). This is a rather tedious but, otherwise, straightforward process. As it turns out, the 10 × 10 matrices M L ,2 and M R,2 have the same exceptional points than their 4×4 first-order counterparts M L ,1 and M R,1 . We defer details to Appendix B.
In Fig. 3 we tune the parameters to be at an exceptional point of the system according to Eq. (33), and plot both steady-state and transient heat currents. For the chosen parameters, LME and GME agree in their steady-state predictions to a very good approximation, see Fig. 3a). However, local and global dynamics do differ significantly at finite time. In contrast, the Redfield equation agrees with the LME at all times, see Fig. 3c). This suggests that, within its error bars, the local approach can be superior to the global one when studying the thermodynamics of multipartite systems with weak internal couplings.

Conclusions
We have analysed one of the central problems in quantum thermodynamics. Namely, the modelling of heat flows across open systems with quantum master equations. We have shown that the two most common approaches-the local and the global master equations-can make very different predictions. Firstly, our results illustrate that the local approach succeeds at capturing dynamical features, in the form of exceptional points, that escape the global master equation. Secondly, we find that, when considering degenerate multipartite open quantum systems with weak internal coupling, the LME also yields much more accurate heat-flow dynamics than the GME, even when both agree in the steady-state.
Furthermore, we have shown that the LME follows directly from the more accurate Redfield equation, and is generally not a weak-coupling limit of the GME. This is always the case for any multipartite weakly-interacting open quantum system that couples to the environment(s) via single-frequency contacts, such as a qubit or a harmonic oscillator. Therefore, for such systems, the LME emerges as an accurate and computationally efficient alternative to the Redfield equation. It proves to be superior to the GME and, in contrast to the Redfield equation, it does guarantee positivity.
These results have profound consequences for quantum thermodynamics. Namely, modelling heat flow in a quantum thermal machine with the local approach instead of the global one, could make a sizeable difference in the predicted heat transfer in, e.g., any thermalising stroke of a finite-time thermodynamic cycle. This could result in a radically different assessment of both performance and power output. However, it is important to remember that the local approach has a limited range of validity. Specifically, it is unsuitable for open systems with strong internal couplings (i.e., large k), and it can lead to unphysical results at odds with thermodynamics. Finding an accurate and scalable master equation for such scenarios still remains an open challenge.
[48] More precisely, the Hamiltonian considered in the analysis of the experiment, i.e., c a a a c + ω h a a a  †  h a a a h + κ (a a a  †  c a a a h + a a a c a a  where H H H SB stands for the dissipative interactions between system (S) and bath (B). Note that we are transferring the magnitude of these interactions into λ, which means that the rescaled H H H SB is now O(1), unlike in the main text (cf., e.g., Eq. (3)).
Our Assuming that λ is small enough so that all O(λ 3 ) terms or smaller can be neglected will be our first key approximation. This will lead to a second-order quantum master equation after taking the time derivative; namely, Notice, however, that this master equation time-non-local and, still, of little practical use. To overcome this problem, we make two additional assumptions: First, we require the initial state to include no correlations between system's and bath's degrees of freedom. Namely Hence, we may replace the upper integration limit in (A.10) by infinity without substantial error This is often referred-to as the Markov approximation and the resulting equation, as (Markovian) Redfield master equation. However, the use of the term 'Markovian' can be problematic since the dynamics generated by this equation is, in general, not even positive [68]. This means that the corresponding dynamical map might not be divisible, and lack of divisibility is often associated with 'non-Markovianity' [73].
Before the final step in our derivation, we must transcribe the shorthand tr BL (t)L(t − s)˜ ˜ ˜ (t) ⊗ B , which gives We may always The explicit form of the interaction-picture system operatorS S S(t) can be found easily by exploiting the decomposition from Eqs. (11) in the main text. Namely, given the properties of the jump operators A A A ω it is easy to see thatS S S(t) = ω A A A ω e −iωt . Putting together all the pieces and the corresponding real and imaginary parts are Γ(ω) = 1 2 γ(ω) + iS(ω). The imaginary part S(ω) is typically ignored. It introduces two effects-a displacement of the energy levels of ; but also, non-trivial dissipative terms. These are, however, typically very small.
Many terms in Eq. (A.15) can be dropped, since they are fast-oscillating and average out to zero over the time-scale defined by the dynamics of (t) [21]. Namely, we can drop all terms for which ω and ω have the same sign (ω × ω > 0), since these oscillate as e ±|ω+ω |t . Eq. (A.13) is often called 'partial Redfield equation' [22,39] It is now time to abandon the interaction picture undoing the corresponding unitary transformation. This gives us In this case, the ordered vector of the covariances is Q := η η η 2 1 , Π Π Π 2 1 , and the constant vector is given by while the matrix of the dynamics is The∆ j have been introduced in Eq. (41) and, analogously, we define here the coefficients In this case, the evaluation of the condition number reveals no exceptional points, just as in the case of the first-order moments. Hence, the discrepancy between the local and global master equation persists at the level of the second-order moments. Now, we present the evolution of the covariances in the case of the Redfield equation. Again, in normal-mode variables, the dynamics is expressed as The matrix of the dynamics takes the form while the constant vector is In these expressions, we have defined the new coefficients Calculating the condition number of the corresponding matrix of eigenvectors, one can readily confirm, once again, the presence of the exceptional points along the exact same 'exceptional lines' (33).
It is worth noting that, since the system considered is Gaussian, these singularities will be present, according to the LME and the Redfield equation, in any n th -order moment. This is so, because any higher-order correlation functions of a Gaussian system can be cast as a combination of first-and second-order moments, and these simultaneously display exceptional points. Conversely, no moments of any order will ever pick up exceptional points according to the global description.

(C.2)
Namely, in the non-resonant case both the eigenvectors and the Bohr frequencies of H H H S depend explicitly on the internal coupling strength k.
Moving now to the Redfield dissipators (A.15) we can write it out in adjoint form as To find the zeroth order term R † (0) , we must simply set k = 0 in P and the normal-mode frequencies.
We thus find that the normal modes rotate back to the local coordinates x x x α and Ω α collapse into the bare frequencies ω α . Importantly, as a result, frequency Ω (0) 1 = ω h will only appear in the hot dissipator R † h and Ω (0) 2 = ω c will be only linked to R † c . Therefore, the double sums in Eq. (C.3) directly transform into the local expression (14) which in the specific case of the system studied is This will generally always be true-regardless of H H H S -provided that the system operator coupling to each local bath enables only one transition at some specific frequency. This is the case, for instance, for all thermal devices that couple to heat baths through 'frequency filters', e.g., harmonic oscillators or transitions involving only two-levels [45].
In the same way, it is easy to show that the relationship between the global and the local dissipators, proposed in [42], remains valid in the out-of-resonance case. The adjoint global dissipator takes the form Again, to find the zeroth order term of the k-expansion of the dissipator, we set k = 0 in the matrix P. Therefore, the global jump operators in Eq.