Pushing the limits of the reaction-coordinate mapping

TL;DR

This paper investigates the reaction-coordinate mapping in quantum thermodynamics, demonstrating that a weak-coupling approximation can accurately reproduce steady states even under strong residual dissipation, but may overestimate heat currents.

Contribution

It reveals that the reaction-coordinate mapping can be effective in the strong-coupling regime, challenging previous assumptions about its limitations, and highlights potential inaccuracies in heat current predictions.

Findings

Weak-coupling treatment reproduces steady states under strong residual dissipation

Master equation captures entire dynamics in overdamped regime

Overestimates heat current despite accurate steady state

Abstract

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.

Figures

Equations

$${\mathbf{H}}_w={\sum }_{\alpha \in \{h,c\}}\left(\frac{1}{2}{\omega }_{\alpha }^2{\mathbf{X}}_{\alpha }^2+\frac{{\mathbf{P}}_{\alpha }^2}{2}\right)+\frac{k}{2}({\mathbf{X}}_h-{\mathbf{X}}_c{)}^2.$$

$${\mathbf{H}}_{\text{diss},\alpha }={\mathbf{X}}_{\alpha }\otimes {\mathbf{B}}_{\alpha }:={\mathbf{X}}_{\alpha }{\sum }_{\mu }{g}_{\mu }^{(\alpha )}{\mathbf{x}}_{\mu }^{(\alpha )}\alpha \in \{h,c\},$$

$$J_{\alpha }(\omega ):=\pi {\sum }_{\mu }\frac{g_{\mu }^{(\alpha )2}}{2{\omega }_{\mu }}\delta (\omega -{\omega }_{\mu })\alpha \in \{h,c\}.$$

$$\mathbf{\delta }{\mathbf{H}}_{w\text{–}\alpha }={\mathbf{X}}_{\alpha }^2{\sum }_{\mu }\frac{g_{\mu }^{(\alpha )2}}{2{\omega }_{\mu }^2}={\mathbf{X}}_{\alpha }^2{\int }_0^{\infty }\frac{\text{d}\omega }{\pi }\frac{J_{\alpha }(\omega )}{\omega }:=\frac{{\delta }_{\alpha }}{2}{\mathbf{X}}_{\alpha }^2,$$

$$\mathbf{H}={\mathbf{H}}_{T_h}+{\mathbf{H}}_{\text{diss},h}+\mathbf{\delta }{\mathbf{H}}_{w\text{–}h}+{\mathbf{H}}_w+\mathbf{\delta }{\mathbf{H}}_{w\text{–}c}+{\mathbf{H}}_{\text{diss},c}+{\mathbf{H}}_{T_c}.$$

$$J_c(\omega )=\frac{\gamma {\lambda }^2\omega }{{\gamma }^2{\omega }^2+({\omega }^2-{\omega }_0^2{)}^2},$$

$$⟨{\mathbf{B}}_{\alpha }(t){\mathbf{B}}_{\alpha }(0)⟩={\int }_0^{\infty }\frac{\text{d}\omega }{\pi }{J}_{\alpha }(\omega )\left(\coth \frac{\omega }{2T_{\alpha }}\cos \omega t-i\sin \omega t\right).$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\mathbf{\varrho }}_w(t)=-i{\text{tr}}_{\bar{w}}[\tilde{\mathbf{H}},\tilde{\mathbf{\rho }}]$$

$$\begin{gathered} \tilde{\mathbf{H}}:={\mathbf{H}}_{T_h}+{\mathbf{H}}_{\text{diss},h}+\mathbf{\delta }{\mathbf{H}}_{w\text{–}h}+{\mathbf{H}}_w+\mathbf{\delta }{\mathbf{H}}_{w\text{–}c} \\ -\lambda {\mathbf{X}}_c{\mathbf{X}}_{\text{RC}}+\frac{1}{2}\left({\omega }_0^2{\mathbf{X}}_{\text{RC}}^2+{\mathbf{P}}_{\text{RC}}^2\right)+\mathbf{\delta }{\mathbf{H}}_{\text{RC–res}} \\ +{\mathbf{X}}_{\text{RC}}{\sum }_{\mu }{\tilde{g}}_{\nu }^{(c)}{\tilde{\mathbf{x}}}_{\nu }^{(c)}+{\sum }_{\nu }{\omega }_{\nu }{\tilde{\mathbf{a}}}_{\nu }^{(c)†}{\tilde{\mathbf{a}}}_{\nu }^{(c)}, \end{gathered}$$

$$\frac{\text{d}}{\text{d}t}{\mathbf{\varrho }}_w(t)=-i{\text{tr}}_{\bar{w}}[\mathbf{H},\mathbf{\rho }],$$

$${\mathbf{\varrho }}_{T_h}\otimes {\mathbf{\varrho }}_w(0)\otimes {\mathbf{\varrho }}_{T_c}={\mathbf{\varrho }}_{T_h}\otimes {\mathbf{\varrho }}_w(0)\otimes {\tilde{\mathbf{\varrho }}}_{\text{RC + res}}(0).$$

$$\boldsymbol{H}_{N}=\frac{1}{2}\,\big{(}\vec{\boldsymbol{X}}^{\mathsf{T}}\mathsf{V}\,\vec{\boldsymbol{X}}+\vec{\boldsymbol{P}}^{\mathsf{T}}\vec{\boldsymbol{P}}\big{)},$$

$${\mathbf{L}}_i({\Omega }_j)=\frac{{\mathsf{P}}_{ij}}{\sqrt{2{\Omega }_j}}{\mathbf{b}}_j,{\mathbf{L}}_i(-{\Omega }_j):={\mathbf{L}}_i({\Omega }_j{)}^{†},$$

$$\frac{\text{d}\boldsymbol{O}(t)}{\text{d}t}=i\left[\boldsymbol{H}_{N},\boldsymbol{O}(t)\right]+\sum\nolimits_{i=1}^{M}\left(\sum\nolimits_{j=1}^{N}\Gamma_{i}(\Omega_{j})\,\big{(}\boldsymbol{L}_{i}(-\Omega_{j})\,\boldsymbol{O}(t)\,\boldsymbol{L}_{i}(\Omega_{j})-\frac{1}{2}\big{\{}\boldsymbol{L}_{i}(-\Omega_{j})\,\boldsymbol{L}_{i}(\Omega_{j}),\boldsymbol{O}(t)\big{\}}_{+}\big{)}\right.\\ \left.+\sum\nolimits_{j=1}^{N}\Gamma_{i}(\Omega_{j})\,e^{-\Omega_{j}/T_{j}}\,\big{(}\boldsymbol{L}_{i}(\Omega_{j})\,\boldsymbol{O}(t)\,\boldsymbol{L}_{i}(-\Omega_{j})-\frac{1}{2}\big{\{}\boldsymbol{L}_{i}(\Omega_{j})\,\boldsymbol{L}_{i}(-\Omega_{j}),\boldsymbol{O}(t)\big{\}}_{+}\big{)}\right),$$

$${\mathbf{H}}_3={\mathbf{H}}_w+\mathbf{\delta }{\mathbf{H}}_{w\text{–}c}-\lambda {\mathbf{X}}_c{\mathbf{X}}_{\text{RC}}+\frac{1}{2}\left({\omega }_0^2{\mathbf{X}}_{\text{RC}}^2+{\mathbf{P}}_{\text{RC}}^2\right)$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\eta }}_j^2⟩$$

$$\frac{\text{d}}{\text{d}t}⟨\{{\mathbf{\eta }}_j,{\mathbf{\pi }}_j{\}}_{+}⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\pi }}_j^2⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\eta }}_j{\mathbf{\eta }}_k⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\eta }}_j{\mathbf{\pi }}_k⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\pi }}_j{\mathbf{\pi }}_k⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\eta }}_j⟩$$

$$\frac{\text{d}}{\text{d}t}⟨{\mathbf{\pi }}_j⟩$$

$$⟨{\mathbf{\eta }}_j^2(\infty )⟩=-\frac{\tilde{\Sigma }({\Omega }_j)}{2\tilde{\Delta }({\Omega }_j){\Omega }_j},$$

$$\displaystyle\big{\langle}\frac{1}{2}\{\boldsymbol{\eta}_{j}(\infty),\boldsymbol{\pi}_{j}(\infty)\}_{+}\big{\rangle}=0,$$

$$⟨{\mathbf{\pi }}_j^2(\infty )⟩=-\frac{\tilde{\Sigma }({\Omega }_j){\Omega }_j}{2\tilde{\Delta }({\Omega }_j)},$$

$$\mathsf{Q}=\left(\begin{array}[]{ccccc}\mathsf{P}_{11}&0&\mathsf{P}_{12}&0&\cdots\\ 0&\mathsf{P}_{11}&0&\mathsf{P}_{12}&\cdots\\ \mathsf{P}_{21}&0&\mathsf{P}_{22}&0&\cdots\\ 0&\mathsf{P}_{21}&0&\mathsf{P}_{22}&\cdots\\ \vdots&\vdots&\vdots&\vdots&\ddots\end{array}\right)\,.$$

$$\begin{gathered} {\dot{\mathcal{Q}}}_{i,\text{me}}^{(N)}={\sum }_{j=1}^N\frac{{\mathsf{P}}_{ij}^2}{2{\Omega }_j}{\Delta }_i({\Omega }_j)\left(\frac{1}{2}{\Omega }_j^2⟨{\mathbf{\eta }}_j^2(\infty )⟩+\frac{1}{2}⟨{\mathbf{\pi }}_j^2(\infty )⟩\right) \\ +\frac{{\mathsf{P}}_{ij}^2}{4}{\Sigma }_i({\Omega }_j). \end{gathered}$$

$$\begin{gathered} \frac{{\text{d}}^2}{\text{d}{t}^2}{\mathbf{X}}_{\alpha }+({\omega }_{\alpha }^2+{\delta }_{\alpha }){\mathbf{X}}_{\alpha }+k({\mathbf{X}}_{\alpha }-{\mathbf{X}}_{\bar{\alpha }}) \\ ={\mathbf{F}}_{\alpha }(t)+{\int }_{t_0}^t\text{d}s{\chi }_{\alpha }(t-s){\mathbf{X}}_{\alpha }(s), \end{gathered}$$

$$\left(\begin{array}[]{cc}-\omega^{2}+\omega_{h}^{2}+\delta_{h}+k-\hat{\chi}_{h}(\omega)&-k\\ -k&-\omega^{2}+\omega_{c}^{2}+\delta_{c}+k-\hat{\chi}_{c}(\omega)\end{array}\right)\left(\begin{array}[]{c}\hat{\boldsymbol{X}}_{h}\\ \hat{\boldsymbol{X}}_{c}\end{array}\right)\coloneqq\mathsf{A}(\omega)\left(\begin{array}[]{c}\hat{\boldsymbol{X}}_{h}\\ \hat{\boldsymbol{X}}_{c}\end{array}\right)=\left(\begin{array}[]{c}\hat{\boldsymbol{F}}_{h}\\ \hat{\boldsymbol{F}}_{c}\end{array}\right)\,.$$

$$\big{\langle}\frac{1}{2}\{\boldsymbol{X}_{\alpha}(t^{\prime}),\boldsymbol{X}_{\beta}(t^{\prime\prime})\}\big{\rangle}=\int_{-\infty}^{\infty}\frac{\text{d}\omega^{\prime}}{2\pi}\int_{-\infty}^{\infty}\frac{\text{d}\omega^{\prime\prime}}{2\pi}\,\sum\nolimits_{\delta\gamma}[\mathsf{A}^{-1}]_{\alpha\gamma}(\omega^{\prime})\,[\mathsf{A}^{-1}]_{\,\beta\delta}(\omega^{\prime\prime})\,\big{\langle}\frac{1}{2}\{\hat{\boldsymbol{F}}_{\gamma}(\omega^{\prime}),\hat{\boldsymbol{F}}_{\delta}(\omega^{\prime\prime})\}\big{\rangle}\,e^{-i\omega^{\prime}t^{\prime}}e^{-i\omega^{\prime\prime}t^{\prime\prime}}$$

$$\text{Im}{\hat{\chi }}_{\alpha }=J_{\alpha }(\omega )\Theta (\omega )-J_{\alpha }(-\omega )\Theta (-\omega )$$

$$\text{Re}{\hat{\chi }}_{\alpha }=\frac{1}{\pi }\text{P}{\int }_{-\infty }^{\infty }\text{d}{\omega }^{\prime }\frac{\text{Im}{\hat{\chi }}_{\alpha }({\omega }^{\prime })}{{\omega }^{\prime }-\omega }$$

$$\displaystyle\big{\langle}\frac{1}{2}\{\hat{\boldsymbol{F}}_{\alpha}(\omega^{\prime}),\hat{\boldsymbol{F}}_{\beta}(\omega^{\prime\prime})\}\big{\rangle}=\frac{1}{2\pi}\,\coth{\Big{(}\frac{\omega^{\prime}}{2T_{\alpha}}\Big{)}}\,\text{Im}\,\hat{\chi}_{\alpha}$$

$$\times {\delta }_{\alpha \beta }\delta ({\omega }^{\prime }+{\omega }^{\prime \prime }),$$

Full text

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 Garg, Onuchic, and Ambegaokar (1985); Hartmann, Goychuk, and Hänggi (2000); Roden et al. (2012), to understanding the thermodynamics of quantum systems Binder et al. (2018), or assisting in the design of nano-structures for quantum-technological applications Alicki et al. (2004); McCutcheon and Nazir (2010); Higgins, Lovett, and Gauger (2013). 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 de Vega and Alonso (2017). Various tools exist to deal with such problems, including exact path-integral methods Feynman and Vernon Jr (2000); Hu, Paz, and Zhang (1992); Tanimura (1990), stochastic Schrödinger equations Stockburger and Grabert (2002); Alonso and de Vega (2005), unitary transformations Wagner (1986); Würger (1998), or Markovian embeddings Garraway (1997); Martinazzo et al. (2011); Woods et al. (2014); Iles-Smith, Lambert, and Nazir (2014). Here, we shall focus on the latter; specifically, on the “reaction-coordinate mapping” Nazir and Schaller (2018).

In a seminal paper by Garg et al. Garg, Onuchic, and Ambegaokar (1985) 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” Hughes, Christ, and Burghardt (2009); Martinazzo et al. (2011); Iles-Smith, Lambert, and Nazir (2014) 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 ‘augmented 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 Martinazzo et al. (2011); Woods et al. (2014), albeit with a much larger system size. Whenever the residual friction (i.e., dissipation strength) is perturbatively small, the problem can be rigorously solved via standard weak-coupling Markovian master equations 111A Fokker–Plank equation may be derived in the opposite large-friction limit Garg, Onuchic, and Ambegaokar (1985); Iles-Smith, Lambert, and Nazir (2014).. 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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Strasberg et al. (2016); Restrepo et al. (2018); Wertnik et al. (2018); Maguire, Iles-Smith, and Nazir (2018); Puebla et al. (2019); Martensen and Schaller (2019); McConnell and Nazir (2019); Lambert et al. (2019) and fermionic Strasberg et al. (2018); Schaller et al. (2018); Restrepo et al. (2019) 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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016), make it particularly appealing as a general-purpose 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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Strasberg et al. (2018), which made us wonder where are its true limitations 222Recently, the RC mapping has been shown to break down at large friction in the zero-temperature limit Lambert et al. (2019). Our calculations here are, however, limited to finite temperatures..

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 González et al. (2017); 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, resulting 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 non-equilibrium 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 Asadian et al. (2013). 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 Caldeira and Leggett (1983); Ford, Lewis, and O’Connell (1988); Weiss (2008), as it is particularly important to use them consistently when performing the reaction-coordinate mapping.

This paper is structured as follows: In Sec. II.1, we introduce our simple model and discuss very briefly the reaction-coordinate mapping. In Sec. II.2, 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 steady-state heat currents. In Sec. II.3 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.1) and dynamical (cf. Sec. III.2) 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

II.1 A two-node non-equilibrium quantum wire

II.1.1 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

[TABLE]

Note that here and in what follows, we set all masses to one. We shall also take $\hbar=k_{B}=1$. The wire is kept out of equilibrium by two linear bosonic baths at temperatures $T_{\alpha}$. Throughout, $\alpha\in\{h,c\}$ stands for ‘hot’ or ‘cold’, i.e., $T_{h}>T_{c}$. Their Hamiltonians can thus be cast as $\boldsymbol{H}_{T_{\alpha}}=\sum\nolimits_{\mu}\omega_{\mu}\,\boldsymbol{a}_{\mu}^{(\alpha)\,\dagger}\boldsymbol{a}^{(\alpha)}_{\mu}$, where $\boldsymbol{a}^{(\alpha)\,\dagger}_{\mu}$ ($\boldsymbol{a}_{\mu}^{(\alpha)}$) is a creation (annihilation) operator of bath $\alpha$ in the collective bosonic environmental mode at frequency $\omega_{\mu}$. In turn, the dissipative interactions between the wire and the baths are

[TABLE]

where the quadratures $\sqrt{2\omega_{\mu}}\,\boldsymbol{x}^{(\alpha)}_{\mu}\coloneqq\boldsymbol{a}_{\mu}^{(\alpha)\,\dagger}+\boldsymbol{a}_{\mu}^{(\alpha)}$ and, as usual, the coupling constants $g_{\mu}^{(\alpha)}$ make up the spectral densities

[TABLE]

Importantly, each system–bath coupling $\boldsymbol{H}_{\text{diss},\alpha}$ requires us to introduce a renormalisation term in the bare Hamiltonian of the wire $\boldsymbol{H}_{w}\mapsto\boldsymbol{H}_{w}+\boldsymbol{\delta}\boldsymbol{H}_{w\text{--}\alpha}$, which compensates for the environmental distortion on the system’s potential Weiss (2008). If we were not to include such terms and let $T_{h}=T_{c}=T$ be arbitrarily large, the exact stationary state would approach $\boldsymbol{\varrho}_{w}(\infty)\sim\exp{[-(\boldsymbol{H}_{w}-\boldsymbol{\delta}\boldsymbol{H}_{w\text{--}h}-\boldsymbol{\delta}\boldsymbol{H}_{w\text{--}c})/T]}$ instead of the classical limit $\boldsymbol{\varrho}_{w}(\infty)\sim\exp{(-\boldsymbol{H}_{w}/T)}$; this should be seen as an important deficiency of the model Caldeira and Leggett (1983). Specifically, these extra terms are

[TABLE]

and the full Hamiltonian of our system is, therefore,

[TABLE]

We take an Ohmic spectrum for the coupling to the ‘hot bath’, i.e., $J_{h}(\omega)=\gamma_{h}\,\omega\,\theta(\omega/\Lambda_{h})$, where $\theta(x)$ is some rapidly decaying function for arguments $x>1$, which places an upper bound on the excitation energies. For practical reasons we choose the algebraic cutoff $\theta(x)=(1+x^{2})^{-1}$, although other choices would not alter our results as long as $\Lambda_{h}$ is large. Such $J_{h}(\omega)$ is referred-to as ‘overdamped’ in the context of energy transfer in molecular systems Iles-Smith et al. (2016). For the coupling of the wire to the cold bath, we take instead the ‘underdamped’ spectrum

[TABLE]

which displays a peak around $\omega_{0}$, whose height and width are essentially controlled by $\lambda$ and $\gamma$, respectively 333Note that in the limit of very large $\gamma$ this becomes $J_{c}(\omega)\sim\frac{\lambda^{2}}{\gamma\,\omega}$.. This is precisely the effective spectral density resulting from the aforementioned ansantz by Garg et al. Garg, Onuchic, and Ambegaokar (1985). The frequency-renormalisation shifts $\delta_{\alpha}$ for these expectral densities are explicitly given by $\delta_{h}=\gamma_{h}\Lambda_{h}$ and $\delta_{c}=\lambda^{2}/\omega_{0}^{2}$.

The decay of the environmental correlation functions $\left\langle\boldsymbol{B}_{\alpha}(t)\,\boldsymbol{B}_{\alpha}(0)\right\rangle$ 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 Breuer and Petruccione (2002)

[TABLE]

While, at finite temperatures, a spectral density like our $J_{h}(\omega)$ 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. However, at sufficiently low temperatures—a regime which we shall not explore here—the bath correlation times can become comparable to the typical system dynamics even for an Ohmic spectral density.

II.1.2 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 $\omega_{0}$—to a residual reservoir with a purely Ohmic spectrum Garg, Onuchic, and Ambegaokar (1985); the coupling between the auxiliary mode and the system being of strength $\lambda$ (see Fig. 1). Put in other words, the dynamics

[TABLE]

generated by

[TABLE]

exactly coincides with that of

[TABLE]

when the coefficients $\{g^{(c)}_{\mu}\}$ in $\boldsymbol{H}_{\text{diss},\,c}$ correspond to Eq. (6) [by virtue of (3)] and the $\{\tilde{g}_{\nu}^{(c)}\}$ in the sixth term on the right-hand side of Eq. (9), to $\tilde{J}_{c}(\omega)=\gamma\,\omega$; technically, some suitable cutoff function $\theta{(\omega/\tilde{\Lambda}_{c})}$ with would be required, the mapping being exact only in the limit $\tilde{\Lambda}_{c}\rightarrow\infty$. Here, $\text{tr}_{\bar{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 $\{\tilde{\boldsymbol{x}}_{\nu}^{(c)}\}$, the creation and annihilation $\{\tilde{\boldsymbol{a}}_{\nu}^{\dagger\,(c)},\tilde{\boldsymbol{a}}_{\nu}^{(c)}\}$ operators in modes at frequency $\omega_{\nu}$; and the joint state of the hot bath, the wire, the reaction coordinate, and the residual reservoir $\tilde{\boldsymbol{\rho}}(t)$. Finally, the newly introduced operators $\boldsymbol{X}_{\text{RC}}$ and $\boldsymbol{P}_{\text{RC}}$ stand for the canonical degrees of freedom of the RC. Note that we have included as well the renormalisation term $\boldsymbol{\delta}\boldsymbol{H}_{\text{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 Garg, Onuchic, and Ambegaokar (1985); Iles-Smith, Lambert, and Nazir (2014); Martinazzo et al. (2011); Strasberg et al. (2016).

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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016). Note that the dynamics generated by Eqs. (8) and (10) only agree if $\boldsymbol{\rho}(0)=\tilde{\boldsymbol{\rho}}(0)$, i.e.,

[TABLE]

In particular, this means that the composite ‘RC + residual reservoir’ should start instead in a joint thermal state at temperature $T_{c}$; that is, $\tilde{\boldsymbol{\varrho}}_{\text{RC + res}}(0)=\boldsymbol{\varrho}_{T_{c}}$, which is not of the form $\boldsymbol{\varrho}_{\text{RC}}(0)\otimes\tilde{\boldsymbol{\varrho}}_{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 Breuer and Petruccione (2002). Luckily, in many cases of practical interest the residual interactions $\tilde{g}_{\mu}^{(c)}$ are sufficiently weak so that the dynamics is faithfully captured under this simple assumption. As we show in Sec. III.2 below, this is indeed the case when working in the overdamped limit. Furthermore, given its uniqueness Subaş ı et al. (2012), 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.2 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 applicability—namely, 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.3, since our $\boldsymbol{H}$ in Eq. (9) is fully linear.

II.2 Markovian master equation and its stationary solution

II.2.1 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 Breuer and Petruccione (2002) in the standard Gorini–Kossakowski–Lindblad–Sudarshan (GKLS) form Gorini, Kossakowski, and Sudarshan (1976); Lindblad (1976). In the present paper we shall only be interested in applying it to a simple 1D chain of three (and, in Sec. III.2, 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 Asadian et al. (2013). 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 Joshi et al. (2014); Levy and Kosloff (2014); Stockburger and Motz (2016); Kołodyński et al. (2018); Maguire, Iles-Smith, and Nazir (2018). 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 Wichterich et al. (2007); Trushechkin and Volovich (2016), or as a limiting case of a discrete collisional process Barra (2015); Barra and Lledó (2018); De Chiara et al. (2018). 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 “plug-and-play” stationary solution below [i.e., Eqs. (19) and (21)] only requires the diagonalisation of the corresponding $N\times N$ interaction matrix.

The Hamiltonian of a general linear network can be cast as

[TABLE]

assuming again that masses are $\mathsf{M}=\mathbbm{1}$. Here, $\vec{\boldsymbol{X}}$ and $\vec{\boldsymbol{P}}$ are $N$-dimensional vectors containing the position and momentum operators of each node, and $\mathsf{V}$ is real and symmetric. Let $\mathsf{P}$ be the orthogonal transformation that brings (12) into the diagonal form \boldsymbol{H}_{N}=\frac{1}{2}\,\big{(}\vec{\boldsymbol{\eta}}^{\mathsf{T}}\mathsf{\Omega}^{2}\,\vec{\boldsymbol{\eta}}+\vec{\boldsymbol{\pi}}^{\mathsf{T}}\vec{\boldsymbol{\pi}}\big{)}, where $\mathsf{\Omega}_{ij}=\Omega_{i}\delta_{ij}>0$ is a diagonal matrix formed of the normal mode frequencies corresponding to the conjugate variables \big{\{}\boldsymbol{\eta}_{i},\boldsymbol{\pi}_{i}\big{\}}_{i\in\{1,\cdots,N\}} (i.e., $\vec{\boldsymbol{\eta}}\coloneqq\mathsf{P}^{\mathsf{T}}\vec{\boldsymbol{X}}$).

The standard derivation of a Born–Markov secular master equation Breuer and Petruccione (2002); Levy and Kosloff (2014); González et al. (2017) now requires to decompose the ‘system–environment’ couplings [in our case, $\boldsymbol{X}_{i}$ for the $M$ nodes coupled to local baths, as per Eq. (2)] as eigen-operators of $\boldsymbol{H}_{N}$. That is $\boldsymbol{X}_{i}=\sum_{j}\boldsymbol{L}_{i}(\Omega_{j})+\boldsymbol{L}_{i}(\Omega_{j})^{\dagger}$ so that $[\boldsymbol{H}_{N},\boldsymbol{L}_{i}(\Omega_{j})]=-\Omega_{j}\boldsymbol{L}_{i}(\Omega_{j})$. These non-Hermitian operators, turn out to be simply

[TABLE]

where \boldsymbol{b}_{j}=\sqrt{\Omega_{j}/2}\,\big{(}\boldsymbol{\eta}_{j}+i\,\boldsymbol{\pi}_{j}/\Omega_{j}\big{)}. With these definitions, the equation of motion for an arbitrary Heisenberg-picture (Hermitian) operator $\boldsymbol{O}(t)$ under the Born–Markov and secular approximations reads Breuer and Petruccione (2002)

[TABLE]

with $\{\cdot,\cdot\}_{+}$ denoting anticommutator and decay rates \Gamma_{i}(\Omega_{j})\coloneqq 2\,J_{i}(\Omega_{j})\,\big{(}1-e^{-\Omega_{j}/T_{i}}\big{)}^{-1}, so that $\Gamma_{i}(-\Omega_{j})/\Gamma_{i}(\Omega_{j})=\exp{(-\Omega_{j}/T_{i})}$, thus reflecting local detailed balance.

The main appeal of Eq. (14) is that it is guaranteed to generate a completely positive and trace-preserving dynamics for the system Gorini, Kossakowski, and Sudarshan (1976); Lindblad (1976), unlike other frequently used weak-coupling master equations Suárez, Silbey, and Oppenheim (1992); Gaspard and Nagaoka (1999). Furthermore, under mild ergodicity assumptions, it admits a unique stationary solution Spohn (1977) which, in the case of a single environmental temperature $T$, is the thermal equilibrium state $\boldsymbol{\varrho}_{N}(t)\propto\exp{(-\boldsymbol{H}_{N}/T)}$ Spohn (1978). Importantly, this means that no renormalisation needs to be done on the Hamiltonian $\boldsymbol{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

[TABLE]

as the system Hamiltonian; i.e, we discard the renormalisation terms $\boldsymbol{\delta}\boldsymbol{H}_{h\text{--}w}$ and $\boldsymbol{\delta}\boldsymbol{H}_{\text{RC--res}}$ in Eq. (9), corresponding to the thermal contact with the hot and the residual environment, respectively.

However, the term $\boldsymbol{\delta}\boldsymbol{H}_{w\text{--}c}$ is—by construction—part of the augmented system after the reaction-coordinate mapping Iles-Smith, Lambert, and Nazir (2014); Strasberg et al. (2016). As we shall see in Sec. III.2 below, disregarding this latter term in the augmented-system Hamiltonian, e.g., on the basis of $\delta_{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.

II.2.2 Equations of motion for the covariances

Applying Eq. (14) to the symmetrised covariances \big{\langle}\frac{1}{2}\{\boldsymbol{r}_{j}(t),\boldsymbol{r}_{k}(t)\}_{+}\big{\rangle}\coloneqq[\mathsf{c}_{\text{me}}^{(N)}]_{jk}(t), where $\vec{\boldsymbol{r}}=(\boldsymbol{\eta}_{1},\boldsymbol{\pi}_{1},\cdots,\boldsymbol{\eta}_{N},\boldsymbol{\pi}_{N})^{\mathsf{T}}$, yields a closed algebra for the ‘covariance matrix’ of the network $\mathsf{c}_{\text{me}}^{(N)}(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.3 below. Specifically, we have

[TABLE]

together with the asymptotically vanishing covariances (for $j\neq k$)

[TABLE]

where $\Sigma_{i}(\Omega_{j})\coloneqq\Gamma_{i}(-\Omega_{j})+\Gamma_{i}(\Omega_{j})$ and $\Delta_{i}(\Omega_{j})\coloneqq\Gamma_{i}(-\Omega_{j})-\Gamma_{i}(\Omega_{j})$. For completeness, the equations of motion for the first-order moments $\langle\boldsymbol{\eta}_{j}\rangle$ and $\langle\boldsymbol{\pi}_{j}\rangle$ are given by

[TABLE]

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 Ferraro, Olivares, and Paris (2005) (that is, $\langle\boldsymbol{r}_{j}(t)\rangle$ and \big{\langle}\frac{1}{2}\{\boldsymbol{r}_{j}(t),\boldsymbol{r}_{k}(t)\}_{+}\big{\rangle}), Eqs. (16)–(18) thus provide a full dynamical description of the problem. Furthermore, since $\langle\boldsymbol{r}_{j}(\infty)\rangle=\langle\boldsymbol{r}_{j}(\infty)\,\boldsymbol{r}_{k}(\infty)\rangle=0$ for $j\neq k$, we can concentrate only in Eqs. (16) as far as the NESS is concerned. Explicitly, this is given by

[TABLE]

where $\tilde{\Sigma}(\Omega_{j})\coloneqq\sum\nolimits_{i=1}^{M}\mathsf{P}_{ij}^{2}/(2\Omega_{j})\,\Sigma_{i}(\Omega_{j})$ and $\tilde{\Delta}(\Omega_{j})\coloneqq\sum\nolimits_{i=1}^{M}\mathsf{P}_{ij}^{2}/(2\Omega_{j})\,\Delta_{i}(\Omega_{j})$. One can then transform $\mathsf{c}^{(N)}(t)$ into the covariance matrix $\mathsf{C}^{(N)}_{\text{me}}(t)$, defined in terms of the original variables $\vec{\boldsymbol{R}}=(\boldsymbol{X}_{1},\boldsymbol{P}_{1},\cdots,\boldsymbol{X}_{N},\boldsymbol{P}_{N})^{\mathsf{T}}$ by means of $\mathsf{C}_{\text{me}}^{(N)}(t)=\mathsf{Q}\,\mathsf{c}_{\text{me}}^{(N)}(t)\,\mathsf{Q}^{\mathsf{T}}$, where

[TABLE]

Importantly, under the secular approximation underpinning the GKLS equation, all the position–momentum covariances $\langle\boldsymbol{X}_{i}(\infty)\,\boldsymbol{P}_{j}(\infty)\rangle$ vanish in steady state. As a result, local current operators defined within the harmonic network would invariably average to zero Wichterich et al. (2007). We shall elaborate more on this in Sec. III.1.2 below. In order to compute the correct stationary heat currents, one can alternatively define the adjoint dissipation super-operators $\mathcal{L}_{i}^{\dagger}$ for each heat bath by rewriting Eq. (14) as $\text{d}\boldsymbol{O}(t)/\text{d}t\coloneqq i[\boldsymbol{H}_{N},\boldsymbol{O}]+\sum_{i=1}^{M}\mathcal{L}_{i}^{\dagger}\,\boldsymbol{O}$. That way, we can cast the steady-state heat current flowing from the $i$th bath into the network as $\dot{\mathcal{Q}}_{i,\text{me}}^{(N)}\coloneqq\langle\mathcal{L}_{i}^{\dagger}\,\boldsymbol{H}_{N}(\infty)\rangle$ Alicki (1979); Kosloff and Levy (2014). In our case, this evaluates to

[TABLE]

In Sec. III.1 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.3). In turn, in Sec. III.2, 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.

II.2.3 A note on the underlying approximations

To conclude this section, let us briefly comment on the approximations underlying the microscopic derivation of Eq. (14) Breuer and Petruccione (2002). First and foremost, it is a second-order perturbative expansion of the exact master equation in the system–environment(s) coupling Gaspard and Nagaoka (1999). 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 $\Omega_{j}$ be well separated as compared to the dissipation rates (i.e., $\min_{j\neq k}\{|\Omega_{j}-\Omega_{k}|,2\Omega_{j}\}\gg\max_{i\in\{1,\cdots,M\}}\,\{\gamma_{i}\}$). Once again, this approximation is incompatible with arbitrarily large dissipation rates $\gamma_{i}$, but may also be easily violated under weak dissipation Wichterich et al. (2007); González et al. (2017). For that reason, the full Redfield equation Suárez, Silbey, and Oppenheim (1992); Gaspard and Nagaoka (1999)—containing all non-secular terms—is often used instead when performing the RC mapping Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Strasberg et al. (2016). 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 Breuer and Petruccione (2002). 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 $\boldsymbol{H}_{N}$ Strasberg et al. (2016). 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 444Indeed, all the standard quantum-thermodynamic arguments based on the contractivity of the dissipative dynamics and uniqueness of the (thermal) fixed point of any of the dissipators of a GKLS equation Spohn (1978); Alicki (1979) hold regardless of whether or not the Lamb shift is included in $\mathcal{L}_{i}$. In particular, one always finds $\mathcal{L}_{i}\,e^{-\boldsymbol{H}_{N}/T_{i}}=0$, i.e., the fixed point of the dynamics is a thermal state with respect to the unshifted Hamiltonian.. 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 Strasberg et al. (2018). Note that this term is not related to the frequency renormalisation discussed in Sec. II.1.1 above.

II.3 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 $\boldsymbol{\rho}(0)=\boldsymbol{\varrho}_{T_{h}}\otimes\boldsymbol{\varrho}_{w}(0)\otimes\boldsymbol{\varrho}_{T_{c}}$, and no restrictions on $\boldsymbol{\varrho}_{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 Riseborough, Hanggi, and Weiss (1985); Ludwig, Hammerer, and Marquardt (2010); Fleming, Roura, and Hu (2011); Correa, Valido, and Alonso (2012); Martinez and Paz (2013); Valido, Ruiz, and Alonso (2015); Freitas and Paz (2014), 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. González et al. (2017), and here we limit ourselves to outline the key steps.

The exact dynamics of the wire obeys the following quantum Langevin equations Ford, Kac, and Mazur (1965); Ford, Lewis, and O’Connell (1988)

[TABLE]

where $\alpha\in\{h,c\}$ and $\bar{c}\coloneqq h$ and $\bar{h}\coloneqq 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 $\chi_{\alpha}(t)$ with a Heavisde step function $\Theta(t)$ [i.e., $\chi_{\alpha}(t)\mapsto\chi_{\alpha}(t)\,\Theta(t)$] Correa et al. (2017). Since we are interested in the steady state of the wire, our aim will be to compute the covariance matrix $\mathsf{C}^{(2)}_{\text{ex}}$ at any finite time $t$ while setting $t_{0}\rightarrow-\infty$.

With this in mind, we can now Fourier-transform Eqs. (22), which yields

[TABLE]

Here, the “hatted” symbols are in the frequency domain, i.e., $\hat{f}(\omega)\coloneqq\int_{-\infty}^{\infty}\text{d}t\,e^{i\omega t}\,f(t)$. Therefore, $(\hat{\boldsymbol{X}}_{h},\hat{\boldsymbol{X}}_{c})^{\mathsf{T}}=\mathsf{A}^{-1}(\omega)\,(\hat{\boldsymbol{F}}_{h},\hat{\boldsymbol{F}}_{c})^{\mathsf{T}}$, so that the objects we wish to compute are

[TABLE]

for $t^{\prime}=t^{\prime\prime}=t$. The position–momentum and momentum–momentum covariances can be obtained by differentiating Eq. (24), which is equivalent to multiplying the integrand by $(-i\omega^{\prime})$ and $(-\omega^{\prime}\omega^{\prime\prime})$, respectively. To carry out the integration in (24) explicitly, we only need the Fourier transform of the dissipation kernels $\hat{\chi}_{\alpha}(\omega)$ and the power spectrum of the environmental forces $\langle\frac{1}{2}\{\hat{\boldsymbol{F}}_{\alpha}(\omega^{\prime}),\hat{\boldsymbol{F}_{\beta}}(\omega^{\prime\prime})\}\rangle$. These are given by Correa et al. (2017)

[TABLE]

where P denotes ‘principal value’, $\delta_{\alpha\beta}$ is a Kronecker delta and $\delta(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}(\omega)=\gamma_{h}\Lambda_{h}^{2}\omega/(\omega^{2}+\Lambda_{h}^{2})$ and $J_{c}(\omega)=\gamma\lambda^{2}\omega/[\gamma^{2}\omega^{2}+(\omega^{2}-\omega_{0}^{2})^{2}]$, which yields

[TABLE]

Note that Eq. (26a) may also be used for the dissipation kernel of the residual bath acting on the augmented system, by merely replacing $\gamma_{h}$ with $\gamma$ and taking a large cutoff.

Summing up, Eqs. (24)–(26) are all we need to fill in the full stationary $4\times 4$ covariance matrix $\mathsf{C}^{(2)}_{\text{ex}}(\infty)$. Note that it is indeed possible to solve the problem not only exactly, but also analytically Correa, Valido, and Alonso (2012); Valido, Correa, and Alonso (2013). In turn, the $6\times 6$ NESS $\mathsf{C}^{(3)}_{\text{ex}}(\infty)$ of the augmented system can be found in a completely analogous way Valido, Correa, and Alonso (2013), by just replacing the ‘vector of forces’ $(\hat{\boldsymbol{F}}_{h},\hat{\boldsymbol{F}}_{c})^{\mathsf{T}}$ by $(\hat{\boldsymbol{F}}_{h},0,\hat{\boldsymbol{F}}_{\text{res}})^{\mathsf{T}}$ and $\mathsf{A}(\omega)$, with

[TABLE]

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 $\alpha$—i.e, $\dot{\mathcal{Q}}^{(N)}_{\alpha,\text{ex}}=i\langle[\boldsymbol{H}_{w\text{--}\alpha},\boldsymbol{H}_{w}]\rangle$—can be cast as Martinez and Paz (2013); Freitas and Paz (2017)

[TABLE]

III Discussion

III.1 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 $\mathsf{C}^{(2)}_{\text{ex}}(\infty)$, as well as that of the augmented (three-node) system, $\mathsf{C}^{(3)}_{\text{ex}}(\infty)$. Alternatively, we can look for the steady state of the augmented system according to the GKLS master equation, i.e., $\mathsf{C}_{\text{me}}^{(3)}$. Benchmarking the RC mapping thus amounts to assessing how “close” is the relevant $4\times 4$ submatrix of $\mathsf{C}_{\text{me}}^{(3)}(\infty)$ to the exact stationary state $\mathsf{C}_{\text{ex}}^{(2)}(\infty)$. We conclude by noting that the covariance dynamics can also be obtained non-perturbatively in the system–baths couplings by means of stochastic propagation and averaging, in linear and weakly non-linear continuous-variable systems Motz, Ankerhold, and Stockburger (2017); Motz et al. (2018).

We thus need to be able to quantify the distance between two covariance matrices $\mathsf{C}_{1}$ and $\mathsf{C}_{2}$. To that end, we resort to the Uhlmann fidelity Uhlmann (1976); Banchi, Braunstein, and Pirandola (2015) $\mathcal{F}(\mathsf{C}_{1},\mathsf{C}_{2})$ which, for arbitrary $N$-mode Gaussian states with vanishing first-order moments, is

[TABLE]

This is a meaningful distance measure, since $\mathcal{F}(\mathsf{C}_{1},\mathsf{C}_{2})=1$ only holds if the states are identical, and $0<\mathcal{F}(\mathsf{C}_{1},\mathsf{C}_{2})\leq 1$.

III.1.1 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 $\mathsf{C}_{\text{me}}^{(3)}(\infty)$ onto the wire degrees of freedom remains nearly identical to the exact stationary state $\mathsf{C}_{\text{ex}}^{(2)}(\infty)$, even at extremely large residual dissipation strengths $\gamma$. In the figure, for instance, the fidelity between the two states falls below $95\%$ only at $\gamma\gtrsim 60$. When it comes to the approximations that justify the GKLS master equation (14) (cf. Sec. II.2.3), this is completely off-limits. Indeed, note that the normal-mode frequencies of the augmented system are, in this example, $\Omega_{1}\simeq 1.31$, $\Omega_{2}\simeq 3.13$, and $\Omega_{3}\simeq 4.02$, which renders the secular approximation problematic already at residual dissipations as small as $\gamma\sim 0.1$. More importantly, $\gamma\simeq 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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Lambert et al. (2019), 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 $\gamma\gtrsim 0.1$, which corresponds to $\log{(\gamma/\gamma_{h})}\gtrsim 5$ [area to the right of the dotted line in Fig. 2(a)]. We can see this in Fig. 2(a), when instead of looking at the reduction of $\mathsf{C}_{\text{me}}^{(3)}(\infty)$ onto the wire, we consider the full augmented system and compare it with the exact three-node solution $\mathsf{C}_{\text{ex}}^{(3)}(\infty)$ (dashed line). Specifically, $\mathcal{F}(\mathsf{C}_{\text{ex}}^{(3)}(\infty),\mathsf{C}_{\text{me}}^{(3)}(\infty))<0.95$ for $\gamma>0.15$, as expected. We are thus not claiming that Markovian master equations in Lindblad 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).

We still need, however, to provide some physical intuition backing this observation. To that end, let us take a detour to comment on recent literature on locality of temperature in quantum many-body lattice systems Ferraro, García-Saez, and Acín (2012); García-Saez, Ferraro, and Acín (2009); Kliesch et al. (2014); Hernández-Santana et al. (2015). It is clear that the reduction of the global thermal state of a large lattice onto a small local subspace can deviate substantially from a local thermal state—this is due to the non-vanishing interactions between the subsystem in question and the rest of the lattice. However, the (non-thermal) state of such sub-lattice may be approximated arbitrarily well as follows: One first envelopes it with a ‘boundary’ or ‘buffer’ region, taken from the surrounding lattice; such augmented system is then set to a thermal state at the global temperature of the full system and then, the auxiliary buffer is traced out Ferraro, García-Saez, and Acín (2012). The result is in good agreement with the local state of interest so long as the boundary is thick enough, relative to some relevant correlation length scale Kliesch et al. (2014); Hernández-Santana et al. (2015). Something similar happens in our example: imposing incorrect (thermal) boundary conditions on an augmented system, we can reproduce the state of the wire faithfully; the techniques only breaks down when the boundary–environment interactions become sufficiently large, so that correlations start to appear between the wire and the residual environment. Making this intuition more precise by studying the correlation sharing structure between wire, RC, and residual bath, goes, however, beyond the scope of the present paper.

III.1.2 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 Kosloff and Levy (2014). 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 $\mathcal{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 $\gamma\sim 60$—where $\mathsf{C}_{\text{ex}}^{(2)}(\infty)$ and the reduction of $\mathsf{C}_{\text{me}}^{(3)}(\infty)$ 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 $\gamma$.

Note that, for us, resorting to the dissipators is indeed the only feasible way to estimate heat currents; $\mathsf{C}_{\text{me}}^{(3)}(\infty)$ is lacking the key covariances $\langle\boldsymbol{X}_{c}\,\boldsymbol{P}_{h}\rangle$ and $\langle\boldsymbol{X}_{c}\,\boldsymbol{P}_{\text{RC}}\rangle$ 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 Wichterich et al. (2007). Alternatively, one could think of waiving the secular approximation to work instead with a Redfield master equation Suárez, Silbey, and Oppenheim (1992); Gaspard and Nagaoka (1999). Although the aforementioned covariances would then cease to be zero, the calculation would continue to yield quantitatively wrong results at very large $\gamma$—this time simply due to the breakdown of the basic weak-coupling assumption. Ultimately, however, the Redfield approach might improve the GKLS results under moderate residual dissipation Strasberg et al. (2016, 2018). Therefore, even in the light of the promising observation made in Sec. III.1.1 above, great care must still be taken when relying on the RC–GKLS mapping to discuss quantum thermodynamics under non-Markovian dissipation.

III.2 Dynamics

One can now ask whether the resilience of the RC–GKLS mapping to strong residual dissipation is exclusively a steady-state feature, or whether it holds throughout the entire dissipative evolution. Unfortunately, we do not have an exact dynamical benchmark—at most, we are able to solve here for the steady state of the exact Eq. (22). We, therefore, chose parameters so that the original two-node problem can be described via a GKLS quantum master equation. We recall, however, that this type of equation can in principle be solved non-perturbatively at finite times with stochastic propagation techniques Motz, Ankerhold, and Stockburger (2017); Motz et al. (2018).

In particular, we scale $\omega_{0}$ and $\lambda$ in the structured spectral density $J_{c}(\omega)$ in Eq. (6) as $\lambda^{2}=\alpha_{1}\,\alpha_{2}\,\gamma$ and $\omega_{0}^{2}=\gamma\,\alpha_{2}$. Taking once again the large friction limit $\gamma\gg 1$ leads to the overdamped spectrum $J_{c}(\omega)\sim\alpha_{1}\,\alpha_{2}\,\omega/(\omega^{2}+\alpha_{2}^{2})$ Strasberg and Esposito (2017).

For our calculations, we will take the numerical values $\alpha_{1}=\gamma_{h}$ and $\alpha_{2}=\Lambda_{h}$. Note that $J_{c}(\omega)\sim\alpha_{1}\,\alpha_{2}\,\omega/(\omega^{2}+\alpha_{2}^{2})$ looks like the Ohmic-algebraic $J_{h}(\omega)$ introduced above, except for a missing factor $\alpha_{2}$ in the numerator. Hence, while $\alpha_{1}$ takes the numerical value of $\gamma_{h}$, it must have units of frequency squared instead of frequency. It is $\alpha_{1}/\alpha_{2}\ll\gamma_{h}$ which plays the role of the dissipation strength in this case. In Fig. 3(a) we plot both the resulting spectral density (solid) along with the Ohmic limiting case of $\gamma\rightarrow\infty$ (dashed). As it can be seen, for our choice of parameters, the corresponding wire–bath coupling ends up being at most $\mathcal{O}(\gamma_{h})\ll 1$, which would justify the weak-coupling approximation and the use of a perturbative 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 $\int_{0}^{t}\text{d}s\,\langle\boldsymbol{B}_{c}(s)\,\boldsymbol{B}_{c}(0)\rangle$, whose saturation time ($\gamma_{h}\tau_{c}\sim 5\times 10^{-4}$) is just below the relevant time scale for the dissipative evolution of the wire ($\gamma_{h}\tau_{w}\sim 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 $\Omega_{1}=0.34$ and $\Omega_{2}=0.97$, while the dissipation rates are both $\mathcal{O}(\gamma_{h})$, which is perturbative.

We 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.1, 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 $\gamma=10^{3}$). As pointed out in Sec. II.1.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 $\langle\boldsymbol{X}_{\text{RC}}^{2}\rangle$ as predicted by the master equation, with its exact asymptotic value. It is remarkable, however, that the covariances for the wire are perfectly reproduced in spite of the extremely large friction $\gamma=10^{3}$. This contrasts with the degradation of fidelity illustrated in Fig. 2(a), and is entirely due to our choice of friction-dependent $\lambda$ and $\omega_{0}$.

Finally, we also take the opportunity here to illustrate the vital importance of the frequency shift $\boldsymbol{\delta}\boldsymbol{H}_{w\text{--}c}$ on the augmented system (cf. Sec. II.1.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 $\omega_{c}\mapsto\omega_{c}^{2}+\lambda^{2}/\omega_{0}^{2}$ 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 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 reaction-coordinate 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 González et al. (2017), 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) non-equilibrium steady state, and the expression for the $M$ stationary heat currents flowing across the network. Note that we have focused exclusively on continuous-variable systems in Gaussian states and hence, extending our conclusions to finite-dimensional or non-linear models would require further work.

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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Strasberg et al. (2016); Newman, Mintert, and Nazir (2017); Strasberg et al. (2018); Wertnik et al. (2018); Schaller et al. (2018); Maguire, Iles-Smith, and Nazir (2018); Restrepo et al. (2018, 2019); McConnell and Nazir (2019); Puebla et al. (2019), our findings suggest that non-Markovian noise featuring broader power spectra—which so far was though to be out of reach for the mapping—may also be modelled in the exact same manner. On the other hand, however, weak-coupling master equations should not be trusted beyond their 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 Strasberg et al. (2016); Newman, Mintert, and Nazir (2017); Restrepo et al. (2018); Wertnik et al. (2018); Strasberg et al. (2018); Schaller et al. (2018); Nazir and Schaller (2018). 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.

We also note that a closely-related systematic technique has been recently put forward to emulate dissipation into structured environments through GKLS-type master equations Tamascelli et al. (2018), which can be used to deal with the strong friction regime. When it comes to extensions of our analysis, it may be possible to improve on the boundary currents by taking the secular approximation back and working with the full Redfield equation Strasberg et al. (2018); Lambert et al. (2019). 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 Iles-Smith, Lambert, and Nazir (2014); Iles-Smith et al. (2016); Strasberg et al. (2016, 2018); Lambert et al. (2019). It is important to bear in mind, however, that Redfield equations may not only violate complete positivity, but even positivity alone Suárez, Silbey, and Oppenheim (1992); Gaspard and Nagaoka (1999), 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.

Acknowledgements

We gratefully acknowledge Janet Anders, James Cresser, Ronnie Kosloff, Neill Lambert, Ahsan Nazir, Philipp Strasberg and Tommaso Tufarelli for useful discussions. This project was funded by the European Research Council (StG GQCOP, Grant No. 637352), the London Mathematical Society (Scheme 3 Grant No. 31826), and the US National Science Foundation (Grant No. NSF PHY1748958). LAC and BM thank the Kavli Institute for Theoretical Physics for their warm hospitality during the program “Thermodynamics of quantum systems: Measurement, engines, and control” and the associated conference.