Nonequilibrium open quantum systems with multiple bosonic and fermionic environments: A hierarchical equations of motion approach
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Nonequilibrium open quantum systems with multiple bosonic and fermionicenvironments: A hierarchical quantum master equation approach
J. B¨atge, Y. Ke, C. Kaspar, and M. Thoss
1, 2 Institute of Physics, University of Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg, Germany EUCOR Centre for Quantum Science and Quantum Computing,University of Freiburg, Hermann-Herder-Str. 3, D-79104 Freiburg, Germany (Dated: February 19, 2021)We present a hierarchical quantum master equation (HQME) approach, which allows a numericallyexact simulation of nonequilibrium transport in general open quantum systems involving multiplebosonic and fermionic environments. The performance of the method is demonstrated for a model ofa nanosystem, which involves interacting electronic and vibrational degrees of freedom and is coupledto fermionic and bosonic baths. The results show the intricate interplay of electronic and vibrationaldegrees of freedom in this nonequilibrium transport scenario for both voltage and thermally driventransport processes. Furthermore, the use of importance criteria to improve the efficiency of themethod is discussed.
I. INTRODUCTION
Open quantum systems, which involve the coupling toone or multiple reservoirs, are ubiquitous in many areas ofphysics, including nanophysics, quantum thermodynam-ics, quantum optics, and quantum information theory.
They are also relevant for a variety of novel technologi-cal developments that involve a nanoscale system coupledto a macroscopic environment, such as quantum informa-tion devices , nanoelectronic setups , and thermoelec-tric generators on the scale of single molecules . Thetheoretical description of open quantum systems requiresin most cases numerical methods. In this context, severalapproximate methods as well as numerically exacttechniques have been developed and applied. Withthe exception of the path-integral approach formulatedby Simine and Segal , to our knowledge, all applicationsof numerically exact approaches focus either on the treat-ment of bosonic or fermionic environments so far. How-ever, especially in the realm of quantum thermodynam-ics and quantum transport the inclusion of both typesof environments in the model is often necessary. To thisend, we extend in this paper the hierarchical quantummaster equation approach (HQME) for an efficient treat-ment of open quantum systems with multiple bosonic andfermionic environments.The HQME method (also called hierarchical equationof motion (HEOM) approach) was originally developedby Tanimura and Kubo for the description of relax-ation dynamics in open quantum systems induced by thecoupling to bosonic environments. Later it was extendedby several groups for the investigation of charge and en-ergy transport caused by fermionic environments. As summarized by Tanimura, so far studies usingthe HQME approach have focused on either bosonic orfermionic environments. In this work, we formulate theHQME method for more general models, depicted in Fig.1, which include multiple bosonic and fermionic environ-ments. This extension allows to study more general sce-narios such as electronic and phononic heat transport FIG. 1. Sketch of an open quantum systems coupled to mul-tiple fermionic and bosonic environments. induced by bias voltages or temperature differences.The outline of this paper is as follows: The generaliza-tion of the model system and the HQME formalism areintroduced in Sec. II. In Sec. III, we discuss the per-formance of the method for examplary models, includingboth voltage and thermally driven charge and heat trans-port. Section IV concludes with a summary. Throughoutthe paper, we use units where ~ = 1 and k B = 1. II. METHOD
We consider an open quantum system described by asubsystem S, which is coupled to both bosonic (B) anda fermionic (F) environments. The total Hamiltonian isgiven by H = H S + H SB + H B + H SF + H F . (1)The fermionic environment H F may consist of multiplemacroscopic noninteracting fermionic baths H F = X α H F α = X α,k ε kα c † kα c kα . (2)Here, c † kα / c kα denote the fermionic creation/annihilationoperator associated with state k of bath α with energy ε kα . Nonequilibrium situations can be introduced by as-signing different chemical potentials µ F α and/or temper-atures T F α to the different baths. Similarly, the bosonicenvironment H B may contain more than one noninter-acting bosonic bath H B = X α H B α = X α,j ω jα b † jα b jα , (3)where b † jα / b jα denote the bosonic creation/annihilationoperator associated with state j of bath α with frequency ω jα . Nonequilibrium situations can be achieved bychoosing different temperatures T B α for different bosonicbaths. For simplicity, we assume vanishing chemical po-tentials µ B α = 0, as it holds true for phonons and pho-tons.For the validity of the subsequent derivation of the hi-erarchical quantum master equation (HQME), we restrictourselves to system-environment coupling terms whichare linear in the bath creation/annihilation operators H SF = X α,k ν kα c † kα W + h.c. , (4a) H SB = V X α,j ξ jα (cid:16) b † jα + b jα (cid:17) . (4b)Both V and W are operators acting only on the subsys-tem degrees of freedom (DOFs). They are not restrictedin their action on the bosonic DOFs, however, W has toact as a generalized annihilaton operator on the fermionicDOFs and overall reduce one fermionic occupation num-ber by one. In contrast, V must not change the fermionicoccupation numbers and may involve sums of even num-bers of fermionic annihilation/creation operators. Theinfluence of the bosonic/fermionic environment on thesystem dynamics is encoded in the so-called spectral den-sity functions, J B α ( ω ) = π X j | ξ jα | δ ( ω − ω jα ) , (5a) J F α ( ε ) =2 π X k | ν kα | δ ( ε − ε kα ) . (5b)With these general specifications of the open quantumsystem, we present the general idea of the derivation ofthe HQME, which is outlined in more detail in Refs. 36and 38. The central quantity of the approach is the re-duced density matrix ρ ( t ) of the subsystem, defined bythe trace over the DOFs of the environment of the totaldensity matrix ̺ ( t ), ρ ( t ) = Tr B+F { ̺ ( t ) } . (6)Employing a bath-interaction picture, we avoid the di-rect appearance of bath related dynamical phases in theequation of motion of the reduced density matrix. In this picture, the time-dependence of the reduced densitymatrix can be written as ρ ( t ) = Tr B+F (cid:8) U ( t, ̺ (0) U † ( t, (cid:9) , (7)with the time-evolution operator U ( t, t ) = T (cid:16) e − i R tt dτH SB ( τ )+ H SF ( τ )+ H S (cid:17) . (8)Here, T denotes the time-ordering operator and the time-dependence of H SB / H SF originates from the chosen bath-interaction picture.Exploiting the Gaussian properties of the nonin-teracting environments, a formally exact path-integralfor the reduced density matrix of the subsystem in-volving a Feynman-Vernon influence functional can bederived. As a consequence of noninteracting Gaussian baths,all information about system-bath couplings is encodedin the time correlation functions of the free baths C B α ( t − τ ) = X j | ξ jα | D b † jα ( t ) b jα ( τ ) + b jα ( t ) b † jα ( τ ) E B α , (9a) C s F α ( t − τ ) = X k | ν kα | (cid:10) c skα ( t ) c ¯ skα ( τ ) (cid:11) F α . (9b)Here, h · i X α = Tr X α { · ρ X α } denotes the expectationvalue with respect to the corresponding free bath ρ X α = exp( − β Xα ( H Xα − µ Xα N Xα ))Tr Xα { exp( − β Xα ( H Xα − µ Xα N Xα )) } and s ∈ {− , + } has beenintroduced with c − kα = c kα and c + kα = c † kα . Further-more, N X α represents the particle number operator ofthe bosonic or fermionc particles and β X α describes theinverse temperature 1 /T X α in the respective bath X α .The time correlation functions are related to the bathspectral densities (cf. Eqs. (5a) and (5b) ) via C B α ( t ) = Z ∞ dω J B α ( ω ) π × (cid:20) coth (cid:18) β B α ω (cid:19) cos( ωt ) − i sin( ωt ) (cid:21) , (10a) C s F α ( t ) = Z ∞−∞ dε J F α ( ε )2 π e siεt e sβ F α ( ε − µ F α ) + 1 . (10b)Within the HQME approach, the correlation functionsare described by sums over exponentials, C B α ( t ) = Λ B α X p η p,α e − γ p,α t , (11a) C s F α ( t ) = Γ F α X q η q,α,s e − γ q,α,s t , (11b)where common strategies to determine the weights η anddecay rates γ include the Matsubara and thePad´e decomposition. Furthermore, Λ B α and Γ F α are coupling strengths of the respective baths and aredefined by the respective spectral density. The exponen-tially decaying time correlation functions ensure a closedset of equations of motion. Under the constraints on the Hamiltonian statedabove, the influence functional factorizes into two partswhich originate from the bosonic and fermionic envi-ronment separately.
The separate influence func-tionals can be obtained as shown in detail in previousworks.
The hierarchical equations of motion areconstructed by consecutive time-derivatives of the influ-ence functional. The combination of the hierarchies orig-inating from the different environments leads to ∂∂t ρ ( m | n ) g | h = − i L S + m X l =1 γ g l + n X l =1 γ h l ! ρ ( m | n ) g | h − X h x A h x ρ ( m | n +1) g | h + x − n X l =1 ( − l C h l ρ ( m | n − g | h − l + X g x B g x ρ ( m +1 | n ) g + x | h + m X l =1 D g l ρ ( m − | n ) g − l | h , (12)with the multi-indeces g = ( p, α ) resp. h = ( q, α, s ),the notation for the multi-index vectors v = v ··· v p , v + x = v ··· v p v x , and v − l = v ··· v l − v l +1 ··· v p , and L S O = [ H S , O ]. To faciliate a unified treatment of thebosonic and fermionic hierarchy, the index managementfor the bosonic part of the hierarchy differs from the con-vention in other publications. Here, ρ (0 | repre-sents the reduced density operator of the subsystem, andthe higher-tier auxiliary density operators (ADOs) ρ ( m | n ) g | h encode the influence of the environment on the subsystemdynamics. Tier by tier the ADOs introduce the effect ofhigher order system-bath correlations. These ADOs alsoinclude the information on time-local transport observ-ables like energy or particle currents. Employingthe generating functional approach such observablesare related to compositions of ADOs.The operators A h x and C h l connect the n th-fermionic-tier ADO to the ( n +1)th- resp. ( n − B g x and D g l connect the m th-bosonic-tier ADO to the ( m +1)th- resp. ( m − A h x ρ ( m | n ) g | h =Γ (cid:16) W s hx ρ ( m | n ) g | h + ( − ( n ) ρ ( m | n ) g | h W s hx (cid:17) , (13a) B g x ρ ( m | n ) g | h =Λ h V, ρ ( m | n ) g | h i , (13b) C h l ρ ( m | n ) g | h = ( − n η h l W ¯ s hl ρ ( m | n ) g | h − η ∗ ¯ h l ρ ( m | n ) g | h W ¯ s hl , (13c) D g l ρ ( m | n ) g | h = η g l V ρ ( m | n ) g | h − η ∗ g l ρ ( m | n ) g | h V, (13d)leading to a two-fold hierarchy of equations of motion.Generalyzing these operators for linear superpositions ofthe system-bath couplings defined in Eqs. (4b) and (4a),the superindices g and h , the operators V and W as well as the coupling strengths Λ and Γ aquire an additionalindex.For applications, only a finite number of poles char-acterizing the baths can be taken into account. Byemploying additional resummation schemes this limita-tion can be overcome. Here, the hierarchy needs to betruncated in a suitable manner. For this purpose pre-vious works have introduced an importance criterion forpure fermionic environments, which estimates the impor-tance of an ADO to the dynamics of the system.
If the importance I assigned to an ADO is below agiven threshold value θ , it is neglected. This reduces thenumerical complexity to a manageable level and allowsfor a systematic convergence. For details of the conver-gence properties of the HQME method, we refer to Refs.41, 64–66. Here, we generalize the criterion to the com-bined hierarchy and assign to each operator ρ ( m | n ) g | h theimportance value I (cid:16) ρ ( m | n ) g | h (cid:17) = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n Y l =1 Γ P a ∈{ ..l } Re [ γ h a ] η h l Re [ γ h l ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) × (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) m Y l =1 Λ P a ∈{ ..l } Re [ γ g a ] η g l Re [ γ g l ] (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) . (14)The importance of an ADO is weighted according to itsinfluence within the fermionic and bosonic hierarchy si-multaneously. Motivated by the factorized structure ofthe HQME, we add a factor considering the bosonic partof the two-fold hierarchy to the importance estimate forthe ADOs in the fermionic hierarchy presented by H¨artle et al. , which stems from the HQME in the stationarycondition. This factor is similarly constructed as forthe fermionic environment. It assesses the importancedue to the respective order in the bosonic part of the hi-erarchy and it reflects the relative importance of an ADOwithin one hierarchical order for the stationary state.
III. ILLUSTRATIVE APPLICATIONS
In this section, we illustrate the performance of theHQME approach for the combination of bosonic andfermionic environments. To this end, we consider a basicmodel for the subsystem involving interacting electronicand vibrational degrees of freedom. The results of thenumerically exact HQME method are compared to thewell-established Born-Markov master equation approach(BMME). Explicit formulas and detailed derivations forthe BMME can be found in Refs. 2, 13, 16, and 67. TheBMME is perturbative in the system-environment cou-pling up to the second, i.e. lowest non-vanishing orderand assumes a separation of relaxation time-scales in thebath and the system.
A. Model
FIG. 2. Sketch of the model consisting of one electronic stateand a single vibrational mode coupled to two bosonic andtwo fermionic baths, describing vibrational heat baths andelectronic leads, respectively.
As sketched in Fig. 2, we consider a model systemthat consists of a single electronic state and a single vi-brational mode, described by the Hamiltonian H S = Ω a † a + ε d † d + λ ( a † + a ) d † d. (15)Here, a † and a are the creation and annihilation oper-ators of the vibrational mode with frequency Ω and d † and d are the creation and annihilation operators of theelectronic state with energy ε . The electronic state iscoupled to the vibrational mode via the coupling con-stant λ . The couplings to the respective environmentshave the form H SF = X α,k (cid:16) ν kα c † kα d + h.c. (cid:17) , (16a) H SB =( a † + a ) X α,j ξ jα (cid:16) b † jα + b jα (cid:17) , (16b)and are specified by the spectral densities J B α ( ω ) =Λ B α ω Ω ω cα ω cα + ω (17a) J F α ( ε ) =Γ F α D α D α + ( ε − µ F α ) . (17b)Here, we include the system frequency Ω in the definitionof the bosonic bath spectral density for two purposes.First, the behavior of the spectral density in the limit ω → ω cα . Sec-ond, the coupling strength Λ B α has the dimension of anenergy and defines a decay time. The Lorentzian sup-pression of high-energy contributions by the cutoff fre-quency and the bandwidth D α is beneficial for the closed Set Γ D α Λ ω cα T Ω ε A 0.005 30 0.005 0.2 0.025 0.2 0.2B 0.05 30 0.0425 0.05 0.025 0.2 0.2TABLE I. Model parameters used in the calculations. Allparameters are given in e V. derivation of the HQME and allows for a straightforwarduse of the Pad´e decomposition scheme to approximatethe correlation functions. Throughout this work, weassume symmetric configurations, i.e. Γ L = Γ R = Γ2 , andΛ L = Λ R = Λ2 .It is noted that this model has been studied exten-sively without bosonic baths using both approximateand numerically exact approaches. Forthe explanation of the transport phenomena, the sys-tem part of the Hamiltonian is often diagonalized by thesmall polaron transformation via H = SHS † with S =exp (cid:0) λ Ω d † d ( a − a † ) (cid:1) . Thereby, the charging and decharg-ing processes of the system induced by the fermionic en-vironment are accompanied by de-/excitations of the vi-brational mode with their transition amplitudes givenby the so-called Franck-Condon matrix elements and the electronic-vibrational interaction term in the sys-tem transforms to a shift of the electronic state energy ε = ε − λ Ω . Applying the small polaron transforma-tion also affects the coupling term to the bosonic envi-ronment. Accordingly, the shift of the electronic stateenergy and the Franck-Condon elements are no longerexact but perturbative up to lowest non-vanishing orderin the system-environment coupling. A non-perturbativegeneralization of the transformation leads to a non-linearsystem-environment coupling. Therefore, we do not usesuch a transformation for the combined fermionic andbosonic HQME considered here.Subsequent to the introduction of the obervables ofinterest, we show results for two parameter regimes spec-ified in more detail in Tab. I. Parameter set A is cho-sen such that the perturbative treatment within theBMME is expected to be valid. In contrast, the system-environment coupling is enhanced as well as the cutofffrequency in the bosonic bath ω cα is reduced in param-eter set B. Both changes weaken the validity of the per-turbative treatment of the model system by the BMME. B. Observables
System observables, such as the population of the elec-tronic level or the vibrational excitation, can be directlyevaluated from the reduced density matrix of the system.For environment-related observables, we employ the gen-erating functional approach to find expressions interms of elements of the hierarchy. For the charge cur-rent from lead F α into the system, we obtain I F α = − e (cid:28) dN F α dt (cid:29) = e Γ F α X h α s h Tr n d ¯ s hα ρ (0 | | h α o . (18)The electronic (vibrational) heat current from lead F α (bath B α ) into the system is given by J F α = − (cid:28) dH F α dt (cid:29) + µ F α (cid:28) dN F α dt (cid:29) = − i Γ F α X h α ω h α Tr n d s hα ρ (0 | | h α o − µ F α e I F α , (19a) J B α = − (cid:28) dH B α dt (cid:29) = i Λ B α X g α ω g α Tr n ( a + a † ) ρ (1 | g α | o + Im { C B α (0) } Tr n ( a + a † ) ρ (0 | o . (19b)It should be noted that there exist different definitionsof the heat current. For example, the expression J B α = − D dH B α dt E − D dH SB α dt E was used in Refs. 61 and 74. Inthe steady state, which is considered in this work, thesedifferent definitions are are equivalent to Eq. (19b), be-cause the term D dH SB α dt E vanishes.It can be shown that expectation values of observables,which are quadratic in system and/or environment anni-hilation or creation operators, are exact for a noninter-acting system ( λ = 0), if all contributions of fermionicand bosonic ADOs up to second tier are taken into ac-count in the calculation. This statement includes thecharge and heat currents and is a generalization of thestatement for fermionic model systems. For the inter-acting system considered here, the tier of ADOs requiredfor convergence is checked in test calculations.
C. Voltage driven charge and heat transport
First, we consider charge and heat transport drivenby an external bias voltage. We start by recapitulatingthe well known charge current-voltage characteristics ofthe model system and briefly discuss the changesinduced by the bosonic environment. To this end, weexpose the model system to a symmetric external biasvoltage Φ, i.e. µ L = − µ R = e Φ2 . As we do not apply atemperature difference and are not considering any otherasymmetries of the contacts, the effect of the bosonicbaths can be described by a single bath with the couplingstrengths Λ = Λ L + Λ R and the environmental tempera-ture T .Figure 3 presents the vibrational excitation h a † a i ,the charge current I , and the differential conductance dId Φ as a function of bias voltage Φ for two differentelectronic-vibrational interaction strengths, λ Ω =1 . . . . . . . . h a † a i BMMEHQME λ Ω =1 . and Λ=0 λ Ω =0 . I [ µ A ] d I d Φ [ µ AV ] Φ [V] . . . . . . . h a † a i BMMEHQME λ Ω =1 . and Λ=0 λ Ω =0 . I [ µ A ] d I d Φ [ µ AV ] Φ [V]
FIG. 3. Vibrational excitation h a † a i , charge current I ,and differential conductance dId Φ as a function of bias voltagecalculated for two different electronic-vibrational interactionstrengths λ . Numerically exact HQME results (solid lines) arecompared to perturbative BMME results (dashed lines). Inpanel (a) and (b), the parameters are chosen according to pa-rameter sets A and B, respectively. The green solid line showsthe HQME results without a bosonic environment (Λ = 0) for λ Ω = 1 . lines) and λ Ω =0 . is slightly visible for weak electron-vibrational interaction ( λ Ω = 0 .
2) and becomes more pro-nounced for stronger interaction strengths ( λ Ω = 1 . e Φ ≈ ε + n Ω) with n ∈ N .At these bias voltages, resonant transport processes be-come active. The reduction of step heights for larger biasvoltages and the according saturation of the charge cur-rent are related to the structure of the Franck-Condonmatrix. Without relaxation of the vibrational mode by abosonic bath (Λ = 0), another prominent phenomenonis the vibrational instability, which describes very highexcitations of the vibrational mode for large biasvoltages and especially for weak electronic-vibrationalinteractions.
As pointed out before, this effectis based on the dominance of electron-hole pair cre-ation processes for the cooling mechanism. Introduc-ing the additional vibrational relaxation by an environ-ment, this dominance is reduced with increasing couplingstrength.
As a result, the vibrational excitation forweak electronic-vibrational interaction is lower than forstronger interactions in the regime of large bias voltages.Furthermore, the vibrational relaxation reduces the pop-ulation of highly vibrational excited states of the system.This affects the charge transport processes for strongelectronic-vibrational interaction in two different ways.For small bias voltages, processes starting with an ex-cited vibration are less likely and the charge current isaccordingly lowered. On the contrary, transitions with ahigh energy transfer into the vibration are favored by theelectronic-vibrational interaction for higher bias voltages.Their likelihood is enhanced by the vibrational environ-ment and consequently the charge current is increased.In addition, it is known that the damping of a bare har-monic oscillator by a bosonic bath leads to a reduction ofthe harmonic oscillator frequency.
The effect of thisfrequency renormalization is reflected in the shift of thepeak positions in the differential conduction to lower biasvoltages with the introduction of the bosonic bath.For weak system-environment coupling, the Born-Markov approximation is expected to be valid. Thisis confirmed by the agreement of the BMME and theHQME results in this parameter regime (cf. Fig. 3 (a)).In contrast, stronger coupling strengths increase the im-portance of co-tunneling and other higher-order pro-cesses, which result in a broadening of the electronic en-ergy level and, for interacting systems, in an additionalshift of the characteristic energies of transport processes(cf. Fig. 3 (b)). As the BMME allows only for sequen-tial tunneling processes, it neglects these effects and de-viations to the numerically exact HQME become morepronounced for stronger coupling strengths.Besides differences in the charge current, the bosonicenvironment also induces a heat current from the elec-tronic into the vibrational degrees of freedom. Figure4 depicts the heat current J vib = − P α J B α from theelectronic environment into the vibrational environmentand its differential behavior as a function of bias volt-age. We observe a step-like increase of the vibrationalexcitation and the heat current with bias voltage, whichis in accordance with the expected excitation of the vi-brational mode and its damping by the vibrational envi- . .
51 0 0 . . . . J v i b [ e V p s ] BMMEHQME λ Ω =1 . λ Ω =0 . d J v i b d Φ [ e p s ] Φ [V] . . . . J v i b [ e V p s ] BMMEHQME λ Ω =1 . λ Ω =0 . d J v i b d Φ [ e p s ] Φ [V]
FIG. 4. Heat current J vib from the electronic into the bosonicenvironment as a function of bias voltage and its differentialbehavior. In panel (a) and (b), the parameters are chosenaccording to parameter sets A and B, respectively. ronment. Analogous to the charge current, the steps inthe heat current are connected to the onset of resonantemission processes that involve more and more vibra-tional quanta. In addition, a decrease of the relative stepheights with decreasing electronic-vibrational interactionis observed, which is caused by the suppression of pro-cesses that involve multiple vibrational quanta. However,the step heights remain almost constant with increasingbias voltage. The latter observation is a result of theinterplay of the reduction of the Franck-Condon matrixelements for processes with a higher energy transfer intothe vibrational degree of freedom and the increase ofthe transfered energy by these processes.Increasing the non-Markovianity of the environmen-tal damping by reducing the cutoff frequency ω c of thebosonic bath and increasing the coupling strengths toboth environments, the BMME results deviate signifi-cantly from the numerically exact HQME results (cf. Fig.4 (b)). Again, we observe that the higher-order processesincluded in the HQME lead to a broadening in the arisingstep structure and a shift of the step positions. D. Thermally driven charge and heat transport
Another important type of transport processes innanostructures is charge and heat transport induced bya temperature difference of the reservoirs. For purelyelectronic transport, such processes result in thermo-electric phenomena. Here, we apply temperature differ-ences ∆ T by heating the left fermionic and bosonic bath T L = T + ∆ T but without introducing an external biasvoltage. The temperature of the right reservoirs remainsconstant.Figure 5 shows the charge current and its differentialbehavior as a function of the applied temperature dif-ference. Overall, the thermally driven charge currentdoes not exhibit such pronounced features as the biasvoltage induced charge current, because the tempera-ture is a smoother driving force. Nevertheless, we ob-serve a suppression of the charge current for strongerelectronic-vibrational interactions, which occurs due tothe Franck-Condon blockade of transport channels. Forweak electronic-vibrational interaction, the thermallydriven charge transport is mainly carried by elastic pro-cesses. The peak in the differential conductance canbe understood by the difference between the Fermi-distribution functions evaluated at the energy of the po-laron shifted electronic state ε . In contrast, many dif-ferent transport channels are contributing to the chargecurrent in the strong electronic-vibrational interactioncase, as discussed for the bias voltage driven transport inthe section above. As these processes have different andhigher threshold energies, their smooth onset induced bythe temperature difference reduces the peak in the dif-ferential behavior of the charge current.For weak-coupling, the BMME and HQME resultsagree very well (cf. Fig. 5 (a)). For stronger coupling,the differences of the charge currents obtained with theHQME and BMME methods (cf. Fig. 5 (b)) are explainedby the energy level broadening. On the one hand, thecharge current for small temperature differences is largerthan predicted by the BMME, because the broadening ofthe characteristic energy of transport processes enablesthe onset of these processes earlier. On the other hand,a part of the transport processes is shifted to higher en-ergies, which get less effective driven by a temperaturedifference. Accordingly, the charge current is reduced forlarge temperature differences by the broadening.Figure 6 presents the heat current leaving the left elec-tronic bath J F L , the heat current leaving the left bosonicbath J B L and their added differential thermal conductiv-ity dJ L d ∆ T = J FL + J BL d ∆ T as a function of the temperature dif-ference ∆ T . In the weak system-environment coupling regime (cf. Fig. 6 (a)), we observe an increase of theheat current out of the hot electronic bath in accordancewith the behavior of the charge current. Due to Franck-Condon blockade of transport channels, the heat currentfor stronger electronic-vibrational interaction is reduced.In comparison to the charge current, the reduction is lesspronounced as the processes with higher thresholds arerealized by electrons that induce a vibrational energytransfer of multiple vibrational quanta. The heat cur-rent out of the hot vibrational bath rises with increas-ing temperature difference without strong dependenceon the electronic-vibrational coupling for a weak system-environment coupling.Stronger system-environment couplings (cf. Fig. 6 (b))affect the electronic heat current analogous to the chargecurrent. As the cooling of the vibration by electrons . . . . I [ µ A ] BMMEHQME λ Ω =1 . λ Ω =0 . d I d ∆ T [ µ A e V ] ∆ T [ e V] .
505 0 0 . . I [ µ A ] BMMEHQME λ Ω =1 . λ Ω =0 . d I d ∆ T [ µ A e V ] ∆ T [ e V] FIG. 5. Charge current I and its differential behavior withrespect to the applied temperature difference. In panel (a)and (b), the parameters are chosen according to parametersets A and B, respectively. . . J F L [ e V p s ] BMMEHQME λ Ω =1 . λ Ω =0 . J B L [ e V p s ] d J L d ∆ T [ µ A e V ] ∆ T [ e V] . . .
205 0 0 . . J F L [ e V p s ] BMMEHQME λ Ω =1 . λ Ω =0 . J B L [ e V p s ] d J L d ∆ T [ µ A e V ] ∆ T [ e V] FIG. 6. Heat currents out of the left contact for the bosonicand fermionic baths and their added differential thermal con-ductivity as a function of temperature difference. In panel(a) and (b), the parameters are chosen according to parame-ter sets A and B, respectively. gets more effective with increasing system-environmentcoupling, the dependence on the electronic-vibrationalcoupling increases as well. E. Computational aspects and numericalconvergence
We finally discuss numerical aspects of the importancecriterion and the numerical convergence of the HQMEmethod. For the computational treatment of the HQMEapproach (cf. Eq. (12)), the number of explicitly con-sidered ADOs is reduced, without loosing any informa-tion, by utilizing the symmetries of the ADOs within one m \ n 0 1 2 3 40 1 26 1 001 19 500 255 7751 4 104 4 004 78 000 1 023 1002 10 260 10 010 195 000 2 557 7503 20 520 20 020 390 000 5 115 5004 35 910 35 035 682 500 8 952 125TABLE II. Total number of independent ADOs of type ρ ( m | n ) g | h in the hierarchy for the bias voltage driven cases. One bosonicbath and two fermionic baths as well as three bosonic andtwelve fermionic Pad´e poles are used. tier, i.e. ρ ( m | n ) g g ··· g m | h ··· h n = ρ ( m | n ) g g g ··· g m | h ··· h n , (20a) ρ ( m | n ) g ··· g m | h h ··· h n = − ρ ( m | n ) g g g ··· g m | h h h ··· h n , (20b) ρ ( m | n ) g ··· g m | h ··· h n =( − m ρ ( m | n ) † g ··· g m | ¯ h n ··· ¯ h . (20c)For an exemplary case, Tab. II shows the number of in-dependent ADOs, which have to be considered in thedifferent tiers of the HQME. It demonstrates the strongincrease of the number of ADOs with increasing order ofthe two-fold hierarchy. Without the application of an im-portance criterion, the large number of ADOs in higherorders renders the numerical simulation intractable, be-cause all of them need to be saved in the memory. Ta-ble III shows that the importance criterion, Eq. (14), re-duces the number of independent ADOs, which have tobe considered, very effectively for the two-fold fermionic-bosonic hierarchy.Next, we consider the consistency of this importancecriterion. To this end, Fig. 7 shows examplarily the effectof using the importance criterion on the convergence ofcharge or heat currents. Specifically, we depict the chargecurrent-voltage and the heat current-voltage character-istics for different threshold values considered. As theimportance estimate for the ADOs within the fermionicpart of the hierarchy is similar to the original definitionby H¨artle et al. , our observations in the bias voltagedriven charge current are also similar. For large biasvoltages, i.e. deep in the resonant transport regime, theresults obtained with a high threshold value show littledeviations to converged results. Approaching the non-resonant regime, e Φ < ε , these deviations become sig- o max \ θ − − − −
02 39 330 1 981 4 957 15 4203 39 383 4 057 21 644 718 4804 39 383 4 223 30 570 19 341 210TABLE III. Number of ADOs, which are considered withinthe HQME for different truncation levels ρ ( m | n ) g | h with m, n ≤ o max and different importance threshold values θ for the pa-rameters from set A (s. Tab. I). . . . . I [ µ A ] θ =10 − θ =10 − θ =10 − θ =10 − J v i b [ e V p s ] bias voltage Φ [V] FIG. 7. Convergence with respect to the used threshold value θ for the charge and heat current as a function of bias voltagefor parameter set B. Included are ADOs ρ ( m | n ) g | h with m, n ≤ nificant. For too large thresholds, we observe unphysicalnegative charge currents or even fail to obtain a station-ary state result. The latter is caused by the occurence ofpositive real parts in the eigenvalues of the propagationmatrix for the system. The influence of these eigenvaluesmight be eliminated by the projection method suggestedby Dunn et al. .The composition of ADOs in the heat current dif-fers from that of the charge current and, thus, the con-vergence behavior with respect to the threshold valuechanges. We observe a deviation for large threshold val-ues over the whole voltage bias window shown. This devi-ation decreases systematically by lowering the thresholdvalue. As observed for the charge current, the conver-gence with respect to the threshold value is slower in thenon-resonant regime.Figure 8 illustrates the convergence of the HQME re-sults for a fixed threshold value of θ = 10 − with increas-ing truncation order o max , which specifies the consideredADOs ρ ( m | n ) g | h m, n ≤ o max . For the parameters consid-ered, the results obtained for the different orders agree inthe non-resonant regime, indicating convergence alreadyfor o max = 2. In the resonant regime, small shifts inthe peak positions in the resonant regime are seen forlower o max , which are more pronounced in the heat cur-rent than in the charge current. These are caused byfrequency corrections introduced by higher-order effectsof the vibrational bath, which require the inclusion ofADOs with at least o max = 3.Finally, we comment on the computational effort re-quired to solve the HQME equations, especially the run-time, using as an example the data in Fig. 3. The datawere obtained by the time-propagation of the HQME intoits stationary state using the Runge-Kutta-Fehlberg algo- . . . . d I d Φ [ µ AV ] d J v i b d Φ [ e p s ] bias voltage Φ [V] FIG. 8. Convergence with respect to the used truncationorder for the charge and heat current as a function of biasvoltage. Included are ADOs ρ ( m | n ) g | h with m, n ≤ , ,
4. Thechosen parameters are ε = 0 . , λ Ω = 1 . , Ω = 0 . , T =0 .
025 eV , Λ = 0 . , Γ = 0 . rithm based on sparse matrix vector multiplication fromthe Intel MKL library. The runtimes for the different pa-rameter regimes needed by a four-threaded execution onfour CPUs of an Intel Xeon E6252 Gold processor are inthe range of five hours for parameter set A and twentyhours for parameter set B. The main reason for the dif-ferent runtimes in the two parameters sets is the differentsystem-environment coupling strengths. The number ofADOs, which have to be considered for converged results,increases strongly for larger coupling strengths. Its effecton the runtime outweighs the effect of reduced propaga-tion times caused by the enhanced damping of the in-ternal system dynamics. We also note that the compu-tational effort does not depend significantly on the biasvoltage and the strength of the electronic-vibrational cou-pling. IV. CONCLUSION
In this work, we have formulated the HQME methodfor open quantum systems involving multiple bosonic andfermionic environments. This extended approach allowsthe study of quantum transport through nanosystems in-duced by bias voltages and/or temperature differences ina numerically exact manner for both types of environ-ments on equal footing. To facilitate efficient simula-tions with the extended HQME method, we have alsoformulated an importance criterion for the elements ofthe hierarchical structure, which reduces the computa-tional effort to a manageable level without compromisingthe ability to control the numerical convergence.0In order to demonstrate the performance of the ex-tended HQME method, we have applied it to a genericinteracting model of electronic and vibrational chargeand heat transport in a molecular junction. The resultsshow the intricate interplay of electronic and vibrationaldegrees of freedom in this nonequilibrium transport sce-nario for both voltage and thermally driven transportprocesses.Applications of the extended HQME method to morecomplex models and/or time-dependent setups are pos-sible and promise to reveal further interesting physicaleffects in nanosystems. The approach may also pro-vide new opportunities to study fundamental questions of quantum thermodynamics in nanosystems beyond lin-ear response and weak coupling in a numerically exactway.
ACKNOWLEDGEMENTS
This work was supported by a research grant of theGerman Research Foundation (DFG). Y.K. thanks theAlexander von Humboldt Foundation for the award of aResearch Fellowship. Furthermore, support by the stateof Baden-W¨urttemberg through bwHPC and the DFGthrough grant no INST 40/575-1 FUGG (JUSTUS 2 clus-ter) is gratefully acknowledged. U. Weiss,
Quantum Dissipative Systems , 2nd ed., Series inmodern condensed matter physics (World Scientific, Sin-gapore, 1999). H.-P. Breuer and F. Petruccione,
The Theory of OpenQuantum Systems (Oxford University Press, Oxford,2007). D. Loss and D. P. DiVincenzo,Phys. Rev. A , 120 (1998). K. D. Petersson, L. W. McFaul, M. D. Schroer,M. Jung, J. M. Taylor, A. A. Houck, and J. R. Petta,Nat. , 380 (2012). P. G. Piva, G. A. DiLabio, J. L. Pitters, J. Zikovsky,M. Rezeq, S. Dogel, W. A. Hofer, and R. A. Wolkow,Nat. , 658 (2005). H. Song, Y. Kim, Y. H. Jang, H. Jeong, M. A. Reed, andT. Lee, Nat. , 1039 (2009). Y. Dubi and M. Di Ventra,Rev. Mod. Phys. , 131 (2011). C. M. Finch, V. M. Garc´ıa-Su´arez, and C. J. Lambert,Phys. Rev. B , 033405 (2009). O. Karlstr¨om, H. Linke, G. Karlstr¨om, and A. Wacker,Phys. Rev. B , 113415 (2011). C. Caroli, R. Combescot, P. Nozieres, and D. Saint-James,J. Phys. C: Solid State Phys. , 21 (1972). P. Hyldgaard, S. Hershfield, J. Davies, and J. Wilkins,Ann. Phys. , 1 (1994). M. Galperin, A. Nitzan, and M. A. Ratner,Phys. Rev. B , 045314 (2006). A. Mitra, I. Aleiner, and A. J. Millis,Phys. Rev. B , 245302 (2004). C. Timm, Phys. Rev. B , 195416 (2008). M. Leijnse and M. R. Wegewijs,Phys. Rev. B , 235424 (2008). R. H¨artle and M. Thoss, Phys. Rev. B , 115414 (2011). H. Ness, S. A. Shevlin, and A. J. Fisher,Phys. Rev. B , 125422 (2001). M. ˇC´ıˇzek, M. Thoss, and W. Domcke,Phys. Rev. B , 125406 (2004). M. Caspary Toroker and U. Peskin,J. Chem. Phys. , 154706 (2007). M. A. Laakso, D. M. Kennes, S. G. Jakobs, and V. Meden,New J. Phys. , 023007 (2014). A. Khedri, V. Meden, and T. A. Costi,Phys. Rev. B , 195156 (2017). A. Khedri, T. A. Costi, and V. Meden,Phys. Rev. B , 195138 (2018). J. Liu and D. Segal, Phys. Rev. B , 155406 (2020). R. H¨utzen, S. Weiss, M. Thorwart, and R. Egger,Phys. Rev. B , 121408(R) (2012). S. Weiss, R. H¨utzen, D. Becker, J. Eckel, R. Egger, andM. Thorwart, Phys. Status Solidi B , 2298 (2013). L. Simine and D. Segal,J. Chem. Phys. , 214111 (2013). L. M¨uhlbacher and E. Rabani,Phys. Rev. Lett. , 176403 (2008). M. Schir´o and M. Fabrizio,Phys. Rev. B , 153302 (2009). J. Klatt, L. M¨uhlbacher, and A. Komnik,Phys. Rev. B , 155306 (2015). H. Wang and M. Thoss,J. Chem. Phys. , 024114 (2009). E. Y. Wilner, H. Wang, M. Thoss, and E. Rabani,Phys. Rev. B , 115145 (2014). H. Wang and M. Thoss,J. Chem. Phys. , 164105 (2016). F. B. Anders and A. Schiller,Phys. Rev. B , 245113 (2006). A. Jovchev and F. B. Anders,Phys. Rev. B , 195112 (2013). G. Cohen and E. Rabani, Phys. Rev. B , 075150 (2011). Y. Tanimura and R. Kubo,J. Phys. Soc. Jpn. , 101 (1989). Y. Tanimura, Phys. Rev. A , 6676 (1990). J. Jin, X. Zheng, and Y. Yan,J. Chem. Phys. , 234703 (2008). X. Zheng, R. Xu, J. Xu, J. Jin, J. Hu, and Y. Yan,Prog. Chem. , 1129 (2012). R. H¨artle, G. Cohen, D. R. Reichman, and A. J. Millis,Phys. Rev. B , 235426 (2013). R. H¨artle, G. Cohen, D. R. Reichman, and A. J. Millis,Phys. Rev. B , 085430 (2015). D. Hou, S. Wang, R. Wang, L. Ye, R. Xu, X. Zheng, andY. Yan, J. Chem. Phys. , 104112 (2015). C. Schinabeck, A. Erpenbeck, R. H¨artle, and M. Thoss,Phys. Rev. B , 201407(R) (2016). L. Z. Ye, X. Zheng, Y. J. Yan, and M. Di Ventra,Phys. Rev. B , 245105 (2016). S. Wenderoth, J. B¨atge, and R. H¨artle,Phys. Rev. B , 121303(R) (2016). C. Schinabeck, R. H¨artle, and M. Thoss,Phys. Rev. B , 235429 (2018). C. Schinabeck and M. Thoss,Phys. Rev. B , 075422 (2020). Y. Tanimura, J. Chem. Phys. , 020901 (2020). R. P. Feynman, Phys. Rev. , 108 (1951). R. Feynman and F. Vernon, Ann. Phys. , 118 (1963). A. Ishizaki and Y. Tanimura,J. Phys. Soc. Jpn. , 3131 (2005). Y. Tanimura, J. Phys. Soc. Jpn. , 082001 (2006). J. M. Moix and J. Cao,J. Chem. Phys. , 134106 (2013). A. Erpenbeck, C. Hertlein, C. Schinabeck, and M. Thoss,J. Chem. Phys. , 064106 (2018). J. Hu, R.-X. Xu, and Y. Yan,J. Chem. Phys. , 101106 (2010). J. Hu, M. Luo, F. Jiang, R.-X. Xu, and Y. Yan,J. Chem. Phys. , 244106 (2011). A. Kato and Y. Tanimura,J. Chem. Phys. , 064107 (2015). A. Erpenbeck and M. Thoss,J. Chem. Phys. , 191101 (2019). For the bosonic hierarchy the index order is not relevant.In contrast the index order is important for the fermionichierarchy. A. Kato and Y. Tanimura,J. Chem. Phys. , 224105 (2016). L. Song and Q. Shi, Phys. Rev. B , 064308 (2017). J. Schwinger, J. Math. Phys. , 407 (1961). Q. Shi, L. Chen, G. Nan, R.-X. Xu, and Y. Yan,The Journal of Chemical Physics , 084105 (2009). B. Popescu, H. Rahman, and U. Kleinekath¨ofer,J. Chem. Phys. , 154103 (2015). I. S. Dunn, R. Tempelaar, and D. R. Reichmann,J. Chem. Phys. , 184109 (2019). H. Rahman and U. Kleinekath¨ofer,J. Chem. Phys. , 244104 (2019). F. Haupt, F. Cavaliere, R. Fazio, and M. Sassetti,Phys. Rev. B , 205328 (2006). J. Koch and F. von Oppen,Phys. Rev. Lett. , 206804 (2005). R. H¨artle, C. Schinabeck, M. Kulkarni,D. Gelbwaser-Klimovsky, M. Thoss, and U. Peskin,Phys. Rev. B , 081404(R) (2018). H. Xie, F. Jiang, H. Tian, X. Zheng, Y. Kwok,S. Chen, C. Yam, Y. Yan, and G. Chen,J. Chem. Phys. , 044113 (2012). J. Koch, F. von Oppen, Y. Oreg, and E. Sela,Phys. Rev. B , 195107 (2004). R. H¨artle, C. Benesch, and M. Thoss,Phys. Rev. Lett. , 146801 (2009). J. Koch, F. von Oppen, and A. V. Andreev,Phys. Rev. B , 205438 (2006). K. A. Velizhanin, H. Wang, and M. Thoss,Chem. Phys. Lett. , 325 (2008). J. Koch, M. Semmelhack, F. von Oppen, and A. Nitzan,Phys. Rev. B , 155306 (2006). R. H¨artle and M. Thoss, Phys. Rev. B , 125419 (2011). H. Grabert, P. Schramm, and G.-L. Ingold,Phys. Rep.168