Energy current and its statistics in the nonequilibrium spin-boson model: Majorana fermion representation
aa r X i v : . [ c ond - m a t . m e s - h a ll ] D ec Energy current and its statistics in the nonequilibrium spin-bosonmodel: Majorana fermion representation
Bijay Kumar Agarwalla and Dvira Segal Chemical Physics Theory Group, Department of Chemistry,and Centre for Quantum Information and Quantum Control,University of Toronto, 80 Saint George St.,Toronto, Ontario, Canada M5S 3H6 (Dated: May 6, 2019)
Abstract
We study the statistics of thermal energy transfer in the nonequilibrium (two-bath) spin-bosonmodel. This quantum many-body impurity system serves as a canonical model for quantum energytransport. Our method makes use of the Majorana fermion representation for the spin operators, incombination with the Keldysh nonequilibrium Green’s function approach. We derive an analyticalexpression for the cumulant generating function of the model in the steady state limit, and showthat it satisfies the Gallavotti-Cohen fluctuation symmetry. We obtain analytical expressions for theheat current and its noise, valid beyond the sequential and the co-tunnelling regimes. Our resultssatisfy the quantum mechanical bound for heat current in interacting nanojunctions. Resultsare compared with other approximate theories, as well as with a non-interacting model, a fullyharmonic thermal junction. . INTRODUCTION The spin-boson (SB) model comprises a two-state system (spin) interacting with a dissi-pative thermal environment, a collection of harmonic modes. It is one of the (conceptually)simplest, yet non-trivial models in the theory of open quantum systems [1, 2]. The modelhas found diverse applications in condensed phases physics, chemical dynamics, and quan-tum optics. In particular, it offers a rich platform for studying complex physical processessuch as dissipative spin-dynamics [1], charge and energy transfer phenomena in condensedphases [1, 3], Kondo physics [2], and decoherence dynamics of superconducting qubits [2, 4].In such applications, the spin system can represent donor-acceptor charge states, a magneticimpurity [1], or a truncated harmonic spectrum, mimicking an anharmonic oscillator [5, 6].The bosonic bath may stand for a collection of lattice phonons, electromagnetic modes,bound electron-hole pairs, and other composite bosonic excitations [1, 2].Beyond questions over quantum decoherence, dissipation, and thermalization, which canbe addressed by the ‘canonical’ SB model, the two-bath, nonequilibrium spin-boson (NESB)model has been put forward as a minimal model for exploring the fundamentals of thermalenergy transfer in anharmonic nano-junctions [5]. When the two reservoirs are maintainedat different temperatures—away from linear response—nonlinear functionality such as thediode effect can develop in the junction [5–9]. More generally, the NESB model serves asa building block for addressing fundamental and practical challenges in thermal conductionin nanoscale gaps [10–14], quantum heat engine operation [14–18], molecular conductionjunctions [12, 19–23] and nano-scale energy conversion devices [21, 24].From the theoretical perspective, the NESB model is an extremely rich platform forstudying nonequilibrium quantum physics. One is interested in studying its transport char-acteristics, including transient dynamics and steady state properties, while covering differ-ent regimes: low-to-high temperatures, weak-to-strong system-bath coupling, adiabatic-to-nonadiabatic spin dynamics, with or without a (magnetic) spin biasing field, from linearresponse to the far from equilibrium regime. This challenge could be tackled by extendingopen quantum system methodologies, previously developed to treat the dissipative dynam-ics of the (traditional) SB model, to treat the more complicated, nonequilibrium, two-bathversion.Among the techniques developed to study the characteristics of the thermal heat current2n the NESB model we recount perturbative quantum master equation tools: Redfield equa-tion [5–7, 9, 25–28], the noninteracting blip approximation (NIBA) [6, 7, 28, 29], as well asKeldysh nonequilibrium Green’s function (NEGF) methods [30, 31]. Computational studieshad further established the non-monotonic behavior of the heat current with the spin-bathcoupling energy, including studies based on the multi-layer multi-configuration Hartree ap-proach [32], the iterative influence functional path integral technique (in the spin-fermionrepresentation) [11, 33], Monte Carlo simulations [10], and the hierarchical equation of mo-tion [34, 35].Beyond the analysis of the thermal conductance or the energy current, in small systemsthe fluctuations of the current are expected to reveal plethora of information, such as currentcorrelations to all orders [36, 37]. In fact, rather than focusing on the thermal conductance,it is more demanding yet highly profitable to pursue the probability distribution P t ( Q ) ofthe transferred energy Q within a certain interval of time t . This measure is also known asfull-counting statistics (FCS) in the context of electron transport. Obtaining the FCS forinteracting systems is a highly desirable, yet formidable task. The FCS of the NESB modelhas been analyzed so far in two different limits: (i) in the sequential tunneling limit i.e.,to the lowest order in the system-bath tunnelling strength, by employing the Redfield-typequantum master equation approach [27, 29], and (ii) in the strong coupling and/or hightemperature regime, following NIBA-type quantum master equations [29, 38]. A theoryinterpolating these two limits was presented in Ref. [28]. However, these studies still missthe low temperature limit [31].In this paper, we use the Schwinger-Keldysh NEGF approach [39, 40] in combinationwith the Majorana fermion representation for the system-spin operators, and obtain thecumulant generating function (CGF) of the NESB model beyond the weak spin-bath couplinglimit. The crucial impetus to introduce the Majorana representation is that in the fermionrepresentation we are able to use Wick’s theorem, thus obtain relevant nonequilibrium spin-spin correlation functions—while including the counting parameter. Our results go beyondthe sequential and co-tunnelling limits, and we are particularly able to capture to all ordersthe low temperature regime. As well, we observe deviations from the weak spin-bath couplinglimit. On the other hand, while our result is valid beyond the strictly weak-coupling limit,it does miss the strong-coupling behavior as received in Refs. [10, 11, 28, 32, 41] using theNIBA approach. The main outcome of our study is an analytic expression for the FCS of the3ESB model, capturing quantum effects, interactions, and far-from-equilibrium function.The paper is organized as follows. We introduce the nonequilibrium spin-boson modeland the Majorana fermion representation in Sec. II. In Sec. III, we present our mainresults for the CGF, followed by a discussion over different limits and numerical examples.We further compare our expressions to previous theories on the NESB model, and to theharmonic oscillator-junction model. We conclude in Sec. V. The derivation of the CGF isexplained in details in the Appendix. II. MODEL
The NESB model comprises a two-state (spin) system coupled to two bosonic reservoirs( ν = L, R ), which are maintained at different temperatures. The generic form of the fullHamiltonian is H = ~ ω σ z + ~ ∆2 σ x + X j,ν ~ ω j,ν b † j,ν b j,ν + σ z X j,ν ~ λ j,ν ( b j,ν + b † j,ν ) . (1)Here, σ i ( i = x, y, z ) are different components of the Pauli matrix, ω and ∆ represents leveldetuning and the hopping between the spin states, respectively. b † j,ν ( b j,ν ) is the creation(annihilation) operator of the j -th phonon mode in the ν -th reservoir. The last term de-scribes the system-bath coupling term with λ j,ν as the coupling strength. For simplicity, wefocus here on the unbiased case with degenerate spin levels ( ω = 0). Performing a unitarytransformation, given by U = √ ( σ x + σ z ), the transformed Hamiltonian reads¯ H = ~ ∆2 σ z + X j,ν ~ ω j,ν b † j,ν b j,ν + σ x X j,ν ~ λ j,ν ( b j,ν + b † j,ν ) . (2)We are interested here in obtaining the steady state energy current and its statistics beyondthe weak system-bath coupling limit. Unlike the Redfield master equation technique, whichcaptures only resonant energy transfer processes due to its underlying weak coupling ap-proximation [6], the Keldysh nonequilibrium Green’s function (NEGF) method offers a wellestablished procedure so as to treat the system-bath interaction in a systematic-perturbativeway [39, 40]. However, the validity of Wick’s theorem is a crucial requirement for practicingthe method. Due to the lack of standard bosonic or fermionic commutation relations for spinoperators, the NEGF approach is in fact unsuitable to be used in the spin representation of4he NESB model. However this problem can be avoided by mapping the impurity spin tofermions, using the Majorana-fermion representation [31, 42, 43].Explicitly, the spin operators can be expressed as ~σ = − i ~η × ~η , i.e ., σ x = − i η y η z , σ y = − iη z η x , σ z = − iη x η y . (3)Majorana fermions satisfy the anti commutation relation, { η α , η β } = 0, for α = β , η α = 1,and unlike the Dirac fermions, they are real η α = η † α . Therefore, these fermions can beconstructed in terms of ordinary Dirac fermions ( f , g ) and their conjugates as η x = ( f + f † ) , η y = i ( f † − f ) , η z = ( g + g † ) . (4)In this context, it is important to introduce the so-called copy-switching operator τ x = − iη x η y η z , (5)in terms of which the Majorana fermions can be expressed as σ α = τ x η α . Note that τ x commutes with all Majorana fermion operators and therefore is a constant of motion. Also, τ x = 1. With the help of this operator, the spin-spin correlator reduces to correlatorinvolving two Majorana fermions h σ α ( t ) σ β ( t ′ ) i = h τ x ( t ) η α ( t ) τ x ( t ′ ) η β ( t ′ ) i = h η α ( t ) η β ( t ′ ) i . (6)In this mixed Majorana-Dirac representation, the full Hamiltonian reads¯ H = ~ ∆2 (cid:0) − f † f (cid:1) + X j,ν ~ ω jν b † j,ν b j,ν + (cid:0) f † − f (cid:1) η z ( B L + B R ) , (7)where B ν ≡ P j ~ λ j,ν ( b j,ν + b † j,ν ) is a ν bath operator coupled to the spin system. Note that inthis representation, the system-bath coupling term is no longer given in a bilinear form. Forlater use, we also identify the components of the Hamiltonian as ¯ H = H S + H L + H R + H SB ,with H S = ~ ∆2 (cid:0) − f † f (cid:1) , H ν = X j ~ ω j,ν b † j,ν b j,ν , H SB = (cid:0) f † − f (cid:1) η z ( B L + B R ) . (8) III. FCS: MAIN RESULTSA. Working expressions for the FCS
The complete information over the energy transport statistics can be obtained from theso-called cumulant generating function, G ( ξ ), for heat exchange. We begin by defining the5nergy current operator as the rate of change of energy in one of the reservoirs, say L , andwrite down the heat current as I L ( t ) = − dH HL ( t ) dt . The operators are written in the Heisenbergpicture, and they evolve with respect to the total Hamiltonian ¯ H in Eq. (7). Therefore, thetotal energy change in the L solid within the time interval t = 0 to t , where t ( t ) is theinitial (final) observation time, is given by the integrated current Q L ( t, t ) = Z tt =0 I L ( t ′ ) dt ′ = H L (0) − H HL ( t ) . (9)Following this definition, we write down the characteristic function Z ( ξ ) based on the two-time measurement protocol [36, 37], Z ( ξ ) = D e iξH L e − iξH HL ( t ) E = D U † ξ/ ( t, U − ξ/ ( t, E , = D T c exp h − i ~ Z c H ξ ( τ ) SB ( τ ) dτ iE , (10)Here, ξ is the “counting-field”, keeping track of the net amount of energy transferred fromthe solid L to the spin. h ... i represents an average with respect to the total density matrix atthe initial time, ρ T (0). We assume a factorized initial state, ρ T (0) = ρ L (0) ⊗ ρ R (0) ⊗ ρ S (0),with reservoirs prepared at a canonical state with inverse temperature β ν = T − ν , ρ ν (0) = e − β ν H ν / Tr ν [ e − β ν H ν ], and an arbitrary state for the spin system ρ S (0). We also use thedefinition, U p ( t, ≡ e ipH L e − i ¯ Ht e − ipH L = e − i ¯ H p t/ ~ , (11)for the counting field-dependent unitary evolution. Here, p = ± ξ/ H p acquires a phase in the system-bath coupling term, modifying onlythe left-bath operators,¯ H p = ~ ∆2 (cid:0) − f † f (cid:1) + X j,ν ~ ω j,ν b † j,ν b j,ν + (cid:0) f † − f (cid:1) η z ( B pL + B R ) (12)Here, B pL = P j ~ λ j,L ( b j,L e − ip ~ ω j,L + b † j,L e ip ~ ω j,L ) is the a bath operator, dressed by the countingfield. In the second line of Eq. (10), the operators are written in the interaction picture withrespect to the non-interacting part of the Hamiltonian H S + H L + H R . T c is the contour-ordered operator which orders operators according to their contour time; earlier contour-timeoperators are placed to the right of later-time terms. In the long time limit, the CGF is6efined as G ( ξ ) ≡ lim t →∞ t ln Z ( ξ ) = lim t →∞ t ∞ X n =1 ( iξ ) n n ! hh Q n ii . (13)Here, hh Q n ii represent cumulants. Specifically, the second cumulant is hh Q ii = h Q i−h Q i .Taking derivatives of the CGF with respect to ξ immediately hands over the current and itshigher order fluctuations, or cumulants. However, instead of working with the CGF directly,one can manipulate the so-called generalized current, defined as I ( ξ ) ≡ ∂ G ( ξ ) ∂ ( iξ ) , (14)by following the nonequilibrium version of Feynman-Hellman theorem first introduced byGogolin et al. [44]— in the context of counting statistics for charge transport. The keyadvantage in treating the generalized current, rather than the CGF, lies in the fact that theproblem can be treated with the diagrammatic NEGF technique, as developed originally—without the counting field [45–48].Using the NEGF with counting fields as developed in [49], an expression for the general-ized energy current can be formally organized as I ( ξ ) = Z ∞−∞ dω π ~ ω h ˜Π
To receive the generalized current, our primary objective is to obtain the components˜Π >xx ( ω ). These terms are obtained using the NEGF method following a first order pertur-bation expansion with respect to the interaction of the bath with the spin. We summarizehere the central results; details are given in the Appendix.The lesser and greater components are obtained to the lowest non-zero order in thenonlinear self-energy. They are given as˜Π
1) + n R ( ω ) ¯ n L ( ω )( e − iξ ~ ω − i . (21)If we eliminate the counting parameter, ξ = 0, Π >xx ( ω ) provides the imaginary componentsof the response function Π Rxx ( ω ), Im [Π Rxx ( ω )] = − (cid:0) Γ L ( ω ) + Γ L ( ω ) (cid:1) ( ω − ∆ ) + ω (cid:2) Γ L ( ω )(1 + 2 n L ( ω )) + Γ R ( ω )(1 + 2 n R ( ω )) (cid:3) , (22)matching the results of Ref. (31).Using these expressions, the CGF for the NESB model, G SB ( ξ ) ≡ R ξ dξ ′ I ( ξ ′ ), is obtainedas G SB ( ξ )= Z ∞−∞ dω π ∆ ω ln n T SB ( ω ; T L , T R ) h n L ( ω )¯ n R ( ω )( e iξ ~ ω −
1) + n R ( ω )¯ n L ( ω )( e − iξ ~ ω − io , (23)8ith the temperature-dependent transmission function T SB ( ω ; T L , T R ) = 4 Γ L ( ω ) Γ R ( ω ) ω (cid:16) ω − ∆ (cid:17) + ω h Γ L ( ω )(1+2 n L ( ω )) + Γ R ( ω )(1+2 n R ( ω )) i . (24)This expression is valid with an arbitrary form for the spectral function Γ ν ( ω ). The CGFfurther satisfies the steady state Gallavotti-Cohen fluctuation symmetry, G ( ξ ) = G ( − ξ + i ( β R − β L )) [51]. Eq. (23) constitutes the main result of our work.The cumulants of the energy flux can be readily obtained by taking derivatives of theCGF with respect to the counting field ξ . For example, the heat current and its noise aregiven by h I i SB ≡ ∂ G SB ( ξ ) ∂ ( iξ ) (cid:12)(cid:12)(cid:12) ξ =0 = Z ∞−∞ dω π ~ ∆ ω T SB ( ω ; T L , T R ) (cid:2) n L ( ω ) − n R ( ω ) (cid:3) , (25) h S i SB ≡ ∂ G SB ( ξ ) ∂ ( iξ ) (cid:12)(cid:12)(cid:12) ξ =0 = Z ∞−∞ dω π ( ~ ∆) n − T SB ( ω ; T L , T R ) (cid:2) n L ( ω ) − n R ( ω ) (cid:3) + T SB ( ω ; T L , T R ) (cid:2) n L ( ω )¯ n R ( ω )+ n R ( ω )¯ n L ( ω ) (cid:3)o . (26)The result for the current agrees with the derivation in Ref. [31]—once we organize ourexpressions, R ∞−∞ dω... → × R ∞ dω ... In the next subsection, we discuss interesting limitsof the general results. C. Special limits
Incoherent sequential tunnelling.
When the system-bath coupling is weak and the reservoirs’temperatures are high, Γ
L,R ≪ ∆ ≤ T L,R , the above generating function reduces to theresult obtained from the Redfield quantum master equation approach [29], when directlyemploying the Born-Markov approximation. We now derive this result. Following Eq. (23),the generalized current can be simplified to I SB ( ξ ) = Z ∞−∞ dω π ∆ M ′ ( ω, ξ )( ω − ∆ ) + ω M ( ω, ξ ) , (27)where M ′ ( ω, ξ ) = ∂M ( ω,ξ ) ∂ ( iξ ) . To the lowest order O (Γ L,R ), working in the limit Γ
L,R ≪ ∆ ≤ T L,R , the poles in the integrand can be approximated by ± n ∆ ± i p M (∆ , ξ ) o . (28)9y employing the residue theorem, the integration in Eq. (27) results in I weakSB ( ξ ) = ∂ √ M ( ξ ) ∂ ( iξ ) and the generating function reduces to G weakSB ( ξ ) = − (cid:16) C (∆) − p M (∆ , ξ ) (cid:17) . (29)This expression matches the result obtained in Ref. [29]. This CGF also respects the fluc-tuation symmetry. It immediately yields the heat current in the weak coupling limit [5] h I i weakSB = ~ ∆ Γ L (∆)Γ R (∆) [ n L (∆) − n R (∆)][Γ L (∆)(1 + 2 n L (∆))] + [Γ R (∆)(1 + 2 n R (∆))] . (30) Co-tunnelling.
At low temperatures, Γ ν ≪ T ν ≤ ∆, the process of sequential tunnelling isexponentially suppressed since incoming phonons are off-resonance—with frequencies belowthe spin energy gap, ω ≪ ∆. The dominant contribution to the current and higher orderfluctuations thus comes from coherent two-phonon co-tunnelling processes. In this limit, thetransmission function of Eq. (24) is given by T coSB ( ω, T L , T R ) ∼ L ( ω )Γ R ( ω ) ω / ∆ ≪
1. Byapproximating ln(1 + x ) ∼ x for small x , we reduce the CGF of Eq. (23) to G coSB ( ξ ) = 2 π Z ω h dω Γ L ( ω )Γ R ( ω )∆ (cid:16) n L ( ω )¯ n R ( ω )( e iξ ~ ω −
1) + n R ( ω )¯ n L ( ω )( e − iξ ~ ω − (cid:17) , (31)with fluctuation symmetry being satisfied. Here, ω h , the upper limit in the integral shouldbe determined by the smaller energy scale, temperature of the cutoff frequency of the baths.The co-tunneling (co) heat current then becomes h I i coSB = 2 π Z ω h dω ~ ω Γ L ( ω )Γ R ( ω )∆ [ n L ( ω ) − n R ( ω )] . (32)This expression was previously achieved in two ways: (i) By using a systematic perturbativetreatment [25], and (ii) working with the so-called Born-Oppenheimer approach for heatexchange [52], by assuming slow bath and a fast (high frequency) impurity. In the case ofan Ohmic bath, Γ ν ( ω ) ∝ ω s with s = 1, the heat current scales as h I i coSB ∝ T L − T R , thusthe thermal conductance scales with T , in agreement with numerically exact simulationson the NESB model [10]. As well, in this low temperature limit the NESB junction behavessimilarly to a fully harmonic junction, as we discuss in Sec. III D.Note that in contrast to the CGF received in Eq. (23) and Eq. (29), the CGF in theco-tunnelling limit is symmetric with respect to Γ L,R ( ω ). Therefore, in this limit the systemdoes not support the thermal rectification effect. Moreover, in this limit the cumulants C n = ∂ n G SB ( ξ ) ∂ ( iξ ) n (cid:12)(cid:12) ξ =0 scale as C n ∝ / ∆ , whereas in the sequential tunneling limit cumulantsgrow as C n ∝ ∆ n . 10 . Comparison between the NESB model and the harmonic oscillator junction In the harmonic oscillator (HO) junction, a single harmonic oscillator of frequency ω ,replaces the spin impurity of the NESB model, Eq. (1). The resulting Hamiltonian is fullyharmonic, and it can be readily solved exactly to yield the CGF [53, 54] G HO ( ξ )= − Z ∞−∞ dω π ln h − T HO ( ω ) (cid:16) n L ( ω )¯ n R ( ω )( e iξ ~ ω −
1) + n R ( ω )¯ n L ( ω )( e − iξ ~ ω − (cid:17)i . (33)Surprisingly, our final expression for the CGF of the NESB model, Eq. (23), is very similarto this expression. The following differences show up: (i) In the HO case the transmissionfunction does not depend on the temperatures of the baths, T HO ( ω ) = 4 Γ L ( ω ) Γ R ( ω ) ω (cid:16) ω − ω (cid:17) + ω (cid:16) Γ L ( ω ) + Γ R ( ω ) (cid:17) . (34)Further, (ii) there is a crucial sign difference in this CGF as compared to G SB ( ξ ) in Eq. (23).This sign difference reflects on the nonlinear nature of the spin. A similar sign-differencebetween harmonic and spin impurity nanojunctions has been observed in vibrationally-assisted electron conducting junctions [21, 24]. The above expression immediately providesthe Landauer expression for the heat current, h I i HO = 14 π Z ∞−∞ dω ~ ω T HO ( ω )[ n L ( ω ) − n R ( ω )] , (35)and the noise h S i HO = 14 π Z ∞−∞ dω ( ~ ω ) (cid:2) T HO ( ω ) ( n L ( ω ) − n R ( ω )) + T HO ( ω ) ( n L ( ω )¯ n R ( ω ) + ¯ n R ( ω ) n L ( ω )) (cid:3) . (36)In the weak coupling limit, the CGF of the HO model reduces to the standard result obtainedby a low order QME [27, 29] G weakHO ( ξ ) = 12 (cid:16) C HO ( ω ) − q C HO ( ω ) − A HO ( ω , ξ ) (cid:17) , (37)with C HO ( ω ) = Γ L ( ω )+Γ R ( ω ) ,A HO ( ω , ξ ) = Γ L ( ω )Γ R ( ω ) (cid:16) n L ( ω ) ¯ n R ( ω )( e iξ ~ ω −
1) + n R ( ω ) ¯ n L ( ω )( e − iξ ~ ω − (cid:17) . h I i weakHO = ~ ω Γ L ( ω )Γ R ( ω )Γ L ( ω ) + Γ R ( ω ) [ n L ( ω ) − n R ( ω )] . (38)The co-tunnelling limit is more subtle, and we exemplify it now when calculating the current.We break the transmission function (34) into two contributions (leaving for a moment thenumerator) T HO ( ω ) = T o ( ω ) + T e ( ω ), T o ( ω ) = ω − ω ( ω − ω ) + [Γ L ( ω ) + Γ R ( ω )] ω T e ( ω ) = ω ( ω − ω ) + [Γ L ( ω ) + Γ R ( ω )] ω . (39)Assuming the hierarchy of energies Γ ν ≪ T ν ≤ ω , we note that the function ω [ n L ( ω ) − n R ( ω )]changes slowly at the vicinity of ω , in the regime where the functions T e,o ( ω ) have significantweight. Therefore, the integral (35) over the odd component (approximately) cancels out,and the current is solely determined by the even term, T e ( ω ) ∼ /ω , to yield h I i coHO = 2 π Z ω h dω ~ ω Γ L ( ω )Γ R ( ω ) ω [ n L ( ω ) − n R ( ω )] . (40)This result reproduces exactly the behavior of the NESB model in the corresponding limit,Eq. (32). This correspondence is not surprising: At low temperatures (smaller than theenergy spacing in the quantum impurity) and at weak system-bath coupling, the NESBand the HO junctions should behave rather similarly. For a comprehensive analysis of theharmonic-mode thermal junction, see Ref. [55]. E. Steady state population and a bound on heat current
Besides transport properties, we use the Majorana formalism and calculate the steadystate population of the ground and excited states in the eigenbasis of the spin. This can beobtained by calculating h σ z i , given as, h σ z i = i Z ∞−∞ dω π h (cid:16) (cid:17) G > Ψ ( ω ) − − (cid:16) − (cid:17) G < Ψ ( ω ) i = − ∆ π Z ∞−∞ dω ω Γ L ( ω ) + Γ R ( ω )( ω − ∆ ) + ω C ( ω ) . (41)The function C ( ω ) is defined in Eq. (21). In the weak coupling limit, we receive the sameresult as obtained in Ref. [6], h σ z i weak = − Γ L (∆) + Γ R (∆)Γ L (∆) [1 + 2 n L (∆)] + Γ R (∆) [1 + 2 n R (∆)] . (42)12he population of the states are p g = (1 − h σ z i ) and p e = (1 + h σ z i ).Recently, a rigorous quantum mechanical bound for the heat current in interacting sys-tems has been derived, valid at the high temperature—yet in the quantum regime [56]. Wenow confirm that the heat current derived in our work, Eq. (25), does not violate the bound.This further affirms the validity and usefulness of our result.In the following analysis we make use of the inequality 0 ≤ [ n L ( ω ) − n R ( ω )] ≤ ( T L − T R ) / ( ~ ω ) for ω > T L > T R . As well, we recall on the positivity of the transmissionfunction T SB ( ω ) >
0. Furthermore, we assume an Ohmic spectral density function for thereservoirs, Γ ν ( ω ) = γ ν ω, ν = L, R (see Ref. [56] for a detailed discussion over differentspectral functions). Putting these pieces together, we conclude that the heat current of Eq.(25) satisfies the following inequality h I i ≤ Z ∞−∞ dω π γ L γ R ∆ ω ( T L − T R )( ω − ∆ ) + ω C ( ω )= 2 π ∆ ( T L − T R ) γ L γ R γ L + γ R Z ∞ dω ω γ L + γ R ( ω − ∆ ) + ω C ( ω )= − ∆( T L − T R ) γ L γ R γ L + γ R h σ z i (43)which precisely matches with the bound organized in Ref. [56] for the NESB model. Weconclude that our expression for the current thus does not violate a fundamental bound,unlike the prediction of the Redfield QME, see Ref. [56]. IV. NUMERICAL RESULTS
In Figs. 1-3, we present simulations demonstrating the behavior of the heat current h I i and the second cumulant h S i , based on Eq. (23), as a function of the system-bath coupling,averaged temperature, and temperature difference. We focus on the following questionsregarding the operation of the NESB nanojunction:(i) How are the current and noise influenced by the system-bath coupling strength? (Fig.1 and 3). (ii) What are the signatures of operation far from equilibrium, as opposed to thelinear response regime? (Fig. 1 and 3) (iii) What is the temperature dependence of the heatcurrent? (Fig. 2) (iv) Thermal diode effect: Can we enhance this effect if we go beyondthe weak spin-bath coupling? (3) (v) What is the relation between the Majorana-basedtreatment and other techniques? (Figs. 1-3).13ig. 1 displays the current and the noise as obtained from Eqs. (26), as well as the weakcoupling (Redfield) limit [11, 29], and the NIBA approximation [11, 41]. We use an Ohmicspectral function for the baths with an exponential cutoff, Γ ν ( ω ) = πα ν ωe − ω/ω c . In accordwith previous results (for the heat current [11]), we find that Redfield equation dramaticallyoverestimates the current and the noise in comparison to the (more accurate) Majoranaand NIBA results. Majorana treatment shows a saturation of the current and its noiseat large α , while under NIBA these quantities quickly decay beyond α ∼ .
15. Since thetemperature is rather high, ∆ = T a , with T a = ( T L + T R ) /
2, we expect the NIBA to berather accurate here [10, 11, 41]. We also confirm in panel (a) that in linear response (LR),the conductance, h I i LR / ∆ T , is proportional to the thermal noise in the junction, in accordwith the Green-Kubo relation, h S i eq = 2 T a h I i LR / ∆ T. (44)Far from equilibrium [see panel (b)], we obviously observe violations of the above relation.However, it is interesting to note that the current and noise still follow a similar functionalform within the three different methods.Fig. 2 displays the temperature dependence of the current and the noise. We study boththe NESB model and a fully harmonic junction, Eq. (35) and (36), and make the followingobservations: (i) Comparing the current in the HO and NESB nanojunctions, anharmonicity,as realized here by the spin, leads to the suppression of the heat current. (ii) At weakcoupling, α ν = 0 .
01, see panels (a1)-(b1), the Majorana and Redfield approaches for theNESB model agree. (iii) At intermediate coupling, α ν = 0 .
2, see panels (a2)-(b2), Redfieldformalism leads to (nonphysical) high currents, even beyond the harmonic limit—at lowtemperatures. (iv) At high temperatures and intermediate coupling, Majorana calculationsshow (a weak) decay of the current with temperature, see panel (a2), an effect expected toshow up in anharmonic nanojunctions [12].Next, we discuss the operation of the NESB as a heat diode, as suggested in Ref. [5].To materialize this effect, it is necessary to (i) include anharmonic interactions, and (ii)introduce a spatial asymmetry [8]. The NESB model naturally includes an anharmonicpotential. We break here the left-right symmetry by using different coupling strengths atthe contacts, α L = α R . In Fig. 3, we analyze the ratio between the forward and backwardcurrents as we switch the temperatures of the two baths, R ≡ |h I ( T L , T R ) i| / |h I ( T R , T L ) i| .14 α ν (a) ∆ T =0.05 T a MajoranaRedfield NIBA sc a l ed c u rr en t and no i s e α ν (b) ∆ T = T a MajoranaRedfield NIBA
FIG. 1. Scaled current h I i / ∆ T (dashed lines) and noise h S i / T a (symbols) for the NESB modelas a function of coupling strength α ν , employing different theoretical schemes: Redfield (blue),Majorana (Red) and NIBA (purple). (a) Results close to equilibrium, ∆ T = 0 . T a . (b) Calcu-lations far-from-equilibrium, ∆ T = T a , demonstrating deviations from the fluctuation-dissipationtheorem. Parameters are ∆ = T a = 1, T L,R = T a ± ∆ T / ω c = 10∆, and α L = α R . We set α L =0.01, 0.2, and modify α R over a broad range of values.Based on Eq. (30), we can readily confirm that under the Redfield formalism the rectifi-cation ratio R does not depend on the absolute value of α (given the linearity of the currentwith α ), only on the ratio α R /α L . In contrast, the Majorana treatment, which goes beyondweak coupling, reveals that the diode effect is enhanced as we increase the coupling strengthitself. This result points out to the crucial role of many-body interactions in realizing thediode function. V. CONCLUSIONS
We have studied the statistics of energy transfer in the nonequilibrium spin-boson model.By combining Majorana fermion representation for the spin operators with the Schwinger-Keldysh Green’s function approach, we were able to derive an analytical expression forthe CGF of the model. This function, which we confirmed here to satisfy the fluctuationsymmetry for heat exchange, hands over the complete information over the energy statisticsin the steady state limit. Our approach goes beyond the weak-coupling (Redfield) and the co-15 −3 T a h I i (a1) 0 1 200.050.1 T a (b1) h S i Majorana NESB Redfield NESB HO junction0 1 200.020.040.06 T a h I i (a2) 0 1 200.51 T a (b2) h S i FIG. 2. Temperature dependence of the heat current and noise for the SB junction (Redfield andMajorana) and the HO model. (a1)-(b1) Weak coupling limit α ν = 0 .
01. (a2)-(b2) Intermediatecoupling, α ν = 0 .
2. Parameters are ∆ = 1, ∆ T = 0 . T a and T L,R = T a ± ∆ T / ω c = 10∆. tunnelling limits. Surprisingly, the CGF of the NESB model has a similar structure as in theharmonic oscillator junction, besides sign differences and the appearance of a temperature-dependent transmission function—in the NESB model. These differences reflect on thenonlinear nature of the spin-boson system.We have presented numerical examples for the heat current and its noise, and comparedour results to previously-developed quantum master equation approaches, namely Redfieldand the NIBA. We have further demonstrated that a heat diode becomes more effective as weincrease the system-bath coupling. Additional improvements to the Majorana formulationpresented here could be made, e.g., by developing a polaron-transformed Majorana fermion-NEGF approach [58]. Future work will be focused on simulating counting statistics in theNESB model beyond perturbative approaches [57].16 α R / α L R e c t i f i c a t i on r a t i o ( a ) Red fi eld α L = 0 . , . α L =0.01Majorana α L =0.2 α R / α L h S i Red fi eld T L < T R Red fi eld T L > T R Majorana T L < T R Majorana T L > T R α L =0.2, | ∆ T | /T a =1( b ) FIG. 3. Thermal diode effect. (a) Rectification ratio, R = h I ( T L = 1 . , T R = 0 . i|h I ( T L =0 . , T R = 1 . i| , as a function of the asymmetry in the system-bath coupling, α R /α L , while fixing α L . (b) Noise h S i for forward and backward operations as a function of the junction asymmetryusing α L = 0 .
2. Parameters are ∆ = 1, ∆ T = T a , ( T L , T R )=(1.5, 0.5) and (0.5, 1.5), ω c = 10∆. ACKNOWLEDGMENTS
The work of DS and BKA was supported by an NSERC Discovery Grant, the CanadaResearch Chair program, and the CQIQC at the University of Toronto.
APPENDIX A: DERIVATION OF THE CUMULANT GENERATING FUNC-TION WITHIN AN NEGF APPROACH
Our goal is to evaluate the generalized current, Eq. (15). It is given in terms of the(dressed) lesser ˜Π
1) + n R ( ω )¯ n L ( ω )( e − iξ ~ ω − (cid:9)i . (A20)with ¯ n ν ( ω ) = 1 + n ν ( ω ). Using Eq. (A4), the lesser and greater components of spin-spincorrelation functions are finally obtained as˜Π
Understanding Quantum Phase Transitions , edited by L. D. Carr (Taylor andFrancis, Boca Raton, 2010).[5] D. Segal and A. Nitzan, Phys. Rev. Lett , 034301 (2005).[6] D. Segal, Phys. Rev. B , 205415 (2006).[7] D. Segal and A. Nitzan, J. Chem. Phys. , 194704 (2005).[8] L.-A. Wu and D. Segal, Phys. Rev. Lett. , 095503 (2009).[9] L.-A. Wu, C. X. Yu, and D. Segal, Phys. Rev. E , 041103 (2009).[10] K. Saito and T. Kato, Phys. Rev. Lett. , 214301 (2013).
11] N. Boudjada and D. Segal, J. Phys. Chem. A, , 11323 (2014).[12] D. Segal and B. K. Agarwalla., Annu. Rev. Phys. Chem. , 185 (2016).[13] J.-S. Wang, J. Wang, J. T. Lu, Eur. Phys. J. B , 381 (2008).[14] N. Li et al. Rev. Mod. Phys. , 1045 (2012).[15] D. Segal and A. Nitzan, Phys. Rev. E , 026109 (2006).[16] D. Segal, Phys. Rev. Lett. , 260601 (2008).[17] J. Ren, P. H¨anggi, and B. Li, Phys. Rev. Lett. , 170601 (2010).[18] T. Chen, B. X. Wang, and J. Ren, Phys. Rev. B , 144303 (2013).[19] L. Simine and D. Segal, Phys. Chem. Chem. Phys. , 13820 (2012).[20] L. Simine and D. Segal, J. Chem. Phys. , 014704 (2014).[21] B. K. Agarwalla, J.-H. Jiang, and D. Segal, Phys. Rev. B , 245418 (2015).[22] J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment ,World Scientific, Singapore , 2010.[23] M. Galperin, M. A. Ratner, A. Nitzan, J. Phys. Condens. Matter , 103201 (2007).[24] B. K. Agarwalla, J.-H. Jiang and D. Segal, Beilstein J. Nanotechnol. , 2129 (2015).[25] T. Ruokola and T. Ojanen, Phys. Rev. B , 045417 (2011).[26] J. Thingna, H. Zhou, J. S. Wang, J. Chem. Phys. , 194101 (2014).[27] J. Ren, P. H¨anggi, and B. Li, Phys. Rev. Lett. , 170601 (2010).[28] C. Wang, R. Jie. J. Cao, Sci. Rep. , 11787 (2015).[29] L. Nicolin and D. Segal, J. Chem. Phys. , 164106 (2011).[30] K. A. Velizhanin, M. Thoss, and H. Wang, J. Chem. Phys. , 084503 (2010).[31] Y. Yang and C. Q. Wu, Euro. Phys. Lett. , 30003 (2014).[32] K. A. Velizhanin, H. Wang, and M. Thoss, Chem. Phys. Lett. , 325 (2008).[33] D. Segal, Phys. Rev. B , 195436 (2013).[34] A. Kato and Y. Tanimura, J. Chem. Phys, , 064107 (2015).[35] J. Cerrillo, M. Buser, T. Brandes, arXiv:1606.05074.[36] M. Esposito, U. Harbola, and S. Mukamel, Rev. Mod. Phys. , 1665 (2009).[37] M. Campisi, P. H¨anggi, and P. Talkner, Rev. Mod. Phys. , 771 (2011).[38] L. Nicolin and D. Segal, Phys. Rev. B , 161414 (2011).[39] J. Rammer, Quantum Field Theory of Non-Equilibrium States (Cambridge University Press,2007).
40] J.-S. Wang, B. K. Agarwalla, H. Li, and J. Thingna, Front. Physics , 673 (2014).[41] D. Segal, Phys. Rev. E , 012148 (2014).[42] A. Shnirman and Y. Makhlin, Phys. Rev. Lett. , 207204 (2003).[43] W. Mao, P. Coleman, C. Hooley, and D. Langreth. Phys. Rev. Lett. , 207203 (2003).[44] A. O. Gogolin and A. Komnik, Phys. Rev. B , 195301 (2006).[45] J. Schwinger, Brownian motion of a quantum oscillator, J. Math. Phys., , 407 (1961).[46] L. P. Kadanoff and G. Baym, Quantum Statistical Mechanics , Benjamin/Cummings, 1962[47] L. V. Keldysh, Diagram technique for nonequilibrium processes, Sov. Phys. JETP, , 1018(1965).[48] P. Danielewicz, Quantum theory of nonequilibrium processes (I) , Ann. Phys., , 239 (1984).[49] H. Li, B. K. Agarwalla, B. Li, and J. -S. Wang, Eur. Phys. J. B , 1 (2013).[50] Y. Meir and N. Wingreen, Phys. Rev. Lett. , 2512 (1992).[51] G. Gallavotti and E. G. D. Cohen, Phys. Rev. Lett. , 2694 (1995).[52] L.-A. Wu and D. Segal, Phys. Rev. E , 051114 (2011).[53] K. Saito and A. Dhar, Phys. Rev. Lett. , 180601 (2007).[54] B. K. Agarwalla, B. Li, and J.-S. Wang, Phys. Rev. E , 051142 (2012).[55] Y. Vinkler-Aviv, A. Schiller, and N. Andrei, Phys. Rev. B , 024307 (2014).[56] E. Taylor and D. Segal, Phys. Rev. Lett. , 220401 (2015).[57] M. Carrega, P. Solinas, A. Braggio, M. Sassetti, and U. Weiss, New J. Phys. , 045030(2015).[58] J. Liu, H. Xu, and C.-Q. Wu, http://dx.doi.org/10.1016/j.chemphys.2016.07.003 Chem. Phys.2016