From Dissipative Dynamics to Studies of Heat Transfer at the

Subscriber access provided by McMaster University Library

Article

From Dissipative Dynamics to Studies of Heat Transfer at the Nanoscale: Analysis of the spin-boson model Nazim Boudjada, and Dvira Segal J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/jp5091685 • Publication Date (Web): 24 Oct 2014 Downloaded from http://pubs.acs.org on October 28, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

From Dissipative Dynamics to Studies of Heat Transfer at the Nanoscale: Analysis of the Spin-Boson Model Nazim Boudjada†,‡ and Dvira Segal∗,‡ Ecole Polytechnique de Montreal, Montreal, Quebec, Canada H3C 3A7, and Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6 E-mail: [email protected]

KEYWORDS: nonelectronic transport; energy transport; heat conduction, quantum dissipation

∗ To

whom correspondence should be addressed Polytechnique de Montreal, Montreal, Quebec, Canada H3C 3A7 ‡ Chemical Physics Theory Group, Department of Chemistry, University of Toronto, 80 Saint George St. Toronto, Ontario, Canada M5S 3H6 † Ecole

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Abstract We study in a unified manner the dissipative dynamics and the transfer of heat in the twobath spin-boson model. We use the Bloch-Redfield (BR) formalism, valid in the very weak system-bath coupling limit, the noninteracting-blip approximation (NIBA), applicable in the non-adiabatic limit, and iterative, numerically-exact path integral tools. These methodologies were originally developed for the description of the dissipative dynamics of a quantum system, and here they are applied to explore the problem of quantum energy transport in a non-equilibrium setting. Specifically, we study the weak-to-intermediate system-bath coupling regime at high temperatures kB T /¯h > ε , with ε as the characteristic frequency of the two-state system. The BR formalism and NIBA can lead to close results for the dynamics of the reduced density matrix (RDM) in a certain range of parameters. However, relatively small deviations in the RDM dynamics propagate into significant qualitative discrepancies in the transport behavior. Similarly, beyond the strict non-adiabatic limit NIBA’s prediction for the heat current is qualitatively incorrect: It fails to capture the turnover behavior of the current with tunneling energy and temperature. Thus, techniques that proved meaningful for describing the RDM dynamics, to some extent even beyond their rigorous range of validity, should be used with great caution in heat transfer calculations, since qualitative-serious failures develop once parameters are mildly stretched beyond the techniques’ working assumptions.

1 Introduction Quantum impurity models, comprising a subsystem in an environment, embody complex processes in condensed phases: electron and exciton transfer in solids, solutions, glasses and biomolecules, 1,2 screening of a magnetic impurity by the Fermi sea electrons, 3 electronic conduction of molecules, 2 and the decoherence behavior of superconducting qubits. 3,4 The dissipative dynamics of impurity models has been explored intensively by time-evolving the subsystem (reduced) density matrix, revealing mechanisms of decoherence and relaxation towards equilibrium. Here, as a case study, we focus on the spin-boson (SB) model with a two-level system (TLS) immersed in a bath of


Page 2 of 46

Page 3 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


harmonic oscillators. 1 Beyond the question of decoherence and dissipation, impurity models can be employed for exploring fundamentals of quantum transport and quantum thermodynamics, when placing the subsystem between two reservoirs maintained e.g. at different chemical potentials or temperatures. In this scenario, the reservoirs exchange charge, spin, or energy carriers through the subsystem, with quantities of interest as (charge, energy, spin) currents in the system, as well as high order cumulants of currents. In the context of nanoscale heat transfer and phononics, 5–7 the “non-equilibrium spin-boson model" (NESB), with a TLS bridging two thermal reservoirs at different temperatures, has been suggested as a toy model for studying the phenomenology of quantum heat transfer in anharmonic junctions, 8,9 see Fig. 1 for a schematic representation. The dissipative dynamics of the SB model has been examined systematically by comparing predictions from different techniques. Results have been organized in several reviews, 1,3 and the problem still provides an active area for exploration, see for example Refs. 10–14 In contrast, the analysis of heat transfer characteristics in the corresponding NESB model is a relatively new problem and a systematic comparison of results from different techniques is still missing. One should note that the computation of transport characteristics in the NESB model (and other non-equilibrium impurity models) relies on a nontrivial extension of open quantum systems methodologies: To calculate the current one needs to follow the dynamics of other operators beyond the reduced density matrix (RDM): two-time correlation functions of subsystem’s operators or expectation values of bath operators. The thermal properties of the NESB nanojunction have been analyzed on the basis of perturbative quantum master equations, 8,9,15–20 Keldysh Green’s function expansions, 21–24 and the noninteracting-blip approximation. 8,25,26 Numerically exact techniques, developed for the study of the (single-bath) SB model, were similarly generalized to explore transport properties: the multilayer multiconfiguration time-dependent Hartree theory, 27 influence functional path integral techniques 28 and Monte-Carlo simulations. 29 Theoretical studies of heat flow in model systems such as the NESB nanojunction are motivated by recent experiments of thermal energy flow across alkane chains, 30–33 proteins 34 and small aro-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

matic molecules. 35–38 These studies aim in exploring the role of vibrational energy flow in e.g. chemical reaction dynamics, protein folding, conformational changes, and molecular electronics. Questions of phononic heat transfer are of great interest in other disciplines. For example, in thermoelectric applications reducing the phononic contribution to the thermal conductivity improves the (heat to work) conversion efficiency; recent experiments reached reduced thermal conductivities in nanocomposites. 39 It is particularly interesting to design a molecular-level or a nanoscale thermal diode, optimally commanding unidirectional energy flow. This would allow control over molecular reactivity, and potentially turn into a building block in phononic (and similarly photonic) devices. Unidirectional heat flow was recently demonstrated in nitrobenzene: 36 Using ultrafast infrared Raman spectroscopy it was shown that energy transfer from the nitro to the phenyl group, or from the nitro to global vibrational modes, was blocked. However, vibrational energy was transferred from the phenyl-localized modes to the nitro modes and to global modes. It is now established that control over quantum energy flow can be achieved by combining many-body interactions with spatial asymmetries. 7–9 Anharmonicity of vibrational modes is thus an essential ingredient for building nontrivial functionalities. While more detailed calculations are imperative to explore particular systems, 40,41 the NESB model with a two-level system, a truncated harmonic vibration, is the simplest-nontrivial model which can allow us to explore the role of anharmonic (many-body) effects in phononic (or photonic) conduction. In this work, we are interested in the problem of quantum heat transfer in anharmonic nanojunctions, particularly when the central object’s coupling energy to the contacts is substantial. Our goal is to examine and compare different techniques, understand their range of validity, and find out when they provide qualitatively correct results in comparison to exact numerical techniques. Focusing on the NESB model, we aim in rectifying the following points: (i) Relation between dissipative dynamics and transport. We study here both the RDM dynamics and the transfer of heat in the NESB model using the weak-coupling (system-bath) BlochRedfield (BR) formalism, the noninteracting-blip approximation (NIBA), valid in the non-adiabatic limit and at high temperatures, and numerically-exact influence functional path integral simula-


Page 4 of 46

Page 5 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


tions. The BR and NIBA techniques provide consistent results when describing the dynamics of the RDM in a certain range of parameters. Does this agreement translate into consistent transport properties? The answer is negative. We show here that even when the BR and NIBA techniques reasonably agree on the RDM dynamics, results significantly deviate when following the heat current behavior: The BR scheme fails in providing the current characteristics, qualitatively and quantitatively, beyond the very weak coupling limit. Similarly, beyond the strict non-adiabatic limit NIBA badly fails in describing transport trends, while it still performs reasonably well in RDM calculations. (ii) Developing approaches for weak-intermediate coupling cases. The BR method administers quantum kinetic equations, and it provides a transparent theory for thermal conduction: a linear enhancement of current with increasing coupling energy to the contacts. Other methods reveal that this trend breaks down immediately beyond the very weak coupling limit. 21,24,25,28,29 However, a careful comparison between different techniques is missing. To study physical situations, e.g., with the molecule moderately or strongly attached to thermal contacts as in Ref., 32,33 it is imperative to develop reliable methodologies that can extend beyond the very weak coupling regime. We play here with four different approaches, BR, 8,9,16,18 NIBA, 8,25 perturbative techniques based on the non-equilibrium Green’s function (NEGF), 21,24 and numerically exact influence functional path-integral simulations. 28 With these tools we study the current characteristics as a function of the contact interaction, as well as the temperature and the frequency of the TLS, and observe a nontrivial non-monotonic performance of the junction, exposing the underlying mechanisms of thermal conduction. The paper is organized as follows. In Sec. 2 we present the model and observables of interest: the reduced density matrix and the heat current, including the linear response coefficient, the thermal conductance. In Sec. 3 we lay down the Bloch-Redfield equations for the RDM and the steady-state current. Sec. 4 presents the corresponding NIBA equations. Sec. 5 describes influence functional path integral approaches. Numerical results for the RDM dynamics and steady-state heat current are included in Sec. 6. In Sec. 7 we summarize our work.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 46

2 Model and Observables of Interest The NESB model includes a two-state system (spin) bridging two bosonic reservoirs (ν = L, R). In the “local" basis (|0i and |1i) the isolated spin Hamiltonian reads

H0 =

h¯ ∆ h¯ ω0 σz + σx , 2 2

(1)

and the total Hamiltonian is given by

H = H0 + ∑ ν ,k

h¯ σz † † λk,ν (bk,ν + bk,ν ) + h¯ ωk bk,ν bk,ν . 2

(2)

The Pauli matrices are defined as σz = |1ih1| − |0ih0|, σx = |0ih1| + |1ih0|, and σy = −i|1ih0| + i|0ih1|, h¯ ω0 is the level detuning (bias), ∆ stands for the tunneling frequency between the spin states, and b†k,ν (bk,ν ) is the creation (annihilation) operator of a boson (e.g. phonon) with a wavenumber k in the ν reservoir. The interaction of the subsystem with the baths can be characterized by a spectral density function, defined as

Jν (ω ) = ∑ λk,2 ν δ (ω − ωk ).

(3)

k

We perform our numerical simulations using an Ohmic form, Jν (ω ) = 2αν ω e−ω /ωc ,

(4)

but the methodologies discussed below can handle other spectral functions. Here αν is a dimensionless interaction parameter between the spin subsystem and the ν reservoir. Below we use the definition α ≡ αL + αR . For simplicity, the cutoff frequency ωc is taken identical in both baths. In the context of electron transfer processes it is useful to define the reorganization energy Erν ≡ d ω Jν (ω )/ω . For Ohmic functions it reduces to Erν = 2αν ωc . R

We principally work here in the so-called non-adiabatic limit of ∆/ωc ≪ 1 and temperatures 6 ACS Paragon Plus Environment

Page 7 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


kB T /¯h∆ & 1. The “non-adiabatic" terminology is delivered from studies of electron transfer reactions in condensed phases, in which the TLS represents electron-donor and acceptor states with a tunneling frequency ∆: “adiabatic processes" ∆ > ωc refer to reactions with fast tunneling electrons relative to the phonon bath. In the opposite non-adiabatic limit ∆ ≪ ωc the characteristic time scale of the bath 1/ωc is short relative to the internal timescale for tunneling. Based on our heat transport results in the non-adiabatic limit we identify four regions: (i) we refer to αν < 0.025 as the very weak coupling regime, (ii) 0.025 < αν < 0.1 corresponds to the weak coupling limit, (iii) 0.1 < αν < 0.5 describes the intermediate regime, and (iv) αν > 0.5 corresponds to the strong coupling limit. The reduced density matrix is organized from the population difference hσz (t)i and the real and imaginary parts of the coherence hσx,y (t)i, hσi (t)i = tr{ρ (t = 0)eiHt/¯h σi e−iHt/¯h },

(5)

with ρ (t = 0) as the initial state of the total density matrix. In what follows we assume a factorized initial condition, ρ (t = 0) = ρL ⊗ ρR ⊗ |1ih1|. The two reservoirs, Hν = ∑k h¯ ωk b†k,ν bk,ν , are separately prepared in a canonical-equilibrium state of temperature Tν = 1/(kB βν ),

ρν =

e−βν Hν , Zν = Trν e−βν Hν . Zν

(6)

At t = 0 the subsystem-bath interaction is turned on and we wait for the (assuming unique) steadystate solution to set in. The heat current can be computed from the transient regime to the steadystate limit; below we focus only on the long-time behavior. It is reached by considering energy leakage at the contacts. For example, at the left contact the heat current operator is defined as dHL −i jˆL ≡ − = [H, HL ]. h¯ dt


(7)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 46

Tracing over all degrees of freedom we write the expectation value

jL ≡ tr[ρ (t = 0) jˆL (t)].

(8)

Operators are written here in the Heisenberg representation. In steady-state, jq = jL = − jR , with the convention of a positive current flowing L to R. In the linear response regime the current is expanded to the lowest order in the temperature difference, jq ∼ κ (TL − TR ), and we obtain the thermal conductance κ from the relation

κ≡

d jq . dTL TL →TR =T

(9)

To practically compute the heat current we should manipulate Eq. (8) into a workable definition. A formally-exact construction has been derived in Ref. 21 from the perturbation expansion of the nonequilibrium Green’s function. This formula, a many-body extension of Landauer’s expression, 42 an analog of the Meir-Wingreen formula for electronic systems, 43 expresses the heat current of the NESB model in correlation functions of the TLS. The linear response limit of this formula was recently examined using Monte-Carlo simulations to explore signatures of Kondo physics in thermal conduction. 29 In the following sections we describe the evaluation of the heat current in the NESB model under different sets of approximations.

Figure 1: Scheme of the non-equilibrium spin-boson model, a minimal picture for studying heat transport in anharmonic nanojunctions. In this work we simulate both the dynamics of the (spin) reduced density matrix ρS and the heat current behavior jq using approximate methods and numerically exact simulation tools.


Page 9 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3 Bloch-Redfield Formalism: Very weak coupling regime 3.1 Hamiltonian The standard Bloch-Redfield equation can be derived from the exact quantum master equation based on the assumption of weak system-bath interactions. 44 It is convenient to develop it in the “energy" basis, the representation in which the spin subsystem is diagonal,

HS =

∑ En|nihn|.

(10)

n=±

The total Hamiltonian should be organized in an additive L-R structure,

H˜ = HS + HL + HR +VL +VR ,

(11)

with the thermal baths denoted each by Hν and the system-bath interaction given in a direct-product form, Vν = Sν ⊗ Bν , Sν =

′ ν ′ |nihn |. ∑′ Sn,n

(12)

n,n

Here Bν and Sν are bath and subsystem operators, respectively. We label the subsystem operator by the index ν = L, R; the impurity may couple to the two baths via distinct operators. Specifically to the NESB model, we diagonalize the isolated TLS of Eq. (1) with a rotation i matrix U = e− 2 θ σy , tan θ = ∆/ω0 . The total Hamiltonian (2) transforms into H˜ ≡ U † HU, with

h¯ ε σ˜ z + ∑ h¯ ωk b†k,ν bk,ν H˜ = 2 ν ,k +

1 (¯hσ˜ z cos θ − h¯ σ˜ x sin θ ) ∑ λk,ν (b†k,ν + bk,ν ), 2 ν ,k


(13)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

where ε =

q

Page 10 of 46

ω02 + ∆2 and σ˜ i are the Pauli matrices in the energy basis, denoted here by |±i. We

now identify the operators of Eq. (12) by Sν = σ˜ z cos θ − σ˜ x sin θ , Bν =

1 h¯ λk,ν (b†k,ν + bk,ν ). ∑ 2 k

(14)

The dynamics of the reduced density matrix under the BR equation has been examined in numerous studies, see for example Refs. 45,46 For completeness, we include in Sec. 3.2 relevant expressions using a notation similar to that employed in Ref. 10 The heat current in the NESB model was only recently derived in a closed form at the level of the BR scheme. 16,18 A workable expression is provided in Sec. 3.3.

3.2 Reduced density matrix In the BR scheme the interaction Vν is treated perturbatively, to the lowest nontrivial order. This results in a master equation for the spin RDM, an integro-differential equation, written here in the local-site basis of Eq. (2), 10,45,46 d hσz (t)i = ∆hσy (t)i, dt Z t Z t d hσy (t)i = ω0 hσx (t)i − ∆hσz (t)i − d τ Gy (τ ) − d τ [Gyy (τ )hσy (t − τ )i + Gyx (τ )hσx (t − τ )i] , dt 0 0 Z t Z t d (15) hσx (t)i = −ω0 hσy (t)i − d τ Gx (τ ) − d τ [Gxx (τ )hσx (t − τ )i − Gyx (τ )hσy (t − τ )i] . dt 0 0 The kernels satisfy

Gx (t) =

∆ ω0 ∆ sin (ε t) ∑ Mν′′ (t), Gy (t) = 2 [1 − cos (ε t)] ∑ Mν′′ (t) ε ε ν ν

Gxx (t) = cos (ε t) ∑ Mν′ (t), Gyy (t) = ν

∆2 + ω02 cos (ε t) Mν′ (t) ∑ 2 ε ν

ω0 sin (ε t) ∑ Mν′ (t), Gyx (t) = ε ν

(16)


Page 11 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and the dissipative terms enclose the correlation function hBν (t)Bν (0)iν = Mν′ (t) − iMν′′ (t), with Mν′ (t) = Mν′′ (t) =

Z ∞ Z0

d ω Jν (ω ) coth(βν h¯ ω /2) cos(ω t),

∞

0

d ω Jν (ω ) sin(ω t).

(17)

Under the Markov approximation, we write a time-local Markovian equation dhσi i/dt = ∑i, j Di, j σ j , see e.g., 10,46,47 and solve it in the long time limit. The equilibrium (eq) solution, in the case of a single bath, satisfies hσz ieq = − ωε0 tanh(β h¯ ε /2), hσx ieq = − ∆ε tanh(β h¯ ε /2), and hσy ieq = 0. In the non-equilibrium two-bath scenario the steady-state (ss) solution of the Markovian equation (reached at t → ∞) satisfies

ω0 JL (ε ) + JR (ε ) , ε JL (ε ) coth(βL h¯ ε /2) + JR (ε ) coth(βR h¯ ε /2) ∆ = hσz iss , ω0

hσz iss = − hσx iss

hσy iss = 0,

(18)

reducing to the correct equilibrium state. These elements build-up the reduced density matrix,

ρSss =





hσx iss − ihσy iss  1 + hσz iss 1 .  2 hσx iss + ihσy iss 1 − hσz iss

We now introduce a short notation, x ≡

JL (ε )+JR (ε ) JL (ε ) coth(βL h¯ ε /2)+JR (ε ) coth(βR h¯ ε /2) ,

(19)

and condense the RDM

into

ρSss

x 1 = Iˆ − 2 ε

∆ ω0 σz + σx , 2 2


(20)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 46

where Iˆ stands for the identity matrix. In the energy basis |±i, ρ˜ Sss = U † ρSssU, resulting in the simple form 1 x ρ˜ Sss = Iˆ − σ˜ z . 2 2

(21)

Explicitly, the population of the two levels in the energy basis follows

p˜ss − =

JL (ε )[1 + nL (ε )] + JR (ε )[1 + nR (ε )] , JL (ε )[1 + 2nL (ε )] + JR (ε )[1 + 2nR (ε )]

ss p˜ss + = 1 − p˜− ,

(22)

with the Bose-Einstein function nν (ε ) = [eβν h¯ ε − 1]−1 . At thermal equilibrium, TL = TR , the spin occupation depends on the temperature of the bath and the spin splitting, but it does not carry information on the coupling strength of the spin to the bath. In contrast, it is significant to observe signatures of the non-equilibrium situation TL 6= TR in the levels’ population, now controlled by the spectral functions of the baths. Thus, out-of-equilibrium the coupling energy to the contacts not only determines the relaxation rate towards steady-state, but it further dictates the steady-state solution. The functional form (22) has been obtained in Refs. 8,25,29 considering the unbiased spin-boson model, ω0 = 0. Thus, the seemingly more complex-biased model does not expose a nontrivial ω0 -dependent controllability. Equations (15) and (20) complete our discussion of the RDM under the Bloch-Redfield formalism: dynamics and steady-state solution. In the next subsection we use these expressions and calculate the steady-state heat current in the NESB model.

3.3 Heat Current A closed expression for the quantum heat current in multi-state nanojunctions has been derived in Ref. 16 under a second order perturbation expansion in the system-bath coupling and the rotating wave approximation. This result was recently extended in Ref. 18 to include anti-rotating wave


Page 13 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


terms, transient effects, and lamb shifts of energies. This derivation is not repeated here; we only review its principles. Briefly, the expectation value of the heat current operator is attained from the definition (7) by time-evolving the density matrix in the energy basis to first order in the systembath interaction term Vν . The overall result is second-order in the interaction parameter since the definition of the heat current itself includes the system-bath interaction operator. We trace the current operator over the baths with a factorized system-bath initial density matrix and take the long time limit. The resulting Bloch-Redfield-type steady-state heat current formula is given by the simple form 18 jν = TrS [ρ˜ Sss Aν ],

(23)

where the relevant RDM is provided by Eq. (21). The matrix Aν depends on the properties of the subsystem; the ν index marks the terminal in which the current is calculated, ν ν Sl, j El, j wνj→l + El,k (wνk→l )∗ . Aνk, j = ∑ Sk,l

(24)

l

Here Ek,l = Ek − El , with the subsystem eigenenergies Eq. (10). The rate constants are given by half-range Fourier transforms of bath correlation functions, wνj→l

1 ∞ = 2 d τ e−iEl, j τ /¯h hBν (τ )B(0)iν h¯ Z 0 Z ∞ 1 ∞ d τ e−iEl, j τ /¯h = d ω Jν (ω ) eiωτ nν (ω ) + e−iωτ (nν (ω ) + 1) , 4 0 0 Z

(25)

where the real part satisfies

ℜ[wνj→l ] =

   nν (El, j /¯h)

if El > E j π Jν (|El, j |/¯h) ×  4  [nν (E j,l /¯h) + 1] if El < E j .


(26)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 46

Since in the energy basis the reduced density matrix of the NESB model is diagonal in steady-state, see Eq. (21), the current in Eqs. (23)-(24) immediately simplifies to ν 2 jν = 2 ∑ |Sk,l | El,k (ρ˜ Sss )k,k ℜ[wνk→l ].

(27)

l,k

Explicitly, for the two-state system this brings out an intuitive structure, ν ν ν ss jν = 2|S−,+ |2 E+,− p˜ss − ℜ[w−→+ ] − p˜+ ℜ[w+→− ] ss = (¯hε ) sin2 (θ )Γν (ε ) p˜ss − nν (ε ) − p˜+ [nν (ε ) + 1] .

(28)

It describes the current, say at the L contact, by the net process of an L-bath induced excitation, multiplied by the ground state population and the energy difference E+,− = E+ − E− , and the relaxation from the excited state, to dispose energy in the L bath. The second line was reached from the definition Γν (ε ) ≡ π2 Jν (ε ), and by identifying the energy difference E+,− = h¯ ε . Employing the steady-state occupations (22) we obtain the closed-form result

jq = (¯hε ) sin2 (θ )

ΓL (ε )ΓR (ε ) [nL (ε ) − nR (ε )]. ΓL (ε )[1 + 2nL (ε )] + ΓR (ε )[1 + 2nR (ε )]

(29)

The heat current can be expanded in orders of ∆T = TL − TR to yield the thermal conductance, a linear response coefficient (9), ΓR (ε )ΓL (ε ) eh¯ εβ (¯hε )2 . κ = sin θ kB T 2 [ΓR (ε ) + ΓL (ε )][1 + 2n(ε )] (eh¯ εβ − 1)2 2

(30)

Here n(ω ) = [eβ h¯ ω − 1]−1 denotes the Bose-Einstein distribution function at the inverse temperature β = 1/(kB T ). Note that ∆2 /ε 2 = sin2 θ . So far our discussion did not assume a particular spectral function. We now employ the Ohmic form (4) with a large cutoff (non-adiabatic limit),


Page 15 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ωc ≫ ε , and receive 29

κ

1 (¯h∆)2 αL αR kB T 2 (αL + αR ) 2 sinh(¯hεβ ) h¯ π ∆2 αL αR β h¯ ε ≪1 . −−−−→ 2 T αL + αR =

πε

(31)

In the classical high temperature limit the thermal conductance is identical for biased (ω0 6= 0) and unbiased (ω0 = 0) models. Furthermore, the prefactor ∆2 evinces on the expected agreement of the BR expression with the NIBA formalism for Ohmic baths, the latter technique strictly holds only in the non-adiabatic regime (∆ ≪ ωc ). The approach discussed here provides the heat current to second order in the system-bath interaction parameter. This scheme can be systematically extended beyond that by working out to higher orders the Dyson expansion of the reduced density matrix equations and the heat current formula. Numerical results to the fourth order in the system-bath coupling parameter are included in Ref., 20 revealing intricate features.

4 Noninteracting-Blip Approximation: Strong Coupling NIBA equations were derived for describing spin polarization dynamics in the SB model based on a path integral influence functional formalism. 1,3 Alternatively, these equations can be recovered by transforming the SB Hamiltonian to the shifted-polaron representation, then working out a second order perturbation theory expansion of the RDM in the dressed tunneling frequency. 45,48 The NIBA scheme is useful in the non-adiabatic limit since it includes only the lowest order (nontrivial) term in a ∆/ωc expansion. 1 For Ohmic baths it provides a good approximation to the polarization dynamics hσz (t)i in the unbiased model at weak damping. It can also faithfully simulate the SB dynamics (polarization and coherences) at strong system-bath interactions, at high temperatures. NIBA fails badly at zero temperature and in biased situations. Particularly, in the SB model it predicts that an infinitesimally small bias drives a strict localization effect in the long time limit,



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 46

hσz i = −1, at zero temperature. In a series of recent studies we had extended NIBA to the out-of-equilibrium regime 8,15,25,26 aiming in simulating heat transport in the NESB nanojunction beyond the BR weak coupling limit. This was achieved by writing down the cumulant generating function of the system, to derive a NIBA-formula for the heat current. 8,25 Recently, we proved that this approximate NIBA expression agrees with numerically-exact simulations of the thermal conductance 29 for Ohmic reservoirs, in the high temperature limit. In Sec. 4.1 we include equations of motion for the spin RDM under NIBA. 1 The NIBA heat current expression was first constructed in Ref.; 8 it was later formally derived from a countingstatistics approach in Ref. 25 In Sec. 4.2 we include relevant expressions.

4.1 Reduced density matrix The exact-formal series expansion of hσz (t)i in ∆ can be organized as a generalized master equation. The dynamics of off-diagonal terms hσx,y (t)i is obtained from the polarization hσz (t)i by exact integral relations, 1 dhσz i = − dt

Z t 0

Ks,z (t − τ )hσz (τ )id τ −

hσx i =

Z t

hσy i =

1 dhσz i . ∆ dt

0

Ks,x (t − τ )d τ −

Z t 0

Z t 0

Ka,z (t − τ )d τ ,

Ka,x (t − τ )hσz (τ )id τ , (32)


Page 17 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The (exact) kernels in Eq. (32) can be truncated to include lowest-nontrivial terms in ∆. This scheme, termed the noninteracting-blip approximation, results in, 1 ′

Ks,z (t) = ∆2 e−Q (t) cos[Q′′ (t)] cos(ω0t), ′

Ka,z (t) = ∆2 e−Q (t) sin[Q′′ (t)] sin(ω0t), ′

Ks,x (t) = ∆e−Q (t) cos[Q′′ (t)] sin(ω0t), ′

Ka,x (t) = ∆e−Q (t) sin[Q′′ (t)] cos(ω0t).

(33)

The function Q(t) = ∑ν Qν (t), Qν (t) = Qν′ (t) + iQν′′ (t) contains real and imaginary components with Z ∞

Jν (ω ) [1 − cos(ω t)][1 + 2nν (ω )], ω2 0 Z ∞ Jν (ω ) ′′ dω Qν (t) = sin(ω t). ω2 0

Q′ν (t) =

dω

In an equilibrium situation, TL = TR , hσx ieq = − ω∆0 tanh

β h¯ ω0 2

and hσz ieq = − tanh

out-of-equilibrium steady-state solution of the Markovian equation satisfies 8,25

hσz iss = [k(−ω0 ) − k(ω0 )]/[k(−ω0 ) + k(ω0 )],

(34)

β h¯ ω0 2

. 3 The

(35)

where the rates are convolutions of the L and R baths-induced rates,

k(ω0 ) =

Z ∞

eiω0t e−QL (t) e−QR (t) dt

−∞

1 = 2π

Z ∞

−∞

kL (ω0 − ω )kR (ω )d ω .

(36)

The “Fermi-Golden-Rule" rate constant kν (ω ) (×∆2 ) was originally derived in the context of reaction rates in donor-acceptor complexes, 2

kν (ω ) =

Z ∞

eiω t e−Qν (t) dt,

−∞


(37)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 46

and it satisfies the detailed balance relation, kν (−ω ) = kν (ω )e−β h¯ ω .

(38)

In the Ohmic case, in the scaling regime, kB T, h¯ ω < h¯ ωc , the ν -bath-induced rate obeys 1 1 kν (ω ) = ωc

h¯ ωc 2π kB T

1−2αν

|Γ(αν + i¯hω /2π kB T )|2 h¯ ω /2kB T , e Γ(2αν )

(39)

with the Gamma function Γ(x). Closed expressions for k(ω ), thus the polarization, are missing in general since the convolution (36) is nontrivial to handle analytically. As we show below, the NIBA heat current expression is built from the Fermi-Golden-Rule rates and the steady-state population of the TLS (in the local basis), pss 1 = (1 + hσz iss )/2, p0 = 1 − p1 , but it does not necessitate the coherences. In equilibrium, the spin population follows eq

p1 =

e−β h¯ ω0 /2 . e−β h¯ ω0 /2 + eβ h¯ ω0 /2

(40)

It is of interest to acquire pss 1,0 approximately-analytically, to study signatures of the non-equilibrium setting on the TLS at strong system-bath couplings.

4.2 Heat current A closed expression for the steady-state heat current under NIBA has been derived by energyunraveling the polarization dynamics, Eq. (32), 25 resulting in 2 Z ∞ h¯ ∆ ss ω d ω [kR (ω )kL (ω0 − ω )pss jq = 1 − kR (−ω )kL (−ω0 + ω )p0 ] . 2 2π −∞

(41)

Considering the unbiased model, the thermal conductance is given by the compact form 26 2 Z ∞ ∆ h¯ 2 κ= ω 2 kR (ω )kL (−ω )d ω . 2 4π kB T 2 −∞


(42)

Page 19 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


In the high T limit we substitute Eq. (39) into Eq. (42) and receive 26

2

κ ≃ A h¯ ≃ A

kB ∆ ωc

∆ ωc 2

Z 1 h¯ ωc 2−2αL −2αR kB T /¯h 2 ω dω kB T 2 kB T 0 h¯ ωc 1−2α , kB T

2

(43)

with A a constant which weakly depends on α through the Gamma function in Eq. (39). We emphasize that under NIBA the heat current, Eq. (41), is only determined by the polarization dynamics and its decay rates to steady-state. This observation explains the satisfactory performance of NIBA (with α ) in heat transfer calculations: 26 NIBA truncation of Ks/a,z carries errors (in the inter-blip correlations of the kernel) only second-order in α . As a result, NIBA dynamics of hσz (t)i is proper for the unbiased and weakly-damped NESB model. In contrast, the kernels Ks/a,x carry first-order errors in α , making them inaccurate even for the unbiased SB model. Nevertheless, fortunately these errors do not propagate into the NIBA heat current formula.

5 Path Integral Simulations: QUAPI and INFPI The numerically exact iterative quasi-adiabatic path-integral (QUAPI) approach has been developed by Makri and Makarov for simulating the reduced density matrix, spin dynamics in the spin-boson model. 49,50 It allows one to treat strong system-bath couplings and to include nonMarkovian effects. Recent efforts directed excitonic energy transfer in biomolecules, treating more complex situations, considering (a single excitation in) multiple sites and prominent vibrations in each site. 51 QUAPI algorithm had been constructed on the grounds of harmonic environments linearly coupled to the subsystem. 49,50 In this case, the effect of the environment on the subsystem’s RDM can be absorbed in an analytic function, the “Feynman-Vernon influence functional". 52 At nonzero temperatures the memory function within the influence functional decays rapidly in time, allowing for its controlled truncation and the development of an iterative time evolution scheme. 49,50



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

This principle has been generalized to construct a QUAPI-type algorithm for anharmonic baths. 53 Furthermore, a related-general approach, the so-called “influence functional path integral" (INFPI) tool has been put forward for treating quantum systems in contact with fermionic reservoirs at finite bias voltages. 54 INFPI can be excersized when an exact analytic form for the influence functional is missing; it is computed numerically using trace identities. Given their conceptual similarity, QUAPI and INFPI were recently combined and applied for the study of charge transport and vibrational excitation and dissipation in donor-acceptor molecular electronic diodes. 55 The standard QUAPI algorithm administers only the dynamics of the reduced density matrix. 49,50 The calculation of other observables, particularly the heat current in the NESB model, necessitates significant technical advances. In contrast, the INFPI algorithm has been formulated for treating a generic impurity Hamiltonian 54,55 and with little effort it can be employed for the study of other observables of quadratic structure, beyond the RDM. For example, INFPI has been applied for the investigation of charge current 54 and equilibration dynamics 56 in the single-impurity Anderson model, More recently, INFPI was adopted to simulate qubit-mediated energy flow between metals. 28 Algorithmic details of QUAPI 49,50 and INFPI 54,55 can be found elsewhere. Here we only highlight the techniques’ working principles and the associated numerical errors. The starting point of these methods involves a Trotter factorization of time evolution operators, e.g., into a free subsystem term and bath-related propagators. Collecting the environmental contributions and tracing over the baths, we reach an influence functional-type expression, see Eq. (44) below. As mentioned above, at finite temperatures and/or a nonzero chemical potential bias bath correlations exponentially decay in time, 57–59 allowing for their truncation beyond a memory time τc . An iterative time evolution scheme can then be constructed by defining an auxiliary quantity (an extension of the observable of interest) on the time-window τc . This time-nonlocal object can be evolved in an iterative manner from the initial condition to the final time t. Simulations with QUAPI and INFPI involve the following numerical errors: (i) Trotter error


Page 20 of 46

Page 21 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


due to the finite time-step adopted in the Trotter breakup δ t, and (ii) an error associated with the truncation of the influence functional to cover a finite time window τc . Furthermore, in INFPI the Fermi sea is discretized. Thus, one should confirm that within the relevant simulation time the number of bath states does not affect results.

5.1 Reduced density matrix The original QUAPI algorithm provided the subsystem’s RDM when coupled to a single harmonic bath. Here QUAPI is (trivially) extended to accommodate two reservoirs of different temperatures. The time evolution of the RDM up to t = N δ t, N is an integer, can be represented in a path integral formulation as, 49,50 − hs+ N |ρS (t)|sN i =

Z

ds± 0

Z

ds± 1 ...

Z

har ± ± ds± (s0 , s1 , ..., s± N ). N−1 I

(44)

Here s± k represent the discrete path of the subsystem on the forward (+) and backward (−) contours (not to be confused with the eigenstates of HS introduced in Sec. 3.1). As an initial condition we assume a product form, ρ (0) = ρL ⊗ ρR ⊗ ρS (0) with the baths separated from the subsystem. The integrand in Eq. (44) is refereed to as an “Influence Functional" (IF); 52 note that in Refs. 49,50,52 it was identified without the free subsystem evolution terms. For a harmonic bath bilinearly coupled to the subsystem the IF is given by an exponential of a quadratic structure, multiplied by free subsystem propagation terms,

I

har

± (s± 0 , ..., sN )

= ×

h

N

k

i

− ν + ν∗ − exp − (s+ k − sk )(ηk,k′ sk′ − ηk,k′ sk′ ) ′ ν k k =0 iH0 δ t − − − + −iH0 δ t + |sN−1 i...hs+ |sN i. hsN |e 0 |ρS (0)|s0 i...hsN−1 |e

∑∑ ∑

(45)

ν depend on the spectral function The free Hamiltonian H0 is defined in Eq. (1), the coefficients ηk,k ′

of the ν bath and its temperature. They were derived in Ref. 49,50 by discretizing the FeynmanVernon IF.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 46

5.2 Heat Current The heat current in the NESB model can be written as a Meir-Wingreen formula combining spinspin correlation functions. 21 Here, rather than generalizing QUAPI to follow correlation functions in a non-equilibrium setup, a nontrivial task, we study the behavior of the NESB model through the related fermionic model, more naturally handled by INFPI. The dissipative dynamics of the spin-boson model with an Ohmic dissipation can be related to the behavior of the spin-fermion (SF) model, considering a fermionic bath with a linear dispersion, as was demonstrated in Refs. 60–62 An important difference that arises in this comparison is that the dimensionless parameter α of the bosonic model takes a different range of values in the fermionic-bath case, as we discuss below. The spin-fermion model of Ref. 60 includes spin-up and spin-down electrons which do not mix, i.e., a spin-flip term is missing in the Hamiltonian. Here, we (conceptually) replace the electron spin degree of freedom by metal leads; electrons are residing in separate leads ν = L, R, and electron hopping between metals is not allowed. With this correspondence, we can directly apply the results of Ref. 60 to our situation, only noting that the two species (spin-up and spin-down or ν = L, R electrons) are prepared in thermal states of different temperatures. We now introduce the SF model, a close variant of the model described in Ref. 60 It comprises the central spin [Eq. (1)] and two metallic leads ν = L, R, 1D electron gases with linear dispersion. The metals are prepared at different temperatures but at the same chemical potential. They are connected indirectly, only through the TLS. Charge transfer processes between the metals are blocked but we allow for excitation energy flow. The total non-equilibrium spin-fermion (NESF) Hamiltonian reads

HSF = H0 + h¯ σz

∑ ′ g p,ν ;p′,ν c†p,ν c p′,ν + ∑ ε p,ν c†p,ν c p,ν .

(46)

ν ,p

ν ,p,p

Here c†p,ν (c p,ν ) represents a fermionic creation (annihilation) operator. The spin polarization couples to intra-bath electron-hole pair generation as in Ref. 28 22 ACS Paragon Plus Environment

Page 23 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


The spin dynamics in the single-bath spin-boson and the spin-fermion model can be written in an equivalent influence functional path integral form. 60–62 Considering here non-equilibrium situations, the NESB and the NESF models, we relate the coupling parameter g p,ν ;p′ ,ν , taken as a constant gν in simulations, to the boson picture via the relation 3,60 2 1 2 αν = atan(πρν (εF )gν ) . 2 π

(47)

Here ρν (εF ) is the density of states at the Fermi energy and αν is a dimensionless dissipation parameter, the prefactor in the bosonic-Ohmic spectral function, see Eq. (4). It is important to point out that the fermionic dissipation factor is bounded, π2 atan[πρν (εF )gν ] ≤ 1, thus we cannot observe the (zero temperature) α = 1 localization effect taking place in the spin-boson model. 3 Furthermore, it becomes exceedingly difficult to simulate strong dissipation scenarios within the SF picture since αν unfavorably translates to large interaction parameters in the fermionic model, forcing us to adopt very small time steps in INFPI simulations so as to maintain a sufficiently small Trotter error. Four comments are now in place: (i) The relation between the single-bath SB and SF models 60 could be immediately extended to associate spin dynamics in the NESB and the NESF models since electron transitions between metals are not allowed, resulting in a multiplicative form for the influence functional. Incorporating electron scatterings between the two metals under a finite voltage-bias would result in a more complicated (non-separable L − R) form for the influence functional. 63 (ii) The relation between the SB and the SF models is valid in the Ohmic case assuming (in the fermionic picture) a constant density of states and a wide conduction band, h¯ ωc ≫ kB T . (iii) The SB and SF models were associated in Refs. 60–62 based on the correspondence of their partition function, or reduced density matrix. The thermal conductance, expressed in terms of an equilibrium spin-spin correlation function, thus exactly relate in these two models. Out-of-equilibrium, it remains a challenge to derive formal influence functional expressions for the expectation values of the energy current operators for the two models, to prove their equivalence. Nonetheless, in Ref. 64



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 46

we demonstrated that in the weak coupling limit spin relaxation (excitation) rates due to electronhole pair generation (annihilation) in a metal of a linear dispersion exactly translate to boson bath induced rates. Thus, in the weak coupling limit the thermal transport properties of the NESF and the NESB models exactly correspond, even far from equilibrium. (iv) A different-common route for exploring the properties of the SB model involves its mapping to either the Kondo Hamiltonian with an anisotropic exchange coupling 65 or the Ising model with long-range interactions. 65 These connections were employed in many studies for exploring the equilibrium properties of the SB model, 3 particularly its thermal conductance. 29 We had recently simulated the energy current characteristics of the unbiased NESF model by adapting the INFPI approach. 28 We discuss next the structure of the current operator in this simulation, for ω0 = 0. First, we transform the Hamiltonian (46) via a unitary transformation H˜ SF = U † HSF U as in Sec. 3.1 into H˜ SF = HS + HF +V,

(48)

where

HS =

h¯ ε σ˜ z , 2

HF = HL + HR , Hν = ∑ ε p,ν c†p,ν c p,ν , p

V = VL +VR , Vν = σ˜ x ∑ h¯ g p,ν ;p′ ν c†p,ν c p′ ,ν

(49)

p,p′

As before, σ˜ x,y,z stand for Pauli matrices in the energy basis. We assume a factorized initial state, ρ˜ (0) = ρ˜ S (0) ⊗ ρL ⊗ ρR , ρ˜ S denotes the reduced density matrix of the subsystem and ρν = e−βν (Hν −µν Nν ) /Trν [e−βν (Hν −µν Nν ) ]. In our simulations below we take µL = µR but assume different (time-zero) temperatures for the thermal baths, TL 6= TR . At t = 0 we put into contact the two Fermi baths through the quantum subsystem, then follow the evolution of the reduced density matrix and the energy current until (quasi) steady-state sets in. Since electron flow is blocked and energy is


Page 25 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


transferred through excitation - de-excitation processes of the TLS, we refer to the energy current here as a “heat current". The heat current in the unbiased (ω0 = 0) model was calculated in the energy basis using the construction 16 i jL = TrS TrF [ jˆL ρ˜ (t)] = − TrS TrF {ρ˜ (t)[HS ,VL ]}, h¯

(50) ˜

˜

obtained from the definition (7) under a steady-state assumption. Here, ρ˜ (t) = e−iHSF t/¯h ρ˜ (0)eiHSF t/¯h , TrF (TrS ) refers to a partial trace over the Fermi-sea electrons (TLS). We identify the current with the contact at which it is evaluated, though jq = jL = − jR is satisfied here. The commutator can be readily performed to yield [HS ,VL ] = i¯h2 ∆σ˜ y ∑ gl,L;l ′ ,L c†l,L cl ′ ,L ,

(51)

jL = (¯h∆)TrS [σ˜ y TrF [AL ρ˜ (t)]] .

(52)

l,l ′

leading to

The bath operator is given by AL ≡ ∑l,l ′ gl,L;l ′ ,L c†l,L cl ′ ,L . We now further define a subsystem operator as ˜

˜

AS (t) ≡ TrF [AL ρ˜ (t)] = TrF [eiHSF t/¯h AL e−iHSF t/¯h ρ˜ (0)],

(53)

and express the current from its matrix elements

jL = (¯h∆)[−i(AS (t))−,+ + i(AS (t))+,− ].

(54)

We use INFPI to time-evolve Eq. (53). 28 These Simulations are compared to BR and NIBA results from Secs. 3.3 and 4.2, respectively. 25 ACS Paragon Plus Environment


α L = α R = 0.025

hσz (t)i

1 0.5

NIBA BR QUAPI

0

(a)

hσy (t)i

hσx (t)i

-0.5 0

(b)

-0.1 -0.2

0.2 0 -0.2 -0.4 -0.6

(c)

0

10

t∆

20

30

Figure 2: Dynamics of the spin subsystem in the NESB model assuming an Ohmic spectral function with ωc = 20∆, and kB TL = kB TR = 2¯h∆, ω0 = 0, αL = αR = 0.025. (a-c) Results for hσx,y,z (t)i are presented in the local basis (|0i and |1i), QUAPI (dotted), NIBA (dashed), BR (dashed-dotted).

hσz (t)i

1

α L = α R = 0.1 (a)

0.5

NIBA BR QUAPI

0

hσx (t)i

0 -0.1

hσy (t)i

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 46

0 -0.1 -0.2 -0.3

(b)

-0.2

(c)

0

10

t∆

20

30

Figure 3: Dynamics of the spin subsystem with the parameters of Fig. 2, αL = αR = 0.1.


Page 27 of 46

hσz (t)i

1

α L = α R = 0.15 (a)

0.5

NIBA BR QUAPI

0

hσx (t)i

0

(b)

-0.1 -0.2 0.1

hσy (t)i

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(c)

0 -0.1 -0.2 0

10

t∆

20

30

Figure 4: Dynamics of the spin subsystem as in Fig. 2, with αL = αR = 0.15.

6 Simulations The dynamics of the reduced density matrix is obtained by solving numerically the BR and NIBA integro-differential equations, Eq. (15) and (32), respectively, using the trapezoidal rule for the inner integral, see Ref. 66 This brute-force approach is adopted with a small time step, δ t = Tp /2000; the period is defined by Tp ≡ 2π /ε . Integrals over frequency are evaluated using the trapezoidal rule, discretized with δ ω ∼ ε /2000, up to the limit 30ωc . QUAPI and INFPI simulations are performed with δ t ∼ Tp /50. Other parameters are kB T /¯h∆ ∼ 1 − 2, ω0 /∆ = 0 − 5, ωc /∆ = 10 − 20 and αν = 10−3 − 0.5. The thermal conductance was calculated by taking a small temperature difference, (TL − TR )/Ta = 0.05, with the averaged temperature Ta ≡ (TL + TR )/2. The averaged temperature is sometimes denoted by T . INFPI simulations better converge for larger temperature differences, kB (TL − TR )/¯h∆ ∼ 1. We use here an Ohmic spectral function, but we are not limited to this form, as other environments can be similarly explored, e.g., a Debye spectral function or a spin bath, mimicked by a harmonic environment with a temperature-dependent spectral function. 67



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 46

6.1 Unbiased model RDM Dynamics. We begin with the unbiased ω0 = 0 model. In Figs. 2-4 we follow the dynamics of the subsystem in the local basis (|0i, |1i) using BR, NIBA and QUAPI, 49,50 extended here to the non-equilibrium two-bath case. At very weak coupling, αν = 0.025, BR and NIBA reasonably agree with QUAPI; at smaller couplings, α < 0.02, the agreement is excellent (not shown). Increasing the coupling to αν = 0.1 − 0.15, translating to a total interaction strength α = 0.2 − 0.3, we note that NIBA and BR reasonably agree with QUAPI over the polarization dynamics, but the real part of off-diagonal elements depart by up to a factor of 1.5 in the steady-state limit. Heat transfer: Coupling to the contacts. Deviations between BR and NIBA in the RDM behavior propagate into significant, qualitative differences in the heat current characteristics. In Fig. 5 we display the thermal conductance of the NESB model and show that within the BR scheme it grows linearly with α , 8 while NIBA demonstrates saturation of the current, then its suppression at large values of α . 26 We also illustrate the behavior of the thermal conductance using the NEGF-Redfield scheme of Ref. 21 In this approach correlation functions in the (exact) heat transfer Meir-Wingreen formula are evaluated at the level of the Redfield theory, which is second order in the system-bath coupling. Another NEGF-based expression for the heat current has been recently proposed in Ref. 24 It was obtained from the diagrammatic expansion, by calculating spin-spin correlation functions via the Majorana-fermion representation of spin operators, truncating diagrams of high order system-bath couplings. This approach results in an elegant expression (ω0 = 0), 2 jq = π

Z ∞ 0

h¯ ω ∆2 ΓR (ω )ΓL (ω )[nL (ω ) − nR (ω )] dω . (ω 2 − ∆2 )2 + [ΓL (ω ) coth(βL h¯ ω /2) + ΓR (ω ) coth(βR h¯ ω /2)]2 ω 2

(55)

This formula resembles the harmonic limit 68 and it nicely interpolates between the very weak coupling equation (29), with a kinetic-type dynamics, and the low-temperature limit, where heat is transferred coherently, not necessarily in-resonance with the central frequency; the integration over frequencies away from ∆ reflects a tunneling, off-resonance, behavior. We find in Fig. 5 that the two NEGF-based approaches perform very well up to αν = 0.1, but


Page 29 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


their predictions beyond that are qualitatively incorrect: These methods fail to provide the decay behavior of the current with α at strong coupling. Nonetheless, even in this regime NEGF-Redfield and the NEGF expression (55) are valuable; their predictions are significantly closer to NIBA than to BR. In Fig. 6 we use INFPI to simulate the heat current in the NESF model. 28 We display the current as a function of the interaction parameter α up to α /2 = 0.1, which we consider as a weakintermediate value. We compare INFPI results to the predictions of BR, NIBA, NEGF-Redfield 21 and NEGF of Ref. 24 We find that BR equation fails beyond α /2 = 0.02, but the other methods excellently agree in this range. Note that in the fermionic language α /2 = 0.1 translates to a phase shift of φν = πρν αν ∼ 0.85. We now explain our classification of system-bath coupling domains based on Figs. 5-6: In the very weak coupling limit αν < 0.025 the Bloch-Redfield treatment is valid within up to 10% deviations from exact results. In the so-called weak coupling regime 0.025 < αν < 0.1 deviations from linearity are apparent, but the current is growing monotonically with α . For stronger couplings, 0.1 < αν < 0.5, the heat current displays a crossover behavior: Fig. 5 shows that beyond αν = 0.15 the thermal conductance drops with increasing couplings to the bath. We refer to this crossover area as the intermediate coupling regime. At even stronger coupling the behavior of jq (α ), or the thermal conductance, is dominated by an exponentially decaying factor of α , see Eq. (43). INFPI results are reported in Fig. 6 only up to α /2 = 0.1 given the computational hardness involved: Simulations are performed from a certain initial condition to (quasi) steady-state, and we need to confirm convergence with respect to the number of electronic states in the fermionic baths, simulation time step and the bath memory time. Increasing the interaction energy requires us to employ more states in the two baths, a finer simulation time step, and a longer bath-decorrelation time, making INFPI computations very costly. Particularly, the nonlinear mapping Eq. (47) unfavorably translates α > 0.3 values to large interaction energies in the fermionic picture. It is useful to comment at this point that a fully harmonic junction is expected to provide larger currents than supported by the (anharmonic) NESB junction. 68 Specifically, under a weak cou-



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

pling approximation the current in harmonic junctions obeys Eq. (29), only missing temperaturedependent distribution functions in the denominator. Beyond the non-adiabatic limit. The non-adiabatic limit describes a junction comprising a lowfrequency (slow) central vibration relative to the cutoff frequency of the bath. In Fig. 7 we abandon the strict non-adiabatic region and explore the dynamics when ∆/ωc = 0.25 with 2.5kB Tν = h¯ ∆. We study the time evolution of the RDM at very weak coupling and find that NIBA fails to reproduce hσx (t)i, deviating by ∼ 40%, while it performs well for the polarization dynamics. We now discuss implications on transport properties. The effect of the spin tunneling frequency ∆ on the steady-state heat current is displayed in Fig. 8, considering a weakly coupled nanojunction. The heat current exhibits a turnover behavior as a function of ∆, explained based on the BR expression (29): The current increases with ∆, the quanta transmitted, but it further requires a sufficiently high thermal occupation factor in matching bath modes. Thus, when the tunneling frequency is high, ∆ > 2kB T /¯h, the current begins to drop due to reduced population of in resonance modes. Comparing BR (dashed) to INFPI (square), we find that the BR scheme reproduces the correct turnover behavior, though the exact position of the maxima is shifted and the magnitude of the current is overestimated at high temperatures. 70 NEGF results are not displayed here; in this weak-coupling example they overlap with BR. In contrast, NIBA equations miss altogether the correct behavior of jq (∆) once ∆/ωc > 0.1 since NIBA only captures non-adiabatic contributions, jq ∝ ∆2 . Moreover, it is important to note that not only does NIBA provide the wrong functional form for jq (∆), it further predicts an erroneous temperaturedependent behavior: While in the non-adiabatic regime at small coupling the current drops here with increasing temperatures, jq ∝ T −1 , a tendency captured by NIBA, the opposite behavior takes place beyond the non-adiabatic regime once the current is controlled by thermal occupation factors in the contacts. Another interesting observation concerns the low temperature kB Ta /¯hωc = 0.07 regime, see Fig. 8(c). In this case INFPI simulations exhibit two peaks in the current jq (∆), around ∆/ωc =0.25, 0.65. The first peak corresponds to absorption and emission processes of a single phonon in either


Page 30 of 46

Page 31 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


the L or R baths, as explained above. In the language of the SF model, these are single electron-hole pair generation or recombination processes. Around ∆/ωc =0.25, the enhancement of the current with an increasing spin frequency is balanced by thermal occupation factors of bath modes. We presume that the second peak corresponds to a similar balance, apparently reflecting two-phonon processes (or two electron-hole pairs) participating in the excitation and relaxation dynamics of the TLS. Additional simulations are required to establish this result. Heat transfer: Temperature dependence. The temperature dependence of the thermal conductance is of particular interest for actual devices; exact Monte-Carlo simulations 29 provided the 2α −1 high temperature behavior of the NESB model κ ∝ kh¯BωTc and the low-temperature scaling

κ ∝ α (T /TK )3 , TK is the Kondo temperature in the system, a function of the microscopic parameters ∆, ωc and α . We had recently proved that the high-temperature scaling is reproduced correctly by the NIBA-heat current formula, see Eq. (43). 26

We examine the temperature dependence of the thermal conductance in Figs. 9 and 10, considering a very weak, weak, and intermediate spin-bath coupling. In Fig. 9 we demonstrate that the BR technique reproduces the correct temperature characteristics at the very weak coupling limit. At low temperatures, an increase of the temperature allows more resonant modes in the bath to be thermally occupied, resulting in an enhancement of the thermal conductance. At high temperatures kB Ta /¯h∆ & 1 the thermal conductance decays with temperature, responding to anharmonicity in the junction. NIBA evidently fails at low temperatures: It misses the low T enhancement of the thermal conductance with increasing T . In Fig. 10 we focus on the weak-to-intermediate coupling regime. We do not display simulations with INFPI here given computational cost, but results from BR and NEGF indicate that a turnover behavior with temperature is again expected to show, though the BR scheme progressively overestimates the conductance as α grows. While NIBA lacks the turnover behavior with temperature, it yields the correct scaling of κ with T in the non-adiabatic limit, see Eq. (43). This scaling is neither reproduced by BR nor by NEGF, see Fig. 10. Complementing this discussion, we note that the NEGF expression (55) reduces to the harmonic result at low temperatures, build-



ing up the scaling κ ∼ T 3 . It is thus crucial to carefully examine parameters of interest to decide which technique is most suitable since approximate approaches, NIBA, BR and NEGF, break down (catastrophically) beyond their strict regime of validity.

0.08

BR NIBA NEGF-Redfield NEGF

κ/(kB ∆)

0.1

κ/(kB ∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 46

0.02 0.01 0 0

0.06

α/2

0.05

0.04 0.02 0 0

0.1

0.2

α/2

0.3

0.4

0.5

Figure 5: Thermal conductance in the unbiased NESB model as a function of α /2 = αL = αR . We use Ohmic spectral functions with ωc = 20∆, ω0 = 0, kB TL = 2.05¯h∆, kB TR = 1.95¯h∆ with BR (dashed), NIBA (▽), NEGF-Redfield from Ref. 21 (+) and the NEGF expression (55) (◦). The inset zooms over the weak coupling regime.

6.2 Biased model RDM dynamics. In Fig. 11 we display the RDM dynamics in the biased model. At weak coupling NIBA misses the correct behavior of the RDM 1 while at strong coupling it reasonably agrees with QUAPI. Thermal conductance. Since the NIBA-heat current expression only depends on the polarization dynamics (steady-state value and the relaxation rates), errors in hσx (t)i do not propagate into the calculation of the heat current, see Fig. 12. The bias affects the conductance in a simple way: In the BR scheme κ ∝ ε / sinh(β h¯ ε ), see Eq. (31). In NIBA calculations we confirmed numerically that κ ∝ ω0 / sinh(β h¯ ω0 ) for α = 0.1 − 0.5. The spin bias thus does not offer a new nontrivial quantum control knob over the heat current, as ω0 does not tangle with the coupling strength α in the non-adiabatic limit. While these results are not encouraging, it is of fundamental interest to explore the behavior of the NESB model with finite detuning using other techniques beyond the BR and NIBA schemes. The NEGF approach of Ref. 21 was developed for the unbiased model, and 32 ACS Paragon Plus Environment

Page 33 of 46

0.04

(a)

jq /(¯h∆2 )

0.08

jq /(¯h∆2 )

0.06

(b)

0.03 0.02 0.01

0.04

0 0

0.02 0 0

0.1

0.2 α/2

0.3

0.05 α/2

0.1

INFPI BR NIBA NEGF-Redfield NEGF

0.4

Figure 6: (a) Heat current in the unbiased NESF model. Simulations were performed in the fermionic picture using electron bands with linear dispersion and a hard cutoff at D/∆ = ±5. The parameters of the fermionic model are matched with the bosonic picture using Eq. (47), αL = αR , ω0 = 0, kB TL = 2¯h∆, kB TR = h¯ ∆. We compare INFPI () to simulations in the bosonic picture, BR (dashed), NIBA (▽), NEGF-Redfield 21 (+) and the NEGF expression (55) 24 (◦). Panel (b) zooms over the small-α regime. INFPI numerical parameters are δ t = 0.1/∆ and Ns = 9, for more details see Ref. 28

α L = α R = 0.025

hσz (t)i

1

(a)

0 NIBA BRE QUAPI

hσx (t)i

-1 0 -0.5 -1

(b)

1 (c)

hσy (t)i

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


0

-1

0

10

t∆

20

30

Figure 7: Dynamics of the unbiased model beyond the non-adiabatic region for ωc = 4∆, kB TL = kB TR = h¯ ∆/2.5, and αL = αR = 0.025. QUAPI (dotted), NIBA (dashed), BR (dashed-dotted).



kB Ta /¯ h ωc = 0.145

0.12

0.095

x 10 (a)

6 Full lines: NIBA Dashed lines: BR Square: INFPI

7

jq /(¯ h ωc2 )

8

0.07 -6

-5

6 jq /(¯ h ωc2 )

5 4

x 10 (b)

4

2

0 0

3

∆/ωc 0.05 -5

3 jq /(¯ h ωc2 )

2 1 0 0

0.2

0.4 ∆/ωc

0.6

x 10

1 0 0

0.8

(c)

2

0.5 ∆/ωc

Figure 8: Heat current as a function of the TLS tunneling frequency. (a) INFPI simulations in the fermionic picture (). The parameters of the fermionic model are matched with the bosonic case using Eq. (47), πρ (εF )g = 0.3, leading to αL = αR = 0.0172. Other parameters are ω0 = 0, kB Ta /¯hωc as indicated in the legend with Ta = (TL + TR )/2 and (TL − TR )/ωc = 0.01. We compare INFPI () results to simulations in the bosonic picture, BR (dashed), NIBA (full). Panel (b) zooms over the non-adiabatic regime. INFPI numerical parameters are δ t = 0.1/∆ and Ns = 9, for more details see Ref. 28 Panel (c) zooms over the low temperature case, kB Ta /¯hωc = 0.07.

0.04 INFPI

BR

NIBA

αν = 0.0172, ωc = 10∆

0.03 κ/(kB ∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 46

0.02

0.01

0 0

0.5

1 kB Ta /¯ h∆

1.5

2

Figure 9: Temperature dependence of the thermal conductance in the unbiased NESB model using a very weak coupling value, αν = 0.0172. BR (dashed), NIBA (⊳), and INFPI ().


Page 35 of 46

BR

NEGF

NIBA αν = 0.15

-1

Ta2α−1

10 κ/(kB ∆)

Ta−1

-2

0

10

10 κ/(kB ∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Ta3 -3

10

αν = 0.05

-5

10

e−¯h∆/k B T a -1

-1

10

0

10 kB Ta /¯ h∆ 0

10

10 kB Ta /¯h∆

Figure 10: Thermal conductance in the unbiased NESB model as a function of temperature in the weak-to-intermediate regime. We use Ohmic spectral functions with ωc = 20∆, αν = 0.15, BR (dashed), NIBA (⊳), and the NEGF expression (55) (◦). NEGF-Redfield formula did not provide physical answers at high temperatures here. Full lines indicate on the scaling relations; NIBA results should retain the asymptotic Ta2α −1 behavior in the scaling limit ∆/ωc ≪ 1 when Ta /ωc ≪ 1. The inset displays a weak coupling example with αν = 0.05. similarly Eq. (55) from Ref. 24 assumes ω0 = 0. The Monte-Carlo study of Ref. 29 again enforced zero bias, as well as early INFPI simulations. 28 Our future work will address this model in more details.

7 Summary We provided a comprehensive analysis of the non-equilibrium spin-boson model, comprising a spin subsystem coupled to two thermal baths of different temperatures. We studied the dynamics of the spin reduced density matrix and the transfer of heat in the model using different techniques: the Bloch-Redfield scheme which is valid in the very weak system-bath coupling limit, the noninteracting-blip approximation, correct in the non-adiabatic (∆ ≪ ωc ) scaling limit, and numerically-exact influence functional path integral simulations. We also compared results to NEGF-based techniques. 21,24 The BR, NIBA and path integral approaches were originally developed for exploring decoherence effects and dissipation dynamics in open quantum systems. Here 35 ACS Paragon Plus Environment


α L = α R = 0.2, ω0 = ∆ hσz (t)i

1

(a)

0.5

NIBA BR QUAPI

0

hσx (t)i

-0.5 0

(b)

-0.1 -0.2

hσy (t)i

0.1

(c)

0 -0.1 -0.2 0

10

t∆

20

30

Figure 11: Dynamics of the biased model ω0 = ∆ in the intermediate coupling regime, αL = αR = 0.2, and ωc = 20∆, kB TL = kB TR = 2¯h∆. QUAPI (dotted), NIBA (dashed), BR (dashed-dotted).

0.03

Top to bottom: ω0 /∆= 0, 1, 2, 5

0.025

BR: Full lines NIBA: dashed-dotted lines

0.02 κ/(kB ∆)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 46

0.015 0.01 0.005 0 0

0.2

0.4

0.6

0.8

1

α/2

Figure 12: Thermal conductance for the biased model assuming Ohmic spectral functions and αL = αR , ωc = 20∆, kB TL = 2.05¯h∆, kB TR = 1.95¯h∆.


Page 37 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


we bring an organized discussion of their extensions to treat transport behavior. Specific observations include: (i) The biased and unbiased NESB models may display similar dynamics within the BR and NIBA approaches, but relatively small deviations in the RDM behavior propagate into strong and qualitative disagreements in the heat current characteristics. (ii) The BR formalism should be used with caution in modeling actual devices since the prediction jq ∝ α fails beyond the very weak coupling limit, providing unphysical-incorrect large conductances. (iii) In the regime of validity for BR and NIBA, the spin bias (detuning) parameter does not offer a new-nontrivial control mean over the heat current; at large detuning the current decays since thermal occupation of high frequency bath modes (above the thermal energy) is reduced. (iv) In the non-adiabatic regime at weak coupling and high temperatures the thermal conductance decreases with increasing temperatures, κ ∝ (TL − TR )/T . This behavior stems from the intrinsic anharmonicity of the junction. This trend is correctly captured by the BR method, NIBA and INFPI. In contrast, beyond the non-adiabatic limit this trend is reversed since thermal occupation factors of bath modes dominate the current, rather than temperature-dependent (anharmonic) scatterings in the junction. As expected, NIBA fails to capture this behavior: It overestimates the current by orders of magnitude and it predicts an enhancement of the current when reducing the temperature, irrespective of the frequency ∆. It is of interest to extend our analysis and further explore the behavior under classical equations of motion, or mixed quantum-classical treatments. 71 Future studies will also examine the timedependent-driven NESB model 72–75 with the objective to understand the role of strong system-bath couplings and quantum coherence in possibly enhancing heat to work conversion efficiency.

Acknowledgement This work was funded by the Natural Sciences and Engineering Research Council of Canada and the Canada Research Chair Program. NB acknowledges support from the CQIQC Undergraduate Summer Research Studentships.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1) Weiss, U. Quantum Dissipative Systems; World Scientific, Singapore, 1999. (2) Nitzan, A. Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems; Oxford, 2006. (3) Leggett, A. J.; Chakravarty, S.; Dorsey, A. T.; Fisher, M. P. A.; Garg, A.; Zwerger, W. Dynamics of the Dissipative Two-State System. Rev. Mod. Phys. 1987, 59, 1-85. (4) Le Hur, K. Understanding Quantum Phase Transitions. Quantum Phase Transitions in SpinBoson Systems: Dissipation and Light Phenomena; ed. Carr, L. D., Taylor and Francis, Boca Raton, 2010. (5) Dhar, A. Heat Transport in Low-Dimensional Systems. Adv. Phys. 2008, 57, 457-537. (6) Wang, J.-S.; Wang, J.; Lu, J. T. Quantum Thermal Transport in Nanostructures. Eur. Phys. J. B 2008, 62, 381-404. (7) Li, N.; Ren, J.; Wang, L.; Zhang, G.; Hänggi, P.; Li, B. Colloquium: Phononics: Manipulating Heat Flow with Electronic Analogs and Beyond. Rev. Mod. Phys. 2012, 84, 1045-1066. (8) Segal, D; Nitzan, A. Spin-Boson Thermal Rectifier. Phys. Rev. Lett. 2005, 94, 034301. (9) Segal, D; Nitzan, A. Heat Rectification in Molecular Junctions. J. Chem. Phys. 2005, 122, 194704. (10) Wang, H.; Thoss, M.; Miller, W. H. Self Consistent Hybrid Approach for Complex Systems; Applications to the Spin-Boson Model with Debye Spectral Density. J. Chem. Phys. 2001, 115, 2991-3005. (11) Wilhelm, F. K.; Kleff, S.; Von Delft, J. The Spin-Boson Model With a Structured Environment: a Comparison of Approaches. Chem. Phys. 2004, 296, 345-353.


Page 38 of 46

Page 39 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(12) Nesi, F.; Paladino, E.; Thorwart, M.; Grifoni, M. Spin-Boson Dynamics: A Unified Approach from Weak to Strong Coupling. Europhys. Lett. 2007, 80, 40005. (13) Nesi, F.; Paladino, , E.; Thorwart, M.; Grifoni, M. Spin-Boson Dynamics Beyond Conventional Perturbation Theories. Phys. Rev. B 2007 76, 155323. (14) Berkelbach, T. C.; Reichman, D. R.; Markland, T. E. Reduced Density Matrix Hybrid Approach: An Efficient and Accurate Method for Adiabatic and non-Adiabatic Quantum Dynamics. J. Chem. Phys. 2012 136, 034113. (15) Segal, D. Heat Flow in Nonlinear Molecular Junctions: Master Equation Analysis. Phys. Rev. B 2006, 73, 205415. (16) Wu, L.-A.; Yu, C. X.; Segal, D. Nonlinear Quantum Heat Transfer in Hybrid Structures: Sufficient Conditions for Thermal Rectification. Phys. Rev. E 2009, 80, 041103. (17) Ruokola, T.; Ojanen, T. Thermal Conductance in a Spin-Boson Model: Cotunneling and Low-Temperature Properties. Phys. Rev. B 2011, 83, 045417. (18) Thingna, J.; Garcia-Palacios, J. L.; Wang, J.-S. Steady-State Thermal Transport in Anharmonic Systems: Application to Molecular Junctions. Phys. Rev. B 2012 85, 195452. (19) Thingna, J. ; Wang, J.-S.; Hänggi, P. Reduced Density Matrix for Non-Equilibrium Steady States: A Modified Redfield Solution Approach. Phys. Rev. E 2013, 88, 052127. (20) Thingna, J.; Zhou, H.; Wang, J.-S. Improved Dyson Series Expansion for Steady-State Quantum Transport Beyond the Weak Coupling Limit - Divergences and Resolution. arXiv:1408.6996. (21) Velizhanin, K. A.; Thoss, M.; Wang, H. Meir-Wingreen Formula for Heat Transport in a Spin-Boson Nanojunction Model. J. Chem. Phys. 2010, 133, 084503.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

˘ Zs ´ Function (22) Wang, J.-S.; Agarwalla, B. K.; Li, H.; Thingna, J. Nonequilibrium GreenâA Method for Quantum Thermal Transport. Front. Phys. 2013. DOI 10.1007/s11467-013-0340x. (23) Vinkler-Aviv, Y.; Schiller, A.; Andrei, N. Single-Molecule-Mediated Heat Current Between an Electronic and a Bosonic Bath. Phys. Rev. B 2014, 89, 024307. (24) Yang, Y.; Wu, C.-Q. Quantum Heat Transport in a Spin-Boson Nanojunction: Coherent and Incoherent Mechanisms. Europhys. Lett. 2014, 107, 30003. (25) Nicolin, L.; Segal, D. Non-equilibrium Spin-Boson Model: Counting Statistics and the Heat Exchange Fluctuation Theorem. J. Chem. Phys. 2011, 135, 164106. (26) Segal, D. Heat transfer in the Spin-Boson Model: A Comparative Study in the Incoherent Tunneling Regime. Phys. Rev. E 2014 90, 012148. (27) Velizhanin, K. A.; Wang, H.; Thoss, M. Heat Transport Through Model Molecular Junctions: A Multilayer Multiconfiguration Time-Dependent Hartree Approach. Chem. Phys. Lett. 2008, 460, 325-330. (28) Segal, D. Qubit-Mediated Energy Transfer Between Thermal Reservoirs: Beyond the Markovian Master Equation. Phys. Rev. B 2013, 87, 195436. (29) Saito, K.; Kato, T. Kondo Signature in Heat Transfer via a Local Two State System. Phys. Rev. Lett. 2013, 111, 214301. (30) Wang, Z.; Carter, J. A.; Lagutchev, A.; Koh, Y. K.; Seong, N.-H.; Cahill, D. G.; Dlott, D. D. Ultrafast Flash Thermal Conductance of Molecular Chains. Science 2007, 317, 787-790. (31) Wang, Z.; Cahill, D. G.; Carter, J. A.; Koh, Y. K.; Lagutchev, A.; Seong, N.-H.; Dlott, D. D. Ultrafast Dynamics of Heat Flow across Molecules. Chem. Phys. 2008, 350, 31-44.


Page 40 of 46

Page 41 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(32) Schwarzer. D.; Kutne, P.; Schroder, J.; Troe, J. Intramolecular Vibrational Energy Redistribution in Bridged azulene-anthracene Compounds: Ballistic Energy Transport Through Molecular Chains. J. Chem. Phys. 2004, 121, 1754-1764. (33) Schwarzer. D.; Hanisch, C; Kutne, P.; Troe, J. Vibrational Energy Transfer in Highly Ex˘ Direct Observation of Energy Flow through cited Bridged Azulene-Aryl Compounds:âAL’ Aliphatic Chains and into the Solvent. J. Phys. Chem. A 2002, 106, 8019-8028. (34) Botan, V.; Backus, E. H. G.; Pfister, R.; Moretto, A.; Crisma, M; Toniolo, C; Nguyen, P. H.; Stock, G.; Hamm, P. Energy Transport in Peptide Helices. Proc. Natl. Acad. Sci. U.S.A. 2007, 104, 12749-12754. (35) Kasyanenko, V. M.; Tesar, S. L.; Rubtsov, G. I.; Burin, A. L.; Rubtsov, I. V. Structure Dependent Energy Transport: Relaxation-Assisted 2DIR Measurements and Theoretical Studies. J. Phys. Chem. B 2011, 115, 11063-11073. (36) Pein, B. C.; Sun, Y.; Dlott, D. D. Unidirectional Vibrational Energy Flow in Nitrobenzene. J. Phys. Chem. A 2013, 117, 6066-6072. (37) Pein, B. C.; Sun, Y.; Dlott, D. D. Controlling Vibrational Energy Flow in Liquid Alkylbenzenes. J. Phys. Chem. B 2013, 117, 10898-10904. (38) Pein, B. C.; Dlott, D. D. Modifying Vibrational Energy Flow in Aromatic Molecules: Effects of Ortho Substitution. J. Phys. Chem. A, 2014, 118, 965-973. (39) Snyder, G. J.; Toberer, E. S. Complex Thermoelectric Materials. Nat. Materials 2008, 7, 105114. (40) Leitner, D. M. Thermal Boundary Conductance and Thermal Rectification in Molecules. J. Phys. Chem. B 2013, 117, 12820-12828.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(41) Xu, Y.; Leitner, D. M. Vibrational Energy Flow through the Green Fluorescent Pro˘ SWater teinâA ¸ Interface: Communication Maps and Thermal Boundary Conductance. J. Phys. Chem. B 2014, 118, 7818-7826. (42) Rego, L. G. C.; Kirczenow, G. Phys. Rev. Lett. 1998, 81, 232-235. (43) Meir, Y.; Wingreen, N. S. Landauer Formula for the Current Through an Interacting Electron Region. Phys. Rev. Lett. 1992, 68, 2512-2516. (44) Redfield, A. G. On the Theory of Relaxation Processes. IBM J. Res. Dev. 1957, 1, 19-31. (45) Aslangul, C.; Pottier, N.; Saint-James, D. Spin-Boson Systems: Equivalence Between the Dilute-Blip and the Born Approximations. J. Physique 1986, 47, 1657-1661. (46) Wertheimer, R.; Silbey, R. On Excitation Transfer and Relaxation Models in LowTemperature Systems. Chem. Phys. Lett. 1980, 75, 243-248. (47) Zarea, M.; Ratner, M. A.; Wasielewski, M. R. Electron Transfer in a Two-Level System Within a Cole-Davidson Vitreous Bath. J. Chem. Phys. 2014, 140, 024110. (48) Dekker, H. Noninteracting-Blip Approximation for a Two-Level System Coupled to a Heat Bath. Phys. Rev. A 1987, 35, 1436-1437. (49) Makri, N.; Makarov, D. E. Tensor Propagator for Iterative Quantum Time Evolution of Reduced Density Matrices. I. Theory. J. Chem. Phys. 1995, 102, 4600-4610. (50) Makri, N.; Makarov, D. E. Tensor Propagator for Iterative Quantum Time Evolution of Reduced Density Matrices. II. Numerical Methodology. J. Chem. Phys. 1995, 102, 4611-4618. (51) Mujica-Martinez, C.; Nalbach, P.; Thorwart, M. Quantification of non-Markovian Effects in the Fenna-Matthews-Olson Complex. Phys. Rev. E 2013, 88, 062719. (52) Feynman, R. P.; Hibbs, A. R. Quantum Mechanics and Path Integrals; McGraw-Hill, NewYork, 1965. 42 ACS Paragon Plus Environment

Page 42 of 46

Page 43 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(53) Makri, N. Iterative Evaluation of the Path Integral for a System Coupled to an Anharmonic Bath. J. Chem. Phys. 1999, 111, 6164-6167. (54) Segal, D.; Millis, A. J.; Reichman, D. R. Numerically Exact Path-Integral Simulation of Nonequilibrium Quantum Transport and Dissipation. Phys. Rev. B 2010, 82, 205323. (55) Simine L.; Segal, D. Path-Integral Simulations with Fermionic and Bosonic Reservoirs: Transport and Dissipation in Molecular Electronic Junctions. J. Chem. Phys. 2013, 138, 214111. (56) Kulkarni, M.; Tiwari, K. L.; Segal, D. Full Density Matrix Dynamics for Large Quantum Systems: Interactions, Decoherence and Inelastic Effects. New J. Phys. 2013, 15 013014. (57) Mitra, A.; Millis, A. J. Spin Dynamics and Violation of the Fluctuation Dissipation Theorem in a Nonequilibrium Ohmic Spin-Boson Model. Phys. Rev. B 2005, 72, 121102(R). (58) Weiss, S.; Eckel, J.; Thorwart, M.; Egger, R. Iterative Real-Time Path Integral Approach to Nonequilibrium Quantum Transport. Phys. Rev. B 2008, 77, 195316. (59) Segal, D.; Reichman, D. R.; Millis, A. J. Nonequilibrium Quantum Dissipation in SpinFermion Systems. Phys. Rev. B 2007, 76, 195316. (60) Chang, L.-D.; Chakravarty, S. Dissipative Dynamics of a Two-State System Coupled to a Heat Bath. Phys. Rev. B 1985, 31, 154-164. (61) Chen, Y.-C. Theory of Quantum Dynamics in Fermionic Environment: An Influence Functional Approach. J. Stat. Phys. 1987, 47, 17-55. (62) Chen, Y.-C. A New Method for Quantum Processes in Fermionic Heat Baths. J. Stat. Phys. 1987, 49, 811-826. (63) Mitra, A.; Millis, A. J. Coulomb Gas on the Keldysh Contour: Anderson-Yuval-Hamann Representation of the Nonequilibrium Two-Level System. Phys. Rev. B 2007, 76, 085342. 43 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(64) Segal, D. Single Mode Heat Rectifier: Controlling Energy Flow Between Electronic Conductors. Phys. Rev. Lett. 2008, 100, 105901. (65) Volker, K. Dynamical Behavior of the Dissipative Two-State System. Phys. Rev. B 1998, 58, 1862-1871. (66) Day, J. T. Note on the Numerical Solution of Integro-Differential Equations. The Computer Journal 1967, 9, 394-395. (67) Makri, N. The Linear Response Approximation and Its Lowest Order Corrections: An Influence Functional Approach. J. Phys. Chem. B 1999, 103, 2823-2829. (68) The heat current through a single mode (frequency ∆) harmonic junction is given by jq = 2 R∞ ¯ωT π 0 h

(ω )[nL (ω ) − nR (ω )]d ω with T (ω ) =

ω 2 ΓL (ω )ΓR (ω ) . 69 (ω 2 −∆2 )2 +[ΓR (ω )+ΓL (ω )]ω 2

(69) Segal, D.; Nitzan, A.; Hänggi, P. Thermal conductance through molecular wires. J. Chem. Phys. 2003, 119, 6840-6855. (70) BR results of Fig. 8 were calculated in the fermionic picture, by following Ref. 28 with the Golden-Rule transition rates evaluated by an explicit integration, to allow the band-edge effects to be properly accounted for. Cutoff effects become more pronounced at high temperatures and large tunneling frequencies, kB T /¯hωc > 0.2, ∆/ωc > 0.1. (71) Wang, J.-S. Quantum Thermal Transport from Classical Molecular Dynamics. Phys. Rev. Lett. 2007, 99, 160601. (72) Grifoni, M.; Winterstetter, M.; Weiss, U. Coherences and Populations in the Driven Damped Two-State System. Phys. Rev. E 1997, 56, 334-345. (73) Ren, J.; Hänggi, P.; Li, B. Berry-Phase-Induced Heat Pumping and Its Impact on the Fluctuation Theorem. Phys. Rev. Lett. 2010, 104, 170601. (74) Chen, T.; Wang, X. B.; Ren, J. Dynamic Control of Quantum Geometric Heat Flux in a Nonequilibrium Spin-Boson Model. Phys. Rev. B 2013, 87, 144303. 44 ACS Paragon Plus Environment

Page 44 of 46

Page 45 of 46

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(75) Uchiyama, C. Non-adiabatic effect on the quantum heat flux control. Phys. Rev. E 2014, 89, 052108.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 46 of 46

jq TL

TR

S

ACS Paragon Plus Environment

From Dissipative Dynamics to Studies of Heat Transfer at the

Recommend Documents