A Reaction coordinate mapping

# Environmental dynamics, correlations, and the emergence of noncanonical equilibrium states in open quantum systems

## Abstract

Quantum systems are invariably open, evolving under surrounding influences rather than in isolation. Standard open quantum system methods eliminate all information on the environmental state to yield a tractable description of the system dynamics. By incorporating a collective coordinate of the environment into the system Hamiltonian, we circumvent this limitation. Our theory provides straightforward access to important environmental properties that would otherwise be obscured, allowing us to quantify the evolving system-environment correlations. As a direct result, we show that the generation of robust system-environment correlations that persist into equilibrium (heralded also by the emergence of non-Gaussian environmental states) renders the canonical system steady-state almost always incorrect. The resulting equilibrium states deviate markedly from those predicted by standard perturbative techniques and are instead fully characterised by thermal states of the mapped system-collective coordinate Hamiltonian. We outline how noncanonical system states could be investigated experimentally to study deviations from canonical thermodynamics, with direct relevance to molecular and solid-state nanosystems.

## I Introduction

A quantum system coupled to its macroscopic environment constitutes a challenging theoretical problem in which the large number of environmental degrees of freedom can lead to both conceptual and practical difficulties Leggett et al. (1987); Weiss (2008); Nitzan (2006); May and Kühn (2011); Breuer and Petruccione (2002). The master equation formalism has been developed to offer a simple and intuitive approach for describing such systems; the complex dynamical evolution of the many-body environment is not tracked explicitly, but instead only its effect on the reduced state of the system of interest is considered, eliminating all information on the environmental state. In most cases one must also rely on a series of assumptions in deriving a tractable master equation Note1 (). These customarily neglect the formation of system-bath correlations and lead to the eventual thermalisation of the system with respect solely to its internal Hamiltonian, resulting in canonical equilibrium states.

Here, by incorporating a collective coordinate of the environment into an effective system Hamiltonian, we develop a master equation formalism that can overcome such restrictions. This enables us to straightforwardly determine key environmental properties as well as track the dynamic generation of correlations between the system and bath. Specifically, we characterise the departure of the environment from its initial Gaussian thermal state due to interactions with the quantum system, and show that the resulting correlations—measured in terms of the mutual information—can have a profound effect on the system dynamics too, even persisting into the steady-state. We demonstrate that system-bath correlations are in fact generated on two distinct timescales, with long-lived correlations leading to a departure of the system steady-state from canonical equilibrium, as would otherwise be expected from a perturbative (Born-Markov) treatment, thus heralding the failure of the accepted statistical mechanics view of thermalisation Breuer and Petruccione (2002). Correctly capturing system-environment correlations is hence shown to be crucial in order to properly describe both the system transient and equilibrium behaviour. As a further, and unique, aspect of our approach we illustrate how noncanonical equilibrium states can still be characterised in terms of thermal states, but now with respect to the effective system-collective coordinate Hamiltonian. This also reveals simple experimental signatures by which deviations from canonical thermodynamics can be observed in real physical systems, for example through measurements of system populations.

## Ii Reaction coordinate mapping and master equation

Our method relies on keeping track of a collective environmental degree of freedom, and to do so we make use of the reaction coordinate mapping Garg et al. (1985); Thoss et al. (2001); Cao and Voth (1997); Hughes et al. (2009); Hartmann et al. (2000); Roden et al. (2012); Martinazzo et al. (2011); Woods et al. (2014); Imamoğlu (1994); Garraway (1997); we take a quantum system coupled to a bosonic environment and map to a model in which a collective mode of the environment, known as the reaction coordinate (RC), is incorporated within an effective system Hamiltonian. We then treat the residual environment within a full second-order Born-Markov master equation formalism. By comparing the RC master equation to the numerically exact hierarchical equations of motion (HEOM) (Tanimura, 1990; Tanimura and Kubo, 1989) we demonstrate essentially perfect agreement in the dynamics across all timescales (see below). Thus, all important system-bath, and indeed intra-bath Bera et al. (2014a, b), correlations are incorporated into the system-RC Hamiltonian in the regimes we study.

Though our approach may be applied quite generally, we shall focus in this work on a two-level system (TLS) described by the spin-boson Hamiltonian Leggett et al. (1987); Weiss (2008); Nitzan (2006); Breuer and Petruccione (2002); May and Kühn (2011); Makri and Makarov (1995); Prior et al. (2010); Alvermann and Fehske (2009); Wang and Thoss (2008); Anders et al. (2007); Winter et al. (2009); Silbey and Harris (1984); McCutcheon and Nazir (2011); Jang et al. (2008); Nazir (2009) (with ):

 H=ϵ2σz+Δ2σx+σz∑kfk(c†k+ck)+∑kνkc†kck. (1)

Here, () are creation (annihilation) operators for bosonic modes of frequency , which couple to the TLS with strength , and () are TLS Pauli operators defined such that . In the absence of the bath, the TLS splitting is determined by the bias and tunneling . The effect of the system-bath interaction can be completely characterised by introducing the spectral density Leggett et al. (1987), .

We now apply a normal mode transformation to Eq. (1) to incorporate the most important environmental degrees of freedom into a new effective system Hamiltonian. We carry out this procedure by first defining a collective coordinate of the environment Garg et al. (1985), the RC, which couples directly to the TLS, and is in turn coupled to a residual harmonic environment, as shown schematically in the upper panel of Fig 1. This leads to a mapped Hamiltonian of the form:

 HRC =ϵ2σz+Δ2σx+λσz(a†+a)+Ωa†a+∑kωkb†kbk +(a†+a)∑kgk(b†k+bk)+(a†+a)2∑kg2kωk, (2)

where the collective coordinate is defined such that:

 λ(^a†+^a)=∑kfk(c†k+ck), (3)

with coupling and frequency  Thoss et al. (2001); Cao and Voth (1997); Hughes et al. (2009). The residual bath, denoted by creation (annihilation) operators (), now couples only to the RC and is characterised by an effective spectral density, . To describe the action of the residual bath on the RC mode we need to relate this spectral density to the original spin-boson . To do so, we follow Refs. Garg et al. (1985); Hughes et al. (2009); Leggett (1984) and replace the TLS with a classical coordinate moving in a potential . As outlined in Appendix A, by considering the Fourier transformed equations of motion for both before and after the mapping, we may relate to the original .

To give a concrete example we now consider a Drude-Lorentz form commonly used to describe molecular systems Gilmore and McKenzie (2005): , with cut-off frequency and coupling strength . The mapping we exploit is not limited to this particular spectral density Woods et al. (2014), however it does allow us to benchmark our results against data attained from the numerically exact HEOM. By taking the RC spectral density to be , we find that the equivalent spin-boson form is given by in the limit . Hence, we recover the Drude-Lorentz spectral density by identifying , , and choosing such that  Thoss et al. (2001).

Diagonalising the first line of Eq. (II) we can account for the original spin-boson interaction term in Eq. (1) to all orders, which ensures that our formalism remains non-perturbative in the system-bath coupling strength . We thus derive a Born-Markov master equation for the reduced state of the composite TLS and RC, , which captures their exact internal dynamics, while treating the coupling to the residual bath to second order (see Appendix B). This results in the Schrödinger picture master equation Breuer and Petruccione (2002); Hausinger and Grifoni (2010):

 ∂ρ(t)∂t=−i[H0,ρ(t)]−[^A,[^χ,ρ(t)]]+[^A,{^Ξ,ρ(t)}], (4)

where we have assumed that only the residual environment remains in thermal equilibrium throughout the evolution. Here, ,

 ^χ= γ∞∫0dτ∞∫0dω ωcos(ωτ)coth(βω2)^A(−τ), (5) ^Ξ= γ∞∫0dτ∞∫0dωcos(ωτ)[H0,^A(−τ)], (6)

with and . To solve our master equation we truncate the RC space as necessary for convergence.

## Iii Benchmarking

In order to demonstrate the validity of the reaction coordinate master equation (RCME), we benchmark its predictions for the TLS population dynamics () against converged data generated using the numerically exact HEOM technique, see Appendix C and Refs. Tanimura, 1990; Tanimura and Kubo, 1989; Ishizaki and Fleming, 2009. To give an illustrative example, in Fig. 1 we have taken parameters that are typical for excitonic energy transfer in molecular systems Ishizaki and Fleming (2009); Prior et al. (2010); Thorwart et al. (2009); Pollock et al. (2013), with a representative value of  cm setting the other energy scales (i.e.  ps at ). The TLS is initialised in state , uncorrelated with both the RC and residual bath, which are taken to be in their respective thermal equilibrium states at a temperature ( K for  cm). We observe essentially perfect agreement between the RCME and the HEOM for a slow environment and for spin-boson coupling strengths [(a)  cm, (b)  cm] encompassing the transition from the coherent to the incoherent regime Note2 (). In contrast, a standard weak coupling approach, treating the interaction term in Eq. (1) to second order, fails even qualitatively to capture the correct system behaviour for any parameters shown. It is thus clear that environmental memory and the generation of correlations with the system—both of which are ignored in the weak coupling calculation—are crucial in order to capture the correct dynamical behaviour in this regime. Moreover, their combined impact on the TLS can be accurately described simply through the TLS-RC coupling. This is somewhat remarkable, given that any non-thermal state effects of the original bath have been reduced solely to the action of a single mode on the TLS.

## Iv Environmental dynamics and correlations

The most important aspect of our formalism is that the inclusion of environmental degrees of freedom into the system Hamiltonian allows us to gain additional insight into the dynamics of both the environmental state and system-environment correlations. We do this by calculating two complementary measures; the RC non-Gaussianity Genoni et al. (2008); Genoni and Paris (2010), which probes the environmental evolution:

 δG[ρRC(t)]=S(ϱ)−S(ρRC(t)), (7)

and the TLS-RC quantum mutual information (QMI) Groisman et al. (2005), characterising the correlations:

 I(ρS:ρRC)=S(ρs)+S(ρRC)−S(ρ). (8)

Here, () is the reduced state of the RC (TLS) and is the von-Neumann entropy. The non-Gaussianity determines the distance from to the nearest Gaussian reference state , and is defined such that iff is Gaussian. The QMI quantifies the total classical and quantum correlations shared between the TLS and RC  Groisman et al. (2005). As detailed in Appendix D, both measures act as rigorous lower bounds for the original spin-boson environment, which enables us to explore properties of the multi-mode bath and its correlations with the system simply through the single mode RC. Furthermore, in the limit that the Born approximation holds between the mapped system and residual environment, then the additive nature of the von-Neumann entropy implies that both these measures become exact for the original spin-boson environment.

Fig. 2 shows the dynamics of the non-Gaussianity and the QMI at both strong and weak system-environment coupling . One of the most striking features of these plots is the presence of two distinct timescales, which is most obvious for weak couplings, but is also present at stronger coupling strengths. At short times, the QMI is oscillatory, an indication of the memory effects implied by system-environment correlations, which also push the RC away from its initial Gaussian state. At longer times, we see that system-environment correlations—and consequently non-Gaussian environmental states—are also generated on a second timescale, and in fact persist into the steady state, with values of both the non-Gaussianity and QMI dependent on the coupling strength as demonstrated in Fig. 3. We therefore find a situation in which the Born-Markov approximation breaks down on all timescales. It describes neither the short time transient dynamics, due to the absence of bath memory effects in the Markov approximation, nor the long time behaviour, due to the generation of significant TLS-bath correlations. This leads to non-Gaussian (and hence non-thermal) environmental states upon tracing out the TLS, neglected in the Born approximation. Our line of enquiry also raises an intriguing question: if correlations can so dramatically affect the Gaussian nature of the environment, how do they impact upon the state of the TLS?

## V Noncanonical equilibrium states

A clue can already be taken from Fig. 3, where we find that the steady-state properties of both the RC and TLS-RC correlations may be described by a thermal state with respect to , i.e.

 ρth=e−βH0Z=1Ze−β(ϵ2σz+Δ2σx+λσz(a†+a)+Ωa†a), (9)

where . Likewise, in Fig. 4 we show that the equilibrium behaviour of the TLS (after tracing out the mode) clearly departs from the canonical statistics expected from a perturbative, weak coupling treatment of the environmental influence, which would be given instead by a thermal state with respect only to the TLS Hamiltonian, .

Fig. 4(a) explores the temperature dependence of the TLS equilibrium state in the energy eigenbasis. To emphasise the departure from canonical statistics we have plotted the population ratio on a logarithmic scale. For a canonical distribution we expect a linear dependence on the inverse temperature, shown by the dash-dot line, and given by , where is the TLS splitting. Crucially, however, the RCME steady state shows a clear deviation from this linear behaviour with decreasing temperature. To provide context, taking our estimate again of  cm relevant to molecular systems, we see that deviations from canonical statistics begin to become apparent around  K, and should thus be observable even at room temperature in such systems. Note, from the inset, that a non-zero level of coherence is also apparent around such temperatures (and lower), again pointing to the noncanonical nature of the TLS equilibrium state. These departures demonstrate that we cannot represent the TLS steady state as a Gibbs distribution over the system eigenstates. They also quantify relatively simple experimental signatures of the breakdown of canonical statistics in open quantum systems. For example, by measuring only the TLS populations over a reasonable temperature range, we may infer the emergence of noncanonical equilibrium states simply by observing a nonlinear temperature dependence as shown.

Finally, in Fig. 4 (b) we examine the TLS steady state as a function of system-bath coupling strength at constant temperature. We again see significant deviations from the (unchanging) canonical thermal state, such as the development of steady state coherences within the TLS eigenstate basis, which are apparent for any non-vanishing coupling strength. Nevertheless, it is evident from all plots in Fig. 4, by the agreement between Eq. (9) and the HEOM, that the noncanonical steady state is still extremely well described by a thermal state across a wide range of parameters, though now in the mapped TLS-RC representation. This behaviour is a clear result of the TLS-environment correlations; the TLS cannot be considered merely as being in a product state with a (thermal) environment. The non-separability of the steady state thus implies that the equilibrium behaviour of the TLS cannot be described as a thermal distribution over its eigenstates either, but rather, one should also consider states of the environment within the description. This is markedly different to standard master equation techniques and statistical mechanics approaches, where steady states are commonly characterised by the bath temperature and system Hamiltonian alone, an artefact of applying the Born approximation.

## Vi Summary

To summarise, by exploiting a collective coordinate mapping, we have derived a master equation valid in the non-adiabatic regime of the spin-boson model. Notably, besides an accurate description of the system dynamics, our approach also allows us to quantify the accumulation of system-environment correlations with time, as well as probe the dynamic evolution of states of the environment. We have shown that properly accounting for the generation of system-environment correlations is essential for describing both the transient dynamics and equilibrium distributions of open systems in this regime. In particular, we have demonstrated that long-lived correlations lead to the emergence of noncanonical system equilibrium states that can be characterised in a simple and intuitive way within our formalism, as thermal states of the system-collective coordinate Hamiltonian. Our approach can be applied to a number of systems of practical relevance. For example, in molecular and solid-state (e.g. superconducting) devices Leggett et al. (1987); Weiss (2008); Breuer and Petruccione (2002); Le Hur et al. (2007); Le Hur (2012); Würger (1998); Mahan (1990); Goldstein et al. (2013); Nitzan (2006); May and Kühn (2011); Garg et al. (1985); Gilmore and McKenzie (2005); Engel et al. (2007); Lee et al. (2007); Collini and Scholes (2009); Santamore et al. (2013); Franco et al. (2013); Peropadre et al. (2013), system-environment coupling can be strong and memory effects important. Taking parameters relevant for energy transfer processes in molecular dimers, we have shown that deviations from canonical statistics could be observable even at room temperature, which raises the intriguing question of the role of noncanonical equilibrium states in larger molecular aggregates Engel et al. (2007); Lee et al. (2007); Collini and Scholes (2009). Finally, a proper understanding of equilibrium states is a vital component in the growing field exploring the thermodynamics of quantum systems; the full implications of noncanonical steady-states in open quantum systems Lee et al. (2012); Martinez and Paz (2013) thus constitutes a fascinating topic for future exploration.

## Vii Acknowledgements

We would like to thank Alex Chin, Jianshu Cao, Javier Cerrillo, Dara McCutcheon, Tom Stace, and Tommaso Tufarelli for interesting discussions. J. I.-S. is supported by the EPSRC and A. N. by Imperial College and The University of Manchester. N. L. acknowledges the hospitality of the Controlled Quantum Dynamics Group at Imperial College.

## Appendix A Reaction coordinate mapping

As in the main text, we start by considering a two level system (TLS) coupled linearly to a harmonic environment described by the spin-boson (SB) Hamiltonian (with ):

 H=ϵ2σz+Δ2σx+σz∑kfk(c†k+ck)+∑kνkc†kck, (10)

where () are the standard TLS Pauli operators in a basis where , and () are creation (annihilation) operators for bosonic modes of frequency , which couple to the TLS with strength . The system-environment coupling is fully specified by the spin-boson spectral density, .

We now apply a normal mode transformation to Eq. (10), leading to a mapped Hamiltonian of the form Martinazzo et al. (2011)

 HRC =ϵ2σz+Δ2σx+λσz(a†+a)+Ωa†a+∑kωkb†kbk +(a†+a)∑kgk(b†k+bk)+(a†+a)2∑kg2kωk, (11)

where we have defined the collective (reaction) coordinate such that

 λ(a†+a)=∑kfk(c†k+ck). (12)

The TLS-reaction coordinate (RC) coupling strength is given by , such that the RC creation and annihilation operators satisfy the canonical commutation relation . This choice of coupling also fixes the RC frequency to be . The RC is now coupled linearly to a residual harmonic environment characterised by the spectral density . There is also a term quadratic in the system operators in Eq. (A), known as the counter term, which is used to renormalise the mode frequency and avoid divergences due to friction Hartmann et al. (2000).

In order to fully specify the mapping described above, we must relate the RC and SB spectral densities. We do this by following the procedure first outlined by Garg et al. Garg et al. (1985) and derive the spectral density, both before and after the mapping, from the classical equations of motion Leggett (1984). Since the spectral density doesn’t contain any information about the system itself, but rather just the coupling strength between the system and environment, it then follows that in Eq. (10) we can swap the TLS for a continuous classical coordinate moving in a potential . This yields a Hamiltonian of the form

 Hq =P2q2+U(q)+q∑k~fk^xk+q2∑k~f2k2ν2k +12∑k(^p2k+ν2k^x2k), (13)

where, for simplicity, we have written Eq. (A) in the position representation, with coordinate and momentum operators of the environment defined as

 ^xk=√12νk(c†k+ck)%and^pk=i√νk2(c†k−ck), (14)

and the coupling between the TLS and the mode of the environment is given by . From this Hamiltonian, we attain a set of classical equations of motion:

 ¨q(t) =−U′(q)−∑k~fkxk(t)−q(t)∑k~f2kν2k, ¨xk(t) =−~fkq(t)−ν2kxk(t). (15)

We can eliminate the bath variables from the equation of motion for the classical coordinate by making use of the Fourier transform, . This leads to an equation of the form , where the Fourier space operator is defined as

 ~K(z) =−z2⎛⎝1+∑k~f2kν2k(ν2k−z2)⎞⎠ =−z2⎛⎜⎝1+∞∫0dνJSB(ν)ν(ν2−z2)⎞⎟⎠. (16)

We solve the integral using the residue theorem, making use of a single pole at , giving . By writing , it follows that Leggett (1984)

 JSB(ω)=1πlimϵ→0+Im[K(ω−iϵ)]. (17)

We now use the same procedure to write in terms of the RC spectral density. Since the mapping at this stage is exact, then the Fourier transformed operator, , will be identical before and after the normal mode transformation. If we once again swap the TLS for a continuous coordinate, and write the RC Hamiltonian given in Eq. (A) in position space, then we have

 Hq =P2q2+U(q)+κq^x+κ22Ω2q2+12(^p2+Ω2^x2) +^x∑k~gk^Xk+^x2∑k~g2k2ω2k+12∑k(^P2k+ωk^X2k), (18)

where we have scaled the coupling strengths such that and , and the position and momentum operators are defined in the usual way. This Hamiltonian leads to classical equations of motion of the form

 ¨q+κ^x+κ2Ω2q=−U′(q), ¨^x+(Ω2+∑k~g2kω2k)^x+κ2q+∑k~gk^Xk=0, ¨^Xα+~gα^x+ω2α^Xα=0. (19)

By moving to Fourier space and eliminating both the RC and environment from the equation of motion for the classical coordinate, we get the expression for the Fourier space operator

 ~K(z)=−z2−κ2Ω2L(z)Ω2+L(z),

where, by using the definition of the RC spectral density, we have

 L(z) =−z2(1+∑k~g2kω2k(ω2k−z2)) =−z2⎛⎜⎝1+4Ω∞∫0JRC(ω)ω(ω2−z2)dω⎞⎟⎠. (20)

Considering the RC to have an Ohmic spectral density of the form , the integral reduces to , in the limit . By plugging this into Eq. (17), we get the condition

 JSB(ω) =1πlimϵ→0+Im[~K(ω−iϵ)] =1π2πΩγκ2ω(Ω2−ω2)2+(2πΩγ)2ω2. (21)

In order for Eq. (A) to be consistent with the original spin-boson spectral density, we identify the relations

 ωc= Ω2πγ, andα= 2λ2πΩ, (22)

which give us

 JSB(ω)=αωcωω2c+ω4ω2cΩ2+ω2. (23)

Finally, by assuming , we get the spin-boson spectral density in Lorentz-Drude form,

 JSB(ω)=αωcωω2+ω2c. (24)

## Appendix B Reaction coordinate master equation

We now wish to use Eq. (A) to derive a master equation that treats the TLS-RC coupling exactly (up to some number of basis states) while the coupling to the residual environment is treated to second order (within a Born-Markov approximation). To this end, we define the interaction Hamiltonian as

 HI=^A⊗^B+~λ ^A2, (25)

where , and . To derive the RC master equation we first move into the interaction picture using the transformation , where , is the TLS-RC Hamiltonian, and . Tracing out the residual environment, we then find a master equation for the reduced TLS-RC density operator , which in the interaction picture may be written as

 ∂~ρ(t)∂t =−i trB[~HI(t),ρ(0)⊗ρB] −t∫0dτ trB[~HI(t),[~HI(t−τ),~ρ(τ)⊗ρB]], =−i ~λ[^A2(t),ρ(0)] −~λ2t∫0dτ[^A2(t),[^A2(t−τ),~ρ(τ)]] −t∫0dτ([^A(t),[^A(t−τ),~ρ(t)]]Γ+(τ) Missing or unrecognized delimiter for \right (27)

Here, we have made use of the Born approximation, that is, we have assumed that the system and residual environment remain in the product state for all time, where . The correlation functions are defined as . In the continuum limit we may write them as

 Γ+(τ)=∞∫0dωJRC(ω)cothβω2cosωτ (28)

and

 Γ−(τ)=i∞∫0dωJRC(ω)sinωτ. (29)

As outlined in the previous section, the RC spectral density takes the form , in the limit that that the cut-off frequency . We can simplify Eq. (B) further by noticing that

 −i~λ[^A2(t),~ρ(t)] =−i~λ[^A2(t),~ρ(0)] −~λ2t∫0dτ[^A2(t),[^A2(t−τ),~ρ(τ)]]. (30)

By substituting this into the above master equation and assuming a Markov limit, that is, we take the time integrals to infinity, we acquire the RC master equation

 ∂~ρ(t)∂t =−i~λ[^A2(t),~ρ(t)] −∞∫0∞∫0dτdωJRC(ω)cothβω2cosωτ ×[^A(t),[^A(t−τ),~ρ(t)]] −i∞∫0∞∫0dτdωJRC(ω)sinωτ ×[^A(t),{^A(t−τ),~ρ(t)}]. (31)

The definition of assumes an infinite cut-off frequency, hence the first and last terms of Eq. (B) are divergent. However, we can eliminate the divergent contributions by integrating the last term by parts, such that

 ∞∫0dτsinωτ^A(t−τ) =−P(^A(t)ω) +∞∫0dτcosωτω∂^A(t−τ)∂τ. (32)

Using the principal value part of this integral we can cancel the counter term, which, after moving back into the Schrödinger picture, gives the master equation

 ∂ρ(t)∂t =−i[H0,ρ(t)] −∞∫0∞∫0dτdωJRC(ω)cosωτcothβω2 ×[^A,[^A(−τ),ρ(t)]] −∞∫0∞∫0dτdωJRC(ω)cosωτω ×[^A,{[^A(−τ),H0],ρ(t)}]. (33)

In order to derive the interaction picture system operators we shall diagonalise the TLS-RC system Hamiltonian, , numerically. Let be an eigenstate of the system Hamiltonian, such that We can now write the position operators in this eigenbasis:

 ^A=∑jkAj,k∣∣φj⟩⟨φk|, (34)

where . In the interaction picture this becomes

 ^A(t)=∑jkAj,keiξjkt∣∣φj⟩⟨φk|, (35)

where is the difference between the and eigenvalues. Using this definition we can include the rates from Eq. (B) into the operators, such that

 ^χ =∞∫0∞∫0dωdτ JRC(ω)cosωτcothβω2^A(−τ) ≈π2∑jkJRC(ξjk)cothβξjk2Ajk∣∣φj⟩⟨φk|, (36) ^Ξ =∞∫0∞∫0dωdτJRC(ω)cosωτω[H0,^A(−τ)] ≈π2∑jkJRC(ξjk)Ajk∣∣φj⟩⟨φk|, (37)

where we have neglected the imaginary Lamb shift terms. We can now write the master equation as

 ∂ρ(t)∂t=−i[H0,ρ(t)]−[^A,[^χ,ρ(t)]]+[^A,{^Ξ,ρ(t)}]. (38)

## Appendix C Heom

The hierarchical equations of motion (HEOM) are a set of time-local equations for the reduced system dynamics, governed by the spin-boson Hamiltonian [Eq. (10)], which capture the bath dynamics and system-bath correlations through a set of auxiliary density matrices. These equations are exact under the assumption of a Lorentz-Drude spectral density, as given in Eq. (24), and an initially separable system-bath state at . Here we assume only a single bath such that the HEOM can be written as

 ˙ρn = −(i[HS,ρn]+K∑m=0nmμm)ρn−iK∑m=0[Q,ρn+m] (39) − iK∑m=0nm(cmQρn−m−c∗mρn−mQ) − Missing or unrecognized delimiter for \right

where is the system-environment coupling operator, , and . The bath correlation functions for the Lorentz-Drude spectral density are

 C=∞∑m=0cmexp(−μmt). (40)

Here , , and the coefficients are

 c0=ωcαH2[cot(βωc/2)−i]/ℏ (41)

and

 cm≥1=2αHωcβμmμ2m−ω2c. (42)

Each density matrix is labelled by an index of positive integers . As here we only have a single bath the integers are defined as . For each “Matsubara” term each index runs from to . The null label defines the system density matrix, and any non-zero label refers to an auxiliary density matrix which encodes the correlations with the bath. The terms in the equation of motion refer to an increase or decrease by of the label index . We take a cut-off in the overall tier and in the Matsubara terms which give convergence in the numerical results.

## Appendix D Quantum mutual information and non-Gaussianity

In this Appendix we show that the mutual information for the TLS and RC acts as a lower bound to the correlations shared between the system and the original multi-mode environment.

Let be the density matrix describing the state of both the system and environment in the original spin-boson representation. We also define the reduced states of the system, , and the environment, . Let be the unitary transformation that maps the spin-boson model to the RC model. Applying this unitary to the reduced state of the system has no effect due to the trace over the environment. However, on the reduced state of the environment this unitary transforms the spin-boson basis to that of the RC and residual bath, that is, .

The quantum mutual information for the system and the spin-boson environment is given by:

 I(ρs:ρE)=S(ρs)+S(ρE)−S(χ),=S(ρs)+S(ρRC+B)−S(~χ), (43)

where is the von-Neumann entropy and is the total state in the RC basis. Here we have used the unitary equivalence of the von-Neumann entropy to write the quantum mutual information in terms of the RC basis.

To proceed we shall make use of the strong subadditivity of the von-Neumann entropy, that is:

 S(~χ)+S(ρRC)≤S(ρS+RC)+S(ρRC+B), (44)

where , with the trace taken over the residual environment. Using this property in conjunction with Eq. (43) gives:

 I(ρs:ρE) ≥S(ρS)+S(ρRC+B)+S(ρRC)−S(ρs+RC) −S(ρRC+B), ≥S(ρs)+S(ρRC)−S(ρS+RC). (45)

Therefore, we have the condition

 I(ρs:ρE)≥I(ρs:ρRC). (46)

Hence, the mutual information between the system and RC acts as a lower bound to the mutual information between the system and spin-boson environment. Furthermore, in the limit that the Born approximation holds between the composite system (TLS and RC) and the residual environment, that is , then the inequality in Eq. (46) becomes an equality due to the additive nature of the von-Neumann entropy.

Similarly, the non-Gaussianity can be shown to be invariant under symplectic transformation (i.e. operators that are quadratic in field operators) and monotonically decreases under partial trace Genoni et al. (2008). This means that the non-Gaussianity of the RC acts as a rigorous lower bound for the non-Gaussianity of the original spin-boson environment, that is:

 δG[ρE]≥δG[ρRC]. (47)

### References

1. A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. Fisher, A. Garg,  and W. Zwerger, Rev. Mod. Phys. 59, 1 (1987).
2. U. Weiss, Quantum Dissipative Systems, 3rd ed. (World Scientific, Singapore, 2008).
3. A. Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer and Reactions in Condensed Molecular Systems (Oxford University Press, Oxford, 2006).
4. V. May and O. Kühn, Charge and Energy Transfer Dynamics in Molecular Systems, 3rd ed. (Wiley, Berlin, 2011).
5. H. P. Breuer and F. Petruccione, The Theory of Open Quantum Systems (Oxford University Press, Oxford, 2002).
6. The most common being the Born-Markov approximations, wherein the bath is assumed to be static, memoryless, and uncorrelated with the system, remaining in thermal equilibrium throughout the dynamical evolution Breuer and Petruccione (2002).
7. A. Garg, J. N. Onuchic,  and V. Ambegaokar, J. Chem. Phys. 83, 4491 (1985).
8. M. Thoss, H. Wang,  and W. H. Miller, J. Chem. Phys. 115, 2991 (2001).
9. J. Cao and G. A. Voth, J. Chem. Phys. 106, 1769 (1997).
10. K. H. Hughes, C. D. Christ,  and I. Burghardt, J. Chem. Phys. 131, 124108 (2009).
11. L. Hartmann, I. Goychuk,  and P. Hänggi, J. Chem. Phys. 113, 11159 (2000).
12. J. Roden, W. T. Strunz, K. B. Whaley,  and A. Eisfeld, The Journal of Chemical Physics 137, 204110 (2012).
13. R. Martinazzo, B. Vacchini, K. H. Hughes,  and I. Burghardt, J. Chem. Phys. 134, 011101 (2011).
14. M. P. Woods, R. Groux, A. W. Chin, S. F. Huelga,  and M. B. Plenio, Journal of Mathematical Physics 55, 032101 (2014).
15. A. Imamoğlu, Phys. Rev. B 50, 3650 (1994).
16. B. M. Garraway, Phys. Rev. B 55, 2290 (1997).
17. Y. Tanimura, Phys. Rev. A. 41, 6676 (1990).
18. Y. Tanimura and R. Kubo, J. Phys. Soc. Jpn. 58, 101 (1989).
19. S. Bera, S. Florens, H. U. Baranger, N. Roch, A. Nazir,  and A. W. Chin, Phys. Rev. B 89, 121108(R) (2014a).
20. S. Bera, A. Nazir, A. W. Chin, H. U. Baranger,  and S. Florens, Phys. Rev. B 90, 075110 (2014b).
21. N. Makri and D. E. Makarov, J. Chem. Phys. 102, 4600 (1995).
22. J. Prior, A. W. Chin, S. F. Huelga,  and M. B. Plenio, Phys. Rev. Lett. 105, 050404 (2010).
23. A. Alvermann and H. Fehske, Phys. Rev. Lett. 102, 150601 (2009).
24. H. Wang and M. Thoss, New J. Phys. 10, 115005 (2008).
25. F. B. Anders, R. Bulla,  and M. Vojta, Phys. Rev. Lett. 98, 210402 (2007).
26. A. Winter, H. Rieger, M. Vojta,  and R. Bulla, Phys. Rev. Lett. 102, 030601 (2009).
27. R. Silbey and R. A. Harris, J. Chem. Phys. 80, 2615 (1984).
28. D. P. S. McCutcheon and A. Nazir, J. Chem. Phys. 135, 114501 (2011).
29. S. Jang, Y.-C. Cheng, D. R. Reichman,  and J. D. Eaves, J. Chem. Phys. 129, 101104 (2008).
30. A. Nazir, Phys. Rev. Lett. 103, 146404 (2009).
31. A. J. Leggett, Phys. Rev. B 30, 1208 (1984).
32. J. Gilmore and R. H. McKenzie, J. Phys.: Condens. Matter 17, 1735 (2005).
33. J. Hausinger and M. Grifoni, Phys. Rev. A 81, 022117 (2010).
34. A. Ishizaki and G. R. Fleming, J. Chem. Phys. 130, 234111 (2009).
35. M. Thorwart, J. Eckel, J. Reina, P. Nalbach,  and S. Weiss, Chem. Phys. Lett. 478, 234 (2009).
36. F. A. Pollock, D. P. S. McCutcheon, B. W. Lovett, E. M. Gauger,  and A. Nazir, New J. Phys. 15, 075018 (2013).
37. We see similarly excellent agreement for the coherences and for as well.
38. M. G. Genoni, M. G. A. Paris,  and K. Banaszek, Phys. Rev. A 78, 060303 (2008).
39. M. G. Genoni and M. G. A. Paris, Phys. Rev. A 82, 052341 (2010).
40. B. Groisman, S. Popescu,  and A. Winter, Phys. Rev. A 72, 032317 (2005).
41. K. Le Hur, P. Doucet-Beaupre,  and W. Hofstetter, Phys. Rev. Lett. 99, 126801 (2007).
42. K. Le Hur, Phys. Rev. B 85, 140506 (2012).
43. A. Würger, Phys. Rev. B. 57, 347 (1998).
44. G. D. Mahan, Many-Particle Physics (Plenum, New York, 1990).
45. M. Goldstein, M. H. Devoret, M. Houzet,  and L. I. Glazman, Phys. Rev. Lett. 110, 017002 (2013).
46. G. Engel, T. Calhoun, E. Read, T. Ahn, T. Mančal, Y. Cheng, R. Blankenship,  and G. Fleming, Nature 446, 782 (2007).
47. H. Lee, Y.-C. Cheng,  and G. R. Fleming, Science 316, 1462 (2007).
48. E. Collini and G. D. Scholes, Science 323, 369 (2009).
49. D. H. Santamore, N. Lambert,  and F. Nori, Phys. Rev. B 87, 075422 (2013).
50. R. L. Franco, B. Bellomo, S. Maniscalco,  and G. Compagno, International Journal of Modern Physics B 27, 1345053 (2013).
51. B. Peropadre, D. Zueco, D. Porras,  and J. J. García-Ripoll, Phys. Rev. Lett. 111, 243602 (2013).
52. C. K. Lee, J. Cao,  and J. Gong, Phys. Rev. E 86, 021109 (2012).
53. E. A. Martinez and J. P. Paz, Phys. Rev. Lett. 110, 130406 (2013).
72369