Effects of electron-phonon interaction on thermal and electrical transport through molecular nano-conductors

The topic of this review is the effects of electron-phonon interaction (EPI) on the transport properties of molecular nano-conductors. A nano-conductor connects to two electron leads and two phonon leads, possibly at different temperatures or chemical potentials. The EPI appears only in the nano-conductor. We focus on its effects on charge and energy transport. We introduce three approaches. For weak EPI, we use the nonequilibrium Green's function method to treat it perturbatively. We derive the expressions for the charge and heat currents. For weak system-lead couplings, we use the quantum master equation approach. In both cases, we use a simple single level model to study the effects of EPI on the system's thermoelectric transport properties. It is also interesting to look at the effect of currents on the dynamics of the phonon system. For this, we derive a semi-classical generalized Langevin equation to describe the nano-conductor's atomic dynamics, taking the nonequilibrium electron system, as well as the rest of the atomic degrees of freedom as effective baths. We show simple applications of this approach to the problem of energy transfer between electrons and phonons.


I. INTRODUCTION
Electron-phonon interaction (EPI) is one of the most important many-body interactions in condensed-matter and molecular systems 1 , responsible for a variety of phenomena, from electrical, thermal conduction, superconductivity to Raman scattering, polaron formation, just to list a few [2][3][4][5][6][7][8] . Its effects on the electrical, thermal, and optical properties of bulk semiconductors and metals have been intensively studied along with the development of many-body theories and experimental techniques. Recent advances in experimental fabrication of meso-and nano-scopic structures have generated tremendous efforts in understanding the effects of EPI on transport properties of reduced-dimensional systems [9][10][11] .
Of special interest are current-induced forces and Joule heating in low-dimensional systems, especially in molecular nano-conductors  . On the one hand, the electrical transport signature of EPI is an invaluable spectroscopic tool to study the structural information of molecular nano-conductors 22,24 . On the other hand, these processes are crucial in maintaining the stability of these conductors 25 , relevant to the continuous scaling down of modern electronic devices. Different theoretical approaches have been developed to study these problems, in many cases separately. Recently, it was realized that non-conservative nature of current-induced forces provides an alternative, deterministic way of energy transfer between electrons and phonons, or more generally atomic motions 15 . It is fundamentally different from the stochastic Joule heating. These advances have motivated the development of methods treating current-induced forces and Joule heating on the same footing [17][18][19]34 .
Equally significantly, there has been an increasing interest in the thermoelectric properties of low dimensional systems [50][51][52][53][54] . A starting point of the theoretical treatment is to ignore the effect of EPI, and study the transport of electrons and phonons separately. But how important the effect of EPI is is a pertinent question, on which much of recent work is devoted to [55][56][57][58][59] . Here, we will look at this problem using the various approaches we have developed.
EPI is a genuine many-body interaction, the exact treatment of which is challenging, if possible at all. One natural approach is to perform perturbation calculation over a certain small parameter. In the most common multi-probe transport setup (see Fig. 1 and Sec. II), this small parameter can be chosen according to the strength of EPI. This strength can be roughly characterized by the ratio between two time scales: the first one corresponds to the phonon period, and the second one corresponds to the electron dwell time 60 in the nano-conductor. If the time electrons spend in the nano-conductor is much shorter than the phonon period, the system is in the weak EPI regime. The small parameter is the EPI matrix. In the other limit, the coupling of the nano-conductor to electrodes is the small parameter, over which one can perform the perturbation expansion. arXiv:1501.06343v2 [cond-mat.mes-hall] 24 Mar 2015 In this review, we summarize our own effort in developing and/or utilizing different theoretical approaches to study the aforementioned problems in different parameter regimes. We discuss some relevant results when possible, but we make no effort on reviewing all of them considering the huge amount of literature. The paper is organized as follows: In Sec. II, we give a brief introduction of the EPI problem starting from the Born-Oppenheimer approximation. We then introduce the system setup and Hamiltonian we use in this paper. In Sec. III, we briefly summarize our use of the nonequilibrium Green's function (NEGF) method to study electron, phonon transport and their interaction perturbatively. We consider several applications of the method. The first one is the effects of EPI on the thermoelectric transport coefficients in a single level model. The second one is the heat transport between electrons and phonons due to EPI. The use of simple models enables us to approach the problems semi-analytically. The last example is a numerical study of the Joule heating and phonondrag effect in carbon nanotubes. In Sec. IV, we consider the case of strong EPI using the quantum master equation (QME) approach. After reviewing the earlier work, the same thermoelectric transport model is re-visited focusing on how the strength of EPI affects the results. In Sec. V, we focus on the current-induced dynamics. Based on the Feynman-Vernon influence functional approach, we derive the semi-classical Langevin equation, taking into account the equilibrium phonon and nonequilibrium electron baths. The final section is our conclusion and remarks.

II. BORN-OPPENHEIMER APPROXIMATION AND ELECTRON-PHONON INTERACTIONS
To discuss the meaning and formulation of the electron and phonon systems and their mutual interaction, we need to start from the Born-Oppenheimer approximation 2,3 . Consider an electron-ion system with a total Hamiltonian H = P i + H e , where P i is the kinetic energy operator for the ions, and H e = P e + U is electron Hamiltonian with kinetic energy of the electrons, P e , and potential energy U = U ee + U ei + U ii , which includes the Coulomb interactions among the electrons and ions. Since the ions are much heavier than the electrons, one can treat the ion kinetic energy term as a small perturbation with the expansion parameter 2 where m e is the mass of an electron and m p mass of an ion (assuming all have the same mass). If the ions are considered infinitely heavy, the ions will not move and the electron wavefunctions satisfy where x represents the set of all coordinates of the electrons, R the positions of all the ions, and α is the electronic state quantum number. The eigen-functions and the eigenvalues depend on R parametrically. We assume an orthonormal set {φ α } that satisfies Eq. (2) has been obtained. To take into account the effects of the ions, we consider a trial full wavefunction in a factored form and consider the variational solution 4 of the full Hamiltonian, min χ Ψ|H|Ψ , subject to the normalization χ|χ = 1. This variational approach is equivalent to omitting the off-diagonal elements (which is the Born-Oppenheimer approximation, see Ref. 2, App. VIII), giving an equation for the ions where · · · means the x-dependence is integrated out but still R-dependent; ∇ R is a multi-dimensional gradient operator with respect to R. Since the left-hand side depends on the electronic quantum number α, the full eigen-energy E and functions also depend on α parametrically, e.g., we may write E β;α . If we assume that the electrons are in its instantaneous ground state, the ions move in a potential surface generated by the electrons. There are no explicit electronphonon interaction (EPI) terms. To account for the EPI, we need to go back to the basis, Eq. (3), and consider the matrix elements αβ|H|α β .
The off-diagonal terms are interpreted as the EPI 5,6 , which are small. If the off-diagonals are omitted, the electrons stay in a given quantum state α. The off-diagonal terms describe the scattering of the electrons to different state α . If ion displacements are small, the most important contribution is from the linear term in the displacement (6) These off-diagonal matrix elements can be used, e.g., in a Fermi-Golden rule calculation of scattering processes. However, the identity (in the sense of effective Hamiltonians) of the electrons and phonons and their mutual interaction are not at all clear. Although EPI plays major role in many physical processes 7 , such as electronic transport and superconductivity, its conceptual foundation is still not very solid. Within the Born-Oppenheimer scheme, it is not clear at all how to transform the original Hamiltonian H into a form of an electron system and independent phonon system and their interaction unambiguously. The problem is related to the fact that in deriving the phonon Hamiltonian (the potential surfaces), the effect of electrons is already used. Thus, putting the electrons back amounts to double counting, see Refs. 4 and 8 for some of the modern treatments.
Instead of pursuing a self-consistent theory of EPI from the Born-Oppenheimer approximation, here in this review, and also in many of the practical applications 61-65 , we adopt a phenomenological point of view, and use the model Hamiltonians as given below in Eqs. (8) and (10). Focusing only the term linear in the displacements away from the equilibrium positions of the ions, we can think of the single electron Hamiltonian H e below having a Rdependence. Taylor expanding it, where |j is the single particle state when ionic system is in equilibrium position R 0 . The extra factor of square root of ion mass m p is because of our convention of displacement variable u. This form of interaction is intuitively understandable and originally proposed by Bloch 66 . In Chap. 4 of Ref. 6, a derivative from Eq. (6) to (7) is given, but the reasoning does not seem to be rigorous. Thus, our starting point of a derivation is a tightbinding Hamiltonian for the electrons, harmonic couplings for the phonons and a standard EPI term. They are taken as given and exact. The charge redistributions and self-consistency for the electrons are not part of the discussion. Symbolically, the total many-body Hamiltonian is given as where the electron part is H 0 e = c † Hc, the phonon part H p = 1 2 (p T p + u T Ku) + V n (u C ). The variable u is mass normalized, u j = √ m j (R j − R 0 j ). Because of this, the conjugate momentum is p =u. V n is the nonlinear force contribution. c is a column vector of the electron annihilation operators, which we can separate into three regions, the left, center, and right, c = (c L , d, c R ) T , T stands for matrix transpose. Similarly u = (u L , u C , u R ) T . Accordingly, the matrices H and K are partitioned into nine regions (submatrices), e.g., such that H 0 . Note that we assume no interaction between the left and right leads (See Ref. 67 for transport when there is a lead-lead coupling). We do a similar partitioning for K using the notation of Ref. 68. The EPI takes the form We assume that the EPI appears only in the central region. A schematic representation of the system setup is shown in Fig. 1.
The separation of the electron and phonon leads makes the theoretical development easier. In reality, they could either be physically separated, or built into one. For example, one electrode could serve both as an electron and a phonon lead, but we assume that we have independent control over their temperatures T α e and T α p , α = L, R.
FIG. 1. Model system considered in this review. The center device, including both electrons, phonons, and their interactions, is coupled with two electron and two phonon leads. Each electron lead is characterized by its chemical potential µα and temperature T α e , and each phonon lead by temperature T α p .

A. Theory
We first consider the case where EPI is weak, so that we can perform a perturbation expansion over the interaction matrix M . In order to do so, we use the NEGF method. Detailed introduction is given in our previous work 32, [68][69][70] . This section can be considered as an application of the general approach developed in Refs. [68][69] to the EPI problem. We use similar notations therein, and only give a brief outline of the approach here.
We denote the electron device Green's function without and with EPI by G 0 and G, the corresponding phonon Green's functions by D 0 and D, and the lead Green's functions without coupling to the center as g α and d α , respectively. The couplings of the device with the leads and that between the electrons and phonons are described by self-energies, with Σ and Π representing that of electron and phonon, respectively. For example, we define the time-ordered electron Green's function including EPI on the Keldysh contour [ Fig. 14 Here, τ /τ is time on the contour, and i/j is index of the electronic states. The contour time order operator T C puts the operators later in the contour to the left. The average · is with respect to the density matrix of the full Hamiltonian. The contour ordered Green's function can be divided into different groups according to the spatial position of i/j, similar to the Hamiltonian. The most interesting one is G C , where i and j are both at the center device region. At the same time, it can be written as a 2 × 2 matrix in time space with G t , Gt, G > , G < the time-ordered, anti-timeordered, greater and lesser Green's functions. The retarded and advanced Green's functions are obtained from them, i.e., G r = G t − G < , and G a = G < − Gt. For the definition and relations among these Green's functions, we refer to the book by Haug and Jauho 71 , and our previous publications 32,68,69 .
To calculate the Green's functions, we use a process of two-step adiabatic switch on. We start from the decoupled system and leads. Each of the electron and phonon leads is at its own equilibrium state, characterized by the temperature T α and/or chemical potential µ α . The corresponding equilibrium Green's functions can thus be defined according to the equilibrium canonical distribution. The initial state of the system is arbitrary and not important in most cases (e.g., for steady state).
At the first step, we switch on the interaction of the center Hamiltonian with the electron and phonon leads. We wait until the electron and phonon subsystem reaches their own nonequilibrium steady state, since the temperature and/or chemical potential of each lead can be different. The two subsystems are quadratic and exactly solvable, and we get the non-interacting center Green's functions G 0 and D 0 from the Dyson equation (we omit the superscript C) Here, we have used a single number to represent the matrix indices and contour time arguments, i.e., G 0 (1, 2) = G 0j 1j2 (τ 1 , τ 2 ). Summation or integration over repeated indices is assumed. g C (d C ) is the center electron (phonon) Green's function without coupling to the L and R leads. The self-energy Σ b = Σ L + Σ R includes contributions from L and R, with Σ α (1, 2) = H Cα g α (1, 2)H αC ; similarly for Π b .
At the second step, we adiabatically switch on the EPI in the center. We perform a perturbation expansion over the interaction Hamiltonian H ep , using Feynman diagramatics. The interacting Green's functions are expressed using similar Dyson equations as Eqs. (13)(14), (4,2). (16) Here, Σ ep and Π ep are electron and phonon self-energies due to EPI. Using Eq. (12), at steady state, we can get the following useful relations in energy/frequency domain We use ε for the energy of electron and ω for the angular frequency of phonon, respectively, and I is the identity matrix. To get an expression for the current, we also need the greater and lesser version of the Green's functions 71 The electrical current (I e ) is expressed as the change rate of the electron number in one of the leads (N α ) times the charge of electron (−e). For example, It can be expressed by the Green's function of the center region and the lead self-energies [71][72][73] , Similarly for the heat current carried by electrons (I h ) and phonons (I p ) We have defined the positive current direction as electrons going from the lead to the center. We dropped the argument of the Green's functions for simplicity. We ignore the spin degrees of freedom, since it is not relevant here. Currents out of the right lead are obtained by replacing index L by R. One can symmetrize the expressions based on energy and charge conservation. The set of coupled equations Eqs. (17)(18)(19)(20)(21)(22) is difficult to solve, due to the many-body EPI. Since the EPI is weak, we consider only the lowest order Feynman diagrams shown in Fig. 2. The expressions for the selfenergies are as follows. The electron Fock self-energies from phonons are The Hartree self-energy does not depend on energy The phonon self-energies from electrons are with ε − = ε− ω. Summation over repeated indices is assumed here. Different charge and energy conserving approximations have been developed in the literature. We will use two of them. In Subsec. III B, we perform an expansion of the current up to the second order in M , following the idea of Ref. 29. In the numerical model calculation in Subsec. III C, we use the self-consistent Born approximation (SCBA), which means we replace G 0 , D 0 by G, D respectively in the above equations.

B. Thermoelectric transport through a single electronic level
We consider a single electronic level H C e = ε 0 d † d, coupled to the left and right electrodes, characterized by the constant level-width broadening Γ α with energy cutoff ε D (see Eq. (40) for the general definition). It interacts with an isolated phonon mode with frequency ω 0 , and H ep = m 0 d † du. In the linear regime, we introduce an infinitesimal change of the chemical potential or temperature at lead L, e.g., µ L = µ + δµ, T L σ = T + δT σ , with σ = e or p, µ and T are the corresponding equilibrium values. We look at the response of the charge and heat current due to this small perturbation. The result, up to the 2nd order in M , is summarized as follows, The linear conductance and the Seebeck coefficient are The coefficients L n are with We have defined andᾱ means the lead different from α. L (1) n is the single electron Landauer result. L (2) n is the quasi-elastic term.
Since we are looking at the linear response regime, f L = f R , we dropped the subscript in Eq. (39). We will also use the Bose-Einstein distribution later When there is no ambiguity, we will also drop the argument T .
In the following, we set the position of the electronic level to ε 0 = 0, and look at the dependence of the conductance on the chemical potential µ. We firstly write down the expressions for the self-energies, and make some observations based on their functional forms.
The Hartree self-energy is real and does not depend on energy 32 At T = 0, we get with Γ = Γ L + Γ R . For large enough ε D , the lower limit term turns to −m 2 0 /(2ω 2 0 ), which is the polaron energy shift. Note here that the 1/2 is due to the fact we use u C k in our definition of H ep (Eq. 10). This is different from the common definition that uses the creation and annihilation operatore a † + a (Eq. 74). We have subtracted the polaron shift term in the following calculation, since it is a constant. After this subtraction, Σ r H is odd in µ with a negative slope near µ = 0. We focus on the µ > 0 regime. It saturates to −m 2 0 /(2ω 2 0 ) for large µ, e.g., µ Γ. At non-zero temperature, the slope and the saturation value change, but the shape of the curve is similar to the T = 0 case. This means that, the Hartree term shifts the electronic level, and reduces the conductance. On the other hand, when µ Γ, the conductance tends to zero, whether we include the EPI or not. Thus, the correction to the conductance due to the Hartree term ∆G eH ≤ 0. It starts from zero at µ = 0, goes back to zero at µ Γ, and it reaches a maximum magnitude at some point in the middle. The described behaviour is schematically shown in Fig. 3 (a). This effectively reduces the broadening of the single level spectral function, also the conductance peak. Since the Seebeck coefficient is related to the logarithmic derivative of the conductance, we expect it to increase the magnitude of the Seebeck coefficient near resonance, and to reduce it off resonance [ Fig. 4].
The imaginary part of the retarded Fock self-energy is 32 It is negative and even in ε. Its role on the differential conductance at the phonon threshold (eV = ω 0 ) has been discussed extensively [74][75][76][77] . The main conclusions are: it reduces the differential conductance at eV = ω 0 for resonant case (µ ∼ 0), where the bare transmission without EPI (T 0 ∼ 1), while it does the opposite for far off resonance case (µ Γ), where T 0 → 0. The transition point between the two opposite behaviors is T 0 = 1/2 if the electronic density of state (DOS) is flat. But, in general, it depends on the system parameters. At nonzero temperature, the sharp threshold broadens out, and the linear conductance is affected: its correction to the conductance ∆G eF I is negative for µ ∼ 0, and positive for µ Γ [ Fig. 3 (c)]. Physically, ImΣ r F gives rise to phonon scattering processes. Its effect can be understood as follows: At T = 0, for small bias (eV < ω 0 ), phonon emission is not possible due to Pauli blocking, while phonon adsorption is not possible due to zero phonon population. So, ImΣ r F does not affect the linear conductance. At high enough temperature, both phonon emission and adsorption are possible even at small bias, due to the broadening of the Fermi distribution, and finite population of phonon modes. The phonon scattering process decreases the conductance on resonance, but increases it far off resonance. As a result, the Seebeck coefficient becomes smaller.
The real part ReΣ r F (ε) is obtained by the Hilbert transform of the imaginary part. At zero temperature, it diverges logarithmically at ε − µ = ± ω 0 77 . Its effects on the conductance and Seebeck coefficient are difficult to analyze. We rely on the numerical result [ Fig. 3 Figure 3 shows the correction to the linear conductance of different self-energy terms as a function of µ. These numerical results confirm our qualitative analysis. By comparing the total conductance at low [ Fig. 3 (e)] and high temperature [ Fig. 3 (f)], we see that, (1) the Hartree term dominates at low temperature, and the G e -µ peak becomes narrower. (2) the Fock term becomes important at high temperature, and the G e -µ peak broadens out. Their effects on the Seebeck coefficient (S) are shown in Fig. 4. At low temperature, when the EPI is included, the magnitude of S gets larger for µ ∼ 0, and smaller for |µ| Γ. At high temperature, the effect of the Fock term results in drop of S. In any case, the correction to S is small for weak EPI. But for the case of strong EPI, the correction could be large (see Subsec. IV C).

C. Heat transport between electrons and phonons
Let us look back at the setup in Fig. 1. We want to study the heat transport between electrons and phonons at finite temperature bias, but zero voltage bias. The simplest setup is that the system couples to one electron and one phonon lead, each at its own temperature, see where, again, summation over repeated indices is assumed. For the ease of analysis, we perform an expansion   of the above expression to 2nd order in M , and it becomes The model system we consider to study the energy transport between electrons and phonons. where and Now the question we ask is whether there is a diode behaviour for the heat transport between electrons and phonons 52 , e.g., Q(∆T ) = Q(−∆T ), with ∆T = T e − T p . This is relevant because in some special situation, e.g., at metal-insulator interface, or insulating molecular junctions, EPI becomes the bottleneck of heat transfer [78][79][80][81][82] . We can define the rectification ratio as If we assume a constant electron DOS, Λ(ω) does not depend on T , we get Q(−∆T ) = −Q(∆T ), and R = 0. The physical reason is that in this case, it is possible to map the electron-hole pair excitation into harmonic oscillators 79,83,84 . Then, it is equivalent to heat transport within a two-terminal harmonic system. We do not expect any rectification effect. To make R = 0, the electronic DOS has to be energy-dependent within the broadening of the Fermi-Dirac distribution given by k B ∆T . This effectively introduces anharmonicity into the system, consistent with previous studies [79][80][81] . We can go one step further, by making a Taylor expansion of the spectral function A(ε) about the Fermi energy, we find that the sign of R is determined by the sign of ∂ 2 A ∂ε 2 (The 1st order term is zero). To check this argument, we calculated the heat current across one-dimensional (1D) metal-insulator junction. The metal side is represented by a 1D tight-binding chain, with hopping element t = −0.1 eV, and onsite energy ε 0 = 0. The insulator side is represented by a 1D harmonic chain with the spring constant K = 0.1 eV/(Å 2 u). The insulator and metal couple through their last two degrees of freedom. Their interaction matrices are Here, m 0 = 0.05 eV/(Å √ u), k represents the phononic degrees of freedom. This means that the system couples Phonon modes can be excited by the mobile electrons due to the EPI effect; i.e., a high bias over the system leads to self-heating (Joule heating). In nanoscale electric devices, the electric current density can be much larger than that in the macroscopic system. The high current density will generate strong Joule heating, which may eventually break the device. In this sense, Joule heating becomes a bottleneck for further increase of the electric current density. Hence, lots of theoretical and experimental efforts have been devoted to understanding the Joule heating phenomenon in the nanoscale electric devices. In experiment, Joule heat can be measured via the thermal-mechanical expansion technique, which records the Joule heat induced temperature rise 85 . The Joule heat can increase the temperature of the molecular junction from room temperature to 463 K, which has been examined through the inelastic electron tunneling spectroscopy 86 . Grosse et al. investigated the nanoscale Joule heating in phase change memory devices 85 . The Joule heating leads to the temperature rise in the phase change memory device, which results in an obvious volume expansion. In another experiment, Joule heating is found to be responsible for the correlated breakdown of nanotube forests 39,41 .
For a system without localized phonon modes, all phonon modes have important contribution to the Joule heat. The Joule heat contributed by these propagating phonon modes have important effects on the electric devices. For example, in graphene transistors, the output electric current will saturate with increasing source-drain voltage 87,88 . This saturated current density can be reduced by 16.5% due to the Joule heating 88 .
The localized phonon modes exist around some defects or nonuniform configurations, such as the free edge, the isotropic doping, interface, etc. This particular type of phonon modes has no direct contribution to the thermal conduction, but localized phonon modes play a particularly important role in the Joule heat phenomenon. They are characteristic for their exponentially decaying vibration displacement; i.e., only a small portion of atoms are involved in the localized vibration. For instance, there are some localized edge phonon modes at graphene nanoribbon's free edge. In these modes, edge atoms vibrate with large amplitude, but the vibrational displacement decays exponentially from the edge into the interior region.
The localized-phonon-mode-induced Joule heat was observed in graphene nanoribbons in experiment, and explained theoretically. Jia, et al. utilized Joule heating to trigger the edge reconstruction at the free edges in the graphene nanoribbons 89 . Engelund, et al. attributed this phenomenon to the Joule heating of the edge phonon modes 90 . There are two conditions for the important Joule heating of the edge phonon modes. First, these localized edge phonon modes can spatially confine the energy at the edges. Second, the electrons interact strongly with the localized edge phonon modes. The mean steadystate occupation of the edge phonon mode can be calculated from the ratio of the current-induced phonon emission rate and damping rate. The effective temperature for the free edge can be extracted by assuming this occupation to be Bose distributed. The effective temperature was found to be as high as 2500 K for bias around 0.55 V. This high effective temperature was proposed to be the origin for the edge reconstruction.
Although Joule heating might be used for selectively bond-breaking 91,92 , its most common outcome is a disaster of device breakdown. The effective temperature is a suitable quantity to describe the Joule heating. In 1998, Todorov studied Joule heating problem in a molecular junction 12 . In his work, the Einstein model is applied to represent the phonon modes in the system, and the electron-electron interaction is ignored as the system size is much smaller than the electron mean free path. Part of the EPI-induced Joule heat will be delivered out of the system by the phonon heat conduction, while the remaining Joule heat gives a high effective temperature. For low ambient temperature, the effective temperature scales with voltage V as T 4 eff ≈ γ 4 V 2 , with γ as an EPIdependent constant. It was shown that the effective temperature can be above 200 K for a very low ambient temperature around 4 K 93 . But at very high bias, the scaling law could differ from this 31 .
There are several experimental approaches to investigate the effective temperature of the electric device induced by Joule heating. The effective temperature can be extracted by measuring some quantities that are temperature-dependent. For example, the breaking force of the single molecular junction is related to the temperature. This force-temperature relationship can be used to estimate the effective temperature 30 . The Raman spectroscopy also depends on the temperature. Hence, it can be used to deduce the effective temperature of Ramanactive phonon modes 35,37,38 .
It has also been shown that EPI has an important effect on the thermal conductance in single-walled carbon nanotubes (SWCNTs) 94 . For them, we apply the Born approximation to consider the EPI effect using the NEGF approach, as the SCBA is computationally more expensive. The phonon thermal current can be calculated by considering the three EPI contributions shown in the Feynman diagrams in Fig. 2. The phonon thermal current flowing from the left lead into the center is given by Eq. (26). The expression for the right lead is analogous. The Joule heat is generated in the system and flows into the leads, so the total Joule heat is the sum of heat currents into both leads, Q = −(I L p + I R p ). The thermal current from Eq. (26) also includes that induced by the temperature gradient, which satisfies I L p = −I R p . Hence, Q gives solely the Joule heat.For metallic SWCNT (10, 10), both electrons and phonons are important heat carriers. The EPI only slightly reduces the electron thermal conductance, but it has a strong effect on the phonon thermal conductance. More specifically, Fig. 7 (a) shows an 'electron-drag' effect on the phonon thermal conductance at 150 K. The phonon thermal conductance becomes negative for high chemical potential value µ > 2.0 eV, which indicates that electrons can help to drag phonons from cold temperature region to the hot temperature region. The 'electron-drag' phenomenon happens at low temperature and high chemical potential, and it does not happen at a higher temperature 300 K as shown in Fig. 7 (b). For semiconductor SWCNT (10, 0), the electronic thermal conductance contributes less than 10% of the total thermal conductance at low bias (e.g. µ = 0.3 eV), while phonons make most significant contribution to the total thermal conductance. Similar 'electron-drag' phe-nomenon also exists in the semiconductor SWCNT (10, 0) at low temperature as shown in Fig. 8 (a).

A. Quantum Master equation formulism
In this section, we introduce the QME approach to consider the case of strong EPI. Before doing that, we should mention that the NEGF method has also been used to treat the strong EPI [95][96][97][98][99][100] . Since the idea behind it is very similar to that of the master equation approach, we choose not to introduce it here.
To simplify the formula, we ignore the coupling of molecular phonon modes to the phonon leads. The model Hamiltonian simplifies to where H S = H C p + H C e + H ep denotes the system Hamiltonian, and V e = V L e + V R e is the system-lead coupling. In the QME formalism we assume the system-lead coupling V e is weak so we can do perturbation on it. We work in the interaction picture with H 0 = H − V e as non-interacting part and V e as the interaction. For simplicity, in this section we use V to represent V e since we don't have V p . The equation of motion for the full density matrix follows the von Neumann equation Here, the subscript I denotes operator in the interaction picture. The time argument in the parentheses means non-interacting evolution O(t) = e iH0t/ Oe −iH0t/ . The above equation can be written in an integral form as One can recursively apply the above equation to get a series expansion of the full density matrix in power of V I . We truncate the series to the second order and differentiate it with respect to time t at both sides of the equation to get the following integro-differential equation We prepare the initial state as a product state of the system and each lead, ρ(t 0 ) =ρ(t 0 ) ⊗ ρ L e ⊗ ρ R e . For the system-lead coupling V we assume it can be written as a product of system operator S and lead operator B as V = α S α ⊗ B α . In such cases we can trace over the lead degrees of freedom to get is the correlation function of the leads. Here we have used the condition that the expectation value of a single lead operator B α is zero. We can now transform back to the Schrödinger picture and extend the initial time t 0 to −∞ to get the QME of Redfield type 101-103 Here we have replacedρ(t 0 ) byρ, which is essential and correct only when one intends to get the 0thorder reduced density matrixρ 0 by solving the above QME 69,104,105 . In the application to the EPI problem, by exact diagonalizing the system Hamiltonian, this Redfield QME can take into account the coherence between electrons and phonons, in contrast to the usual rate equation approach 106,107 . We write the above equation in the eigenbasis of the system Hamiltonian H S to obtain 102,108 where the relaxation tensor reads 109 The transition coefficients are given by where ∆ kj = E k − E j are the energy spacings of the system Hamiltonian.
Since we are only interested in the steady state, we impose the condition dρ/dt = 0 at t = 0 and solve the above equation order by order with respect to V . One can find that all the off-diagonal elements of the 0thorder reduced density matrix vanish in steady state and the diagonal elements can be evaluated via the matrix equation 109 i R ii nnρ (0) ii = 0, together with the constraint of Tr[ρ (0) ] = 1.
For the calculation of currents, we go through a similar derivation as the QME. The electronic current operator J e and heat current operator J h can be written in the form J e(h) = α S α ⊗ B α e(h) . The expectation value of currents can be calculated according to I e(h) = Tr[ρ I (t)J I (t)]. Since we are interested in the lowest order of current, we truncate Eq. (57) to the lowest order and plug in to get By taking the trace over the leads and transforming back to the Schrödinger picture one obtains the current at t = 0 as where C αβ e(h) (t) = B α (t)B β e(h) (0) is the correlation function between the lead operators occurring in the systemlead coupling Hamiltonian and the current operator definition. For the same reason as the derivation of master equation, hereρ(t 0 ) needs to be replaced byρ to get correct steady state results. Since we are calculating lowest order of current, we can use the 0th order reduced density matrixρ (0) . Written in the eigenbasis of system Hamiltonian, the above equation becomes I e(h) = Tr ρ (0) I r e(h) with the reduced current operator defined as where the transition coefficients are .
Up to now the QME formalism is general and not restricted to any specific form of system or leads Hamiltonians. For the application to transport problems with EPI as concerned in this review, the system-coupling Hamiltonian is considered as a tunneling Hamiltonian 72 . In such case the system operator S, leads operator B and B can be specified The infinite nature of the leads can be specified by defining a continuous spectra function for the leads Throughout this section we use a wide band spectra function for the electronic leads with Lorentzian cut-off as The non-vanishing correlation functions can be evaluated via and C 12 (t) = C 12 L (t) + C 12 R (t), C 12 e (t) = −eC 12 L (t), C 21 e (t) = eC 21 L (t). We note that here the upper index 1 or 2 refers to the two components of B and B e/h given above. For the system Hamiltonian, we focus only on the single electronic level coupled to a single phonon mode representing the center of mass of the molecule. In such case, the system Hamiltonian will reduce to with ω 0 the angular frequency of the phonon mode and λ denotes the EPI strength. This type of Hamiltonian has been well-studied in the context of molecular junction 32,33,55,97,99,110-119 . The QME formalism treats the nonlinearity of EPI exactly. Therefore in this section we will focus on strong EPI regime with emphasis on the EPI strength dependence of the electron/heat current. In the following we will discuss the effect of EPI on the electronic transport properties, including the phonon sidebands, negative differential resistance, thermoelectric properties and local heating effects.

B. Phonon sidebands and negative differential resistance
One of the earliest findings of the vibrational effects on the electronic transport through a molecular quantum dot is the appearance of the phonon sidebands in the I − V characteristics. When electrons transport through a single electronic level, the differential conductance (dI/dV ) will manifest a peak at resonant level when plotted against the voltage bias (V ). However, when the electronic level is interacting with a vibrational mode, replica side peaks will appear at the side of the resonant peaks. These side peaks are called phonon sidebands. A simple reason for the appearance of phonon sidebands is due to the fact that the electrons can emit or absorb phonons when they pass through the molecule. Therefore, the distance between each adjacent peaks is always equal to a single phonon energy. The phonon sidebands attract wide attentions in molecular junction systems. Experimentally, phonon sidebands were found in 1980s 120 and then were utilized to identify vibrational modes in molecular junctions [121][122][123][124] and quantum wires 24,[27][28][29]125 . Theoretically, at the beginning the sidebands were investigated by using scattering theory, which gives the transmission probability T (ε, ε ) for an electron to passing through an EPI system. The electron is coming from vacuum at energy ε and leaving at energy ε 126 . The scattering theory predicts side peaks in the transmission probability, which qualitatively justifies the phonon sidebands in molecular junction systems. However, prediction of phonon sidebands in the lead-molecule-lead junctions is a much more difficult task. A simple generalization to take the Fermi-Dirac statistics nature of the electron leads into account is to weight the exact transmission probability with the Fermi-Dirac distribution of each lead, i.e., by multiplying the transmission probability T (ε, ε ) by a factor f L (ε)[1 − f R (ε )] as a new transmission probability for electron going from the left lead to the right lead through the nano-conductor. This approximation is called single particle approximation (SPA). Plenty of earlier work is in this framework and it predicts Lorentzian type of phonon sidebands with the same width. But obviously this method assumes each electron transports independently through the junction, where the many-body effects are ignored. As a result, it overestimates the currents and it is not able to predict the quantized conductance e 2 /h either.
Based on the NEGF technique, more rigorous methods merged in dealing with the nonlinearity in EPI, such as the Green's function equation of motion method (EOM) 95 , SCBA [96][97][98]127 and nearest neighbor crossing approximation (NNCA) 100 . All of the above are Green's function based formalisms with different kinds of approximations. In general these approaches predict that the phonon side peaks are much sharper than the SPA approach. This sharpness is closely related to the Pauli exclusion, which the SPA approach failed to take into account. Other than Green's function based methods, another approach is to use rate equation of electron occupation probability in the molecule, via calculating the transition coefficients of the electron to tunnel from the molecule to each lead and vice versa. This method assumes the transport is an electron tunneling process and the electron will lose its phase information when it resides in the molecule. Therefore it will be valid when the molecule-lead coupling is weak and the coherence of the electron and phonon in the molecule can be neglected. For all these formalisms, we would like to point out that one should take care of the phonon distribution. Treating the phonon at equilibrium distribution at a fixed temperature could be valid when the EPI strength is much weaker than the coupling strength between the phonon and its environment. However, when the environmental influence is weak, one should consider phonons in nonequilibrium states. This nonequilibrium treatment of phonon distribution can have pronounced effect on I − V characteristics due to the fact that the current induced vibrational excitation can be significant 127,128 . Besides the peak distances and peak width discussed above, other aspects characterizing the sidebands include the weights of the zero-phonon band and the number of peaks. The investigation of these properties mainly focuses on the effects of EPI strength, Fermi energy of the molecule, chemical potential and temperature of the leads. In general, the higher order peaks will be suppressed by the Frank-Condon factor 77,128,129 . In the framework of NEGF, Chen et al. found that the weight from zero-phonon band will decrease monotonically with increase of EPI strength and temperature while the weights of higher order sidebands will increase and then decrease 97 . The chemical potentials of the leads will influence the presence of the sidebands at both sides of the zero phonon peaks [ Fig. 9, panel (a) and (b)]. If one keeps the chemical potential of one lead fixed and increases the chemical potential of the other lead, the phonon sidebands will appear only at one side of the 0th order peak [Fig 9, panel (c)]. However, if one fixes the Fermi-level of the molecule and changes the chemical potentials of both leads, phonon sidebands will appear at both sides 97,98 . We note that, in this section, Fermi-level of the molecule is defined as (µ L + µ R )/2.
In the framework of QME formalism described earlier, which is exact in the weak system-lead coupling limit, we find phonon sidebands for the single electronic level interacting with a single phonon mode. In this case the zero-phonon peak will occur at the renormalized resonant level due to polaron shift (ε 0 − λ 2 / ω 0 ) and the sidebands will appear at every distance of ω 0 . The peaks will appear at each side of zero-phonon peak under symmetric change of the lead chemical potentials while only appear at one side if we fix chemical potential of one lead. The EPI strength will not only shift the peaks, but also modulate the weights of the peaks. We find that upon increasing either EPI strength or temperature, the weight of the zero-phonon peak will decrease, which is consistent with the previous work 97 . However, interesting phenomena happen when the renormalized level of the quantum dot gets close to the Fermi-level of the dot, [ε 0 − λ 2 /ω 0 ≈ (µ L + µ R )/2], i.e., additional peaks appear at each side. Those peaks have distance ω 0 with the zero-order peak at the opposite side. However, the two major peaks really merge together, and those additional peaks disappear again. In the case of asymmetric change of lead chemical potential (right most panel of Fig. 9 (c)), we found peaks appear at both sides when the zero-phonon peaks merge together.
Another interesting perspective of the I − V characteristics is the phenomenon of the negative differential resistance (NDR), where the current decreases with the increase of voltage bias. The NEGF formalism predicts that NDR is impossible in ballistic electronic transport, but it will emerge in the presence of EPI. NDR has been both theoretically investigated 76,128,130 and experimentally measured 131 . An important reason of NDR is due to the redistribution of the molecular states. As discussed in the previous section, the zero-phonon band carries major portion of electronic current. The probability for the molecule to be in that state will be related to the chemical potential of the leads. If one increases the bias via increasing the chemical potential of one lead, one actually lifts the Fermi-level of the molecule as well. If one brings the Fermi-level of molecule far away from the eigenenergy of zero phonon state, the probability of the molecule to be in that state will decrease and thus the current will decrease. Based on this analysis one can draw several immediate conclusions: 1. If one increases the bias simultaneously for both leads in pace and keeps the Fermi-level of the molecule fixed, there will be no NDR. 2. If one treats the phonon fixed in the equilibrium distribution, NDR will not appear 128 . 3. If the chemical-potentialvarying lead couples stronger to the molecule than the chemical-potential-fixed lead, the redistribution will be more sensitive, thus the NDR will be enhanced. Figure 10 shows the NDR predicted in the QME formalism. We find that the NDR appears in the symmetric system-lead coupling. Moreover, NDR will be enhanced if the chemical-potential-varying lead couples stronger to the molecule than the other lead, but it will disappear in the other way around. We also find that NDR is most pronounced in the moderate EPI regime, while it is less significant in both weak and strong EPI regimes 128 .

C. Vibrational effect on thermoelectricity
In this part, we study the effects of EPI on the thermoelectric properties of the nano-conductor. We first look at the thermoelectric current, which is the electronic current induced by a temperature difference between the leads. Thermoelectric current exhibits quite different features comparing with voltage-bias current. It will increase monotonically and smoothly with the increasing of the temperature bias. Therefore, there will be no phonon sidebands. This is expected, because under temperature bias, the tunneling channel of the electron is always restricted to the molecular state that is close to the chemical potential of the leads. The increasing of temperature will only excite more conducting electrons, but not be able to extend extra tunneling channels. Therefore, there will be no sudden change of thermoelectric current. Due to the same reason, there will be no NDR effect as the increasing of the temperature bias will always make more electrons to be involved in tunneling. The restriction on tunneling channels also makes the thermoelectric current much smaller than the voltage-bias current.
The dependence of electronic conductance on the chemical potentials of the leads has been discussed in the Subsec. III B. The sign of the Seebeck coefficient indicates that the currents will change direction with different chemical potential. Another interesting perspective is to find the dependence of currents on EPI strength. Figure 11 shows the plot of voltage-bias current [panel (a)] and thermoelectric current [panel (b)] with the EPI strength λ. For voltage-bias current, we find that the current decays with EPI strength in general. More precisely, the maximum current one can achieve via adjusting the ε 0 is decreasing monotonically with EPI strength. This is due to the Frank-Condon blockage of the current. However, for each ε 0 we can find the enhancement of current due to EPI, and that is mainly because of the polaron shift which can bring the electron resonance level into the conduction band from outside. For the thermoelectric current, the profile is quite different. We see that the current can change sign with the increase of EPI strength when ε 0 is higher than the Fermi level, which indicates that the EPI can switch the charge carriers of the quantum dot between electrons and holes 59 . For each ε 0 there exist two optimized values of EPI strength such that the thermoelectric current maximizes. This optimized λ shifts left with decrease of ε 0 until disappear one by one at the λ = 0 end.
One important quantity to describe the efficiency of the thermoelectric material is the figure of merits Z e T . It is related to the electronic conductance G e , Seebeck coefficient S, thermal conductance G h via the formula Z e T = G e S 2 T /G h . Here we use the notation Z e T to denote the figure of merits of the system by ignoring the thermal conductance due to phonons. So G h here only takes account the thermal conductance due to electrons. The effect of phonons on the figure of merit Z e T is rather complicated, which is closely related to the electron Fermi energy 55,118,119,132 , the phonon energy 113 , the temperature and chemical potentials of the leads 118,119,132,133 . Figure 12 shows the dependence of the electronic conductance, thermal conductance, Seebeck coefficient and figure of merit on the electron energy ε 0 and EPI strength λ. The major effect of EPI is on the thermal conductance of the molecular junction 113 .
The EPI can open extra channels, from which high energy electrons can tunnel from hotter lead to colder lead, while low energy electrons tunnel in the opposite direction. Therefore, the thermal conductance is enhanced while the electronic conductance is not affected too much. Due to the increase of the thermal conductance, Z e T will be suppressed drastically 113,119 . Though the figure of merits Z e T will be reduced quickly under the influence of phonon scattering in weak EPI regime, it will gradually saturate at strong EPI regime. We also would like to point out that when the electron energy is close to the Fermi energy of the quantum dot, the Seebeck coefficient and hence the Z e T will become very small. However, the phonon scattering can renormalized the electron energy via polaron shift and thus the Seebeck coefficient can be enhanced. As a result, EPI can enhance Z e T in this particular parameter regime.

D. Local heating
Local heating is an important phenomenon in molecular junction, not only due to its own importance affecting the stability of the system, but also due to its close relation to the phonon sidebands 127 , NDR and thermoelectric effect 56,133 . The distribution of the phonon states in current-carrying system can be far away from equilibrium 127 , or in some cases, may even lead to phonon instability 17,43,134,135 . Phonons can be excited significantly by the voltage bias 128 , and in turn affect the I − V characteristics. At each peak of the phonon sidebands, one can actually find a vibrational excitation event 128,136 . Previous study also found that the local heating can enhance the thermoelectric efficiency 56,133 . However, in most cases local heating is not preferred because it will affect the stability of the system 94 and introduce noise to the measurements [137][138][139][140] . Therefore, lot of effort has been put into cooling the system using electronic current, such as using superconducting single electron transistor 141 , or double quantum dots 140 .
Here we study the effects of the electronic current on the phonon mode of the nano-conductor. To investigate the local heating, effective temperatures are usually defined in various ways 9,32,51,[142][143][144] . In this part, since we only have one single phonon mode, instead of defining the effective temperature, we use average phonon number to characterize the local heating effect. The way to specify the electronic current induced heating is to compare the nonequilibrium phonon number n neq with equilibrium phonon number n eq . When the molecule is weakly coupled to the leads, the molecular states statistics will follow canonical distributionρ = e −βH S /Tr(e −βH S ) and thus the equilibrium phonon number can be calculated exactly as 59 The first term is the Bose-Einstein distribution function, the second term is a correction due to the polaron energy shift. The nonequilibrium phonon number can be calculated from the nonequilibrium reduced density matrix obtained from the QME. Figure 13 shows the difference of phonon numbers under voltage bias (top) and temperature bias (down). For voltage bias, ∆n is always positive which indicates that the system is always heated up. The yellow regime where ∆n ≈ 0 is the regime where the electronic current vanishes. When there is electronic current passing through, the heating effect is generally more pronounced in stronger EPI regime. However, for the temperature bias, we find both heating (∆n > 0) and cooling regimes (∆n < 0). Therefore, the local heating effect is not only related to the magnitude of electronic current, but also related to the energy each electron carries when it tunnels into the molecule. For the thermoelectric current, low energy electron can tunnel from the cooler lead to the molecule, absorb a phonon and tunnel to the hotter lead, resulting in cooling of the molecule. Such process is impossible in the voltage-bias case in the present setup as the electron is always flowing from the higher chemical potential side to the lower chemical potential side.

V. CURRENT-INDUCED SEMI-CLASSICAL LANGEVIN DYNAMICS
In the previous two sections, we mainly look at the effect of phonons on the electric and thermoelectric transport properties of electrons, which is also the focus of most published work. But, to study the current-induced forces, and their effect on atomic dynamics, we need to turn around. In this section, starting from the total Hamiltonian H tot , we derive a semi-classical Langevin equation to describe the atomic dynamics of the system, coupled to both phonon and (nonequilibrium) electron leads. The Langevin equation applies to the weak EPI regime, since similar to Sec. III our expansion parameter is the interaction matrix M . However, the derivation is not limited to the form of H tot . Following the same procedure we discuss below, we can also do an adiabatic expansion of the electron influence functional over the velocity of the ions. In that case, the Langevin equation applies to slow ions, but the EPI could be of arbitrary magnitude [17][18][19]145 .
The advantage of the Langevin equation approach is that, we can easily include the anharmonic phonon interaction, as in other molecular dynamics method. The anharmonic interaction is crucial in dealing with currentinduced dynamics. This is because the phonon modes that interact with electrons are normally high-frequency ones, while the low frequency modes conduct heat from the system to surrounding electrodes more effectively. The energy transfer from high to low frequency modes is possible only when anharmonic interaction is included. Although possible, it is not a trivial task to incorporate anharmonic interactions in the NEGF or QME approach.

A. Initial states and reduced density matrices
We assume at a remote past t 0 and earlier time, the central system is decoupled with the leads and there is also no EPI, so that the electrons and phonons are also decoupled. The density matrix is assumed to be product of known equilibrium states. For example, the left lead is specified by Exactly what to take for the center is not important as for steady state with t 0 → −∞, the results do not depend on it (except maybe in very subtle cases).
The density matrix (of the whole system) is governed by the von Neumann equation, and formally we can write where assuming t > t 0 . For the other case of t < t 0 , the timeorder operator T should be replaced by the anti-time order operator. We are interested only in the center, so the leads degrees of freedom will be traced out. For notational simplicity, we assume only one left lead. The result for two or more leads is trivially generalized. We defineρ where the first reduced density matrix only eliminates the lead, while the last one eliminates the electrons as well, leaving only an effective density matrix for the phonons. The procedure to eliminate the lead for both the electrons and phonons follows the standard method of Feynman and Vernon 146 , except that we need to be careful for the electrons which are fermions. Since there is no coupling between electrons and phonons in the leads, the phonon and electron degrees of freedom can be done separately (the initial density matrix is a product of the two).

B. Influence functional for phonons
The matrix elements of the density operator ρ(t) is taken in the basis of the coordinates u. Following the standard treatment 103,147 , the density matrix is then given by The path C 2 consists of two segments (see Fig. 14(a)) running from time t with coordinate variable u back to t 0 at variable u 0 , and then second segment running from To eliminate the lead, the phonon Lagrangian is split into L = L C + L L − V where V is the coupling potential energy term between the lead and center, given by Taking the trace of Eq. (82) means that we identify u L and (u ) L as the same variable u L and integrate it out. As far as variable u L is concerned, the path is of Keldysh type, i.e., running from t 0 to t from above and then back from t to t 0 , as shown in Fig. 14(b). The reduced density matrix is theñ with the influence functional The above expression can be further simplified. Firstly, the initial distribution of the lead is in thermal equilibrium, ρ L 0 ∝ e −β L H L . Secondly, we can rewrite the path integral formula back into the operator form; thus, we obtain where Z L is the canonical partition function, T C is contour order operator, the subscript eq.L stands for equilibrium average with respect to the left lead. One more transformation can be made to simplify it further. The time-dependence (or rather, contour time τ dependence) is understood to be in the Heisenberg picture governed by the Hamiltonian H L + V (τ ) which is different in the forward and backward direction. We can work in the interaction picture, thus eliminate the explicit H L from the formula. The resulting equation is This is the same equation [Eq. (5)] given in Ref. 148.
The influence functional can be calculated explicitly since the interaction V is a quadratic form, and the contour operator naturally leads to contour ordered Green's functions and Wick's theorem is valid. We define the contour ordered phonon Green's function of the lead as Expanding the exponential (or a cumulant expansion, expanding and taking logarithm), we get The first order term (and all odd order in V terms) vanishes, because u L = 0. We have defined the ubiquitous contour ordered self-energy due to the lead as It can be shown that the last line in the above derivation is an exact result. It is instructive to clarify and compare with other notations used in the literature. A commonly used notation is (e.g., Ref. 103) Using the rules C dτ → σ t t0 σdt, and Π(τ, τ ) → Π σ,σ (t, t ), u C (τ ) → u σ (t) where σ = + or − for the upper and lower branch, the relations among the various self-energies, and symmetry relation Π > ij (t, t ) = Π < ji (t , t), we can rewrite Eq. (88) in the form of Eq. (90). By comparison, we find α(t, t ) ≡ L(t, t ) is the notation used by Schmid 149 .

C. Influence functional for electrons
The derivation of the influence functional for the electrons is similar except that we have to deal with grassmann integrals [150][151][152][153] . We follow the approach of Weinberg 154 . To do the trace over the leads we need a specific representations for the operator ρ. For the electrons, we use the coherent state characterized by grassmann numbers such that The hat denotes operator, without hat, it is a grassmann number. The state |c has explicit form given as The orthogonality is in the form Similar states for the creation operatorĉ † j can be constructed with eigenvalue (a grassmann number)c, given the following results (similar to inner product of eigen states of u and its conjugate momentum p and completeness of the eigenstates).
The · · · is an extra + or − sign factor which we'll not keep track. The tilde over the product sign means that order is in exactly the opposite canonical order [e.g., that of Eq. (95)]. With the above very sketching outline, the fermion evolution operator can be represented as a path integral of the form with the action The lead influence functional can be then obtained with the same procedure as that for the phonons. The result involves an integral kernel which is exactly the contour ordered self-energy of the lead Σ L (τ, τ ).

D. Reduced density matrix for phonon in center
Putting things together, the reduced density matrix, when the leads are eliminated, has the form The ordinary number (column vector) u and grassmann numberc, c involve only the degrees of freedom of the center. For notational simplicity, we have dropped the superscript C. Note that the electron terms do not have the characteristic factor 1/2 asc and c are independent variables. The dependence onc and c is a bi-linear form, thus the path integral over them can be done analytically. This gives the reduced density matrix of the phonon only, as The influence functional to the phonons due to electrons is given by (110) Interpreting the τ in the above as Keldysh variable defined on C has a problem. As agreed, the contour is supposed to be on C 2 with the t 0 end connected with the initial distribution of the electrons ρ 0 at the center. However, if we assume that in the limit t 0 → −∞ the results should not depend on the distribution of the center, we can ignore this initial distribution and it is completely fixed by the lead. But we cannot give a mathematically sound justification here.
Similar to that in Ref. 154 for the field theory of quantum electrodynamics, we want to put the influence functional in an exponential form. This can be done using the formula, Det(A) = e TrlnA , and the expansion of the function ln(1 + x) = x − x 2 /2 + · · · for small x, given where G −1 0 and y are matrices indexed by lattice sites j as well as contour time τ . And, if the proper metric for a discretization of the time is chosen so that det(G −1 0 ) can be meaningful, we can identify G 0 as the electron contour ordered Green's function when there is no EPI defined in Sec. III. Since det(G −1 0 ) is independent of u, the effective action for the phonon is only determined by the exponential factor, which is a polynomial (functional) in u. With some caveat regarding the initial distribution, Eqs. (109) and (111) offer a formally exact solution to the problem.

E. Semi-classical approximation
If we ignore the linear term in u which produces a constant force, the effect of which is to shift the equilibrium positions, and also neglect higher order contributions, we end up with a quadratic form for the effective action dτ u(τ ) T Π tot (τ, τ )u(τ ). (117) with Π tot = Π L + Π R + Π ep as in Eq. (20). We have also included the right phonon lead. A generalized Langevin equation can be derived from the above action 18,149 where F n = −∂V n /∂u, Π r tot is the retarded total selfenergy, and the noises satisfy We note that the effect of the electron leads to the phonons has exactly the same form as that of the phonon leads. The self-energy consists of a sum of contributions of the two sources.

F. Applications
Before discussing the applications, we note that similar generalized Langevin equation as Eq. (118) can be derived by doing an adiabatic expansion over the momenta of the ions to the 2nd order [17][18][19]155 .
Its most important feature is the inclusion of the quantum nature of the electron and phonon leads. For example, the zero point motions of atoms are correctly taken into account, and proved critical in determining the thermal and structural properties of materials made from light elements [158][159][160][161][162][163]170,171 .
Including the correct Bose distribution of the phonons opens a way to study the quantum ballistic phonon transport by doing classical molecular dynamics. In Refs. 156 and 157 the transition from ballistic to diffusive phonon thermal transport is studied using this approach. In Refs. 17 and 34, including the nonequilibrium electrons, current-induced dynamics have been studied. Several interesting effects have been predicted or confirmed, and their effects on the stability of the system are studied. For example, it has been shown that (1) the current-induced forces are not conservative 15 , (2) the atoms feel an effective magnetic force, originating from the Berry phase of the electrons 17 . Moreover, the power of the Langevin approach is to be able to include the anharmonic phonon-phonon interactions classically, and treat the EPI quantum-mechanically. This enables one to study the energy transport between different phonon modes, between electrons and phonons at the same time. The exploration of its power and potential is still under way.
As an example, we consider the heat generation in a 4 × 2 graphene armchair ribbon, see Ref. 157 for the definition of structure parameters. The electronic structure and EPI matrix are obtained from a combined SIESTA 172 + TranSIESTA 173 + Inelastica 64 calculation, while the Brenner potential is used for the inter-atomic interaction. To reduce the simulation time, we have ignored the energy-dependence of the electronic structure. As a result, the electronic friction becomes time-local. The Langevin equation becomes (after an integration by part) (121) where F C is the force from the second-generation Brenner potential, dΓ(t)/dt = Π r (t) is the phonon retarded self-energy due to two leads, V is applied bias voltage. The eV ξ − u term gives a nonconservative force as ξ − is antisymmetric. ξ = ξ L + ξ R is the noise due to left and right phonon leads, while f is the noise due to electron bath. The expression for the electronic friction and noise correlation is the same as Eqs. (56)(57)(58)(59)(60)(61)(62)(63) in Ref. 18. Further implementation details can be found in Ref. 157. The phonon heat current is calculated using (122) where α = L, R. In steady state, the energy flow balances, and the heat generation is calculated according to Q = −I L p − I R p . The rate of heat generation is plotted in Fig. 15. The result is for the same configuration as shown in Ref. 157 of Fig. 4. Each data point takes about 4 days on a typical Opteron CPU. The error bars are quite small. The heat generation at zero bias should be zero. However, we get a small value. This has to do with the cut-off used in the noise for the electrons. We have used an abrupt cut-off for the spectrum at ω = 1.29 eV. The calculation demonstrates the feasibility of computing the Joule heating current, intrinsically a quantum effect at nanoscale, by classical molecular dynamics.

VI. CONCLUSIONS
In summary, using a tight-binding-like Hamiltonian for the EPI, Eq. (8), we have introduced three different approaches to study the effect of EPI in different parameter regimes. We focused on the electronic, phononic, and thermoelectric transport properties of nano-conductors, in a general multi-probe setup. For each approach, we started with the theoretical derivation of the main equations. This was followed by applications in models or simple systems, mainly for illustration purpose.
Applications of these approaches to more interesting problems are straightforward. Examples of such prob-lems are: (1) application of the NEGF and QME approach to nonequilibrium thermoelectric transport to study the thermoelectric efficiency at finite power output, (2) application of the QME approach to look at system where both EPI and electron-electron interaction are important, (3) combining the generalized Langevin equation with first-principles or tight-binding electronic structure package to study current-induced dynamics in realistic nano-conductors, especially to explore how the electron-dissipated heat is transferred in and out of the nano-conductor. With these available tools, more interesting and important systems can be investigated.