Current-induced forces and hot-spots in biased nano-junctions
Jing-Tao Lü, Rasmus B. Christensen, Jian-Sheng Wang, Per Hedegård, Mads Brandbyge
CCurrent-induced forces and hot-spots in biased nano-junctions
Jing-Tao L¨u,
1, 2, 3, ∗ Rasmus B. Christensen, Jian-Sheng Wang, Per Hedeg˚ard, and Mads Brandbyge School of Physics and Wuhan National High Magnetic Field Center,Huazhong University of Science and Technology, Wuhan, China Department of Micro- and Nanotechnology, Technical University of Denmark, Kongens Lyngby, Denmark Niels Bohr Institute, Nano-Science Center, University of Copenhagen, Copenhagen, Denmark Department of Physics, National University of Singapore, 117551 Singapore, Republic of Singapore Center for Nanostructured Graphene (CNG), Department of Micro- and Nanotechnology,Technical University of Denmark, Kongens Lyngby, Denmark
We investigate theoretically the interplay of current-induced forces (CIF), Joule heating, andheat transport inside a current-carrying nano-conductor. We find that the CIF, due to the electron-phonon coherence, can control the spatial heat dissipation in the conductor. This yields a significantasymmetric concentration of excess heating (hot-spot) even for a symmetric conductor. Whencoupled to the electrode phonons, CIF drive different phonon heat flux into the two electrodes.First-principles calculations on realistic biased nano-junctions illustrate the importance of the effect.
PACS numbers: 85.75.-d, 85.65.+h,75.75.+a,73.63.Fg
Introduction–
Current-induced forces and Joule heat-ing both originate from the coupling between electronsand phonons[1], one of the most fundamental many-bodyinteractions responsible for a wide range of phenomena inmolecular and condensed-matter physics. Their vital rolein maintaining the electronic device stability is furtherpromoted at nanoscale. Our understanding of the twoclosely related effects, especially their interplay in nano-and atomic-conductors is still under development[2–10].Several forces, present only in the nonequilibrium situa-tion, have been discovered theoretically. Among themare the non-conservative (NC) “wind force”, and theBerry-phase (BP) induced pseudo-magnetic force. Dif-ferent from the stochastic Joule heating[11–23], the NCand BP forces can generate deterministic energy andmomentum transfer between the current-carrying elec-trons and the vibrations in the conductor[2–6]. In care-fully designed devices, this effect may be used to driveatomic motors[2, 7]. Meanwhile, it can also impactthe stability of the device[3, 24, 25]. To this end, thevibrational/phononic[1] heat transport and heat distri-bution in the presence of current flow becomes an urgentproblem to investigate.The electrode phonons play an important role asheat sinks for the locally dissipated Joule heat in theconductor[15]. However, the effects on the heat trans-port of the deterministic CIF, and the momentum trans-fer from the current has so far not been explored.To address this question, we go beyond the previoustreatments[3, 26] considering localized vibrations in theconductor, and include coupling to the phonons in theelectrodes[27]. Employing the semi-classical generalizedLangevin equation(SGLE),[26, 28–30], we find that, inaddition to energy transfer, the CIF also influence howthe excess vibrational energy is distributed in the junc-tion and transported to the electrodes. Using first-principles calculations, we demonstrate how symmetric current-carrying nano-junctions typically possess a sig-nificant asymmetric excess heat distribution with heat accumulation at hot-spots in the junction. At the sametime the phonon heat flow to the two electrodes differs.This behavior is governed by the phases of the electronand phonon wavefunctions, and is a result of electron-hole pair symmetry breaking in the electronic structure.It will have important implications, and should be takeninto account when considering junction disruption athigh bias[24, 31].
Method –
In the SGLE approach we adopt thetwo-probe transport setup, where a “bottleneck” nano-junction (system) is connected to left( L ) and right( R )electrodes. We consider the case where the system regionis characterized by a significant current density and de-viation from equilibrium. The current-carrying electronsare treated as a nonequilibrium bath, coupling linearlywith the system displacement, while the remaining atomsin L and R form two phonon baths interacting with thesystem also via a linear coupling. The electron-phonon(e-ph) coupling Hamiltonian can be written as H eph = (cid:88) i,j,k M kij ( c † i c j + h.c. )ˆ u k . (1)Here, ˆ u k = √ m k ˆ x k is the mass-normalized displacementaway from the equilibrium position of the k -th atomicdegrees of freedom, with m k the mass, and ˆ x k the dis-placement operator from equilibrium position; c † i ( c j ) isthe electron creation(annihilation) operator for the i -( j -)th electronic state in the junction. The coupling matrix, M kij , is local in real space, non-zero in the system andneglected in L, R . We treat the e-ph interaction pertur-batively using the electron and phonon states obtainedfrom the Born-Oppenheimer approximation. In order tofocus on the effect of CIF, we will ignore the change ofHamiltonian due to the applied voltage.The SGLE describing the dynamics of the systematoms reads,¨ U ( t ) − F ( U ( t )) = − (cid:90) t Π r ( t − t (cid:48) ) U ( t (cid:48) ) dt (cid:48) + f ( t ) , (2) a r X i v : . [ c ond - m a t . m e s - h a ll ] F e b where, U is a vector composed of the mass-normalizeddisplacements of the system, and F ( U ( t )) is the force vec-tor from the potential of the isolated system. We adoptthe harmonic approximation, F ( U ( t )) = − KU ( t ), with K being the dynamical matrix. The effect of all bath de-grees of freedom is hidden in the terms on the right-handside of the SGLE. Each of them contains separate contri-butions from the L , R phonons, and the electron bath ( e ),such that Π r = Π rL + Π rR + Π re and f = f L + f R + f e . Thephonon self-energy Π r describes the time-delayed backac-tion of the bath on the system due to its motion[3, 26, 28–30]. The second quantum term f ( t ) is a random force(noise) due to the thermal, or current-induced fluctu- ation of the bath variables. It is characterized by thecorrelation matrix (cid:104) f α ( t ) f Tα ( t (cid:48) ) (cid:105) = S α ( t − t (cid:48) ), with α = L, R, e . The two phonon baths ( L and R ) are assumedto be in thermal equilibrium. Their noise correlation S L/R is related to the Π rL/R through the fluctuation-dissipation theorem, S L/R ( ω ) = ( n B ( ω, T ) + )Γ L/R ( ω )with Γ L/R ( ω ) = − rL/R ( ω ), n B the Bose distributionfunction (using atomic units, (cid:126) = 1). Due to the electri-cal current, the electronic bath is not in equilibrium. Wedefine the coupling-weighted electron-hole pair density ofstates as,[3, 26]Λ αβkl ( ω ) = 2 (cid:88) m,n (cid:104) ψ m | M k | ψ n (cid:105)(cid:104) ψ n | M l | ψ m (cid:105) ( n F ( ε n − µ α ) − n F ( ε m − µ β )) δ ( ε n − ε m − ω ) , (3)with n F the Fermi-Dirac distribution, and ψ n the elec-tron scattering state originating from the n -th channel ofelectrode α when there is no e-ph interaction. The noisecorrelation and the backaction term of the electron bathcan now be written as, S e ( ω ) = − π (cid:88) αβ (cid:20) n B ( ω − ( µ α − µ β )) + 12 (cid:21) Λ αβ ( ω ) , (4)Π re ( ω ) = −
12 ( H{ Γ e ( ω (cid:48) ) } ( ω ) + i Γ e ( ω )) , (5)Γ e ( ω ) = − π (cid:88) αβ Λ αβ ( ω ) , (6)where H{ A } is the Hilbert transform of A .In the absence of electrical current, the electrons serveas an equilibrium thermal bath, similar to phonons. How-ever, in the presence of current, the term ( ∼ ImΛ
RLkl , k (cid:54) = l ) becomes important. It may coherently couple two vi-brational modes ( kl ) inside the system leading to non-zero NC and BP forces. In Eq. (3) we observe thatthese effects depend on the phase of the electronic wave-function, and thus the direction of electronic current.Furthermore, the coherent coupling breaks time-reversalsymmetry of the noise correlation function, S e ( t − t (cid:48) ) (cid:54) = S e ( t (cid:48) − t ). Hereafter, we denote these forces by asym-metric CIF , and focus on their role for the excess heatdistribution and heat transport in the junction.We will consider the case where all baths are at thesame temperature ( T ), and the electron bath is subjectto a nonzero voltage bias ( eV = µ L − µ R ). To look at theexcess heating, we calculate the kinetic energy of atom n from its local displacement correlation function, and obtain E n = (cid:88) σ = x,y,z (cid:90) + ∞ ω diag { D r SD a } n,σ ( ω ) dω π . (7)Here D r ( D a ) is the eV -dependent phonon retarded (ad-vanced) Green’s function, S is the sum of noise correla-tion function from all the baths, and diag { A } n,σ meansthe diagonal matrix element of A , corresponding to the n -th atom’s σ degrees of freedom.To study heat transport, we calculate the phonon heatcurrent flowing into the bath L as the product of thevelocity of the system degrees of freedom, and the forceexerted on them by bath L . Applying time average, usingthe solution of the SGLE, we arrive at a Landauer-likeexpression (Sec. I, Supplemental Materials (SM)) J L = − (cid:90) + ∞−∞ ω tr (cid:2) Γ L ( ω ) D r ( ω )Λ RL ( ω ) D a ( ω ) (cid:3) × ( n B ( ω + eV ) − n B ( ω )) dω. (8)Defining the time-reversed phonon spectral function fromthe left bath ˜ A L = D a Γ L D r , and similarly A e = D r Λ RL D a , we can write the trace in Eq. (8) in differ-ent formstr[Γ L D r Λ RL D a ] = tr[Γ L A e ] = tr[Λ RL ˜ A L ] . (9)Equations (8) is analogous to the Landauer ornon-equilibrium Green’s function formula for elec-tron/phonon transport. In our present case the energycurrent is driven by a non-thermal electron bath withthe bias showing up in the Bose distributions and in thecoupling function, Λ RL , between phonons and electricalcurrent. The two forms in Eq. (9) emphasize two as-pects of the problem. In the first version emphasis is onthe coupling, Γ L of the system vibrations as describedby A e , to the phonons of the leads. This is a generalformula, which does not explicitly depend on the situa-tion we are considering here, namely that the source ofenergy is the non-equilibrium electron bath. This aspectis emphasized in the second version. Here the couplingto the electrical current, Λ RL is made explicit, and thecomplete phonon system including the coupling to leadsare in the function ˜ A L . In both forms the asymmetricCIF show up in the different versions of the A functions.The forces are responsible for the build up of vibrationalenergy inside the junction, a fact that is present in thetwo phonon Green’s functions D r and D a . Apart fromthis effect the non-equilibrium nature of the electron sys-tem shows up in the explicit factor Λ RL in the secondversion of Eq. (9). This will develop an imaginary partwhich is not present in equilibrium.Applying these formulas to a minimal model, in Sec. IIof the SM, we have shown analytically that the asymmet-ric CIF, especially the NC force, generate an asymmetricphonon heat flow and energy distribution, even for a left-right symmetric system. First Principles calculations–
Next we turn to numer-ical calculation for two concrete nano-junctions. Weuse SIESTA/TRANSIESTA[32, 33] to calculate the elec-tronic transport, vibrational modes, e-ph coupling em-ploying Ref. 34, and coupling to electrode phonons us-ing Ref. 35, with similar parameters. The effect of cur-rent on the stability of gold single atomic junctions hasbeen studied for more than a decade[31, 36]. Here wefirst consider a symmetric single atom gold chain be-tween two Au(100) electrodes(Fig. 2 inset).[37, 38] Wehave previously[39] studied the asymmetric forces in thissystem neglecting the coupling to electrode phonons.Figure 1 shows the average excess kinetic energy(∆ E n = E n ( eV ) − E n (0))[40–45] of atoms along the chainfor three different Fermi level E F . The structure is al-most mirror symmetric. When we turn off the asym-metric CIF (ImΛ RL = 0) as in previous studies[14, 46],the heating profile follows this symmetry. However, oncewe include them, the kinetic energy of one side becomesmany times higher than that of the other. Meanwhile,the total kinetic energy stored in the system increase sig-nificantly. Further analysis shows that both effects aredue to the NC force (Fig. 2 in SM).We now turn to the phonon heat current calculatedusing Eq. (8), shown in Fig 2 (a). The inclusion of theasymmetric CIF drives much larger heat current into the L bath. Intuitively, this is due to the asymmetric energyaccumulation induced by the NC force, e.g., modifying D r /D a in Eqs. (8-9). However, there is another contri-bution at low bias. Ignoring the bias-induced change of˜ A L , we get opposite heat flow into L and R ( J L = − J R )due to tr[ImΛ RL Im ˜ A L ]. This term drives asymmetricheat flow even in the linear response regime, contributingwith a correction to the thermoelectric Peltier coefficient(Sec. I(A) of SM). In the next section, we will show thatit can be understood as asymmetric excitation of left-and right-travelling phonon waves. (a) (c) (e) e ( ) ( m e V ) (b) t i c E n e r g y ( (d) (f) x c e ss K i n e t E FIG. 1: Excess kinetic energy of each atom in a gold chain(inset of Fig. 2(a)) at V = 1 . T = 300 K, with (bottom)and without (top) the asymmetric CIF. The total energy dif-ference between the two cases is due to the non-conservativeforce contribution. The blue dots and the colored plot of eachatom are from the full calculation. The asymmetric heatingis qualitatively reproduced by only considering electron cou-pling with vibrational modes (1) and (2) in the inset of Fig. 2(a), as shown by red triangles. (a)-(b) E F = − . E F = 0, and (e)-(f) E F = 0 . From Fig. 1 (b)-(d) and 2 (b), we see that the positionof E F is controlling the direction and magnitude of theasymmetry. According to the analysis in Sec. IV of SM,this could be due to the phase change of the electronicwavefunction with E F . Thus we expect that the direc-tion of electron flow is essential in the description of theatomic dynamics in the junction, as indicated in recentexperiments[8].The second system we consider is an armchairgraphene nanoribbon (a-GNR) with partial hydrogenpassivation, shown in Fig. 3 (a). This example is in-spired by experiments showing current-induced edge-reconstructions in graphene[47] where the physical mech-anism was attributed to Joule heating[48]. In Fig. 3 (a),the four pairs of unpassivated carbon dimers give rise tolocalized high-frequency vibrations interacting stronglywith electrical current. Consequently, the excess en-ergy is mainly stored in the dimers and nearby atoms(Fig. 3 (b),(d)), consistent with the experimental find-ings in Ref. 47. Including the asymmetric CIF leads tosymmetry breaking of the heating profile along the cur-rent direction. Contrary to experiments on the gold chain E F may in this case be tuned by gating. We predict theresulting hot-spot to move from “down-stream” to “up- e (a) e (1)(1)(2) (b) FIG. 2: (a) Bias dependence of the phonon heat current, go-ing into the left and right phonon baths. Solid lines includethe asymmetric CIF ( ∼ ImΛ RL ), dashed lines do not, andthe dash-dotted lines ignore the change of phonon spectral( D r /D a ) due to NC and BP forces. In the inset, we showthe two vibrational modes that couple most strongly with theelectrical current, with vibrational energy at (1) 19 and (2)18 meV. (b) Phonon heat current going into the left (red, cir-cle) and right (black, square) baths at V = 1V, for differentFermi levels to illustrate the importance of the phase of theelectron wavefunctions. stream” w.r.t. the electron current when tuning from E F = 1 . E F = − . E F can be under-stood as follows (Sec. III of SM). For a mirror-symmetricsystem with electron-hole symmetry, the asymmetricheating and heat flow is absent. When E F crosses theelectron-hole symmetric point, the dominant current-carriers contributing to inelastic transport change fromelectrons to holes, or vice versa. Thus, the hot-spot e (a) (b) (d) m e V ) c E n e r gy ( m (c) (e) ce ss K i n e ti ( ) ( ) E x c FIG. 3: (a) Structure of a partially passivated armchairgraphene ribbon considered. The two sides of the ribbon ishydrogen passivated except in the device region, enclosed bythe solid lines. (b)-(c) The excess kinetic energy of each atomwithout and with the asymmetric CIF, at V = 0 . T = 300K, E F = 1 . E F = − . moves from one side to the other. Interestingly enough,similar effect in micrometer scale has been observed ex-perimentally in graphene transistors[49, 50] and elec-trodes of molecular junctions[23]. Here we show that itis equally important at atomic scale, and related to theasymmetric CIF. Scattering analysis –
The asymmetric heating andphonon heat flow at low bias can be qualitatively un-derstood from the momentum transfer between electronsand phonons. To show this, we consider a simple 1Dmodel with a local e-ph interaction which involve the dis-placement of the n - and n + 1-th atoms (junction) (Sec.IV of SM), H eph = (cid:88) j ∈{ n,n +1 } − m ˆ u j ( c † j c j +1 − c † j c j − + h.c. ) . (10)For eV >
0, the important process is the inelastic elec-tronic transition from the filled, left scattering states withmomentum k L to the empty, right states with k R . It isstraightforward to show that the emission probability ofa right-travelling phonon with momentum q is differentfrom that of a left-travelling mode, − q , due to the differ-ence in matrix elements for the processes,∆ M LR = | M qLR | − | M − qLR | ∼ sin( q ) sin( k L − k R ) . (11)Consequently, the left- and right-travelling steady statephonon populations become different, resulting in asym-metric heat flow.In conclusion, we have presented a theory showing thatCIF in nano-junctions lead to asymmetric distributionsand transport of the excess heat. We derived a Landauer-like formula for the excess heat transport. Employingfirst-principles calculations, we demonstrate that the size of the asymmetry can be crucial for current-induced pro-cesses at the atomic scale.We thank T. N. Todorov, D. Dundas, and T.Markussen for discussions and the Danish Center forScientific Computing (DCSC) for computer resources.This work is supported by the Lundbeck Foundation(R49-A5454), National Natural Science Foundation ofChina (Grants No. 11304107, 61371015), and the Fun-damental Research Funds for the Central Universities(HUST:2013TS032). ∗ Electronic address: [email protected][1] We use phonons and vibrations interchangably, although,strictly speaking, phonons are defined only in systemswith translational invariance.[2] D. Dundas, E. J. McEniry, and T. N. Todorov, NatureNanotech. , 99 (2009).[3] J. T. L¨u, M. Brandbyge, and P. Hedeg˚ard, Nano Lett. , 1657 (2010).[4] N. Bode, S. V. Kusminskiy, R. Egger, and F. von Oppen,Phys. Rev. Lett. , 036804 (2011).[5] T. N. Todorov, D. Dundas, A. T. Paxton, and A. P. Hors-field, Beilstein Journal of Nanotechnology , 727 (2011).[6] I. A. Pshenichnyuk and M. ˇC´ıˇzek, Phys. Rev. B ,165446 (2011).[7] R. Bustos-Mar´un, G. Refael, and F. von Oppen, Phys.Rev. Lett. , 060802 (2013).[8] C. Schirm, M. Matt, F. Pauly, J. C. Cuevas, P. Nielaba,and E. Scheer, Nature Nanotech. , 645 (2013).[9] P. J. Wheeler, R. Chen, and D. Natelson, Phys. Rev. B , 155411 (2013).[10] B. Cunningham, T. N. Todorov, and D. Dundas, Phys.Rev. B , 115430 (2014).[11] N. J. Tao, Nature Nanotech. , 173 (2006).[12] M. Galperin, M. A. Ratner, and A. Nitzan, J.Phys.:Condens. Matter , 103201 (2007).[13] M. Galperin, M. A. Ratner, A. Nitzan, and A. Troisi,Science , 1056 (2008).[14] Z. Huang, F. Chen, R. D’Agosta, P. A. Bennett,M. Di Ventra, and N. Tao, Nature Nanotech. , 698(2007).[15] M. Tsutsui, M. Taniguchi, and T. Kawai, Nano Lett. ,3293 (2008).[16] Y. Asai, Phys. Rev. B , 045434 (2008).[17] R. H. M. Smit, Y. Noat, C. Untiedt, N. D. Lang, M. C.van Hemert, and J. M. van Ruitenbeek, Nature , 906(2002).[18] W. Y. Wang, T. Lee, I. Kretzschmar, and M. A. Reed,Nano Lett. , 643 (2004).[19] J. G. Kushmerick, J. Lazorcik, C. H. Patterson,R. Shashidhar, D. S. Seferos, and G. C. Bazan, NanoLett. , 639 (2004).[20] Z. Ioffe, T. Shamai, A. Ophir, G. Noy, I. Yutsis, K. Kfir,O. Cheshnovsky, and Y. Selzer, Nature Nanotech. , 727(2008).[21] D. R. Ward, D. A. Corley, J. M. Tour, and D. Natelson,Nature Nanotech. , 33 (2011).[22] K. Kaasbjerg, T. c. v. Novotn´y, and A. Nitzan, Phys. Rev. B , 201405 (2013).[23] W. Lee, K. Kim, W. Jeong, L. A. Zotti, F. Pauly, J. C.Cuevas, and P. Reddy, Nature , 209 (2013).[24] R. H. M. Smit, C. Untiedt, and J. M. van Ruitenbeek,Nanotechnology , S472 (2004).[25] G. Schulze, K. J. Franke, A. Gagliardi, G. Romano,C. S. Lin, A. L. Rosa, T. A. Niehaus, T. Frauenheim,A. Di Carlo, A. Pecchia, et al., Phys. Rev. Lett. ,136801 (2008).[26] J.-T. L¨u, M. Brandbyge, P. Hedeg˚ard, T. N. Todorov,and D. Dundas, Phys. Rev. B , 245444 (2012).[27] J.-S. Wang, Phys. Rev. Lett. , 160601 (2007).[28] R. P. Feynman and F. L. Vernon, Ann. Phys. , 118(1963).[29] A. Caldeira and A. Leggett, Physica A , 587 (1983).[30] A. Schmid, J. Low Temp. Phys. , 609 (1982).[31] Y. Oshima and Y. Kurui, Phys. Rev. B , 081404(2013).[32] J. Soler, E. Artacho, J. Gale, A. Garcia, J. Junquera,P. Ordejon, and D. Sanchez-Portal, J. Phys.:Condens.Matter , 2745 (2002).[33] M. Brandbyge, J. L. Mozos, P. Ordejon, J. Taylor, andK. Stokbro, Phys. Rev. B , 165401 (2002).[34] T. Frederiksen, M. Paulsson, M. Brandbyge, and A.-P.Jauho, Phys. Rev. B , 205413 (2007).[35] M. Engelund, M. Brandbyge, and A. P. Jauho, Phys.Rev. B , 045427 (pages 11) (2009).[36] H. Yasuda and A. Sakai, Phys. Rev. B , 1069 (1997).[37] H. Ohnishi, Y. Kondo, and K. Takayanagi, Nature ,780 (1998).[38] A. I. Yanson, G. R. Bollinger, H. E. van den Brom,N. Agra¨ıt, and J. M. van Ruitenbeek, Nature , 783(1998).[39] J. T. L¨u, P. Hedeg˚ard, and M. Brandbyge, Phys. Rev.Lett. , 046801 (2011).[40] Y. Dubi and M. Di Ventra, Phys. Rev. B , 115415(2009).[41] Y. Dubi and M. Di Ventra, Phys. Rev. E , 042101(2009).[42] P. A. Jacquet, J Stat Phys , 709 (2009).[43] P. A. Jacquet and C.-A. Pillet, Phys. Rev. B , 125120(2012).[44] J. P. Bergfield, S. M. Story, R. C. Stafford, and C. A.Stafford, ACS Nano , 4429 (2013).[45] Another way of quantifying the heating is to use the localtemperature defined in some way. We tried to use themethod in Refs. [40-44]. The result is shown in Fig. 3 of the SM. The overall heating profile agrees with Fig. 1.[46] T. Frederiksen, M. Brandbyge, N. Lorente, and A.-P.Jauho, Phys. Rev. Lett. , 256601 (2004).[47] X. Jia, M. Hofmann, V. Meunier, B. G. Sumpter,J. Campos-Delgado, J. M. Romo-Herrera, H. Son, Y.-P. Hsieh, A. Reina, J. Kong, et al., Science , 1701(2009).[48] M. Engelund, J. A. F¨urst, A. P. Jauho, and M. Brand-byge, Phys. Rev. Lett. , 036807 (2010).[49] M. Freitag, H.-Y. Chiu, M. Steiner, V. Perebeinos,and P. Avouris, Nature Nanotechnology , 497 (2010),1004.0369.[50] M.-H. Bae, Z.-Y. Ong, D. Estrada, and E. Pop, NanoLetters , 4787 (2010), 1004.0287. SUPPLEMENTAL MATERIALSI. DERIVATION OF THE PHONON HEAT CURRENT EQ. (8)
We start from the semi-classical generalized Langevin equation (SGLE) (Eq. (2) in the main text). To study theenergy transport, we look at the energy increase of the system per unit time˙ E S ( t ) = ddt (cid:18)
12 ˙ U T ˙ U + 12 U T KU (cid:19) = − ˙ U T (cid:32)(cid:88) α (cid:90) t −∞ Π rα ( t − t (cid:48) ) U ( t (cid:48) ) dt (cid:48) − f α ( t ) (cid:33) , α = L, R, e. (12)Note that the system includes only the atomic degrees of freedom. We can define the energy current flowing into thebath α from the system J α ( t ) ≡ ˙ U T (cid:18)(cid:90) t −∞ Π rα ( t − t (cid:48) ) U ( t (cid:48) ) dt (cid:48) − f α ( t ) (cid:19) . (13)At steady state we have − ˙ E s ≡ J e + J L + J R ≡ J e + J ph = 0 . (14)We can write the expression for the average energy current in the frequency domain, J α ≡ lim T → + ∞ T (cid:90) T (cid:28) ˙ U T ( t ) (cid:18)(cid:90) t Π rα ( t − t (cid:48) ) U ( t (cid:48) ) dt (cid:48) − f α (cid:19)(cid:29) dt = lim T → + ∞ T (cid:90) dω π (cid:104) ˙ U † ( ω ) (Π rα ( ω ) U ( ω ) − f α ( ω )) (cid:105) (15)Now we use the solution of the Langevin equation U ( ω ) = − D r ( ω ) f ( ω ) , (16) D r ( ω ) = (cid:2) ω − K − Π r ( ω ) (cid:3) − , (17)Π r ( ω ) = Π L ( ω ) + Π R ( ω ) + Π e ( ω ) , (18)and the noise correlation function (cid:104) f α ( ω ) f α ( ω (cid:48) ) (cid:105) = δ ( ω + ω (cid:48) ) S α ( ω ) , (19) S ( ω ) = S L ( ω ) + S R ( ω ) + S e ( ω ) , (20)to get ( (cid:126) = 1) J α = i (cid:90) + ∞−∞ dω π ω Tr [Π rα ( ω ) D r ( ω ) S ( ω ) D a ( ω ) + S α ( ω ) D a ( ω )] (21)= i (cid:90) + ∞ dω π ω Tr [Π rα ( ω ) D r ( ω ) S ( ω ) D a ( ω ) + S α ( ω ) D a ( ω ) (22) − Π rα ( − ω ) D r ( − ω ) S ( − ω ) D a ( − ω ) − S α ( − ω ) D a ( − ω )] . (23)The two phonon baths ( L and R ) are assumed to be in thermal equilibrium. Their noise correlation S L/R isrelated to the Π rL/R through the fluctuation-dissipation theorem, S L/R ( ω ) = ( n B ( ω, T ) + )Γ L/R ( ω ) with Γ L/R ( ω ) = − rL/R ( ω ), n B the Bose distribution function (using atomic units, (cid:126) = 1). The noise correlation of the electronbath is given by Eqs. (3)-(4) in the main text. Using the following properties( D r ) † ( ω ) = D a ( ω ) , (Π r ) † ( ω ) = Π a ( ω ) , Γ( ω ) = i (Π r ( ω ) − Π a ( ω )) , (24) S † ( ω ) = S ( ω ) , S ( − ω ) = S ∗ ( ω ) , D r ( − ω ) = ( D r ) ∗ ( ω ) , D a ( − ω ) = ( D a ) ∗ ( ω ) , (25)and taking transpose of Eq. (23), we get a compact form J α = (cid:90) + ∞ dω π ω Tr [Γ α ( ω ) D r ( ω ) S ¯ α ( ω ) D a ( ω ) − S α ( ω ) D r ( ω )Γ ¯ α ( ω ) D a ( ω )] . (26)This result has a clear physical meaning. Here, Γ α characterizes coupling of the α bath to the system, and S ¯ α represents the energy source from all other baths. The first term in the trace represents energy flow into bath α fromother baths; while the second one represents the opposite process. A. Current-induced phonon heat transport
Now suppose all the baths are at the same temperature ( T ), but the electron bath is subject to a nonzero bias( eV ). The energy current injecting into the phonon bath ( L ) is J L = (cid:90) + ∞ dω π ω Tr [Γ L ( ω ) D r ( ω ) S ¯ α ( ω ) D a ( ω ) − S L ( ω ) D r ( ω )Γ ¯ α ( ω ) D a ( ω )] (27)= (cid:90) + ∞ dω π ω Tr [Γ L ( ω ) D r ( ω ) S e ( ω ) D a ( ω ) − S L ( ω ) D r ( ω )Γ e ( ω ) D a ( ω )] . (28)To go from Eq. (27) to (28), we notice that the energy flow from L to R is the same as that from R to L , since theyare at the same temperature. Thus, the only energy source is the electron bath. Using Eqs. (3-6) in the main text,the heat current now reads J L = − (cid:88) α (cid:54) = β (cid:90) + ∞ dωω Tr (cid:2) Λ αβ ( ω ) D a ( ω )Γ L ( ω ) D r ( ω ) (cid:3) ( n B ( ω − ( µ α − µ β )) − n B ( ω )) (29)= − (cid:90) + ∞ dωω Tr (cid:2) Λ RL ( ω ) D a ( ω )Γ L ( ω ) D r ( ω ) (cid:3) ( n B ( ω + eV ) − n B ( ω )) . (30)Define the time-reversed phonon spectral function ˜ A L ( ω ) = D a ( ω )Γ L ( ω ) D r ( ω ), we can write it in other equivalentforms J L = − (cid:90) + ∞−∞ dωω Tr (cid:104) Λ LR ( ω ) ˜ A L ( ω ) (cid:105) ( n B ( ω − eV ) − n B ( ω )) (31)= − (cid:90) + ∞−∞ dωω Tr (cid:104) Λ RL ( ω ) ˜ A L ( ω ) (cid:105) ( n B ( ω + eV ) − n B ( ω )) . (32)Similar equation holds for J R J R = − (cid:90) + ∞−∞ dω ω Tr (cid:104) Λ RL ( ω ) ˜ A R ( ω ) (cid:105) ( n B ( ω + eV ) − n B ( ω )) (33)= − (cid:90) + ∞−∞ dω ω Tr (cid:104) Λ LR ( ω ) ˜ A R ( ω ) (cid:105) ( n B ( ω − eV ) − n B ( ω )) . (34)Let’s look at the low bias situation. We ignore the change of ˜ A L/R , and replace it with ˜ A L , the counterpart of ˜ A L without coupling to electrons. The asymmetric current-induced forces ( ∼ ImΛ RL ) drive a heat current J L,p = 12 (cid:90) + ∞−∞ dω ω Tr (cid:104) ImΛ RL ( ω )Im ˜ A L ( ω ) (cid:105) (cid:18) coth (cid:18) ω + eV k B T (cid:19) − coth (cid:18) ω k B T (cid:19)(cid:19) , (35)The expression for J R,p is obtained by replacing ˜ A L with ˜ A R . From Im ˜ A L + Im ˜ A R = 0, we get J L,p = − J R,p . Thatis, the heat flowing into bath L and R is opposite. This makes J L (cid:54) = J R , even for a symmetric structure. Furthermore,in the linear response regime, considering thermoelectric transport, from Eq. (35) we get a correction to the Peltiercoefficient due to electron-phonon interaction: The applied bias drives a phonon heat current from one phonon bathto the other.From the derivation of Eq. (35), and (21), we observe that, the first term with ω coth(( ω + eV ) / (2 k B T )) in J L,p is contributed by the fluctuating force in the SGLE, while the second term with ω coth( ω/ (2 k B T )) is from thedeterministic NC force. If the bias | eV | is much higher than the phonon frequency, the contribution from NC forcedominates. This can be seen from the symmetry of the functions, as follows: for high enough bias, ω coth(( ω + eV ) / (2 k B T )) is close to be odd in ω , e.g., ignoring ω in the coth function. But ω coth( ω/ (2 k B T )) is even in ω .Meanwhile, the trace in Eq. (35) can be approximated by an even function for small ω . Thus, the contribution of theNC force dominates. The above analysis based on Eq. (35) is correct to the 2nd order in M . Going beyond the 2ndorder, we notice that in Eqs. (31-34), the deterministic NC and BP force modifies the phonon spectral function ˜ A L/R ,while the fluctuating force has no effect on it. Altogether, we conclude that, the asymmetric noise has a negligiblecontribution to the asymmetric heat flow.
II. MINIMAL MODEL
We now consider a minimal model with two atomic vibrations. In additional to electrons, they couple symmetricallyto the left and right phonon bath, respectively. This gives rise to lifetime broadening of γ e and γ ph , respectively. Thephonon Green’s function is written as D r ( ω ) = 1 N (cid:18) Ω − ω − a − ibω − ω + a + ibω Ω (cid:19) . (36)Here, Ω = ω − ω + iγ t ω , γ t = γ e + γ ph , N = Ω − ( ω + a + ibω )( ω − a − ibω ), a and b are due to NC and BPforces, respectively. Finally, ω is the atomic vibration frequency, and ω characterizes the coupling between the twosites. We have ignored a term ∼ − iγ (cid:48) e ω in the off-diagonals of D r ( ω ). The advanced Green’s function is D a = ( D r ) † .We also have ˜Π L ( ω ) = − iγ ph ω (cid:18) (cid:19) . (37)From these, we get the time-reversed phonon left spectral function˜ A L ( ω ) = 2 ωγ ph | N | (cid:18) | Ω | − ( a + ibω + ω )Ω ∗ − ( a − ibω + ω )Ω ω + a + b ω + 2 aω (cid:19) , (38) A. Heat current
To calculate the heat current, we assumeΛ αβ ( ω ) ≈ − ω − ( µ α − µ β )) (cid:18) λ αβ λ αβ + iλ αβ λ αβ − iλ αβ λ αβ (cid:19) . (39)This means we ignore the energy dependence of the electronic properties within the bias window. We can now evaluatethe trace in Eq. (32),Tr (cid:104) ImΛ RL Im ˜ A L (cid:105) = − | N | ω ( ω + eV ) γ ph (cid:20) ω c eV λ ( ω − ω ) + eV λ γ t − λ ω γ t (cid:124) (cid:123)(cid:122) (cid:125)(cid:21) , (40)Tr (cid:104) ReΛ RL Re ˜ A L (cid:105) = − | N | ω ( ω + eV ) γ ph (cid:20) | Ω | λ + λ (cid:18) ω + eV λ (cid:18) ω ω c (cid:19)(cid:19) − λ ω ( ω − ω ) + 2 eV λ (cid:20) λ (cid:18) ω ω c γ t + ( ω − ω ) (cid:19) − λ ω (cid:21)(cid:124) (cid:123)(cid:122) (cid:125) . (41)We have drooped the RL superscript in λ s for notational simplicity, and used a = − eV λ RL , b = a/ω c . Here ω c is onthe order of the electron bandwidth. Substituting back into Eq. (32), we find that those terms in the curly bracketsof Eqs. (40) and (41), due to the asymmetric current-induced forces ( ∼ λ ), induce asymmetric heat flow (odd in eV )to the left and right phonon bath.0 B. Average kinetic energy
We now calculate the average kinetic energy difference between the two atomic sites. If we take a general noisecorrelation for the electronic bath S e ( ω, eV ) = (cid:88) α,β = { L,R } g αβ ( ω ) (cid:18) λ αβ λ αβ + iλ αβ λ αβ − iλ αβ λ αβ (cid:19) , (42)with g αβ ( ω ) = 2 π ( ω − ( µ α − µ β )) coth (cid:18) ω − ( µ α − µ β )2 k B T (cid:19) . (43)Using Eq. (6) in the main text, we get∆ E e = (cid:88) α,β = L,R (cid:90) ω g αβ ( ω ) | N | (cid:104) λ αβ ( a ( ω − ω ) + bω γ t ) − aλ αβ ω − λ αβ γ t ωω (cid:105) dω π . (44)We look at the nonequilibrium contribution first S none ( ω ) = S e ( ω, eV ) − S e ( ω, eV (cid:29) ω , similar arguments toSec. I A show that the main contribution comes from the real part of the two terms with α (cid:54) = β , and the asymmetricnoise is negligible, S none ( ω ) ∼ π | eV | (cid:18) λ RL λ RL λ RL λ RL (cid:19) , eV (cid:29) ω . (45)Consequently, we get ∆ E non ≈ (cid:90) πω eV λ RL | eV || N | (cid:20) λ RL (cid:18) ω − ω + ω γ t ω c (cid:19) − λ RL ω (cid:21) dω π . (46)The BP force contribution is negligible if ω c is the largest energy scale of the problem. If we further ignore λ RL to beconsistent with Eq. (36), we get ∆ E non ≈ − πeV | eV | λ RL λ RL ω (cid:90) ω | N | dω π . (47)For the equilibrium part, including contribution from phonon baths, we get∆ E equ ≈ − eV λ RL ( γ ph + γ e ) ω (cid:90) ω | N | coth (cid:18) ω k B T (cid:19) dω π , (48)The total difference ∆ E e = ∆ E non + ∆ E equ . We see that ∆ E e = 0 if λ RL = 0. Thus, the asymmetric current-induced forces generate asymmetric energy distribution, with the NC force contributes predominantly. The asymmetryis enhanced by coupling to phonon baths ( γ ph in Eq. (48)). III. ELECTRON-HOLE SYMMETRY
Assuming symmetrical voltage drop across the conductor, we define the zero energy as the equilibrium Fermi level.The left and right chemical potential are at eV / − eV /
2, respectively. The Λ-function now readsΛ αβkl ( ω, eV ) = 2 (cid:88) m,n (cid:104) ψ m | M k | ψ n (cid:105)(cid:104) ψ n | M l | ψ m (cid:105) [ n F ( ε n − s α eV / − n F ( ε m − s β eV / δ ( ε n − ε m − ω ) , (49)where we have written explicitly its eV dependence, and s L = 1, s R = −
1. It has the following properties:Λ αβkl ( ω, eV ) = Λ αβlk ∗ ( ω, eV ) , Λ αβkl ( ω, eV ) = − Λ βαlk ( − ω, eV ) . (50)For the convenience of further analysis, we now use A α ( ε ) = 2 π (cid:88) n | ψ n (cid:105) δ ( ε − ε n ) (cid:104) ψ n | , (51)1to write it asΛ αβkl ( ω, eV ) = 2 (cid:90) dε π (cid:90) dε (cid:48) π Tr (cid:2) M k A α ( ε ) M l A β ( ε (cid:48) ) (cid:3) [ n F ( ε − s α eV / − n F ( ε (cid:48) − s β eV / δ ( ε − ε (cid:48) − ω ) , (52)Note that the spectral function A α ( ε ) is Hermitian. If we use a real-space basis set, it is a complex matrix. We definethe system has electron-hole symmetry ifRe A α ( ε ) = Re A α ( − ε ) , or A α ( ε ) = A ∗ α ( − ε ) , (53)The two conditions are equivalent since Re A α ( ε ) is related to Im A α ( ε ) by Hilbert transform, which changes theirsymmetry with respect to ε . Using Eq. (53) in (52), together with Eq. (50), we find thatΛ LR ( ω, eV ) = Λ RL ( ω, − eV ) . (54)Here, we have further assumed that the electron-phonon interaction matrix is real. This is a reasonable assumption, ifwe ignore the bias-dependence of the electronic Hamiltonian, and consider Cartesian phonon index, without externalmagnetic field. Substituting it into Eqs. (31-32) and (33-34), we find that J α ( eV ) = J α ( − eV ) , α = L, R. (55)We reach the conclusion that the heat flow into L and R are the same if the system has electron-hole symmetry andthere is a symmetrical voltage drop across the conductor. IV. SCATTERING ANALYSIS
The asymmetric heating and heat flow at low bias can be qualitatively understood as momentum transfer betweenelectrons and phonons. To show this, we consider a simple one-dimensional (1D) model. The electronic subsystemis described by a nearest neighbour tight-binding Hamiltonian, and the phonon subsystem by a harmonic oscillatormodel. To simplify the analysis, we assume that the tight-binding hopping parameter and the spring constant betweenall the nearest sites are the same. But the analysis can be easily extended to more general case, where our conclusionin this section still holds. The electron and phonon states are described by scattering waves originating from L and R . We introduce a local e-ph interaction on two atomic sites n and n + 1 (junction), that is, the displacement of the n - and n + 1-th atoms modifies the electronic hopping elements nearby linearly, e.g., for phonon mode q , M q ∼ − e iq − e iq − e iq − e iq . (56)For positive bias eV >
0, the main process contributing to phonon emission is the inelastic electronic transitionfrom the filled, left scattering states ψ L to the empty, right states ψ R . The transition rate is proportional to themodulus square of the matrix element, | M qLR | = | M − qRL | ∼ cos
12 ( q − k L + k R ) . (57)The emission probability of a right-travelling phonon mode q ( q >
0) is different from that of a left-travelling mode, − q . The difference is ∆ M LR = | M qLR | − | M − qLR | ∼ sin( q ) sin( k L − k R ) , (58)and as a result, the left and right-travelling steady state phonon populations becomes different. The difference changessign upon changing the current direction which reveal the importance of electron momentum.Next, we use the retarded phonon Green’s function to consider the response, | r (cid:105) , of the phonon system to theasymmetric excitation, | s (cid:105) ∼ ( · · · e iq · · · ) T , at n and n + 1, and find |(cid:104) m | r (cid:105)| ∼ |(cid:104) m | D r | s (cid:105)| ∼ (cid:40) cos q + | q | , m (cid:28) n cos q −| q | , m (cid:29) n (59)2 b iqa be iqa a (a) n n+ be − iqa a (b)(b) n n+ FIG. 4: The motion of atom n and n + 1 has a phase shift of ± q . Within the configuration space of u n and u n +1 , the twosituations correspond to elliptical motion in opposite directions. The excitation probabilities of these two modes differs. Thisresults in (1) an asymmetric heat flow to the left and right phonon bath, (2) polarization of the motion within configurationspace ( u n , u n +1 ). where obviously the response differs at the left and right side of the perturbation (Fig. 4).We conclude that the applied bias breaks the population balance between left and right electron scattering states.Consequently, electrons excite the left and right travelling phonon states differently resulting in transfer of both energyand momentum to the phonons. The momentum transfer generates a different phonon energy flux to the left andright for the spatially symmetric system under bias. A schematic diagram of these processes are shown in Fig. 4. Ifwe turn on the e-ph interaction at all the sites, the interaction matrix becomes M qLR ∼ δ ( k L − k R − q − N π ). Theasymmetric phonon excitation reduces to the rule of crystal momentum conservation in the periodic structure.To make connection with the current-induced NC and BP force, in Fig. 4 we illustrate the orbital of the two phononexcitation within the configuration space of ( u n , u n +1 ). They are elliptical and related by time-reversal. From thispoint of view, the current-induced NC and BP forces polarize the atomic orbital motion, and generate a net angularmomentum. The heat flow into the two electrodes becomes different due to this elliptical polarization.Finally, it is instructive to compare the scattering analysis against the Langevin approach. In fact, one can showthat Im (cid:104) ψ L | M n | ψ R (cid:105)(cid:104) ψ R | M n +1 | ψ L (cid:105) ∼ sin( k L − k R ) , (60)and (Im ˜ A L ) n +1 ,n ( ω q ) ∼ sin( q ) . (61)So, comparing Eqs. (35) and (58), we can see that the asymmetric heat flow can indeed be understood as a result ofasymmetric excitation of left- and right-travelling phonon waves. V. SUPPORTING FIGURES FROM THE FIRST-PRINCIPLES CALCULATION (a) V ) (c) (e) e E n e r g y ( m e (b) ss K i n e t i c E (d) (f) E x e s FIG. 5: Excess kinetic energy averaged over atoms in each column at V = 1 . T = 300 K, similar to Fig. 2 in the maintext. The top part (a), (c), (e) shows results without the asymmetric current-induced forces, whiled the bottom part (b), (d),(f) shows results that include only the BP and asymmetric fluctuating force. This shows the contribution from the BP andfluctuating force is negligible. Buttiker ProbeButtiker Probe (b’) (d’) (f’) K ) T ( x FIG. 6: Another way of characterizing heating in the chain is to use the B¨uttiker probe (Refs. [40-44] of the main text) to‘measure’ the temperature of each atom. (b’), (d’) and (f’) show the ‘measured’ temperature of each atom using this methodwhen including all the forces. The overall heating profile agree with Fig. 2(b), (d), (f) in the main text. e (a) (b) (d) y ( m e V ) n e ti c E n e r gy ( ) ( ) E x ce ss K i n (c) (e) E FIG. 7: Excess kinetic energy averaged over atoms in each column at V = 0 . T = 300 K, similar to Fig. 3 in the main text.The top rows are results without the asymmetric current-induced forces, while the results in the bottom row include the BPand asymmetric fluctuating force. Again, their contribution to the asymmetric heating is negligible. -2 -1 0 1 200.511.522.5 T e Energy (eV)(a) T p h Energy (eV)(b) -2 -1 0 1 201020 D O S e Energy (eV)(c) D O S p h Energy (eV)(d)
FIG. 8: Additional information on the Graphene nanoribbon calculation. (a) Electronic transmission ( T e ). (b) Phononictransmission ( T ph ). (c) Electronic density of states ( DOS e ) projected to the device region. (d) Phononic density of states( DOS phph