Pushing the limits of the reaction-coordinate mapping

The reaction-coordinate mapping is a useful technique to study complex quantum dissipative dynamics into structured environments. In essence, it aims to mimic the original problem by means of an 'augmented system', which includes a suitably chosen collective environmental coordinate---the 'reaction coordinate'. This composite then couples to a simpler 'residual reservoir' with short-lived correlations. If, in addition, the residual coupling is weak, a simple quantum master equation can be rigorously applied to the augmented system, and the solution of the original problem just follows from tracing out the reaction coordinate. But, what if the residual dissipation is strong? Here we consider an exactly solvable model for heat transport---a two-node linear"quantum wire"connecting two baths at different temperatures. We allow for a structured spectral density at the interface with one of the reservoirs and perform the reaction-coordinate mapping, writing a perturbative master equation for the augmented system. We find that: (a) strikingly, the stationary state of the original problem can be reproduced accurately by a weak-coupling treatment even when the residual dissipation on the augmented system is very strong; (b) the agreement holds throughout the entire dynamics under large residual dissipation in the overdamped regime; (c) and that such master equation can grossly overestimate the stationary heat current across the wire, even when its non-equilibrium steady state is captured faithfully. These observations can be crucial when using the reaction-coordinate mapping to study the largely unexplored strong-coupling regime in quantum thermodynamics.


I. INTRODUCTION
Understanding the dynamics of open quantum systems in structured environments is central to nearly all aspects of quantum research-from modelling the chemistry of biomolecules [1][2][3] , to understanding the thermodynamics of quantum systems 4 , or assisting in the design of nanostructures for quantum-technological applications [5][6][7] .Unfortunately, treating open systems in complex environments is extremely challenging, the main reason being the absence of a clear-cut timescale separation between system and evironmental dynamics 8 .Various tools exist to deal with such problems, including exact path-integral methods [9][10][11] , stochastic Schrödinger equations 12,13 , unitary transformations 14,15 , or Markovian embeddings [16][17][18][19] .Here, we shall focus on the latter; specifically, on the "reaction-coordinate mapping" 20 .
In a seminal paper by Garg et al. 1 a very simple ansatz was put forward for the structure of the environment modulating the rate of an electron-transfer process in a biomolecule.Essentially, it assumes that a distinct collective environmental coordinate-the reaction coordinate (RC)-couples strongly to the donor-acceptor system, which can be thought-of as a two-level spin.In this construction, the combined effect of all other environmental degrees of freedom would merely cause semiclassical friction on the spin-RC composite.It is then possible to view the spin as an open system and work out its dissipative dynamics via, e.g., exact path-integral methods.
Interestingly, the ansatz can be "turned on its head" 17,19,21 and viewed as a Markovian embedding technique.Namely, an arbitrarily complicated environment may be iteratively decomposed by first, extracting a collective environmental coordinate and working out the coupling of the resulting 'aug-mented system' to the remaining 'residual environment'.By repeating this procedure sufficiently many times, one ends up with an open-system model with the simplest friction-like Ohmic dissipation 17,18 , albeit with a much larger system size.Whenever the residual friction is perturbatively small, the problem can be rigorously solved via standard weak-coupling Markovian master equations 22 .This provides a simple route to tackle otherwise intractable open quantum systems, especially when a single iteration of the procedure suffices for the problem at hand.The reaction-coordinate mapping has been applied extensively to open quantum systems strongly coupled to both bosonic 19,[23][24][25][26][27][28][29][30][31] and fermionic [32][33][34] reservoirs.Its relative ease of use and the neat physical picture that emerges from it, in terms of, e.g., system-environment correlation-sharing structure 19,23 , make it particularly appealing as a generalpurpose open-system tool.Unfortunately, relying on perturbative master equations imposes a priori severe limitations on the parameter ranges in which the method can be used.Intriguingly, however, it has resisted benchmarking at finite temperatures over a wide friction range 19,23,32 , which made us wonder where are its true limitations 35 .
In this paper, we set out precisely to "push" the method to the limit, by deliberately taking the forbidden large friction limit in a minimal heat-transport setup.Our biggest advantage is that we work with an exactly solvable model 36 ; we can thus always benchmark the accuracy of the mapping without having to approximate the exact dynamics numerically.Under steady-state conditions, we find that the RC mapping does work accurately even under extremely large friction, in spite of the fact that the underlying master equation breaks down.We also find that overdamped dynamics, re-sulting in strong residual friction, is accurately captured by this method.Importantly, however, when the residual friction is strong and one relies on weak-coupling master equations to compute heat (or particle) currents across the nonequilibrium open system of interest, the results can be completely flawed and yet, appear physically consistent.This observation can have important consequences when using the reaction-coordinate mapping to explore the thermodynamics of strongly coupled nanoscale open systems; verifying that the method approximates the state of an open system correctly is certainly not enough to trust it with the calculation of quantum-thermodynamic variables.
As a by-product of our master-equation analysis of the augmented system subjected to friction, we derive here a (global) Born-Markov secular quantum master equation for a general linear network of harmonic nodes coupled to arbitrarily many equilibrium environments.This generalises the customarily used local master equations applied to quantum transport problems through weakly interacting networks 37 .We also write the ensuing non-equilibrim steady state, and explicit formulas for the corresponding stationary heat currents.Finally, we discuss the dos and don'ts of the often confusing Hamiltonian frequency-renormalisation counter-terms that appear in quantum Brownian motion [38][39][40] , as it is particularly important to use them consistently when performing the reactioncoordinate mapping.
This paper is structured as follows: In Sec.II A, we introduce our simple model and discuss very briefly the reactioncoordinate mapping.In Sec.II B, we provide the general quantum master equation that we shall later apply on our augmented system.Rather than reproducing the standard textbook derivation from the microscopic system-bath(s) model, we limit ourselves to provide here the key steps, and write down instead the full equations of motion explicitly, along with their stationary solutions, and the corresponding steadystate heat currents.In Sec.II C we outline the exact solution of both our original problem and that of the augmented system undergoing (arbitrarily strong) friction.We then proceed to discuss the steady-state (cf.Sec.III A) and dynamical (cf.Sec.III B) benchmarks to the reaction coordinate mapping, commenting both on the approximation to the state of the system and to the stationary heat currents flowing across it.Finally, in Sec.IV, we wrap up and draw our conclusions.

II. THE MODEL AND ITS SOLUTION
A. A two-node non-equilibrium quantum wire

Full Hamiltonian
As already advanced, our model consists of a two-node chain (or "quantum wire") of harmonic oscillators with a linear spring-like coupling of strength k (see Fig. 1), that is Note that here and in what follows, we set all masses to one.We shall also take = k B = 1.The wire is kept out of equilibrium by two linear bosonic baths at temperatures T α .Throughout, α ∈ {h, c} stands for 'hot' or 'cold', i.e., T h > T c .Their Hamiltonians can thus be cast as H H H T α = µ ω µ a a a (α) † µ a a a (α) µ , where a a a (α) † µ (a a a (α) µ ) is a creation (annihilation) operator of bath α in the collective bosonic environmental mode at frequency ω µ .In turn, the dissipative interactions between the wire and the baths are where the quadratures 2ω µ x x x (α) µ a a a (α) † µ + a a a (α) µ and, as usual, the coupling constants g (α) µ make up the spectral densities Importantly, each system-bath coupling H H H diss,α requires us to introduce a renormalisation term in the bare Hamiltonian of the wire H H H w → H H H w + δ δ δH H H w-α , which compensates for the environmental distortion on the system's potential 40 .If we were not to include such terms and let T h = T c = T be arbitrarily large, the exact stationary state would approach instead of the classical limit w (∞) ∼ exp (−H H H w /T ); this should be seen as an important deficiency of the model 38 .Specifically, these extra terms are and the full Hamiltonian of our system is, therefore, We take an Ohmic spectrum for the coupling to the 'hot bath', i.e., J h (ω) = γ h ω θ(ω/Λ h ), where θ(x) is some rapidly decaying function for arguments x > 1, which places an upper bound on the excitation energies.In practice, we choose the algebraic cutoff θ(x) = (1 + x 2 ) −1 , although other choices would not alter our results as long as Λ h is large.Such J h (ω) is referred-to as 'overdamped' in the context of energy transfer in molecular systems 23 .For the coupling of the wire to the cold bath, we take instead the 'underdamped' spectrum which displays a peak around ω 0 , whose height and width are essentially controlled by λ and γ, respectively 41 .The frequency-renormalisation shifts δ α corresponding to these expectral densities are explicitly given by δ h = γ h Λ h and δ c = λ 2 /ω 2 0 .The decay of the environmental correlation functions B B B α (t) B B B α (0) gives an idea of the bath's memory time, and to which extent a simple Markovian relaxation process can be a good approximation to the actual dynamics.Specifically 42 FIG. 1. Sketch of the non-equilibrium quantum wire with nodes at frequencies ω h and ω c and internal coupling k.The dissipative interaction between node ω h and the corresponding (hot) bath, at T h , is characterized by an Ohmic spectral density, e.g., J h (ω) ∼ γ h ω.As a result, the corresponding environmental correlation time is short.Furthermore, the dissipation strength γ h is assumed to be perturbatively weak.On the contrary, the (cold) bath at T c features long-lived correlations due to the structured spectral density J c (ω) = γλ 2 ω/ γ 2 ω 2 + ω 2 − ω 2 0 2 .The resulting dynamics can be mimicked exactly by coupling a reaction coordinate at frequency ω 0 to the system with strength λ.This composite makes up the augmented system which, in turn, couples to a residual reservoir-customarily assumed to be in equilibrium also at T c -via Jc (ω) ∼ γ ω.This guarantees that the residual environmental correlations for the augmented system are short-lived.However, if the 'friction' coefficient γ in J c (ω) is large, so is the residual dissipation strength.Crucially, this clashes with the weak-coupling approximation which underpins any perturbative quantum master equation that could be written for the augmented three-node system.
While a spectral density like our J h (ω) typically leads to very short correlation times, consistent with the Markovian approximation, a spectrum such as (6) can give rise to very long-lived correlations and thus, to a much more complex dynamics.

The reaction-coordinate mapping in a nutshell
To circumvent this problem one may try to exploit the fact that Eq. ( 6) is the effective spectral density for a system which couples indirectly-namely, through a bosonic mode, or reaction coordinate, of frequency ω 0 -to a residual reservoir with a purely Ohmic spectrum 1 ; the coupling between the auxiliary mode and the system being of strength λ (see Fig. 1).Put in other words, the dynamics generated by exactly coincides with that of when the coefficients {g (c) µ } in H H H diss, c correspond to Eq. (6) [by virtue of (3)] and the {g (c) ν } in the sixth term on the right-hand side of Eq. ( 9), to Jc (ω) = γ ω; technically, some suitable cutoff function θ(ω/ Λc ) with would be required, the mapping being exact only in the limit Λc → ∞.Here, tr w amounts to tracing over all degrees of freedom except for the wire.The boldface symbols with tilde correspond to operators completely or partly supported in the residual reservoir; in our case, the quadratures {x x x (c) ν }, the creation and annihilation {ã a a † (c) ν , ã a a (c) ν } operators in modes at frequency ω ν ; and the joint state of the hot bath, the wire, the reaction coordinate, and the residual reservoir ρ ρ ρ(t).Finally, the newly introduced operators X X X RC and P P P RC stand for the canonical degrees of freedom of the RC.Note that we have included as well the renormalisation term δ δ δH H H RC-res arising from the coupling between the RC and the residual bath [cf.Eq. ( 4)].Accessible and rigorous derivations of the equivalence between Eqs. ( 10) and ( 8) can be readily found in the literature 1,17,19,24 .
There is, however, an important caveat regarding the initial condition for the augmented system.It is common practice to assume that the residual reservoir is in equilibrium at temperature T c , just like the original physical bath (see Fig. 1); and to initialise the auxiliary RC in a thermal state at T c , uncorrelated from the rest 19,23 .Note that the dynamics generated by Eqs. ( 8) and (10) In particular, this means that the composite 'RC + residual reservoir' should start instead in a joint thermal state at temperature T c ; that is, ˜ RC + res (0) = T c , which is not of the form RC (0) ⊗ ˜ T c .Hence, there could be large initial correlations between the RC and the residual reservoir, especially at low T c .Importantly, the absence of correlations with the environment is central to the derivation of the most common quantum master equations 42 .Luckily, in many cases of practical interest the residual interactions g(c) µ are sufficiently weak so that the dynamics is faithfully captured under this simple assumption.As we show in Sec.III B below, this is indeed the case when working in the overdamped limit.Furthermore, given its uniqueness 43 , the non-equilibrium steady state (NESS) of our linear wire is always correctly reproduced by the augmented system, regardless of the initial condition for the RC.
Before moving on, let us briefly recapitulate: Our original problem consists of two interacting oscillators locally coupled to two heat baths.The coupling to one of them is of the form (6), which complicates the analysis as it is likely to produce non-Markovian dissipation (i.e., with long memory times).Luckily this precise dissipative dynamics can be exactly mimicked by replacing the problematic thermal contact with one auxiliary oscillator undergoing purely Markovian dissipation.In a suitable parameter range, this 'augmented' three-oscillator model can thus be tackled via a standard master equation (as we do in Sec.II B below), which would allow us to recover the original dynamics by just tracing out the auxiliary coordinate.The "twist" of this paper is that we push such master equation far beyond its range of applicabilitynamely, we allow for a very strong residual dissipation on the augmented system-and benchmark its prediction for the steady state of the wire against the exact stationary solution of the problem.This can always be obtained with the methods outlined in Sec.II C, since our H H H in Eq. ( 9) is fully linear.

The (global) GKLS master equation
We will now outline the derivation of the adjoint quantum master equation for an arbitrary linear network of N harmonic nodes, locally coupled to M baths.This is a Born-Markov secular master equation 42 in the standard Gorini-Kossakowski-Lindblad-Sudarshan (GKLS) form 44,45 .In the present paper we shall only be interested in applying it to a simple 1D chain of three (and, in Sec.III B, also two) harmonic oscillators with heat baths coupled at both ends.Nonetheless, the general equation is of independent interest, as it can be applied to many problems in quantum transport.
It is important to stress that we treat dissipation globally, as opposed to the widespread 'local' or 'additive' approach 37 .That is, we acknowledge that even if each bath couples locally to one node of the network, the ensuing dissipation affects the system as a whole, due to the internal interactions.Indeed, the local approach is known to lead to severe physical inconsistencies 27,[46][47][48][49] .Rigorously, such local equations are only acceptable when understood as either the lowest-order term in a perturbative expansion of a global master equation in the internal coupling strength 50,51 , or as a limiting case of a discrete collisional process [52][53][54] .In any case, addressing dissipation locally is often the only practical way forward in large interacting non-linear open systems-exact diagonalisation of the full many-body Hamiltonian is, otherwise, required.Remarkably, finding, e.g., the NESS, which sets the transport properties of any interacting linear network, with the "plugand-play" stationary solution below [i.e., Eqs. ( 19) and ( 21)] only requires the diagonalisation of the corresponding N × N interaction matrix.
The Hamiltonian of a general linear network can be cast as assuming again that masses are M = 1.Here, X X X and P P P are Ndimensional vectors containing the position and momentum operators of each node, and V is real and symmetric.Let P be the orthogonal transformation that brings (12) into the diagonal form where Ω i j = Ω i δ i j > 0 is a diagonal matrix formed of the normal mode frequencies corresponding to the conjugate variables η η η i , π π π i i∈{1,••• ,N} (i.e., η η η P T X X X).
The standard derivation of a Born-Markov secular master equation 36,42,47 now requires to decompose the 'systemenvironment' couplings [in our case, X X X i for the M nodes coupled to local baths, as per Eq. ( 2)] as eigen-operators of ).These non-Hermitian operators, turn out to be simply where b b b j = Ω j /2 η η η j + i π π π j /Ω j .With these definitions, the equation of motion for an arbitrary Heisenberg-picture (Hermitian) operator O O O(t) under the Born-Markov and secular approximations reads 42 dO O O(t) with {•, •} + denoting anticommutator and decay rates The main appeal of Eq. ( 14) is that it is guaranteed to generate a completely positive and trace-preserving dynamics for the system 44,45 , unlike other frequently used weak-coupling master equations 55,56 .Furthermore, under mild ergodicity assumptions, it admits a unique stationary solution 57 which, in the case of a single environmental temperature T , is the thermal equilibrium state N (t) ∝ exp (−H H H N /T ) 58 .Importantly, this means that no renormalisation needs to be done on the Hamiltonian H H H N to recover the correct equilibrium state in the high-temperature limit.For that reason, when applying Eq. ( 14) to the three-node augmented system, we take as the system Hamiltonian; i.e, we discard the renormalisation terms δ δ δH H H h-w and δ δ δH H H RC-res in Eq. ( 9), corresponding to the thermal contact with the hot and the residual environment, respectively.
However, the term δ δ δH H H w-c is-by construction-part of the augmented system after the reaction-coordinate mapping 19,24 .As we shall see in Sec.III B below, disregarding this latter term in the augmented-system Hamiltonian, e.g., on the basis of δ c being small, can yield the wrong dynamics for the wire at intermediate times, even if the short-time evolution and the steady state are reproduced accurately.

Equations of motion for the covariances
Applying Eq. ( 14) to the symmetrised covariances 1 2 {r r r j (t), r r r k (t)} + [c (N) me ] jk (t), where r r r = (η η η 1 , π π π 1 , • • • , η η η N , π π π N ) T , yields a closed algebra for the 'covariance matrix' of the network c (N) me (t), where the sub-index 'me' stands for 'master equation' and allows to differentiate it from the 'ex' (for 'exact') covariance matrix, that we will compute in Sec.II C below.Specifically, we have together with the asymptotically vanishing covariances (for j k) where For completeness, the equations of motion for the first-order moments η η η j and π π π j are given by Since our Hamiltonian (12) is quadratic in position and momenta, any Gaussian initial state of the network will remain Gaussian at all times.In turn, given that Gaussian states are fully characterised by their first-and second-order moments 59 (that is, r r r j (t) and 1 2 {r r r j (t), r r r k (t)} + ), Eqs. ( 16)-( 18) thus provide a full dynamical description of the problem.Furthermore, since r r r j (∞) = r r r j (∞) r r r k (∞) = 0 for j k, we can concentrate only in Eqs.(16) as far as the NESS is concerned.
Explicitly, this is given by where Σ(Ω j ) We can also define the adjoint dissipation super-operators L † i for each heat bath by rewriting Eq. ( 14) That way, we can cast the stationary heat current flowing from the ith bath into the network as Q(N) i,me 60,61 .In our case, this evaluates to In Sec.III A below, we shall apply the general equations ( 19) and ( 21) to the simple three-oscillator chain making up the augmented system for our quantum wire (cf.Fig. 1), and compare them with the exact stationary state and heat currents (see Sec. II C).In turn, in Sec.III B, we compare the reduced dynamics of the augmented system with the time-evolution of the two-node wire in a parameter regime where Eqs. ( 16)-( 18) are also directly applicable to the original problem.

A note on the underlying approximations
To conclude this section, let us briefly comment on the approximations underlying the microscopic derivation of Eq. ( 14) 42 .First and foremost, it is a second-order perturbative expansion of the exact master equation in the systemenvironment(s) coupling 56 .Therefore, it is only meaningful under the assumption of weak dissipation.In addition, the Markov approximation has been performed by neglecting any memory effects in the dissipative process, since environmental correlations are assumed to be very short-lived.Note that it may well be the case that environmental correlations are indeed short while the dissipation is strong; recall that the bath memory time is essentially determined by the "shape" of the spectral density [cf.Eq. ( 7)].In such situation, the Markov approximation would be valid, but the weak coupling assumption would be violated.
The completely positive GKLS form ( 14) is attained after performing the secular approximation which, in our case requires that all normal-mode frequencies Ω j be well separated as compared to the dissipation rates (i.e., min Once again, this approximation is incompatible with arbitrarily large dissipation rates γ i , but may also be easily violated under weak dissipation 36,50 .For that reason, the full Redfield equation 55,56 -containing all non-secular terms-is often used instead when performing the RC mapping 19,23,24 .As we will see in Sec.III below, even if both the weak coupling and the secular approximation are violated on the augmented system, the two-node reduction of the resulting state may still provide an excellent approximation to the exact steady state of the wire. As a final remark, notice that Eq. ( 14) does not include the so-called Lamb shift term 42 .This is a Hamiltonian-like contribution to the master equation, dissipative in origin.The Lamb shift is often neglected for being a 'small' contribution when compared with the bare Hamiltonian H H H N 24 .It is safe to say that, when working with a GKLS quantum master equation, the Lamb shift is entirely irrelevant for the thermodynamics of steady-state energy-conversion processes 62 .Interestingly, however, when the Redfield equation is used instead (e.g., due to the inadequacy of the secular approximation), the Lamb shift can have noticeable effects 32 .Note that this term is not related to the frequency renormalisation discussed in Sec.II A 1 above.

C. Exact stationary solution
The stationary state of our two-node wire can be obtained exactly, with no other assumptions than a factorised initial state of the form ρ ρ ρ(0) = T h ⊗ w (0) ⊗ T c , and no restrictions on w (0).Importantly, the problem can be solved analytically regardless of the spectral densities at the boundaries.These linear open systems have been extensively studied in the literature [63][64][65][66][67][68][69] , as they are among the few which admit an exact solution under strong dissipation.Full details about the calculation of the steady state and stationary heat currents for the Hamiltonian in Eq. ( 5) were given by González et al. 36 , and here we limit ourselves to outline the key steps.
The exact dynamics of the wire obeys the following quantum Langevin equations 39,70 where α ∈ {h, c} and c h and h c.As we can see, the coherent evolution of the two coupled (and renormalised) oscillators is affected by environmental driving and dissipation (terms on the right-hand side).Importantly, the upper limit of the integral can be extended to infinity by supplementing the dissipation kernel χ α (t) with a Heavisde step function Θ(t) [i.e., χ α (t) → χ α (t) Θ(t)] 71 .Since we are interested in the steady state of the wire, our aim will be to compute the covariance matrix C (2) ex at any finite time t while setting t 0 → −∞.With this in mind, we can now Fourier-transform Eqs.(22), which yields Here, the "hatted" symbols are in the frequency domain, i.e., f (ω) so that the objects we wish to compute are for t = t = t.The position-momentum and momentummomentum covariances can be obtained by differentiating Eq. ( 24), which is equivalent to multiplying the integrand by (−iω ) and (−ω ω ), respectively.To carry out the integration in ( 24) explicitly, we only need the Fourier transform of the dissipation kernels χα (ω) and the power spectrum of the environmental forces 1 2 { F F F α (ω ), F F F β (ω )} .These are given by 71 Im χα where P denotes 'principal value', δ αβ is a Kronecker delta and δ(x) is a Dirac delta.The integration in Eq. (25b) can be readily performed for the overdamped and underdamped spectral densities of interest, i.e., J h (ω Note that Eq. (26a) may also be used for the dissipation kernel of the residual bath acting on the augmented system, by merely replacing γ h with γ and taking a large cutoff.Summing up, Eqs. ( 24)-( 26) are all we need to fill in the full stationary 4 × 4 covariance matrix C (2) ex (∞).Although analytic solutions are possible 66,72 , we proceed numerically.In turn, the 6×6 NESS C (3) ex (∞) of the augmented system can be found in a completely analogous way 72 , by just replacing the 'vector of forces' ( F To conclude this section, let us introduce the exact stationary heat currents, for comparison with Eq. (21).A direct calculation shows that the change in the energy of our wire (or the augmented system) due to dissipative interactions with bath III. DISCUSSION

A. Steady state and stationary heat currents
We are now ready to put the reaction-coordinate mapping to the test.Using Eqs. ( 24)-( 26), we can compute the exact stationary covariance matrix of the original (two-node wire) problem C (2) ex (∞), as well as that of the augmented (threenode) system, C (3) ex (∞).Alternatively, we can look for the steady state of the augmented system according to the GKLS master equation, i.e., C (3) me .Benchmarking the RC mapping thus amounts to assessing how "close" is the relevant 4 × 4 submatrix of C (3) me (∞) to the exact stationary state C (2) ex (∞).We thus need to be able to quantify the distance between two covariance matrices C 1 and C 2 .To that end, we resort to the Uhlmann fidelity 74,75 F (C 1 , C 2 ) which, for arbitrary Nmode Gaussian states with vanishing first-order moments, is This is a meaningful distance measure, since F (C 1 , C 2 ) = 1 only holds if the states are identical, and 0 < F (C 1 , C 2 ) ≤ 1. ex (∞) and the relevant two-node reduction of C (3) me (∞).This is achieved simply by eliminating rows and columns related to the reaction-coordinate variables X X X RC and P P P RC from the 6×6 matrix C (3) me (∞).The abscissa corresponds to the log of the friction coefficient γ of the underdamped spectral density in Eq. ( 6), at the interface between the wire and the cold bath (normalised by the dissipation strength into the hot bath γ h ).Only in the shaded-grey area the fidelity drops below 95%.(dashed) Fidelity between C (3) ex (∞) and C (3) me (∞) as a function of the normalised friction γ/γ h .(b) Steady-state heat currents coming from the hot (red) and cold (blue) baths into the wire as a function of γ/γ h for the same parameters as in panel (a).The solid lines correspond to the exact calculation from Eq. (28a), while the dashed ones result from applying the GKLS master equation to the augmented system, as per Eq. ( 21).The inset is a zoom into the low friction limit and the shaded-grey area is the same as in (a).In both panels, the wire is characterised by ω h = 1, ω c = 3, and k = 0.8; the overdamped Ohmic spectral density at the hot interface, by γ h = 10 −3 and Λ h = 10 3 ; the underdamped spectral density at the cold interface, by λ = 0.9 and ω 0 = 4; and the baths, by temperatures T h = 3.3 and T c = 1.2.

Steady states
In Fig. 2(a) we illustrate our steady-state benchmark for the RC mapping (solid line).Strikingly, we find that the reduction of C (3) me (∞) onto the wire degrees of freedom remains nearly identical to the exact stationary state C (2) ex (∞), even at extremely large residual dissipation strengths γ.In the figure, for instance, the fidelity between the two states falls below 95% only at γ 60.When it comes to the approximations that justify the GKLS master equation ( 14) (cf.Sec.II B 3), this is completely off-limits.Indeed, note that the normalmode frequencies of the augmented system are, in this example, Ω 1 1.31, Ω 2 3.13, and Ω 3 4.02, which renders the secular approximation problematic already at residual dissipations as small as γ ∼ 0.1.More importantly, γ 60 can by no means be considered small and hence, a perturbative expansion of the generator of the dissipative dynamics is out of the question.Our extensive numerics show that this surprising observation is not due to a lucky parameter choice but rather, a generic feature.It is also consistent with the excellent agreement previously reported in other (non-linear) models 19,23,31 , between the reduction of the master-equation-propagated augmented system and the numerical solution to the original problem.
As surprising as this observation may seem, there is nothing contradictory in it-indeed, the GKLS master equation does break down for γ 0.1, which corresponds to log (γ/γ h ) 5. We can see this in Fig. 2(a), when instead of looking at the reduction of C (3) me (∞) onto the wire, we consider the full augmented system and compare it with the exact three-node solution C (3) ex (∞) (dashed line).Specifically, F (C (3) ex (∞), C (3) me (∞)) < 0.95 for γ > 0.15, as expected.We are thus not claiming that Markovian master equations in Lidblad form are generally valid for strong coupling situations.What we find is that non-equilibrium energy transfer processes through open quantum systems in complex environments can be captured faithfully over a much wider parameter range than previously thought, by combining the RC mapping with a GKLS master equation (RC-GKLS mapping).

Steady-state heat currents
Besides faithfully reproducing the NESS of an open quantum system, one would also like to learn about the stationary heat currents that it supports, especially when viewing it as a 'continuous thermal device' for quantum thermodynamics 61 .To do so from the RC-mapped picture, we need to gauge the energy per unit time crossing the boundary between either bath and the augmented system; this can only be achieved by using to the corresponding GKLS dissipators L i [cf.Eq. ( 21)].Under strong coupling, however, these are certainly not valid generators of the dissipative dynamics.A priori, one should thus expect a substantial mismatch between the GKLS stationary heat currents and their exact values in this regime.In Fig. 2(b) we can indeed see that for γ ∼ 60-where C (2) ex (∞) and the reduction of C (3) me (∞) differ only by 5%-the master equation overestimates the heat currents by an order of magnitude, and fails to capture, even qualitatively, their behaviour for larger friction γ.
Note that, for us, resorting to the dissipators is indeed the only feasible way to estimate heat currents; C (3) me (∞) is lacking the key covariances X X X c P P P h and X X X c P P P RC needed to evaluate the dissipative change in the energy of the heat baths [cf.Eq. (28a)].In fact, this has been criticised as one of the most unsatisfactory features of GKLS-type quantum master equations 50 .Alternatively, one could think of waiving the secular approximation to work instead with a Redfield master equation 55,56 .Although the aforementioned covariances would then cease to be zero, the calculation would continue to yield quantitatively wrong results at very large γthis time simply due to the breakdown of the basic weakcoupling assumption.Ultimately, however, the Redfield approach might improve the GKLS results under moderate resid- ual dissipation 24,32 .Therefore, even in the light of the promising observation made in Sec.III A 1 above, great care must still be taken when relying on the RC-GKLS mapping to discuss quantum thermodynamics under non-Markovian dissipation.

B. Dynamics
One can now ask whether the resilience of the RC-GKLS mapping to strong residual dissipation is exclusively a steadystate feature, or whether it holds throughout the entire dissipative evolution.Unfortunately, we do not have any exact dynamical benchmark-at most, we are able to solve for the steady state of the exact Eq. (22).To circumvent this problem, we chose parameters so that the original two-node problem can be described via a GKLS quantum master equation.
In particular, we scale ω 0 and λ in the structured spectral density J c (ω) in Eq. ( 6) as λ 2 = α 1 α 2 γ and ω 2 0 = γ α 2 .Taking once again the large friction limit γ 1 leads to the overdamped spectrum J c (ω) 76 , where we pick α 1 = γ h and α 2 = Λ h , with appropriate dimensions.In Fig. 3(a) we plot both the resulting spectral density (solid) along with the Ohmic limiting case of γ → ∞ (dashed).As it can be seen, for our choice of parameters, the corresponding wire-bath coupling ends up being O(γ h ) 1 which would justify the weak-coupling approximation and the use of a per-turbative master equation.
The next step towards a GKLS equation is to certify the validity of the Markov approximation: we must ensure that the decay of the bath correlation functions computed in Eq. ( 7) is sufficiently fast when compared to the dynamics of the wire.In Fig. 3(b) we plot the integrated correlation t 0 ds B B B c (s) B B B c (0) , whose saturation time (γ h τ c ∼ 5 × 10 −4 ) is just below the relevant time scale for the dissipative evolution of the wire (γ h τ w ∼ 10 −3 ) [compare with Fig. 4(a)].We thus confidently say that the Markov approximation holds.For the parameters chosen, the secular approximation is also not a problem (cf.caption of Fig. 4).Namely, the normal-mode frequencies of the wire are Ω 1 = 0.34 and Ω 2 = 0.97, while the dissipation rates are both O(γ h ), which is perturbative.
thus take the time evolution of the two-node wire according to the master equation ( 14), as valid approximation to the exact dissipative dynamics, and a good benchmark for the RC mapping.Just like we did in Sec.III A, we also apply a GKLS master equation to the resulting three-node augmented system; again in spite of the fact that it is totally unjustified (the residual dissipation is γ = 10 3 ).As pointed out in Sec.II A 2, initially, we assume no correlations between the reaction coordinate, the wire, and either of the two baths, and initialise the RC in a thermal state at the original temperature of the cold bath.
Our results are plotted in Fig. 4. As we can see, the RC-GKLS mapping (open dots) accurately approximates the dynamics of the covariances of the wire (solid black line), and it does so during the entire evolution.However, as expected from the results in Fig. 2(a), it fails to capture the covariances of the reaction coordinate itself.We show this in Fig. 4(d) by comparing the stationary value of X X X 2 RC as predicted by the master equation, with its exact asymptotic value.
Finally, we also take the opportunity here to illustrate the vital importance of the frequency shift δ δ δH H H w-c on the augmented system (cf.Sec.II A 1).Note that before the mapping we do not include any shifts in the Hamiltonian of the wire, since we are tackling the original problem via a master equation.However, for the mapping to be an identity, the frequency of the 'cold oscillator' must be shifted as in ω c → ω 2 c + λ 2 /ω 2 0 when applying the master equation to the augmented system.For our choice of parameters this means tuning it from 0.5 to 0.501, which might seem totally negligible.Indeed, the short-time dynamics [cf.Fig. 4(a)] and the stationary state [cf.Fig. 4(c)] remain virtually unaffected when the shift is not taken into account (dashed grey lines).At intermediate times, however, the effects of the shift become evident, as shown in Fig. 4(b)-neglecting it does cause the RC-GKLS mapping to break down.

IV. CONCLUSIONS
We have benchmarked the reaction-coordinate mapping in an exactly solvable linear model, consisting of a two-node chain of harmonic oscillators.These are individually coupled to two baths at different temperatures and thus, support a steady-state heat current.The mapping takes this setup into a h at different stages of the dynamics, according to a master equation applied directly on the two-node wire (solid black), on the three-node augmented system (open circles), and on an augmented system whose cold frequency has not been suitably shifted as per the RC mapping (grey dashed).(d) (solid black) exact steady-state value of X X X 2 RC superimposed to the asymptotic value of this covariance according to the master equation, acting on a shifted (open circles) and unshifted (dashed grey) augmented system.Here, ω h = 0.1, ω c = 0.5, k = 0.4, T h = 0.6, T c = 0.5, and the rest of parameters are the same as in Fig. 3. three-oscillator augmented system, which is also exactly solvable.The idea, however, is to tackle the augmented system via a weak-coupling Markovian master equation.What we found can be summarised as follows: • The reduction of the stationary state of the augmented system onto the degrees of freedom of the two-node wire-according to the master equation-resembles very closely the exact steady state.This can be so, even in regimes of parameters for which the approximations underpinning the master equation break down; specifically, the secular approximation and even the basic weak-coupling assumption.
• Even when the stationary state of the wire is captured faithfully by the master-equation approach, the joint state of all three nodes of the augmented system can differ very substantially from the exact solution of the augmented problem.This happens whenever the underlying approximations cease to be justified.
• More importantly, the non-equilibrium steady state of the wire may be accurately reproduced by the master equation acting on the augmented system and yet, the stationary heat currents obtained from it can be quantitatively and even qualitatively wrong.
• At least in the overdamped limit, the reactioncoordinate mapping succeeds in approximating the state of the wire not only asymptotically, but throughout the entire dissipative dynamics.
In addition, we discussed the subtleties surrounding the frequency renormalisation shifts appearing as a result of the system-environment(s) coupling, and illustrated the importance of using them consistently.We also presented in full detail a consistent Markovian master equation in GKLS form that generalises previous results 36 , and can be directly applied to an arbitrary network of of N harmonic oscillators locally connected to M heat baths at different temperatures.We explicitly provided the corresponding (Gaussian) nonequilibrium steady state, and the expression for the M stationary heat currents flowing across the network.
Our results have two important consequences when dealing with virtually intractable problems involving nano-and micro-scale systems in non-Markovian baths, such as biological environments.On the one hand, they raise hopes of relying on the combination of 'reaction-coordinate mapping' and 'weak-coupling master equations' beyond the strict range of applicability of the latter.Although the mapping had been successfully applied to open-systems strongly coupled to highly structured environments 19,[23][24][25][26][27][28]30,[32][33][34]77 , our findings suggest that non-Markovian noise featuring broader power spectra may also be modelled in the exact same manner. On the othe hand, however, weak-coupling master equations should not be trusted beyond their strict range of applicability when calculating boundary heat currents-even if these appear to be thermodynamically consistent, they may be serious overestimations.It is pertinent to keep this it in mind when using the reaction-coordinate mapping to extend quantum thermodynamics into the strong coupling regime, an interesting line which currently attracts an increasing attention 20,[24][25][26]32,33,77 .Put simply, being able to replicate accurately the exact numerical propagation of an open system with the reaction-coordinate technique does not guarantee that the boundary heat (or particle) currents calculated from the corresponding master equation are equally accurate.This is our main message.
It may be possible to improve on the boundary currents by taking the secular approximation back and working with the full Redfield equation 31,32 .It would thus be interesting to generalise Eqs. ( 16)-( 19) and ( 21) to allow for non-secular contributions, and benchmarking those instead.After all, as already mentioned the reaction-coordinate mapping is often combined with Redfield rather than GKLS quantum master equations 19,23,24,31,32 .It is important to bear in mind, however, that Redfield equations may not only violate complete positivity, but even positivity alone 55,56 , which seriously compromises the consistency of any quantum-thermodynamic variables derived from it.This generalisation lies, however, beyond the scope of this paper, and will be tackled elsewhere.

2 FIG. 2 .
FIG. 2. (a)(solid) Uhlmann fidelity between C(2) ex (∞) and the relevant two-node reduction of C(3) me (∞).This is achieved simply by eliminating rows and columns related to the reaction-coordinate variables X X X RC and P P P RC from the 6×6 matrix C(3) me (∞).The abscissa corresponds to the log of the friction coefficient γ of the underdamped spectral density in Eq. (6), at the interface between the wire and the cold bath (normalised by the dissipation strength into the hot bath γ h ).Only in the shaded-grey area the fidelity drops below 95%.(dashed) Fidelity between C(3) ex (∞) and C(3) me (∞) as a function of the normalised friction γ/γ h .(b) Steady-state heat currents coming from the hot (red) and cold (blue) baths into the wire as a function of γ/γ h for the same parameters as in panel (a).The solid lines correspond to the exact calculation from Eq. (28a), while the dashed ones result from applying the GKLS master equation to the augmented system, as per Eq.(21).The inset is a zoom into the low friction limit and the shaded-grey area is the same as in (a).In both panels, the wire is characterised by ω h = 1, ω c = 3, and k = 0.8; the overdamped Ohmic spectral density at the hot interface, by γ h = 10 −3 and Λ h = 10 3 ; the underdamped spectral density at the cold interface, by λ = 0.9 and ω 0 = 4; and the baths, by temperatures T h = 3.3 and T c = 1.2.