Testing the validity of the local and global GKLS master equations on an exactly solvable model

When deriving a master equation for a multipartite weakly-interacting open quantum systems, dissipation is often addressed \textit{locally} on each component, i.e. ignoring the coherent couplings, which are later added `by hand'. Although simple, the resulting local master equation (LME) is known to be thermodynamically inconsistent. Otherwise, one may always obtain a consistent \textit{global} master equation (GME) by working on the energy basis of the full interacting Hamiltonian. Here, we consider a two-node `quantum wire' connected to two heat baths. The stationary solution of the LME and GME are obtained and benchmarked against the exact result. Importantly, in our model, the validity of the GME is constrained by the underlying secular approximation. Whenever this breaks down (for resonant weakly-coupled nodes), we observe that the LME, in spite of being thermodynamically flawed: (a) predicts the correct steady state, (b) yields the exact asymptotic heat currents, and (c) reliably reflects the correlations between the nodes. In contrast, the GME fails at all three tasks. Nonetheless, as the inter-node coupling grows, the LME breaks down whilst the GME becomes correct. Hence, the global and local approach may be viewed as \textit{complementary} tools, best suited to different parameter regimes.

Interestingly, one may use Eq. (1) to model a continuous (quantum) thermodynamic cycle [7,8]. By coupling the open system (i.e. the working substance) to various heat baths at different temperatures and possibly also to a periodic external drive, a stationary non-equilibrium state builds up. The direction of the corresponding steady-state heat currents may be controlled by suitably engineering the spectrum of the working substance. Hence, we can speak of 'quantum heat engines' or 'quantum compression/absorption refrigerators' [9], which have attracted a lot of attention in recent years [10][11][12][13][14].
In such quantum heat devices, the stationary incoming heat currents {Q α } and the power output −P are defined as [7] where ρ ρ ρ ∞ is the steady state of the working substance, and D α denotes the GKLS dissipation super-operator associated with bath α.
Owing to properties (i) and (ii) above, the stationary heat currentsQ α satisfy the relation ∑ αQα /T α ≤ 0, which is the Clausius inequality. In other words, addressing the dynamics of quantum heat devices with GKLS quantum master equations ensures thermodynamic consistency.
When modelling open quantum systems made up of multiple weakly-interacting parts coupled to local environments, it is commonplace to build the corresponding master equation by simply adding the local dissipators for the relaxation of each individual component (ignoring their coherent interactions). That is, for a multipartite system with Hamiltonian H H H = ∑ j h h h j + kV V V , where V V V contains all the internal couplings (of strength k), one would write 1 Although Eq. (3) is in GKLS form, property (ii) ceases to hold, as the dissipators D (k=0) α do not match the Hamiltonian H H H, but rather the non-interacting ∑ j h h h j . Consequently, describing heat transport with the local master equation (3) may lead to thermodynamic inconsistencies: Heat could, for instance, flow against the temperature gradient [15], or non-vanishing steadystate heat currents could be present even if all reservoirs are set to the same temperature [16].
These observations, strongly advise to follow the standard procedure to consistently obtain the correct global dissipators D α [3]. However, doing so may become particularly challenging when dealing with large systems, e.g. long harmonic or spin chains. Moreover, in such cases the capital assumption that the dissipation time scale is by far the largest in the problem is likely to break down as the spectrum of the system becomes denser [17]; Eq. (1) would then lack a microscopic justification. These difficulties explain the popularity of simple approaches based on weak internal coupling approximations such as Eq. (3) [17,18]. In this paper we wish to put such local approaches to the test.
In particular, we choose an exactly solvable model consisting of a two-node harmonic chain weakly coupled on both edges to two heat baths at different temperatures. Our system is set up so that, when the inter-node coupling strength becomes comparable or smaller than the node-baths dissipative couplings, the secular approximation underlying Eq. (1) may break down. This allows us to gauge to which extent the local master equation (LME) remains an accurate description. Interestingly, we find that the local approach yields an excellent approximation to the steady state, the stationary heat currents, and the asymptotic quantum and classical correlations, in the regime of parameters in which the global master equation (GME) fails even qualitatively. More generally, it follows that heat conduction through arbitrarily large harmonic chains can be correctly modelled within the local approach always provided that the internal couplings are sufficiently weak. The present work thus adds to the efforts of Refs. [15][16][17][18][19][20][21][22][23][24] to clarify the dos and don'ts of modelling heat transport through multipartite open quantum systems. This paper is structured as follows: In Sec. II A we outline the steps of the microscopic derivation of the GKLS quantum master equation. We then proceed to derive and solve such an equation for our specific model in Sec II B. The alternative local master equation is obtained in Sec. II C. Before proceeding to benchmark both approaches, in Sec. III we sketch how the exact steady-state solution of the system may be obtained by solving the quantum Langevin equations. We then devote Sec. IV to present and discuss our results. Finally, in Sec. V we summarize and draw our conclusions.

II. DERIVING MARKOVIAN MASTER EQUATIONS
A. The model, the Markovian master equation and its steady state We will consider a two-node 'quantum wire' (see Fig. 1) consisting of mechanically-coupled harmonic oscillators with bare frequencies ω c and ω h and coupling strength k > 0. Each node will be weakly connected to a bosonic bath, i.e. an infinite collection of uncoupled harmonic modes in thermal equilibrium (at temperatures T c < T h ). The total Hamiltonian may be cast as where the masses of the nodes have been set to m c = m h = 1, and the constants g α,µ stand for the coupling strength between node α and each of the environmental modes (α, µ). Also, in all what follows we shall seth and the Boltzmann constant k B to 1. We will refer to the first three terms in the right-hand side of Eq. (4) as the free ( as opposed to the last term H H H int , which describes the system-baths interaction. For later convenience, we shall also introduce the notation B B B α := ∑ µ g α,µ x x x α,µ . We will group the system-baths cupling constants in the spectral density functions defined as J α (ω) := π ∑ µ g 2 α,µ 2m µ ω µ δ (ω −ω µ ). In particular, we will choose 1D baths with the Ohmic form where Λ is a high-frequency cutoff (max{ω c , ω h } Λ) and the parameter λ captures the dissipation strength. Note that the bath operators B B B α are thus O(λ ).
For completeness, we will now briefly sketch a simple procedure to obtain the standard second-order Markovian generator for the reduced dynamics of the system (see Ref. [3] for full details). Let us take the Liouville-von Neumann equation in the interaction picture whereH H H int (t) := e iH H H 0 t H H H int e −iH H H 0 t t and [·, ·] stands for a commutator. Formally integrating Eq. (6) and assuming that the initial condition is such that tr{ρ ρ ρ(0)H H H int } = 0 yields the following equation of motion for the system: Here,σ σ σ := tr Bρ ρ ρ and tr B {· · · } denotes trace over the baths. We will now assume that the dissipation strength λ is so weak that when starting from a factorized initial condition ρ ρ ρ 0 := σ σ σ (0) ⊗ τ τ τ c ⊗ τ τ τ h the propagated stateρ ρ ρ(t) σ σ σ (t) ⊗ τ τ τ c ⊗ τ τ τ h remains approximately factorized at all times. τ τ τ α∈{c,h} are thermal states for the hot and cold bath. We will also replaceσ σ σ (s) inside the integral in Eq. (7) byσ σ σ (t), thus making it time-local. The change of variables s → t − s yields the Redfield equation [25,26] Notice that the resulting state ρ ρ ρ(t) does keep a memory of the initial condition ρ ρ ρ(0) and hence, Eq. (8) is non-Markovian. However, provided that the integrand above decays sufficiently fast, one might set t → ∞ in the upper limit of integration, which is referred-to as Born-Markov approximation. This approximation is justified whenever λ 2 min{T, Λ}. A further step still remains to be undertaken in order to bring the resulting 'Markovian Redfield' master equation to the canonical GKLS form-the secular approximation. Let us examineH H H int (t) more closely. One may always decompose where γ α (ω) = 2Re ∞ 0 ds e iω s tr B {B B B α (t)B B B α (t − s)}. Note that we are completely ignoring Im ∞ 0 ds e iω s tr B {B B B α (t)B B B α (t − s)}, which would eventually lead to a mere shift (of order λ 2 ) on the energy levels of the Hamiltonian (Lamb shift) [3]. The secular approximation consists in time-averaging Eq. (9) over an interval of the order of the dissipation time T D ∼ λ −2 . All terms with ω = ω above can then be discarded provided that they oscillate fast as compared with T D . Returning to the Schrödinger picture, finally leaves us with the GKLS quantum master equation where {·, ·} + stands for anti-commutator. In the next section, we will concentrate in obtaining the specific form of the operators L L L ω α for the Hamiltonian in Eq. (4). Because Eq. (4) is overall quadratic in position and momenta, the steady state will be Gaussian and thus, fully characterized by its first and second order moments [27]. In fact, one can easily see that X X X α = P P P α = 0, where · denotes stationary average. As a result, the steady state will be specified only by the 4 × 4 covariance matrix, with elements [Γ] kl := 1 2 {R R R k ,R R R l } + , with R R R = (X X X c ,P P P c ,X X X h ,P P P h ) T .
Since we wish to calculate the covariances [Γ] kl rather than the state σ σ σ , it will be more convenient to work with the adjoint master equation [3] which, for an arbitrary system observable in the Heisenberg B. The global master equation The first step to derive consistent L L L ω α operators will be to rotate H H H S into its normal modes. These are where the angle ϑ is and, in turn, The corresponding normal-mode frequencies write as After this transformation, the Hamiltonian (4) rewrites as where Π Π Π s = dη η η s / dt. By writing η η η s = (a a a s + a a a † s )/ √ 2Ω s (with a a a s and a a a † s being annihilation and creation operators on mode Ω s ) one can see that X X X α = L L L Looking back to the right-hand side of Eq. (9), we see that there are 16 terms associated with 5 different open decay channels, oscillating as e i(ω −ω)t at frequencies |ω − ω| = {0, 2Ω + , 2Ω − , Ω + + Ω − , Ω + − Ω − }. Provided that the nodes are sufficiently detuned, i.e. δ ω λ 2 , the secular approximation is guaranteed to be valid for any value of the coupling k. However, if ω h − ω c became comparable or smaller than the dissipation strength λ 2 , there would be no justification to discard the non-secular terms oscillating at Ω + − Ω − when k becomes very small. Indeed, defining R ± := 2k ± 4k 2 + δ 4 ω one may write as From Eq. (17) it is clear that for the secular approximation to be valid the dissipation rate must be such that which, in the limit of resonant nodes simplifies to λ 2 k/ω c . Hence, we can anticipate that Eq. (10) will fail to describe nearly resonant weakly coupled nodes, which is precisely the regime in which we shall focus our analysis. The only additional ingredient required to build Eq. (11) are the decay rates γ α (±Ω ± ). A direct calculation leads to where n α (ω) := (e ω/T α − 1) −1 is the bosonic occupation number for frequency ω at temperature T α . Note that γ α (−ω) = exp (−ω/T α )γ α (ω). Combining all the above, and after tedious but straightforward algebra, we can obtain a closed set of equations of motion for the covariances 2 η η η 2 ± , Π Π Π 2 ± , and {η η η ± ,Π Π Π ± } + . Note that · denotes here instantaneous average.
where the following notations have been introduced Below, it will be convenient to break down each of these coefficients into the sum of its two constituent terms, as −Ω + = cos 2 ϑ γ c (−Ω + )/(2Ω + ), will also be employed later on. The stationary solution to Eq. (20) is simply η η η 2 ± = −Σ ± /(2∆ ± Ω ± ), Π Π Π ± = −Ω ± Σ ± /(2∆ ±), and {η η η ± ,Π Π Π ± } + = 0, so that the non-zero elements of the asymptotic covariance matrix in the original quadratures reads where Finally, the steadystate heat currents can be written aṡ Using Eqs. (21) we can cast (23) asQ from where it is clear that T h > T c entailsQ G h = −Q G c > 0; that is, heat always flows from the hotter bath into the colder one. 2 It is indeed enough to consider the equations of motion for the mode occupation numbers a a a † ± a a a ± (which are decoupled), to fully solve the dynamics.

C. The local master equation
Recall from Sec. I that, while the local master equation looks formally identical to the GME, the choice of operators L L L ω α in the local approach is not consistent with the Hamiltonian H H H S . As already advanced and provided that the coupling k is weak, one could derive two independent local dissipators D (k=0) α , acting on the cold and hot nodes separately, to then construct an approximate equation of motion such as dσ σ σ / dt −i[H H H S ,σ σ σ ] + ∑ α∈{c,h} D (k=0) α σ σ σ , as an alternative to Eqs. (20). One might argue that this is a convenient strategy whenever finding all energy eigenstates of the full interacting Hamiltonian is hard, as these are required to write the decomposition {L L L ω α } of the system operator coupled to each bath [18]. In our simple example, however, the local approach leads to a more complicated dynamics than the global one-all 10 independent covariances are needed in order to obtain a closed set of equations of motion.
Specifically, within the local approach one decomposes X X From (25), the equations of motion for the elements of the corresponding covariance matrix Γ L are found to be {X X X h ,P P P h } + + {X X X c ,P P P c } + − ν 2 c X X X c P P P h − ν 2 h X X X h P P P c + 2ω α , and the angled brackets denote again instantaneous average. The stationary solution of Eq. (26) is cumbersome but the steady-state heat currents can be compactly cast aṡ As anticipated above and unlike Eq. (20), Eq. (26) does not necessarily yield a thermodynamically consistent steady state: One could even encounter striking situations for whichQ L h = −Q L c < 0 for T h > T c orQ L α = 0 for T h = T c , as illustrated in [15,16].

D. Comment on the general validity of the local approach for modelling heat transport under weak internal coupling
In spite of its thermodynamic inconsistencies, as it was pointed out in Ref. [18] the LME (25) can be formally understood as the lowest order in the perturbative expansion D α = D (0) . The LME (25) would therefore be correct up to O(λ 2 k) and any thermodynamic inconsistency encountered should fall within this 'error bar'.
Note that the GME is itself a perturbative master equation which neglects corrections of order O(λ 3 ) and below [28]. However, it is guaranteed to give rise to thermodynamically consistent steady-state heat currents [7,10], as it enjoys the GKLS form (cf. Sec. I). Interestingly, it is the secular approximation which endows the GME with thermodynamic consistency: The Markovian Redfield equation (9), i.e. the previous step in the derivation of Eq. (10), is known to break positivity [28] and caution must be exercised when using it [29].
Coming back to our problem of describing heat transport in the limit of quasi-resonant weakly-coupled nodes, notice that the secular approximation is not problematic when invoked in the derivation of Eq. (25). Indeed, the operators L L L ω α may be expanded as a power series in k in Eq. (9). At the zeroth order in k, each heat bath would contribute to the right-hand side of Eq. (9) with one non-oscillatory secular term and two fast-rotating non-secular terms at frequencies ±2ω α . These may be safely averaged out provided that ω α λ 2 and regardless of the detuning between the nodes. Consequently, and unlike Eq. (20), the LME should correctly describe the stationary properties of our system when k/ω c λ 2 .
More generally, one can claim that energy transport through an arbitrarily long harmonic chain is correctly captured by a LME within its range of validity; that is, whenever the inter-node couplings are weak. The claim can be made extensive to heat fluxes on spin chains, which were already addressed in Ref. [17] via a perturbative master equation relying on 'weak internal couplings', precisely in order to bypass the problems created in the GME by the secular approximation.
Finally, let us note that a natural alternative to the LME in our problem would be to incorporate the problematic decay channel of frequency Ω + − Ω − into the GME, thus arriving to a partial Markovian Redfield equation (cf. Appendix). However, scaling up the system in the number of nodes would quickly render this approach too involved to be practical.

III. EXACT NON-EQUILIBRIUM STEADY STATE
The steady state for an all-linear model can also be found exactly by solving the corresponding quantum Langevin equations [30][31][32]. In this section, we will limit ourselves to outline the procedure to calculate the stationary covariances for our particular problem, while full details on its application to similar settings can be found in e.g. Refs. [33][34][35][36][37].
To begin with, we must mention that the bare frequencies of the nodes need to be shifted so as to compensate for the distortion caused by the system-baths interaction. This eventually allows to recover the correct high temperature limit [31]. Hence, in the reminder of this section, we shall make the replacement ω 2 α →ω 2 α , whereω 2 α := ω 2 α + ∑ µ g 2 α,µ /(m α,µ ω 2 α,µ ). For our choice of spectral density (5), the shift amounts simply to π −1 ∞ 0 dω J(ω)/ω = λ 2 Λ. Starting from Eq. (4), one may write the Heisenberg equations of motion for all degrees of freedom. Formally solving for x x x α,µ and inserting the result into the equations for X X X α yields the quantum Langevin equations These are the equations of motion for two coupled harmonic oscillators, each of which is perturbed by the noise F F F α (t) and relaxes according to the dissipation kernel χ α (t). Specifically, these are defined as The only assumption that we will make in order to find the exact steady state is, once again, that system and baths are initialized in the factorized initial condition ρ ρ ρ 0 = σ σ σ (t 0 ) ⊗τ τ τ c ⊗τ τ τ h . We shall also take t 0 → −∞. Let us first concentrate on the (stationary) covariance 1 2 {X X X α (t),X X X α (t)} + , which may be written in terms of the Fourier transformX X X α (ω) := ∞ −∞ dt X X X α (t)e iωt as In turn,X X X(ω) α can be directly found after Fourier-transforming Eqs. (28), which yields From Eq. (29a), one can show that 1 2 where the Dirac delta δ (·) is not to be confused with the Kronecker delta δ α,α . Consequently, the integral in Eq. (30) for e.g. α = α = c writes as where we are exploiting the fact that our J(ω) is an odd function. The position-momentum and momentum-momentum covariances are readily obtained as e.g. 1 2 {P P P α (t ),X X X α (t )} + = 1 In order to calculateχ α (ω) it is useful to note that Imχ α (ω) = J(ω)Θ (ω) − J(−ω)Θ (−ω), and that Reχ α (ω) and Imχ α (ω) are related through the Kramers-Kronig relation where P indicates Cauchy principal value. For our choice of spectral densityχ h (ω) =χ c (ω) = λ 2 Λ 2 /(Λ − iω) which, combined with Eqs. (5), (31), and (32), allows us to compute all the elements of the exact steady-state covariance matrix Γ. Finally, following Refs. [38,39], we can cast the exact steady-state heat currents aṡ Both the steady state covariances and the corresponding heat currents can be seen to perfectly coincide with those obtained from the Markovian Redfield equation derived in the Appendix, always provided that the Born-Markov approximation holds.

A. Steady state and stationary heat currents
In this section we will compare the steady states and the stationary heat currents predicted by the global, local, and exact approaches. We shall be especially interested in setting up the wire with quasi-resonant nodes (δ ω λ 2 ) so as to confirm our intuition that the LME can succeed in describing the system when the GME breaks down (cf. Sec. II D).
On the contrary, if the detuning is set to δ ω λ 2 , the steady state of the GME can be seen to disagree with the exact solution when the inter-node coupling k/ω c approaches or falls below the dissipation strength λ 2 [cf. Fig. 2(b)]. Recall that this is entirely due to elimination of the non-secular decay channel at frequency Ω + − Ω − (cf. Sec. II B). Importantly, the LME is still valid so long as k/ω c ω c , regardless of the breakdown of the secular approximation. Eventually, as k decreases further, the nodes effectively decouple, and the GME correctly predicts a steady state made up of two uncorrelated thermal modes.
One can also make use of Eqs. (23), (27), and (34) to compare the steady-state heat currents. Once again, under large detuning δ ω , both the local and global approach are in good agreement with the exact solution (vanishingly small heat currents), except for when the inter-node coupling becomes comparable to the node frequencies, which invalidates the LME. Interestingly, in Fig. 1(c) we can see that the local approach does indeed violate the second law of thermodynamics by predicting heat transport against the temperature gradient (i.e.Q h = −Q c < 0) for any k [15]. The magnitude of this violation, however, loosely falls within the 'error bars' O(λ 2 k) [18] of the LME.
On the other hand, Fig. 1(d) shows again a situation in which δ ω λ 2 . Remarkably, we observe that the global approach largely overestimates the magnitude of the steady-state heat currents where F(Γ, Γ G ) falls below 1 (i.e. in the grey area). The LME, however, yields a quantitatively good estimate in all the range of parameters for which it is valid.
We have thus illustrated that the breakdown of the secular approximation may render the predictions of the global master equation qualitatively wrong, while the local approach, in spite of its thermodynamic inconsistency, proves to be an accurate working tool within its range of applicability.

B. Steady-state correlations
As we shall now see, the GME also fails qualitatively in assessing the node-node correlations (both classical and quantum) when the secular approximation breaks down. This is not the case for the LME.  2. (color online) (top row) Uhlmann fidelity F between the exact steady state Γ and the approximations Γ G and Γ L calculated wihtin the global (solid) and local (dashed) approach, as a function of the coupling k at fixed dissipation strength λ 2 = 10 −3 . In (a) frequencies and temperatures were set to ω h = 2, ω c = 1, T h = 3, and T c = 2, so that δ ω λ 2 and the secular approximation is justified. Hence, the global GKLS equation is in perfect agreement with the exact result. The LME starts to break down around k ∼ 0.1, i.e. when the internode coupling becomes comparable to the node frequencies. In (b) the nodes are quasi-resonant (ω c = 1 and δ 2 ω = 2 × 10 −6 ), while the temperatures are the same as in (a). Due to the breakdown of the secular approximation, the global GKLS equation becomes unreliable. The shaded grey area corresponds to 1 − F(Γ, Γ G ) ≥ 10 −4 . In contrast, the LME remains accurate in that regime of parameters. (bottom row) Stationary incoming heat currents from the hot (red) and cold (blue) baths, as given by the global (thin solid), local (dashed), and exact (thick transparent) approaches. The parameters in (c) are the same as in (a). As it can be seen, the LME violates the second law of thermodynamics predicting reversed heat currents for all k. Finally, the parameters in (d) are the same as in (b). In the shaded grey region, in which the secular approximation breaks down, the GME greatly overestimates the magnitude of the steady-state heat currents, while the LME perfectly follows the exact result. For all four plots Λ = 10 3 . Recall that we work in units of m c = m h =h = k B = 1.
We measure the total correlations between the 'cold' and 'hot' nodes of the wire by means of the quantum mutual information I(σ σ σ c h ) := S(σ σ σ c ) + S(σ σ σ h ) − S(σ σ σ c h ), where S(ρ ρ ρ) = − tr {ρ ρ ρ logρ ρ ρ} is the von Neumann entropy and σ σ σ α := trᾱ σ σ σ c h stands for the reduced state of node α (the subindices ' c' and 'h' are added to the starionary state of the wire to emphasize its bipartite nature). The von Neumann entropy of an n-mode Gaussian state can be written as [42] S(Γ) = ∑ n j=1 where the ν j are the n symplectic eigenvalues of the generic 2n × 2n covariance matrix Γ. These can be obtained form the spectrum {±iν 1 , · · · , ±iν n } of J −1 Γ. For Γ to be physical, the symplectic spectrum must satisfy ν j ≥ 1 2 . In our case, note that e.g. the single-mode covariance matrix Γ c results from retaining only the first two rows and columns of the two-mode covariance matrix of the full system, i.e. those related to the 'cold quadratures' {x x x c , p p p c }.
As we can see from Fig. 2(a) the inter-node correlations can be both overestimated and underestimated by the global master equation whenever the secular approximation fails. In contrast, the LME assesses I faithfully. Note from Eq. (22) that the stationary covariances x x x c p p p h and p p p c x x x h (i.e. Γ 14 and Γ 23 ) are neglected in the global approach. Indeed, it is easy to see from the corresponding Markovian Redfield equation (cf. Appendix) that these covariances are related to the excluded non-secular term at frequency Ω + − Ω − . From Fig. 2(b) we observe that the deficit in total quantum correlations predicted by the GME around k/ω c λ 2 in Fig. 2(a) 2(b), that the peak in the total correlations at lower k is explained by the fact that the GME overestimates x x x c x x x h ; once again, within the region in which the secular approximation breaks down. In (c) we plot the steady state inter-node entanglement, as quantified by the logarithmic negativity E N within the global (solid black), local (dashed black), and exact (thick transparent blue) approaches. The observation of non-vanishing asymptotic entanglement requires large ratios ω/T α and very large coupling strengths k. Unfortunately, this prevents entanglement from being observed in the problematic region k/ω c λ 2 . Interestingly, the LME predicts a saturation of E N for large k, which is anyway far beyond its range of applicability. In (c) ω c = ω h = 10, T c = 1, T h = 2, λ 2 = 10 −3 , and Λ = 10 3 .
It is possible to split the total correlations into a quantum and a classical share (blue dotted and red dashed lines in Fig. 2(a), respectively). We will say that a bipartite quantum state ρ ρ ρ AB has quantum correlations with respect to B if there exists no local measurement on B that leaves the marginal of A unperturbed. This notion of quantumness of correlations is captured by the discord Q ← (ρ ρ ρ AB ) [43,44]. Given a complete set of projectors {Π Π Π B j } on B, ρ ρ ρ A| j := tr B {Π Π Π B j ρ ρ ρ AB Π Π Π B j } denotes the post-measurement marginal of A conditioned on the outcome j, occuring with probability p j = tr{Π Π Π B j ρ ρ ρ AB }. Note that discord is not symmetric, i.e. the quantumness of correlations as revealed by measurements on B need not coincide with the quantumness of correlations as revealed by measurements on A.
Note as well that, due to the explicit minimization over all local measurements on B, the evaluation of Q ← is often very challenging. Luckily, restricting the optimization to the set of Gaussian positive operator valued measurements makes it possible to obtain a closed formula for two-mode Gaussian states (see Ref. [45][46][47] for full details). The difference between the total correlations and the quantum discord is referred-to as classical correlations C ← (ρ ρ ρ AB ) := I(ρ ρ ρ AB ) − Q ← (ρ ρ ρ AB ). As shown in Fig. 2(a), both quantum and classical correlations behave very similarly to the mutual information within the global approach. This is not the case, however, for the LME [cf. inset in Fig. 2(a)]: at large coupling strengths (i.e. beyond its range of validity) the local approach may overestimate the amount of quantum correlations present between the nodes, although at sufficiently large couplings, all correlations are largely underestimated.
Finally, we may want to look at the inter-node entanglement [48]. Entanglement is a somewhat stronger form of quantum correlation since a state can display non-zero discord and yet be unentangled, but not the other way around. In the case of two-mode Gaussian states, quantum entanglement can be gauged by the logarithmic negativity E N , which writes as [49,50] whereν j are the symplectic eigenvalues of the partially-transposed covariance matrixΓ. This is obtained from Γ by simply changing the sign of all covariances involving e.g. the momentum p p p c and either of the 'hot' quadratures. The buildup of steady-state entanglement requires much larger inter-node coupling k and large ratios ω α /T α as shown in Fig. 2(c). While there is no reason for the GME not to accurately capture the entanglement as k → ∞, the LME wrongly predicts a saturation in the stationary logarithmic negativity in that limit. One can obtain the correct scaling of entanglement at strong coupling from the GME which, for resonant nodes, simplifies to We have studied a simple model for heat transport between two heat baths at different temperatures when weakly connected through a two-node quantum wire. Due to the weak dissipative wire-baths coupling, it is possible to address the problem via second-order Markovian quantum master equations. In particular, we consistently derived the GKLS master equation via a global treatment of dissipation, and found its steady state, the stationary heat currents through the wire, and the asymptotic internode quantum and classical correlations. For comparison, we adopted the popular local approach, which addresses dissipation on each node individually (i.e. ignoring the effects of the inter-node coupling). Since our model is linear, its steady state can be obtained exactly by resorting to quantum Langevin equations. This provided us with means to quantitatively compare the performance of the global and the local approaches.
As expected, we found that the local approach is only valid when the internal coupling between the nodes of the wire is weak. Furthermore, as previously noted, we observed that the local approach does break the second law of thermodynamics [15], although any violations can be bounded with suitably-defined error bars within its range of applicability [18].
Interestingly, our setup allows us to consider very weak internal couplings, comparable with the dissipation strength. In this regime, the crucial secular approximation breaks down if, in addition, the nodes are nearly resonant. As a result, the predictions of the global master GME become qualitatively wrong-the magnitude of the stationary heat currents is largely overestimated, and key features of the correlation-sharing structure are not captured by the GME. On the contrary, the LME does accurately describe the stationary properties of the wire. This agrees with previous observations on the complementarity of GME and LME when describing dynamics [19]. More generally, the usage of the local approach in the treatment of heat transport through arbitrarily long harmonic or spin chains [17] may be justified provided that the internal couplings are weak enough, and always keeping in mind that the predictions of the LME should by accompanied by the corresponding error estimates [18].
In spite of these encouraging observations, the local approach should not be used lightly, especially in quantum thermodynamics. Even though the LME may be an excellent working tool that even outperforms the canonical global GKLS master equation in certain regimes, it might as well lead to qualitatively wrong conclusions, a priori within its range of applicability. For instance, it has been shown that a local modelling of quantum thermodynamic cycles completely fails to account for heat leaks and internal dissipation effects [51,52] that can become dominant in the operation of the device in question. As a result, e.g. intrinsically irreversible models may be wrongly classified as endoreversible. This is a reminder that perturbative equations of motion for open quantum systems must always be handled with care.
NOTE ADDED: During the preparation of this manuscript we became aware of the related work by Patrick P. Hofer et al. [53], where local and global approach are compared in a quantum heat engine model. In order to compensate for the deficiencies of the GME one may simply take into consideration the problematic non-secular term corresponding to the Ω + − Ω − channel. Eqs. (9) and (11)  where the operators L L L ω α are those defined in Sec. II B. In principle, a full set of 10 dynamical variables would be necessary to obtain all steady-state covariances. We shall choose D D D ±± := i(a a a † ± a a a † ± − a a a ± a a a ± ), S S S ±± := a a a † ± a a a † ± + a a a ± a a a ± , D D D +− := i(a a a † + a a a † − − a a a + a a a − ), S S S +− := a a a † + a a a † − + a a a + a a a − , d d d +− := i(a a a † + a a a − − a a a + a a a † − ), s s s +− := a a a † + a a a − + a a a + a a a † − , and n n n ± := a a a † ± a a a ± . As it turns out, the stationary averages of the first six variables vanish (i.e. D D D ±± = S S S ±± = D D D +− = S S S +− = 0), so that we are left with only four relevant observables. The corresponding equations of motion write as d y y y/ dt = B y y y + b, where y y y = (n n n + ,n n n − ,d d d +− ,s s s +− ) T , the non-zero elements of b are given by and the coefficients of the matrix B read