# Electron-vibration effects on the thermoelectric efficiency of molecular
junctions

C. A. Perroni, D. Ninno, and V. Cataudella  CNR-SPIN and Dipartimento di
Fisica, Universita’ degli Studi di Napoli Federico II ,
Complesso Universitario Monte S. Angelo, Via Cintia, I-80126 Napoli, Italy

###### Abstract

The thermoelectric properties of a molecular junction model, appropriate for
large molecules such as fullerenes, are studied within a non-equilibrium
adiabatic approach in the linear regime at room temperature. A self-consistent
calculation is implemented for electron and phonon thermal conductance showing
that both increase with the inclusion of the electron-vibration coupling.
Moreover, we show that the deviations from the Wiedemann-Franz law are
progressively reduced upon increasing the interaction between electronic and
vibrational degrees of freedom. Consequently, the junction thermoelectric
efficiency is substantially reduced by the electron-vibration coupling. Even
so, for realistic parameters values, the thermoelectric figure of merit can
still have peaks of the order of unity. Finally, in the off-resonant
electronic regime, our results are compared with those of an approach which is
exact for low molecular electron densities. We give evidence that in this case
additional quantum effects, not included in the first part of this work, do
not affect significantly the junction thermoelectric properties in any
temperature regime.

##  I Introduction

In solid state systems, a voltage induces a temperature gradient and vice
versa. These phenomena are known as thermoelectric effects. nolas  In some
semiconducting materials, nolas ; ioffe  these thermoelectric effects can be
strong enough to allow either the fabrication of devices converting wasted
heat into electrical energy or the realization of solid-state coolers. A key
requirement to improve the energy conversion efficiency is to increase the
electrical conductance ($G$) and the Seebeck coefficient ($S$) reducing the
electronic ($G_{K}^{el}$) and lattice ($G_{K}^{ph}$) contributions to the
thermal conductance $G_{K}$. A value of the dimensionless figure of merit
$ZT=GS^{2}T/G_{K}$ of the order of $1$ is considered a fundamental
prerequisite of a useful thermoelectric devices. nolas ; shakouri  A clear
limitation of the thermoelectric technology is that three mutually
contraindicated properties of the same material have to be optimized. In
metals, for instance, $ZT$ is typically limited by the Wiedemann-Franz law,
stating that the ratio $G_{K}^{el}/(GT)$ is a constant (the Lorenz number)
independent of the metal specificities.

In order to optimize the thermoelectric efficiency, the possibility of
controlling materials at the nanoscale has been advanced. shakouri ; dressel1
; dressel2 ; koumoto  For example, large values of $ZT$ can be obtained in
semiconducting nanowires with highly peaked densities of states. dressel3  It
has been predicted that also molecular devices can be efficient for conversion
of heat into electric energy. arad ; dubi  The improvement of thermoelectric
efficiency derives from the discreteness of energy levels that leads to the
violation of the Wiedemann-Franz law. Mahan  Therefore, the emerging field of
molecular thermoelectrics can be very interesting for both the basic physics
and applications and it is mainly for this reason that the subject has
attracted a lot of attention in recent years. majum1 ; baheti ; majum2 ; finch
; murphy ; galperin ; koch ; leijnse

Single-molecule measurements have focused on the Seebeck coefficient whose
sign directly provides the sign of charge carriers involved into the transport
mechanisms. majum1 ; majum2 ; datta ; cuevas ; datta1  Most measurements have
taken into account simple small molecules where the transport is dominated by
the HOMO (highest occupied molecular orbital). majum1 ; baheti  Recently, more
controllable alignment between Fermi level and molecular orbitals (whose
energy separation is still of the order of $0.5$ eV) has been achieved with
larger molecules whose transport is dominated by the LUMO (lowest unoccupied
molecular orbital). Measurements of conductance and Seebeck coefficient in
junctions based on $C_{60}$ molecules have been performed considering three
different metallic electrodes (platinum, gold, and silver).majum2  A very high
magnitude of single-molecule thermopower has been reported ($S$ of the order
of $-30$ $\mu V$/K). However, the application of a gate voltage in these kinds
of measurements remains elusive. Moreover, heat transport in molecular devices
remain poorly characterized owing to experimental challenges. dubi ; wang1 ;
wang2  Recently, heat transport and dissipation have been investigated in
junctions with benzene-like molecules, cuevas1  even if the study has been
limited to the elastic transport regime.

Non-interacting models of molecular junctions, using a Landauer-type approach,
datta1  are generally not accurate since intramolecular interactions typically
constitute the largest energy scales of the problem. The electron-vibration
coupling, indeed, significantly affects the transport characteristics of
molecular devices. cuevas ; galperin1  In particular, either the molecule
center of mass oscillations Park  or thermally induced acoustic phonons Qin
can be the source of coupling between electronic and vibrational degrees of
freedom. The electron-vibration coupling has been studied in a fully out-of-
equilibrium linear response regime with different theoretical tools, ranging
from rate equation siddi ; mitra ; dong  (in the regime of weak tunnel
coupling between molecule and electrodes) to non equilibrium Green function
formalism fang ; galperin2 ; wrudzinski  (for perturbative and intermediate to
strong strength of interaction). In particular, in devices with large
molecules or carbon nanotube quantum dots, the low energy of the relevant
vibrational degrees of freedom has been exploited to devise a non-equilibrium
adiabatic approach. This method is semiclassical for the vibrational dynamics
but it is valid for arbitrary strength of electron-vibration coupling. mozy ;
pistol ; hussein ; alberto ; alberto1

The thermopower $S$ has been analyzed in a single level molecule with
vibrational coupling by means of the master equation approach finding that it
is sensitive to the interplay between electrons and molecular vibrations. koch
The effect of electron-vibration coupling on the thermoelectric properties of
few level molecules or dots has been studied by the rate equation leijnse  and
Green function galperin ; yang ; zianni ; ren ; tagani  formalism. However, in
these papers, the phonon thermal contribution $G_{K}^{ph}$ to the figure of
merit $ZT$ has not been calculated or discussed. Recently, this contribution
to the thermoelectric efficiency has been investigated in a molecular junction
where a benzene molecule is directly connected to platinum electrodes. hsu
However, the electron-vibration inelastic effects on the thermoelectric
properties have been evaluated only at a perturbative level.

In this paper, we have studied the thermoelectric properties of a molecular
junction with electron-vibration coupling within the linear response regime
focusing on the phonon thermal contribution $G_{K}^{ph}$ to the figure of
merit $ZT$ at room temperature. The non-equilibrium adiabatic approach has
been used to solve the junction model which takes into account the interplay
between the low frequency center of mass oscillation of the molecule and the
electronic degrees of freedom. Parameters appropriate for junctions based on
$C_{60}$ molecules connected between different metallic leads have been
considered. We have found that the semiclassical $G_{K}^{ph}$ typically
overcomes the electronic thermal conductance $G_{K}^{el}$, and it gets
enhanced with increasing the electron-vibration coupling. Moreover, the
increase of the electron-vibration coupling makes the ratio $G_{K}^{el}/GT$
closer to the Lorentz number. Actually, the figure of merit $ZT$ can be
substantially reduced by these effects, even if, for realistic parameters of
the model, it can still have peaks of the order of unity. Finally, in the off-
resonant regime, where the thermoelectric properties show peak values, we have
compared the results of the adiabatic approach with those of a fully quantum
formalism which is exact for low electron level density. We have stressed
that, in the regime of low temperatures, the additional quantum effects in
$G_{K}^{ph}$ not included in the adiabatic approach poorly affect the
thermoelectric efficiency.

The paper is organized as follows. In Sec. II, the model of molecular junction
is proposed. In Sec. III, the adiabatic approach is quickly explained. In Sec.
IV, the results within the adiabatic approach are discussed. In Sec. V, the
comparison of the results within the adiabatic approach with those of a method
exact in the regime of low level occupation is performed. Two Appendices close
the paper: Appendix A, where the derivation of the Langevin equation for the
center of mass oscillator is reported, and Appendix B, where some results
about the oscillator damping rate induced by the electron-vibration coupling
and position distribution function are commented.

##  II Molecular junction model

Figure 1: (Color online) Sketch of the molecular junction studied in this
work. The curved lines between dots (lead atoms) depict charge electron
hoppings in the lead bulks ($t^{\prime}$) and between lead and molecule ($t$).
The broken lines between dots (lead atoms) depict springs in the lead bulks
(with elastic constant $k^{\prime}$) and between lead and molecule (with
elastic constant $k$), therefore, they denote the phonon part. The Left Lead
and the Right Lead are kept at chemical potential $\mu_{L}$, temperature
$T_{L}$ and chemical potential $\mu_{R}$, temperature $T_{R}$, respectively.

The molecular junction model includes both electric and elastic coupling
between leads and molecule, whereas the molecule is described as a single
electron level coupled to the molecule center of mass vibration (see Fig. 1
for a scheme of the device). The Hamiltonian $\hat{H}$ of the junction is,
then, given by

$\hat{H}={\hat{H}}_{el}+{\hat{H}}_{ph}+{\hat{H}}_{int},$ (1)

where the Hamiltonian ${\hat{H}}_{el}$ (${\hat{H}}_{phon}$) describes
electronic (vibrational) degrees of freedom of both leads and molecule and
${\hat{H}}_{int}$ describes the coupling between electronic and vibrational
degrees of freedom of the molecule. We assume that the electronic and
vibrational degrees of freedom in metallic leads are not interacting cuevas ;
Haug , and, therefore, the electron-vibration coupling is assumed effective
only on the molecule.

The molecule is modeled as a single electronic level locally interacting with
a single vibrational mode (see Fig. 1). This means that the focus is on a
molecular electronic orbital which is sufficiently separated in energy from
other orbitals. For example, in a free $C_{60}$ molecule, the LUMO energy
differs from the HOMO energy for more than 1 eV. Even when the degeneracy of
the LUMO is removed by the contact with Ag, the splitting gives rise to levels
which are separated by an energy of the order of $0.5$ eV. Lu  In this last
situation, the low energy level of $C_{60}$ gets aligned with the Fermi level
of Ag majum2 ; Lu  increasing the thermopower. Finally, the choice of a single
level model is motivated by the fact that a system with these electronic
features can improve the thermoelectric performances. Mahan

Since sizable figures of merit ZT are obtained for temperatures much higher
than Kondo temperature of molecular junctions, cuevas  and large Seebeck
coefficients are found when the molecular orbital is poorly populated, we
focus our analysis in a regime where the effects of Coulomb local interactions
are negligible. Hence, the electronic Hamiltonian ${\hat{H}}_{el}$ of Eq. (1)
is assumed spinless:

$\hat{H}_{el}=\epsilon{\hat{d}^{{\dagger}}}{\hat{d}}+\sum_{q,\alpha}\varepsilon_{q,\alpha}{\hat{c}^{{\dagger}}_{q,\alpha}}{\hat{c}_{q,\alpha}}+\sum_{q,\alpha}\left[V_{q,\alpha}{\hat{c}^{{\dagger}}_{q,\alpha}}{\hat{d}}+h.c.\right],$
(2)

where the molecular electronic level has energy $\epsilon$ and
${\hat{d}^{{\dagger}}}({\hat{d}})$ are creation (annihilation) operators on
the molecule. The presence of a gate in the junction can be simply simulated
by changing the value of the local energy $\epsilon$. cuevas  The operators
${\hat{c}^{{\dagger}}_{q,\alpha}}({\hat{c}}_{q,\alpha})$ create (annihilate)
electrons with momentum $q$ and energy
$\varepsilon_{q,\alpha}=\xi_{q,\alpha}-\mu_{\alpha}$ in the left ($\alpha=L$)
or right ($\alpha=R$) free metallic leads. The difference of the electronic
chemical potentials in the leads, $\mu_{L}$ and $\mu_{R}$ respectively,
provides the bias voltage $V_{bias}$ applied to the junction:
$\mu_{L}-\mu_{R}=eV_{bias}$, with $e$ electron charge. The left and right
leads will be considered as thermostats in equilibrium at the temperatures
$T_{L}$ and $T_{R}$, respectively, with temperature difference $\Delta
T=T_{L}-T_{R}$ (see Fig. 1). Therefore, the left and right electron leads are
characterized by the free Fermi distribution functions $f_{L}(\omega)$ and
$f_{R}(\omega)$, respectively. The electronic tunneling between the molecular
dot and a state $q$ in the lead $\alpha$ has the amplitude $V_{q,\alpha}$. As
usual for metallic leads, the density of states $\rho_{q,\alpha}$ is assumed
flat about the small energy range relevant for the molecular orbital, making
valid the wide-band limit: $\rho_{q,\alpha}\mapsto\rho_{\alpha}$,
$V_{q,\alpha}\mapsto V_{\alpha}$. Therefore, the full hybridization width of
the molecular orbital is $\hbar\Gamma=\sum_{\alpha}\hbar\Gamma_{\alpha}$, with
$\hbar$ Planck constant and the tunneling rate
$\Gamma_{\alpha}=2\pi\rho_{\alpha}|V_{\alpha}|^{2}/\hbar$. In junctions with
$C_{60}$ molecules, $\hbar\Gamma$ has been estimated to be of the order of
$20$ meV. natelson ; mra  In the following, we consider the symmetric
configuration: $\Gamma_{L}=\Gamma_{R}=\Gamma/2$; $\mu_{L}=\mu+eV_{bias}/2$,
$\mu_{R}=\mu-eV_{bias}/2$, with $\mu$ average chemical potential;
$T_{L}=T+\Delta T/2$, $T_{R}=T-\Delta T/2$, with $T$ average temperature.

In analogy with the electronic model, we consider only one relevant
vibrational mode for the molecule. We will focus on the center of mass mode,
which is expected to have the lowest frequency for large molecules. Moreover,
this mode should have a frequency smaller than the Debye frequency of the
metallic leads impeding the formation of localized vibrational states on the
molecule. Therefore, this mode can favor the exchange of heat between the
leads. Besides, in $C_{60}$ molecules, experimental results provide compelling
evidence for a coupling between electron dynamics and the center of mass
motion. Park

In Eq.(1), the Hamiltonian ${\hat{H}}_{phon}$ describes the vibrations of the
molecule center of mass mode, the free phonon modes of the leads, and the
coupling between them:

${\hat{H}}_{ph}={\hat{H}}_{cm}+\sum_{q,\alpha}\hbar\omega_{q,\alpha}{\hat{a}^{{\dagger}}_{q,\alpha}}{\hat{a}_{q,\alpha}}+\sum_{q,\alpha}\left(C_{q,\alpha}{\hat{a}_{q,\alpha}}+h.c.\right){\hat{x}}.$
(3)

The center of mass hamiltonian ${\hat{H}}_{cm}$ is

${\hat{H}}_{cm}={{\hat{p}}^{2}\over 2m}+\frac{k{\hat{x}}^{2}}{2},$ (4)

where ${\hat{p}}$ and ${\hat{x}}$ are the center of mass momentum and position
operator, respectively, $m$ is the total large mass, $k$ is the effective
spring constant, and the low frequency is $\omega_{0}=\sqrt{k/m}$. For
$C_{60}$ junctions, $\hbar\omega_{0}$ has been estimated to be of the order of
$5$ meV. Park  Therefore, for a large molecule such as $C_{60}$, the adiabatic
regime, $\omega_{0}<<\Gamma$, is valid.

In Eq.(3), the operators
${\hat{a}^{{\dagger}}_{q,\alpha}}({\hat{a}}_{q,\alpha})$ create (annihilate)
phonons with momentum $q$ and frequency $\omega_{q,\alpha}$ in the lead
$\alpha$. As shown in Fig. 1, the left and right phonon leads will be
considered as thermostats in equilibrium at the temperatures $T_{L}$ and
$T_{R}$, respectively, which we assume to be the same as that of the electron
leads. In order to investigate heat exchange with the molecule, the leads
phonon spectrum is assumed to be acoustic. For silver (atomic number Z=47),
gold (Z=79), and platinum (Z=78) leads considered in experimental
measurements, majum2  the Debye frequency is such that $\hbar\omega_{D}$ is of
the order of $18.5$ meV, $15.1$ meV, and $20.7$ meV, respectively. kittel
Therefore, $\hbar\omega_{D}\simeq 15-20$ meV for these leads. In any case, as
for any large molecule, the center of mass mode is such that
$\omega_{0}<<\omega_{D}$.

In Eq.(3), the coupling between the center of mass position and a phonon $q$
in the lead $\alpha$ is given by the elastic constant $C_{q,\alpha}$. In order
to characterize this interaction, one introduces the spectral density
$J(\omega)$:

$J(\omega)=\frac{\pi}{2}\sum_{q,\alpha}\frac{C_{q,\alpha}^{2}}{M\omega_{q,\alpha}}\delta(\omega-\omega_{q,\alpha})=m\omega\tilde{\gamma}(\omega),$
(5)

with $M$ atomic mass of the leads and $\tilde{\gamma}(\omega)$ complex
frequency dependent memory-friction kernel of the oscillator.weiss  In the
regime $\omega_{0}<<\omega_{D}$, $\tilde{\gamma}(\omega)$ can be approximated
as real and independent of frequency, providing the damping rate $\gamma$:
$\tilde{\gamma}(\omega)\simeq\gamma$. weiss  In analogy with the electronic
model, we consider the symmetric configuration:
$\gamma_{L}=\gamma_{R}=\gamma/2$. If the center of mass is elastically coupled
with a neighbor atom of the leads by a spring with constant $k$, one gets
$\gamma\simeq 16k^{2}/(mM\omega_{D}^{3})$. Taking the mass $m$ of the $C_{60}$
molecule and the atomic mass $M$ of Ag, Au, and Pt, $\hbar\gamma$ is of the
order of $7.68$ meV, $7.74$ meV, and $2.98$ meV, respectively. The smallest
value of coupling to phonon baths is due to the largest Debye frequency
$\omega_{D}$ of platinum. In any case, $\hbar\gamma\simeq 3-8$ meV for these
metals, therefore $\omega_{0}$ is of the same order of $\gamma$.

Finally, in Eq.(1), the interaction term ${\hat{H}}_{int}$ derives from
electrostatic interactions since it is induced by the charges injected by the
leads onto the molecule. When the center of mass displacement is not large, it
has been shown that this interaction is provided by a linear coupling between
the electron occupation on the molecule,
${\hat{n}}={\hat{d}^{{\dagger}}}{\hat{d}}$, and the $\hat{x}$ operator of the
oscillator:

${\hat{H}}_{int}=\lambda{\hat{x}}{\hat{n}},$ (6)

where $\lambda$ is the electron-oscillator coupling strength. alberto1  This
coupling gives rise to a model known in the literature as Anderson-Holstein
model. cuevas  In the following, the electron-vibration interaction will be
often described in terms of the coupling energy $E_{P}=\lambda^{2}/(2k)$. As
reported in experimental measurements, Park  the effects of $E_{P}$ are not
negligible in junctions with $C_{60}$ molecules.

##  III Adiabatic approach

In this paper, we focus on the dynamics of the molecule center of mass which
is generally slow for large molecules ($\hbar\omega_{0}\simeq 5$ meV for
junctions with $C_{60}$ molecules). Moreover, the interest is for the
thermoelectric properties close to ambient temperature $T_{A}$. Therefore, in
the following, we will study the system under arbitrary electron-vibration
coupling in the range $\hbar\omega_{0}\ll k_{B}T$ ($k_{B}$ Boltzmann constant)
assuming that the dynamics of the center of mass can be treated as classical.

###  III.1 Electron dynamics dependent on oscillator parameters

Since the oscillator is assumed to be classic, the electronic dynamics is
equivalent to that of a time dependent level with energy
$E_{0}(t)=\epsilon+\lambda x(t)$, where $x(t)$ is the oscillator position.
Using the Keldysh formalism, Haug ; alberto  one can exactly solve the Dyson
and Keldysh equations for the molecular Green functions integrating out
exactly the lead electronic degrees of freedom.

Exploiting the fact that the full hybridization width, $\hbar\Gamma$, between
the molecule and the leads, is much larger than the oscillator energy
$\hbar\omega_{0}$, one can determine the adiabatic expansion of the electronic
Green’s functions considering the explicit dependence of the electronic
quantities on the oscillator position $x$ and making expansions on the
oscillator velocity $v=p/m$. alberto ; alberto1 ; alberto2 ; perroni ;
perroni1  For instance, at the lowest order, the molecular level spectral
function is position dependent and given by

$A(\omega,x)=\frac{\hbar\Gamma}{(\hbar\omega-\epsilon-\lambda
x)^{2}+(\hbar\Gamma)^{2}/4}.$ (7)

Adiabatic expansions can be easily evaluated for all the quantities which can
be expressed as a function of the level Green functions. In Appendix A, the
expansion of the level occupation $N(x,v)$ is required in order to make a
self-consistent calculation of the electron and oscillator dynamics.

###  III.2 Dynamics of the center of mass oscillator

In this subsection, we analyze the slow dynamics of the molecule center of
mass. The effects of the electron degrees of freedom and phonon baths give
rise to a generalized Langevin equation for the the center of mass oscillator.
We are mostly interested on the junction transport properties close to room
temperature, therefore in the regime
$\hbar\omega_{0}\ll\hbar\omega_{D}\simeq\hbar\Gamma<k_{B}T$.

The generalized Langevin equation for the oscillator dynamics is derived in
Appendix A. The resulting motion equation

$m\frac{dv}{dt}=F_{det}(x,v)+\xi(x,t)$ (8)

has the deterministic force $F_{det}(x,v)$ and the position dependent
fluctuating force $\xi(x,t)$. The deterministic force

$F_{det}(x,v)=F_{gen}(x)-A_{eff}(x)v,$ (9)

can be decomposed into a generalized force $F_{gen}(x)$

$F_{gen}(x)=-kx+F_{\lambda}(x),$ (10)

where $F_{\lambda}(x)$ is due to the electron-vibration coupling, and a
dissipative force with an effective position dependent positive definite term
$A_{eff}(x)$

$A_{eff}(x)=A_{\lambda}(x)+m\gamma,$ (11)

with $A_{\lambda}(x)$ due to the electron-vibration interaction. The
fluctuating force $\xi(x,t)$ in Eq.(8) is such that

$\langle\xi(x,t)\rangle=0,\;\;\;\;\langle\xi(x,t)\xi(x,t^{\prime})\rangle=D_{eff}(x)\delta(t-t^{\prime}),$

where the effective position dependent noise term $D_{eff}(x)$ is

$D_{eff}(x)=D_{\lambda}(x)+k_{B}(T_{L}+T_{R})m\gamma,$ (12)

with $D_{\lambda}(x)$ determined by the electron-vibration coupling. In
equilibrium conditions at temperature $T=T_{\alpha}$ and chemical potential
$\mu$ ($V_{bias}=0$ and $\Delta T=0$), one gets a generalized fluctuation-
dissipation relation $D_{eff}(x)=2k_{B}TA_{eff}(x)$, since this equation is
verified for each fixed position $x$.

We have numerically solved the Langevin equation under generic non-equilibrium
conditions using a generalized Runge-Kutta algorithm. alberto ; honey ; honey1
As a result of the numerical calculations, the oscillator distribution
function $Q(x,v)$ and the reduced position distribution function $P(x)$ can be
determined giving insights on the oscillator positions relevant for the
dynamics.

In Appendix B, we thoroughly discuss the features of the parameter
$A_{\lambda}(x)/m$ in the linear response regime. The peak value of
$A_{\lambda}(x)/m$ is always smaller than the values of $\gamma$ considered in
this paper even for strong electron-vibration coupling $E_{P}$. Actually, the
quantity $A_{\lambda}(x)/m$ will not strongly affect the phonon conductance
$G_{K}^{ph}$ for realistic $E_{P}$ couplings. Typically, the effects due to
the electron-vibration coupling on the oscillator dynamics do not represent a
large perturbation with respect to those induced by the coupling to phonon
leads. Obviously, as discussed in Appendix B, the effects of the electron-
vibration coupling depend on the occupation of the electronic level. In the
regime of low occupation, the dynamics of the oscillator is poorly influenced
by these effects even for strong $E_{P}$.

###  III.3 Electron and oscillator quantities

In the adiabatic regime, once solved the Langevin equation, one can derive the
behavior of a electronic observable $O(x,v)$ dependent on oscillator
parameters. For example, we calculate a central quantity of our work, the
electronic spectral function $A(\omega)$, making the average of the function
$A(\omega,x)$ in Eq.(7):

$A(\omega)=\int_{-\infty}^{+\infty}dxP(x)A(\omega,x).$ (13)

Within this approach, the features of the electronic system can be fully
characterized in general non-equilibrium conditions. We point out that, within
the adiabatic regime, charge and electronic energy currents are conserved on
the stationary non-equilibrium states. In this paper, we will focus on the
regime of linear response around the average chemical potential $\mu$ and the
temperature T ($\Delta T\rightarrow 0$, $V_{bias}\rightarrow 0$). We will
evaluate the electronic conductance $G$

$G=\left(\frac{e^{2}}{\hbar}\right)\left(\frac{\hbar\Gamma}{4}\right)\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}A(\omega)\left[-\frac{\partial
f(\omega)}{\partial(\hbar\omega)}\right],$ (14)

where $f(\omega)=1/(\exp{[\beta(\hbar\omega-\mu)]}+1)$ is the free Fermi
distribution corresponding to the chemical potential $\mu$ and the temperature
$T$, with $\beta=1/k_{B}T$. Then, we will calculate the Seebeck coefficient
$S=-G_{S}/G$, with

$G_{S}=\left(\frac{e}{\hbar}\right)\left(\frac{\hbar\Gamma}{4T}\right)\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}(\hbar\omega)A(\omega)\left[-\frac{\partial
f(\omega)}{\partial(\hbar\omega)}\right].$ (15)

Finally, we will determine the electron thermal conductance
$G_{K}^{el}=G_{Q}+TG_{S}S$, with

$G_{Q}=\left(\frac{1}{\hbar
T}\right)\left(\frac{\hbar\Gamma}{4}\right)\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}(\hbar\omega)^{2}A(\omega)\left[-\frac{\partial
f(\omega)}{\partial(\hbar\omega)}\right].$ (16)

Static quantities relative to the center of mass oscillator can be derived
starting from the distribution functions $Q(x,v)$ and $P(x)$. However, dynamic
quantities relative to the oscillator have to be inferred directly from the
numerical stochastic dynamics. As discussed in Appendix A, the phonon energy
currents are conserved for generic non-equilibrium states:
$J_{R}^{ph}=-J_{L}^{ph}$, with $J_{\alpha}^{ph}$ current from the $\alpha$
phonon lead. Therefore, within the adiabatic approach used in this paper, the
total energy conservation is satisfied since the energy currents are
separately conserved for the electron and vibrational channels. In this paper,
we will focus on the phonon thermal conductance $G_{K}^{ph}$ calculated within
the linear response regime wang3  around temperature T as

$G_{K}^{ph}=\lim_{\Delta T\rightarrow 0}\frac{(J_{L}^{ph}-J_{R}^{ph})}{2\Delta
T}.$ (17)

The total thermal conductance is $G_{K}=G_{K}^{el}+G_{K}^{ph}$ making feasible
the evaluation of the figure of merit $ZT$. When the coupling of the center of
mass mode to the metallic leads is absent ($\gamma=0$), $G_{K}=G_{K}^{el}$, so
that $ZT=ZT^{el}$, which can be used to characterize the electronic
thermoelectric properties.

In this paper, we will consider parameters appropriate to junctions with
$C_{60}$ molecules. Therefore, $\hbar\Gamma\simeq 20$ meV will be the energy
unit. As a consequence, $\Gamma$ will be the frequency unit. We will assume
$\omega_{0}=0.25\Gamma$ and vary $\gamma$ from $0.15\Gamma$ to $0.40\Gamma$
(simulating the effects of different metallic leads). We will measure lengths
in units of $\lambda/k$, times in units of $1/\Gamma$, temperatures in units
of $\hbar\Gamma/k_{B}$ (ambient temperature $T_{A}\simeq 1.25$ in these
units). Finally, we will assume the average chemical potential $\mu=0$.

##  IV Results within the adiabatic approach

In this paper, we will discuss linear response transport properties in
different conditions trying to clarify the role of the electron-vibration
coupling. In particular, we will focus on the phonon energy transmission and
on the electronic level variation with respect to the leads chemical
potential. These variations can be controlled, in our model, changing the
molecule energy level $\epsilon$ from $\epsilon=0$, that coincide with the
lead chemical potential (resonant case), to a very different value
$|\epsilon|>>0$ (off-resonant case).

Since our aim is to clarify the effects due to the electron-vibration
coupling, in the next subsection, we start discussing, for comparison, the
electronic response functions in the absence of electron-vibration interaction
($E_{P}=0$) and coupling to phonon leads ($\gamma=0$). In the second
subsection, we will discuss the effects of electron-vibration coupling in the
absence of coupling to phonon leads, then, in the final subsection, the
features with full coupling.

###  IV.1 Results in the absence of electron-vibration interaction and
coupling to phonon leads

As shown in the panel (d) of Fig. 2, the thermoelectric effect is maximum at
specific off-resonant conditions where $ZT^{el}$ can reach values of the order
of unity or even larger, depending on temperature. In fact, we have analyzed
the behavior of transport properties for different temperatures finding that
$ZT^{el}$ gets larger with increasing temperature. These behaviors result from
the relevant role played by the Seebeck coefficient for these regimes of
parameters.

Figure 2: (Color online) Panel (a): Electron conductance G (in units of
$e^{2}/h$) as a function of the level energy (in units of $\hbar\Gamma$) for
different temperatures $T$ (in units of $\hbar\Gamma/k_{B}$). Panel (b):
Seebeck coefficient (in units of $k_{B}/e$) as a function of the level energy
(in units of $\hbar\Gamma$) for different temperatures $T$ (in units of
$\hbar\Gamma/k_{B}$). Panel (c): Electron thermal conductance $G_{K}^{el}$ (in
units of $k_{B}\Gamma$) as a function of the level energy (in units of
$\hbar\Gamma$) for different temperatures $T$ (in units of
$\hbar\Gamma/k_{B}$). For comparison, the thermal conductance quantum
$g_{0}(T)=\pi^{2}k_{B}^{2}T/(3h)$ at $T=0.05\hbar\Gamma/k_{B}$ is shown by
means of the magenta double dot-dash line. Panel (d): Electronic dimensionless
figure of merit $ZT^{el}$ as a function of the level energy (in units of
$\hbar\Gamma$) for different temperatures $T$ (in units of
$\hbar\Gamma/k_{B}$). In all the plots, electron-oscillator coupling $E_{P}=0$
and oscillator damping rate $\gamma=0$ (absence of coupling to phonon leads).

In the off-resonant regime, the charge conductance is expected to be small. As
shown in the panel (a) of Fig. 2, for $|\epsilon|>>0$, the conductance is much
smaller than the conductance quantum $e^{2}/h$ ($e^{2}/h$ is about $3.87\times
10^{-5}$ S). In particular, for $\epsilon=20$, $G$ is of the order of
$10^{-3}$ $e^{2}/h$ in agreement with the order of magnitude of experimental
data in $C_{60}$. majum2  On the other hand, in the resonant case
$\epsilon=0$, the conductance reaches the quantum at low temperatures and,
then, rapidly decreases to values of the order of $0.1$ $e^{2}/h$ with
increasing temperature.

As shown in the panel (c) of Fig. 2, for both resonant and off-resonant
conditions, $G_{K}^{el}$ is small when measured in units of $k_{B}\Gamma$
($k_{B}\Gamma$ is about $419.8$ pW/K for $\hbar\Gamma\simeq 20$ meV). In the
resonant case, $G_{K}^{el}$ decreases with increasing temperature and it is of
the order of $0.01-0.02k_{B}\Gamma\simeq 4-8$ pW/K for different temperatures,
therefore it is smaller than the molecular conductance (per chain) of the
order of $50$ pW/K measured in hydrocarbon molecules anchored to a gold
substrate. wang2  Next, we will stress that the electron-vibration coupling
will induce an increase of the thermal conductivity providing results
compatible with experimental estimates. At low temperatures
($T=0.05\hbar\Gamma/k_{B}$ is about $12$ K), in the resonant case, the result
practically coincides with the thermal conductance quantum
$g_{0}(T)=\pi^{2}k_{B}^{2}T/(3h)$ at that temperature ($g_{0}(T)\simeq
9.456\times 10^{-13}(W/K^{2})T)$. jezouin  In the off resonant regime
$|\epsilon|>>0$, $G_{K}^{el}$ gets larger with increasing temperature, but, in
the investigated temperature range, it is always smaller than
$0.01k_{B}\Gamma$.

As shown in the panel (b) of Fig. 2, the Seebeck coefficient $S$ shows large
variations that are responsible for the behavior of $ZT^{el}$ shown in panel
(d) of Fig.(2). At low temperatures and in quasi-resonant regime
($\epsilon\simeq 0$), $S$ is very small whereas, in the off-resonant condition
and at high temperatures, the thermopower can be very large. When $\epsilon$
is positive (n-type behavior), $S$ is negative. In particular, for
$\epsilon=20$, $S$ is about $-0.45k_{B}/e\simeq-38.5\mu V/K$ in agreement with
the magnitude of experimental data in $C_{60}$. majum2  On the other hand,
when $\epsilon$ is negative (p-type behavior), $S$ is positive. Therefore, as
expected, the sign of $S$ is sensitive to the charges responsible for the
transport. In the off-resonant case, the peak values of $S$ are quite large
since they are of the order of a few $k_{B}/e$ ($k_{B}/e$ is about $86$ $\mu$
eV/K).

Summarizing, in the off-resonant condition and at high temperatures, the
reduction of the conductance $G$ is fully compensated by the strong increase
of the Seebeck coefficient $S$. Moreover, for $\gamma=0$, the thermal
conductance is small. Therefore, as shown in the panel (d) of Fig. 2,
$ZT^{el}$ can acquire values larger than $1$. Finally, we stress that the peak
values of $ZT^{el}$ at room temperature are almost coincident with maxima and
minima of the Seebeck coefficient $S$.

###  IV.2 Results in the absence of coupling to phonon leads

Figure 3: (Color online) Spectral function (in units of $1/\hbar\Gamma$) as a
function of the frequency $\omega$ (in units of $\Gamma$) for different values
of $E_{P}$ (in units of $\hbar\Gamma$) at $\epsilon=1\hbar\Gamma$ (Upper
Panel) and $\epsilon=10\hbar\Gamma$ (Lower Panel). In all the plots,
$T=1.25\hbar\Gamma\ k_{B}$ (close to room temperature),
$\omega_{0}=0.25\Gamma$, and oscillator damping rate $\gamma=0$ (absence of
coupling to phonon leads).

In this subsection, we focus on the effects of the electron-vibration coupling
on the electronic response functions still in the absence of coupling to
phonon leads ($\gamma=0$). We focus on a temperature close to room temperature
($T=1.25$).

First, we analyze the behavior of the spectral function as a function of the
electron-vibration coupling. In the upper panel of Fig. 3, we show the
spectral function for different values of the electron-vibration coupling in
the quasi-resonant case $\epsilon=1$. We point out that there is a strong
transfer of spectral weight toward low frequency with increasing $E_{P}$. In
particular, at $E_{P}=1$, the spectral function is peaked on the chemical
potential $\mu=0$. Indeed, when $\epsilon=E_{P}$, one gets the level density
$n=0.5$, that is, the half-filling condition occurs. Moreover, at $E_{P}=2$,
the peak is much smaller in frequency, and the spectral function acquires
asymmetric features. As shown in the lower panel of Fig. 3, a different
behavior takes place in the off-resonant regime ($\epsilon=10$ is considered
in the figure). For the considered values of $E_{P}$, the spectral function
gets enlarged, but its peak position is quite rigid. Therefore, the behavior
of the spectral function is different in the resonant and off-resonant regime
for fixed values of $E_{P}$.

Figure 4: (Color online) Panel (a): Electron conductance G (in units of
$e^{2}/h$) as a function of the level energy (in units of $\hbar\Gamma$) for
different values of $E_{P}$ (in units of $\hbar\Gamma$). Panel (b): Seebeck
coefficient (in units of $k_{B}/e$) as a function of the level energy (in
units of $\hbar\Gamma$) for different values of $E_{P}$ (in units of
$\hbar\Gamma$). Panel (c): Electron thermal conductance $G_{K}^{el}$ (in units
of $k_{B}\Gamma$) as a function of the level energy (in units of
$\hbar\Gamma$) for different values of $E_{P}$ (in units of $\hbar\Gamma$).
Panel (d): Electronic dimensionless figure of merit $ZT^{el}$ as a function of
the level energy (in units of $\hbar\Gamma$) for different values of $E_{P}$
(in units of $\hbar\Gamma$). In all the plots, $T=1.25\hbar\Gamma\ k_{B}$
(close to room temperature), and oscillator damping rate $\gamma=0$ (absence
of coupling to phonon leads).

The features of the spectral function affect the behavior of the thermoeletric
properties. As shown in Fig. 4, the most relevant effect of the coupling
$E_{P}$ on the conductance $G$ (panel (a)) and the Seebeck coefficient $S$
(panel (b)) is to shift the curves and reduce the magnitude of the response
function. The shift of the conductance peak and of the value where the Seebeck
coefficient vanishes is given by $E_{P}$ ($n=0.5$ for $\epsilon=E_{P}$). We
point out that, at fixed level energy, unlike the conductance $G$, the Seebeck
coefficient shows a large sensitivity to the change of the coupling $E_{P}$.
For example, this occurs for energies close to the minimum and the maximum.
For larger values of $\epsilon$, there is an inversion in the behavior of $S$
with increasing the coupling $E_{P}$.

A different behavior is shown by the thermal conductance $G_{K}^{el}$ (panel
(c) of Fig. 4). Indeed, $G_{K}^{el}$ can be enhanced with increasing the
electron-oscillator coupling $E_{P}$. For example, at the resonance and
$E_{P}=1$, $G_{K}^{el}$ is of the order of $0.05k_{B}\Gamma\simeq 20$ pW/K for
different temperatures, therefore it gets closer to the conductance of the
order of $50$ pW/K measured in hydrocarbon molecules. wang2  Next, we will see
that the total conductance, including also the phonon contribution, is even
closer to these experimental data. Actually, within the adiabatic approach,
more energetic channels open with increasing the electron-vibration coupling
since the effective level is renormalized by the oscillator dynamics that
becomes less localized. The reduction of $G$ and $S$ combined with the
enhancement of $G_{K}^{el}$ leads to a sensible reduction of the figure of
merit $ZT^{el}$ with increasing $E_{P}$. For $E_{P}=0$, the peak value of
$ZT^{el}$ is around $3$, while, for $E_{P}=1$, the peak value is smaller than
2. Therefore, even if one neglects the role of phonon thermal conductance, the
electron-vibration coupling is able to induce an important reduction of the
figure of merit.

###  IV.3 Results in the case of full coupling

Figure 5: (Color online) Upper Panel: Phonon thermal conductance $G_{K}^{ph}$
(in units of $k_{B}\Gamma$) as a function of electron-vibration coupling
$E_{P}$ (in units of $\hbar\Gamma$) for different values of level energy
$\epsilon$ (in units of $\hbar\Gamma$). Lower Panel: Phonon thermal
conductance $G_{K}^{ph}$ (in units of $k_{B}\Gamma$) as a function of level
energy $\epsilon$ (in units of $\hbar\Gamma$) for different values of
electron-vibration coupling $E_{P}$ (in units of $\hbar\Gamma$). In all the
plots, $T=1.25\hbar\Gamma/k_{B}$ (close to room temperature), the oscillator
damping rate $\gamma=0.15\Gamma$, and $\omega_{0}=0.25\Gamma$.

Finally, we analyze the case when also the coupling of the center of mass mode
to the metallic leads is present. First, as shown in Fig. 5, we focus on the
phonon thermal conductance $G_{K}^{ph}$. We start from the minimum value of
$\gamma$ ($\gamma\simeq 0.15\Gamma$) considered in this work. We find that,
for weak electron-vibration coupling $E_{P}$ (see upper panel of Fig. 5) or in
the off-resonant regime $|\epsilon|\gg 0$ (see lower panel of Fig. 5),
$G_{K}^{ph}$ reaches its lowest value that is close to $0.04k_{B}\Gamma\simeq
16$ pW/K, the numerical value obtained when electron-vibration effects are
neglected (see the analytical estimate of $G_{K}^{ph}$ given in Eq.(40) of
Appendix A). This value corresponds only to the contribution given by the
phonon leads. We point out that this asymptotic value of $G_{K}^{ph}$ is
always larger than the values of $G_{K}^{el}$ shown in Fig. 2 (corresponding
to $E_{P}=0$). Therefore, $G_{K}^{ph}$ plays the major role in determining the
total thermal conductance $G_{K}$ for weak electron-vibration coupling.

In the lower panel of Fig. 5, we show that $G_{K}^{ph}$ always gets larger
with increasing the electron-vibration coupling $E_{P}$. Actually, the
electron-oscillator coupling gives rise to an additional damping rate whose
effect is to enhance the thermal conductivity $G_{K}^{ph}$ (see also Appendix
A). In the quasi-resonant regime, the increase of $G_{K}^{ph}$ can be also
favored by a softening of the oscillator frequency alberto . Moreover, we
notice that the conductances $G_{K}^{el}$ and $G_{K}^{ph}$ tend to acquire
similar values with increasing the electron-vibration coupling. For example,
at the resonance and $E_{P}=1.0$, $G_{K}^{el}$ and $G_{K}^{ph}$ are both of
the order of $0.05k_{B}\Gamma\simeq 20$ pW/K, therefore the total conductance
is of the order of $0.1k_{B}\Gamma\simeq 40$ pW/K, a value perfectly
compatible with the conductance of the order of $50$ pW/K measured in
hydrocarbon molecules. wang2  Finally, if one considers larger values of
$\gamma$ (for example $\gamma\simeq 0.4\Gamma$), $G_{K}^{ph}$ plays an even
more important role in $G_{K}$.

As discussed above and in Appendix B, the effects of the electron-vibration
coupling on the oscillator dynamics depends not only on the strength of the
coupling $E_{P}$, but also on the occupation of the electronic level.
Actually, the behavior of $G_{K}^{ph}$ is strongly dependent on the value of
level energy $\epsilon$. As shown in the upper panel of Fig. 5, in the quasi-
resonant case ($\epsilon=2$), the increase of $G_{K}^{ph}$ as a function of
the electron-vibration coupling $E_{P}$ is marked. Actually, for $E_{P}=2$,
the value of $G_{K}^{ph}$ is doubled. On the other hand, in the off-resonant
regime of low level occupation, the dynamics of the oscillator is poorly
influenced by the electron-vibration effects, even if $E_{P}$ is not small.
Finally, in the lower panel of Fig. 5, we have analyzed the behavior of
$G_{K}^{ph}$ as a function of the level energy $\epsilon$ for different values
of $E_{P}$. As expected, $G_{K}^{ph}$ shows the largest deviations from the
asymptotic value in the quasi-resonant case. We point out that the peak value
is practically coincident with the value of $E_{P}$, therefore, $G_{K}^{ph}$
is strongly sensitive to the renormalizations of the electron level induced by
the electron-vibration coupling.

Figure 6: (Color online) Dimensionless figure of merit $ZT$ as a function of
level energy $\epsilon$ (in units of $\hbar\Gamma$) for different values of
electron-vibration coupling $E_{P}$ (in units of $\hbar\Gamma$) at
$\gamma=0.15\Gamma$ (Upper Panel) and $\gamma=0.40\Gamma$ (Lower Panel). In
all the plots, $T=1.25\hbar\Gamma/k_{B}$ (close to room temperature), and
$\omega_{0}=0.25\Gamma$.

In Fig.6, we focus on the figure of merit $ZT$ for different values of
electron-vibration coupling $E_{P}$ at $\gamma=0.15\Gamma$ (upper panel) and
$\gamma=0.40\Gamma$ (lower panel). From the comparison with the results shown
in Figs.2 and 4, we stress that the role of the phonon thermal conductance
$G_{K}^{ph}$ is important in inducing a suppression of $ZT$. For example,
already at $E_{P}=0$, the peak value of $ZT$ is decreased by a factor of $3$
for $\gamma=0.15\Gamma$ and a factor of $5$ for $\gamma=0.40\Gamma$.
Obviously, the electron-vibration coupling provides an additional decrease.
However, in the intermediate coupling regime ($E_{P}=0.5$ in Fig.6), the
reduction of $ZT$ due to the electron-vibration coupling is not strong. Only
in the strong coupling regime, $ZT$ acquires peak values less than half unity.
Finally, we stress that, in any coupling regime, the peak values of $ZT$ are
always linked to the maxima and the minima of the Seebeck coefficient.
Summarizing, the cooperative effects of phonon leads and electron-vibration
coupling on the molecule are able to weaken the thermoelectric performance of
this kind of device. However, within a realistic weak to intermediate
electron-vibration coupling regime, the figure of merit $ZT$ is still of the
order of unity, making these devices useful for thermoelectric applications.

Figure 7: (Color online) Ratio $G_{K}^{el}/(GT)$ in units of the Lorenz number
$L=\pi^{2}k_{B}^{2}/(3e^{2})$ (Upper Panel) and dimensionless figure of merit
$ZT$ (Lower Panel) as a function of temperature (in units of
$\hbar\Gamma/k_{B}$) for different values of electron-vibration coupling
$E_{P}$ (in units of $\hbar\Gamma$). In all the plots, $\gamma=0.15\Gamma$ and
$\omega_{0}=0.25\Gamma$.

As shown in Fig.2, the temperature plays an important role in enhancing the
Seebeck coefficient $S$, and, consequently, $ZT$. Therefore, in Fig.7, we
study the temperature behavior of the response functions. In the upper panel
of Fig.7, we report the ratio $G_{K}^{el}/(GT)$ in units of the Lorenz number
$L=\pi^{2}k_{B}^{2}/(3e^{2})$ (L is about $2.44\times 10^{-8}$ W $\Omega
K^{-2}$). When this ratio is one (dot line in the upper panel of Fig.7), the
Wiedemann-Franz law is satisfied. We stress that, in this system, this law is
followed only at very low temperatures. murphy  Indeed, with increasing
temperature, already at $E_{P}=0$, the Wiedemann-Franz law is violated in our
system due to a peaked density of states. We point out that, as reported in
the upper panel of Fig.7, the violation of this law as a function of
temperature is reduced with increasing the coupling $E_{P}$.

It has been suggested that the violation of the Wiedemann-Franz law favors the
increase of the figure of merit $ZT$ in system with discrete density of
states. Mahan  Actually, as shown in the lower panel of Fig.7, the largest
$ZT$ are present for $E_{P}=0$. Besides, the figure of merit $ZT$ shows a
maximum around $1.5$ close to a temperature of the order of $2T_{A}$. This
structure is maintained with increasing the electron-vibration coupling. In
particular, for the intermediate coupling $E_{P}=0.5$, the maximum value for
$ZT$ is around $1$. Therefore, the thermoelectric performances of the junction
can be still optimized by varying the temperature.

##  V Comparison with an approach exact in the regime of low level occupation

The regime of low level occupation is the most important to get large Seebeck
coefficients $S$ and, consequently, the figure of merit. In the off-resonant
regime $|\epsilon|\gg 0$, the electronic Green functions of the model can be
exactly calculated if, as assumed in this paper, the wide band limit is used
for the leads. Mahan1 ; almb  This calculation takes fully into account the
quantum nature of the oscillator, so that it is valid also at very low
temperatures ($T\leq\omega_{0}$). We stress that, within the off-resonant
regime, the oscillator dynamics is very weakly perturbed by the effects of the
electron-vibration coupling, but it is very sensitive to the coupling to
phonon leads.

We focus on the retarted Green function in real time $G_{R}(t)$ at the average
temperature $T$

$G_{R}(t)=G_{R}^{(0)}(t)\exp{\left[-\frac{1}{2}E_{P}^{2}\phi(t)\right]},$ (18)

which is non-perturbative in the electron-vibration coupling, as it is in the
adiabatic approach. Indeed, in Eq.(18), $G_{R}^{(0)}(t)$ is the Green function
without electron-vibration effects

$G_{R}^{(0)}(t)=-\frac{i}{\hbar}\theta(t)\exp{\left[-it\left(\frac{\epsilon}{\hbar}-i\Gamma/2\right)\right]},$
(19)

with $\theta(t)$ Heaviside function, and $\phi(t)$ is the phonon term

$\phi(t)=i\int_{0}^{t}dt_{1}\int_{0}^{t}dt_{2}D(t_{1}-t_{2}),$ (20)

with $D(t)$ phonon Green function. Due to the coupling of the oscillator to
phonon leads, $D(t)$ is not the free Green function, but contains the damping
$\gamma$. At finite temperature $T$, $D(t)$ is

$D(t)=\theta(t)D_{1}(t)+\theta(-t)D_{2}(t)$ (21)

where $D_{1}(t)$ is

$D_{1}(t)=-\frac{i}{\hbar}\left[(N_{0}+1)\exp{(-i\tilde{\omega}_{0}t)}+N_{0}\exp{(i\tilde{\omega}_{1}t)}\right]$
(22)

and $D_{2}(t)$

$D_{2}(t)=-\frac{i}{\hbar}\left[(N_{0}+1)\exp{(i\tilde{\omega}_{0}t)}+N_{0}\exp{(-i\tilde{\omega}_{1}t)}\right],$
(23)

with $N_{0}$ Bose distribution function at temperature $T$,
$\tilde{\omega}_{0}=\omega_{0}-i\gamma/2$, and
$\tilde{\omega}_{1}=\omega_{0}+i\gamma/2$. Upon taking the Fourier transform
of $G_{R}(t)$ in Eq.(18), we calculate the spectral function. Using this
spectral function, we can evaluate the electron transport properties through
Eqs.(14),(15), and (16) in analogy with the adiabatic approach.

Figure 8: (Color online) Spectral function as a function of the frequency (in
units of $\Gamma$) for adiabatic approach (solid line) and low density
approach (dashed line for $\gamma=0$, dash-dotted line for $\gamma=0.15$) at
temperature $T=1.25\hbar\Gamma/k_{B}$ (Upper Panel) and temperature
$T=0.05\hbar\Gamma/k_{B}$ (Lower Panel). In all the plots,
$E_{P}=0.5\hbar\Gamma$, $\epsilon=5\hbar\Gamma$, and $\omega_{0}=0.25\Gamma$.

The spectral function calculated here will be compared with that obtained
within the adiabatic approach in Eq.(13). First, in the upper panel of Fig. 8,
we consider the spectral functions at $T=1.25\hbar\Gamma/k_{B}$, which is
close to room temperature. Moreover, we consider the off-resonant regime
$\epsilon=5$, which is very close to the minimum of the Seebeck coefficient.
The spectral weight up to $0$ (position of the chemical position) indicates
that the level occupation $n$ is small (less than $0.1$). The agreement
between the spectral functions calculated within the two approaches is
excellent. The peak positions for both approaches are at $\omega=\epsilon$ and
the widths of the curves perfectly match. The role of $\gamma$ is not
relevant, since, in any case, it is much smaller than $\Gamma$. Obviously,
with decreasing the temperature, the two approaches tend to differ. In the
lower panel of Fig.8, we have considered the worst case of very low
temperature ($T=0.05\hbar\Gamma/k_{B}$). We point out that the agreement
between the two approaches is still good. Indeed, the approach exact at low
molecule occupation slightly favors a small transfer of spectral weight at
high frequency. In any case, strong similarities in the spectral function will
give rise to analogous behaviors of electron transport properties within the
two approaches. wingreen

In the regime of low molecule occupation, the phonon thermal conductance
$G_{K}^{ph}$ at the average temperature $T$ can be calculated neglecting all
the renormalization effects due to electron-vibration coupling but retaining
all the quantum contributions wang4 ; pana :

$G_{K}^{ph}=\frac{k_{B}\gamma^{2}}{2}\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{\omega^{2}}{(\omega_{0}^{2}-\omega^{2})^{2}+(\omega\gamma)^{2}}F(\omega),$
(24)

with $F(\omega)$ given by

$F(\omega)=\left(\frac{\hbar\omega}{k_{B}T}\right)^{2}\frac{\exp{\left(\frac{\hbar\omega}{k_{B}T}\right)}}{\left[\exp{\left(\frac{\hbar\omega}{k_{B}T}\right)}-1\right]^{2}}.$
(25)

Figure 9: (Color online) Figure of merit $ZT$ as a function of temperature (in
units of $\hbar\Gamma/k_{B}$) for different values of parameters and
approaches. Solid line: low density quantum approach for $\gamma=0.15\Gamma$;
Dash line : adiabatic approach for $\gamma=0.15\Gamma$; Dash-dot line: low
density quantum approach for $\gamma=0$; Dash-double dot line : adiabatic
approach for $\gamma=0$. In the figure, $E_{P}=0.5\hbar\Gamma$, and
$\epsilon=5\hbar\Gamma$.

The quantity $RG_{K}^{ph}$ will indicate the ratio between the conductance in
Eq.(24) and that calculated in Eq. (40) of Appendix A based on the adiabatic
semiclassical approach. In the off-resonant case corresponding to
$\epsilon=5\hbar\Gamma$, the effects due to electron-vibration coupling
$E_{P}=0.5\hbar\Gamma$ on the oscillator dynamics are negligible. Due to the
quantum effects, $F(\omega)$ is smaller than unity implying that the ratio
$RG_{K}^{ph}$ shares the same behavior. In the inset of Fig. 9, we report this
ratio showing that it goes rapidly to 1 with increasing temperature. Actually,
at $T=\hbar\Gamma/k_{B}$, this ratio is already around $0.95$. Quantum
corrections are relevant only for temperatures of the order of
$\hbar\omega_{0}/k_{B}$ (much smaller than $\hbar\Gamma/k_{B}$).

Finally, we have compared the figure of merit calculated within the adiabatic
semiclassical approach with that determined in this section which retains
quantum corrections. We still focus on the off-resonant regime corresponding
to $\epsilon=5$, neglecting the effects of the electron-vibration coupling on
the phonon thermal conductance (we use Eq. (24) within the low density quantum
approach). At very low temperatures, the figure of merit is definitely small.
As reported in Fig. 9, $ZT$ rapidly gets larger with increasing temperature.
First, we have compared $ZT$ for the two approaches in the low temperature
range at $\gamma=0$. We find that, even in this regime, the agreement between
the two approaches is good. $ZT$ calculated in the adiabatic approach is
slightly larger than that obtained in the low density approach. Then, we have
considered $ZT$ for the two approaches at $\gamma=0.15\Gamma$. Therefore, we
include the contribution from the phonon thermal conductance $G_{K}^{ph}$.
Since, as shown in the inset of Fig. 9, $G_{K}^{ph}$ with quantum terms in
Eq.(24) is smaller than the adiabatic semiclassical quantity, the $ZT$
calculated within the two approaches perfectly match. Therefore, the
thermoelectric properties discussed in this paper are consistently described
also in the low temperature regime. The comparison between two approaches
cannot be extended in the regime of high temperatures since the level density
$n$ increases making the quantum low density approach less valid. For
$\epsilon>5$, the range of temperature for the comparison increases finding a
perfect agreement between the two approaches.

##  VI Conclusions

In this paper, the thermoelectric properties of a molecular junction have been
studied within the linear response regime at room temperature. In particular,
we have analyzed the role played by the phonon thermal contribution
$G_{K}^{ph}$ on the figure of merit $ZT$ in the presence of electron-vibration
coupling. The interplay between the low frequency center of mass oscillation
of the molecule and the electronic degrees of freedom has been investigated
using a non-equilibrium adiabatic approach. Parameters appropriate to $C_{60}$
molecules connected with different metallic leads have been considered. The
semiclassical $G_{K}^{ph}$ is typically of the same order of or larger than
electronic thermal conductance $G_{K}^{el}$. Both conductances are affected by
the changes in the occupation of electron levels, and they get larger with
increasing the electron-vibration coupling. Moreover, deviations from the
Wiedemann-Franz law are progressively reduced with increasing the electron-
vibration coupling. Therefore, the figure of merit $ZT$ depends appreciably on
the behavior of $G_{K}^{ph}$ and electron-vibration coupling. Indeed, for
realistic parameters of the model, $ZT$ can be substantially reduced, but it
can still have peaks of the order of unity with enhancements due to
temperature increase. Finally, we have compared the results of the adiabatic
approach with those of a formalism which is exact for low electron level
density. We have pointed out that the additional quantum effects included in
$G_{K}^{ph}$ poorly influence the thermoelectric properties in any regime of
temperatures.

The nanoscopic junction investigated in this paper is advantageous compared to
bulk or other low-dimensional structures in providing a mechanism to keep the
phonon thermal conduction low. Actually, the enhancement of $G_{K}^{ph}$ due
to the electron-vibration coupling at most provides a factor of two to a value
that is small compared to bulk conductances. The phonon thermal conductance
depends not only on the properties of the metallic leads, but also of the
tunneling barriers, shakouri  therefore, the phonon conduction can be made
negligible selecting barrier materials with low lattice thermal conductivity
(in our model this could correspond to a strong reduction of damping rate
$\gamma$). In any case, in this paper, we have pointed out that, even if one
neglects the contribution from phonon thermal conductance, the electron-
vibration coupling ($E_{P}$ in our model) is able to reduce the figure of
merit. In order to improve the thermoelectric efficiency, molecules and
metallic leads (which screen not only electron-electron but also electron-
vibration interactions on the molecule) have to be selected to ensure a weak
coupling between electronic and vibrational degrees of freedom.

Focus of the paper has been on the off-resonant electronic regime, where the
thermoelectric properties show peak values, close to room temperature.
Obviously, the effects of Coulomb local interactions are expected to be
negligible within these conditions. It would be interesting to extend the
analysis to all the electron density regimes including the electron-electron
interaction. Work in this direction is in progress. Finally, we point out that
the electron-vibration interaction investigated in this paper is linked to the
charge density injected by the external leads onto the molecule. Another
possible source of coupling could come from the renormalization of the lead-
molecule hopping integral induced by the center of mass movement. koch  Due to
the large mass of the molecules investigated in this paper, we expect that the
coupling through electron level density plays a major role.

##  Appendix A Generalized Langevin equation for the center of mass oscillator

In this Appendix, we report the derivation of the Langevin equation for the
molecule center of mass oscillator.

The resulting Langevin equation for the oscillator dynamics

$m\frac{dv}{dt}=\xi(x,t)+F_{det}(x,v)$ (26)

has the position dependent fluctuating force term $\xi(x,t)$ and the
deterministic force $F_{det}(x,v)$

$F_{det}(x,v)=-kx+F^{el}(x,v)+F_{L}^{ph}(x,v)+F_{R}^{ph}(x,v),$ (27)

where $F^{el}(x,v)$ is the force due to the effect of all the electronic
degrees of freedom, and $F_{\alpha}^{ph}(x,v)$ is the force due to the
coupling to the $\alpha$ lead phonon degrees of freedom.

Within the adiabatic regime, in Eq.(27), $F^{el}(x,v)=-\lambda N(x,v)$, where
$N(x,v)$ is the electronic level occupation with an explicit dependence on the
oscillator parameters. In order to make a self-consistent calculation, one
needs the adiabatic expansion of the level occupation $N(x,v)$:

$N(x,v)\simeq N^{(0)}(x)+vM^{(1)}(x),$ (28)

where $N^{(0)}(x)$ is the zero order ”static” term

$N^{(0)}(x)=\frac{1}{2}\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}[\sum_{\alpha}f_{\alpha}(\omega)]A(\omega,x),$
(29)

with $A(\omega,x)$ given in Eq.(7),
$f_{\alpha}(\omega)=1/(\exp{[\beta_{\alpha}(\hbar\omega-\mu_{\alpha})]}+1)$
free Fermi distribution of the $\alpha$ lead corresponding to the chemical
potential $\mu_{\alpha}$ and the temperature $T_{\alpha}$
($\beta_{\alpha}=1/k_{B}T_{\alpha}$), and $M^{(1)}(x)$

$M^{(1)}(x)=-\hbar\lambda\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}\frac{\left[\sum_{\alpha}f_{\alpha}(\omega)\right](\hbar\Gamma)^{2}I(\omega,x)}{\left[I^{2}(\omega,x)+(\hbar\Gamma)^{2}/4\right]^{3}},$
(30)

with $I(\omega,x)=\hbar\omega-\epsilon-\lambda x$, is proportional to the
first order ”dynamic” contribution through the velocity $v$ and it is
sensitive to charge fluctuations. Therefore, using Eqs.(29) and (30), one gets
$F^{el}(x,v)=-\lambda N^{(0)}(x)-A_{\lambda}(x)v$, with
$A_{\lambda}(x)=\lambda M^{(1)}(x)$ positive definite position dependent
dissipative term. In the regime investigated in this work, in Eq.(27), one has
$F_{\alpha}^{ph}(x,v)=-m\gamma_{\alpha}v$, with $\gamma_{\alpha}=\gamma/2$.
Summarizing, the deterministic force $F_{det}(x,v)$

$F_{det}(x,v)=F_{gen}(x)-A_{eff}(x)v,$ (31)

consists of a generalized force $F_{gen}(x)$

$F_{gen}(x)=-kx-\lambda N^{(0)}(x),$ (32)

and an effective position dependent damping term $A_{eff}(x)$

$A_{eff}(x)=A_{\lambda}(x)+m\gamma.$ (33)

In the adiabatic regime, exploiting the effect of the electron and phonon
environment on the slow center of mass motion, the fluctuating force
$\xi(x,t)$ in Eq.(26) is composed of three independent terms

$\xi(x,t)=\xi^{el}(x,t)+\xi_{L}^{ph}(t)+\xi_{R}^{ph}(t),$ (34)

where $\xi^{el}(x,t)$ is due to the electronic degrees of freedom such that

$\langle\xi^{el}(x,t)\rangle=0,\;\;\;\;\langle\xi^{el}(x,t)\xi^{el}(x,t^{\prime})\rangle=D_{\lambda}(x)\delta(t-t^{\prime}),$

with

$D_{\lambda}(x)=\hbar\lambda^{2}\frac{(\hbar\Gamma)^{2}}{4}\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi}\frac{\sum_{\alpha,\eta}f_{\alpha}(\omega)\left[1-f_{\eta}(\omega)\right]}{\left[I^{2}(\omega,x)+(\hbar\Gamma)^{2}/4\right]^{2}},$
(35)

$\alpha,\eta=L,R$, and $\xi_{\alpha}^{ph}(t)$ is due to the $\alpha$ phonon
lead such that

$\langle\xi_{\alpha}^{ph}(t)\rangle=0,\;\;\;\;\langle\xi_{\alpha}^{ph}(t)\xi_{\alpha}^{ph}(t^{\prime})\rangle=2K_{B}T_{\alpha}m\gamma_{\alpha}\delta(t-t^{\prime}).$

Combining the three terms, one gets a fluctuating force $\xi(x,t)$ such that

$\langle\xi(x,t)\rangle=0,\;\;\;\;\langle\xi(x,t)\xi(x,t^{\prime})\rangle=D_{eff}(x)\delta(t-t^{\prime}),$

where the effective position dependent noise term $D_{eff}(x)$ is

$D_{eff}(x)=D_{\lambda}(x)+K_{B}(T_{L}+T_{R})m\gamma,$ (36)

From the solution of the Langevin equation, the oscillator distribution
function $Q(x,v)$ and the reduced position distribution function $P(x)$ are
determined. In equilibrium conditions at temperature $T$, one has
$P_{eq}(x)=C\exp{[-\beta V_{gen}(x)]}$, with C normalization constant,
$\beta=1/(K_{B}T)$, and $V_{gen}(x)$ potential energy derived by the
generalized force $F_{gen}(x)$ in Eq.(32).

In order to calculate the thermal conductance, one can determine the
vibrational energy currents directly from the derivative of the oscillator
energy. pana  The oscillator is directly in contact with phonon leads, but
only indirectly with electron leads since there is the molecular level (see
Fig.1). In analogy with the terms in the deterministic (see Eq.(27)) and
fluctuating (see Eq.(A)) force acting on the oscillator, the total energy
current J involving the oscillator is composed of three terms brand :

$J=J^{el}_{\lambda}+J_{L}^{ph}+J_{R}^{ph},$ (37)

where $J^{el}_{\lambda}$ originates from the electron level and depends on the
electron-vibration coupling

$J^{el}_{\lambda}=\langle v\left[\xi^{el}(x,t)-A_{\lambda}(x)v\right]\rangle,$
(38)

and $J_{\alpha}^{ph}$ comes from the $\alpha$ phonon lead

$J_{\alpha}^{ph}=\langle
v\left[\xi_{\alpha}^{ph}(t)-m\gamma_{\alpha}v\right]\rangle.$ (39)

These quantities have to be evaluated along the dynamics. Once the stationary
state is reached, the energy conservation requires that the total energy
current J vanishes. Within the numerical simulations, we have not only found
the total energy conservation, but that also $J^{el}_{\lambda}$ vanishes.
Actually, the current mediated by the molecular electron level is not
effective on the stationary state. Therefore, as emphasized in the main text,
we have numerically calculated the phonon thermal conductance in the linear
regime. wang3 ; wang4

###  A.1 Weak electron-vibration coupling

In the regime of weak electron-vibration coupling, the Langevin equation is
simplified since the dissipative and fluctuating terms do not depend on the
position x. Actually, in Eq.(31), $A_{eff}(x)$ is replaced by $A_{eff}$ since
one gets $A_{\lambda}=\lambda M^{(1)}(x=0)$, with $M^{(1)}(x)$ given in
Eq.(30). Therefore, within the weak coupling, there is a simple damping rate
$\gamma_{\lambda}=A_{\lambda}/m$. Analogously, in Eq.(A), $D_{eff}(x)$ is
substituted by $D_{eff}=D_{\lambda}(x=0)$ in Eq.(35).

Within the weak-coupling regime, the Langevin equation is linear, therefore it
can be analitically solved. The distribution functions $Q(x,v)$ and P(x,v) are
Gaussian at and out of equilibrium. In particular, we have focused on the
energy conservation finding analitically that not only $J=0$, but also
$J^{el}_{\lambda}=0$ and $J_{R}^{ph}=-J_{L}^{ph}$. The quantity $G_{K}^{ph}$
has been explicitly evaluated as

$G_{K}^{ph}=\frac{K_{B}\gamma(\gamma+\gamma_{\lambda})}{2}\int_{-\infty}^{+\infty}\frac{d\omega}{2\pi}\frac{\omega^{2}}{(\omega_{0}^{2}-\omega^{2})^{2}+\omega^{2}(\gamma+\gamma_{\lambda})^{2}}.$
(40)

Therefore, there is a renormalization of $\gamma$ due to $\gamma_{\lambda}$,
but not everywhere. We stress that, as expected, the direct link of
$G_{K}^{ph}$ to $\gamma$ is still present. Moreover, as emphasized in the next
Appendix, the order of magnitude of $\gamma_{\lambda}$ is always smaller than
the values of $\gamma$ considered in this work within the weak coupling
regime. Finally, we point out that the conductance $G_{K}^{ph}$ with
$\gamma_{\lambda}=0$ represents a refencence value even when the electron-
vibration coupling is strong, but the level occupation is very low. Actually,
even in this case, the oscillator dynamics is not perturbed by the effects of
electron-vibration interaction.

##  Appendix B Electron-vibration damping rate and oscillator position
distribution

In this Appendix, we thoroughly discuss the features of the parameter
$A_{\lambda}(x)/m$ in the linear response regime at the average temperature
$T$. Moreover, in the same regime, we compare its $x$ dependence with that of
the oscillator position distribution $P(x)$. The static distribution $P(x)$ is
essentially the equilibrium distribution, so that it does not depend on the
values of $A_{\lambda}(x)/m$. However, a comparison of $x$ dependence between
$A_{\lambda}(x)/m$ and $P(x)$ will clarify the conditions under which the
electron-vibration interaction can affect the dynamics of the center of mass
oscillator. We will use the same parameters and units of the main text.
Therefore $\hbar\Gamma\simeq 20$ meV will be the energy unit and $\Gamma$ the
frequency unit. We will assume $\omega_{0}=0.25\Gamma$ and the average
chemical potential $\mu=0$. We will measure lengths in units of $\lambda/k$,
times in units of $1/\Gamma$, temperatures in units of $\hbar\Gamma/K_{B}$. In
this Appendix, we will consider the oscillator properties for $\gamma=0$ since
we are interested on the effects induced by the electron-vibration coupling.

Figure 10: (Color online) Upper panel: Oscillator position distribution $P(x)$
(in units of $k/\lambda$) as a function of position x (in units of
$\lambda/k$) for different values of level energy $\epsilon$ (in units of
$\hbar\Gamma$) at temperature $T=1.25\hbar\Gamma/k_{B}$ (close to room
temperature) and electron-vibration coupling $E_{P}=2\hbar\Gamma$. Lower
panel: Electron-vibration damping rate $A_{\lambda}(x)/m$ (in units of
$\Gamma$) as a function of position x (in units of $\lambda/k$) for different
values of level energy $\epsilon$ (in units of $\hbar\Gamma$) at
$T=1.25\hbar\Gamma/k_{B}$ (close to room temperature) and
$E_{P}=2\hbar\Gamma$.

As reported in the lower panel of Fig. 10, the peak values of
$A_{\lambda}(x)/m$ are always smaller than the values of $\gamma$ considered
in this paper ($\gamma\simeq 0.15-0.40\hbar\Gamma$) even for the intermediate
to strong value $E_{P}=2\hbar\Gamma$ of the electron-vibration coupling. This
means that the effects due to the electron-vibration coupling on the
oscillator dynamics do not typically represent a large perturbation with
respect to those induced by the coupling to phonon leads. Obviously, as
reported in the figure, the effects of the electron-vibration coupling depends
on the occupation of the electronic level. Indeed, the peak of
$A_{\lambda}(x)/m$ largely shifts passing from the quasi-resonant case
($\epsilon=2\hbar\Gamma$) to the off-resonant condition
($\epsilon=20\hbar\Gamma$). In order to better quantify the effects of
electron-vibration coupling on the oscillator dynamics, in the upper panel of
Fig. 10, we report the oscillator position distribution $P(x)$ with varying
the level energy $\epsilon$. In the quasi-resonant case
($\epsilon=2\hbar\Gamma$), the peak positions of $A_{\lambda}(x)/m$ and $P(x)$
are almost coincindent. By the way, within the units considered in this paper,
the peak position of the distribution $P(x)$ is about $-2nE_{P}/\hbar\Gamma$,
with $n$ density of the electron level. Therefore, with increasing $\epsilon$,
the density $n$ is strongly reduced, so that the peak position of $P(x)$
quickly goes to zero. In the off-resonant cases ($\epsilon=10\hbar\Gamma$ and
$\epsilon=20\hbar\Gamma$), in the linear regime considered in this Appendix,
the distribution $P(x)$ is practically the Gaussian of the free harmonic
oscillator at temperature T. On the other hand, in the off-resonant cases, the
peak positions of $A_{\lambda}(x)/m$ go toward the direction opposite to the
peaks of $P(x)$. In the regime of low occupation, the dynamics of the
oscillator is not influenced by these effects for the value
$E_{P}=2\hbar\Gamma$, which is not negligible.

Figure 11: (Color online) Upper panel: Oscillator position distribution $P(x)$
(in units of $k/\lambda$) as a function of position x (in units of
$\lambda/k$) for different values of electron-vibration coupling $E_{P}$ (in
units of $\hbar\Gamma$) at temperature $T=1.25\hbar\Gamma/k_{B}$ (close to
room temperature) and level energy $\epsilon=10\hbar\Gamma$. Lower panel:
Electron-vibration damping rate $A_{\lambda}(x)/m$ (in units of $\Gamma$) as a
function of position x (in units of $\lambda/k$) for different values of
electron-vibration coupling $E_{P}$ (in units of $\hbar\Gamma$) at temperature
$T=1.25\hbar\Gamma/k_{B}$) (close to room temperature) and level energy
$\epsilon=10\hbar\Gamma$.

In order to understand the interplay between the changes of $\epsilon$ and
$E_{P}$, in Fig. 11, we have analyzed the behavior of $P(x)$ (upper panel) and
$A_{\lambda}(x)/m$ (lower panel) for different values of $E_{P}$. We have
considered the off-resonant case $\epsilon=10\hbar\Gamma$. The distribution
$P(x)$ is practically the Gaussian of the free harmonic oscillator at
temperature T for different values of $E_{P}$ (it is so, since there is a
change in the position unit with varying $E_{P}$). On the other hand, the peak
positions of $A_{\lambda}(x)/m$ shift toward zero with increasing $E_{P}$.
Only for the large coupling $E_{P}=4\hbar\Gamma$, the superposition between
$P(x)$ and $A_{\lambda}(x)/m$ is not negligible.

###### Acknowledgements.

The authors would like to thank financial support from the project GREEN
(PON02-00029-2791179) funded by Ministero dell’Istruzione, dell’Università e
della Ricerca.

## References

  * (1) G. S. Nolas, J. Sharp and J. Goldsmid, Thermoelectrics: Basic Principles and New Materials Developments (Springer, 2010).
  * (2) A. F. Ioffe, Semiconductor Thermoelements and Thermoelectric Cooling (Infosearch Limited, London, 1957).
  * (3) A. Shakouri, Annu. Rev. Mater. Res. 41, 399 (2011).
  * (4) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 12727 (1993).
  * (5) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 16631 (1993).
  * (6) K. Koumoto and T. Mori, Thermoelectric Nanomaterials - Materials Design and Applications (Springer-Verlag, Berlin Heidelberg, 2013).
  * (7) Y.-M. Lin, X. Sun, and M. S. Dresselhaus, Phys. Rev. B 62, 4610 (2000).
  * (8) S. V. Aradhya and L. Venkataraman, Nat. Nanotechnol. 8, 399 (2013).
  * (9) Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 (2011).
  * (10) G. D. Mahan and J. Sofo, Proc. Natl. Acad. Sci. USA 93, 7436 (1996).
  * (11) P. Reddy, S.-Y. Jang, R. A. Segalman, and A. Majumdar, Science 315, 1568 (2007).
  * (12) K. Baheti, J. A. Malen, P. Doak, P. Reddy, S.-Y. Jang, T. D. Tilley, A. Majumdar, and R. A. Segalman, Nano Lett. 8, 715 (2008).
  * (13) S. K. Yee, J. A. Malen, A. Majumdar, and R. A. Segalman, Nano Lett. 11, 4089 (2011).
  * (14) C. M. Finch, V. M. Garcia-Suarez, and C. J. Lambert, Phys. Rev. B 79, 033405 (2009).
  * (15) P. Murphy, S. Mukerjee, and J. Moore, Phys. Rev. B 78, 161406(R) (2008).
  * (16) M. Galperin, A. Nitzan, and M. A. Ratner, Mol. Phys. 106, 397 (2008).
  * (17) J. Koch, F. von Oppen, Y. Oreg, and E. Sela, Phys. Rev. B 70, 195107 (2004).
  * (18) M. Leijnse, M. R. Wegewijs, and K. Flensberg, Phys. Rev. B 82, 045412 (2010).
  * (19) M. Paulsson and S. Datta, Phys. Rev. B 67, 241403(R) (2003).
  * (20) J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment (World Scientific Publishing Company, 1st edition, 2010).
  * (21) S. Datta, Lessons from Nanoelectronics: A New Perspective on Transport (World Scientific Publishing Company, 2012).
  * (22) R. Y. Wang, R. A. Segalman, and A. Majumdar, Appl. Phys. Lett. 89, 173113 (2006).
  * (23) Z. Wang, J. A. Carter, A. Lagutchev, Y. K. Koh, N.-H. Seong, D. G. Cahill, and D. D. Dlott, Science 317, 787 (2007).
  * (24) W. Lee, K. Kim, W. Jeong, L. A. Zotti, F. Pauly, J. C. Cuevas, and P. Reddy, Nature 498, 209 (2013).
  * (25) M. Galperin, A. Nitzan, and M. A. Ratner, J. Phys.: Condens. Matter 19, 103201 (2007).
  * (26) H. Park, J. Park, A. K. L. Lim, E. H. Anderson, A. P. Alivisatos, and P. L. McEuen, Nature 407, 57 (2000).
  * (27) H. Qin, A. W. Holleitner, K. Eberl, and R. H. Blick, Phys. Rev. B 64, 241302(R) (2001).
  * (28) L. Siddiqui, A. W. Ghosh, and S. Datta, Phys. Rev. B 76, 085433 (2007).
  * (29) A. Mitra, I. Aleiner, and A. J. Millis, Phys. Rev. B 69, 245302 (2004).
  * (30) B. Dong, X. L. Lei and N. J. M. Horing, Appl. Phys. Lett. 90, 242101 (2007).
  * (31) T.-F. Fang, Q.-F. Sun, and H.-G. Luo, Phys. Rev. B 84, 155417 (2011).
  * (32) M. Galperin, A. Nitzan, and M. A. Ratner, Phys. Rev. B 78, 125320 (2008).
  * (33) W. Rudzinski, J. Phys.: Condens. Matter 20, 275214 (2008).
  * (34) D. Mozyrsky, M. B. Hastings, and I. Martin, Phys. Rev. B 73, 035104 (2006).
  * (35) F. Pistolesi, Ya. M. Blanter, and I. Martin, Phys. Rev. B 78, 085127 (2008).
  * (36) R. Hussein, A. Metelmann, P. Zedler, and T. Brandes, Phys. Rev. B 82, 165406 (2010).
  * (37) A. Nocera, C.A. Perroni, V. Marigliano Ramaglia, and V. Cataudella, Phys. Rev. B 83, 115420 (2011).
  * (38) A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, and V. Cataudella, Phys. Rev. B 86, 035420 (2012).
  * (39) K.-H. Yang, Y.-L. Zhao, Y.-J. Wu, Y.-P. Wu, Phys. Lett. A 374, 2874 (2010).
  * (40) X. Zianni, Phys. Rev. B 82, 165302 (2010).
  * (41) J. Ren, J.-X. Zhu, J. E. Gubernatis, C. Wang, and B. Li, Phys. Rev. B 85, 155443 (2012).
  * (42) M. Bagheri Tagani and H. Rahimpour Soleimani, Physica B, 413, 86 (2013).
  * (43) B. C. Hsu, C.-W. Chiang, and Y.-C. Chen, Nanotechnology 23, 275401 (2012).
  * (44) H. Haug and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, (Springer, Berlin, 2008).
  * (45) X. Lu, M. Grobis, K. H. Khoo, S. G. Louie, and M. F. Crommie, Phys. Rev. Lett. 90, 096802 (2003).
  * (46) L. H. Yu and D. Natelson, Nano Lett. 4, 79 (2005).
  * (47) J. Mravlje, A. Ramsak, and T. Rejec, Phys. Rev. B 74, 205320 (2006).
  * (48) C. Kittel, Introduction to Solid State Physics (John Wiley & Sons., 2004, 8th ed.).
  * (49) U. Weiss, Quantum Dissipative Systems, (World Scientific Publishing Company, 2008, 3rd edition).
  * (50) A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, G. Cantele, and V. Cataudella, Phys. Rev. B 87, 155435 (2013).
  * (51) C. A. Perroni, A. Nocera, and V. Cataudella, Europhys. Lett. 103, 58001 (2013).
  * (52) C. A. Perroni, F. Romeo, A. Nocera, V. Marigliano Ramaglia, R. Citro, and V. Cataudella, arXiv:1307.6834.
  * (53) R. L. Honeycutt, Phys. Rev. A 45, 600 (1992).
  * (54) R. L. Honeycutt, Phys. Rev. A 45, 604 (1992).
  * (55) J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007).
  * (56) S. Jezouin, F. D. Parmentier, A. Anthore, U. Gennser, A. Cavanna, Y. Jin, and F. Pierre, Science 342, 601 (2013).
  * (57) G. D. Mahan, Many-Particle Physics (Springer, 3rd edition, 2000)
  * (58) C. O. Almbladh and P. Minnhagen, Phys. Rev. B 17, 929 (1978).
  * (59) N. S. Wingreen, K. W. Jacobsen, J. W. Wilkins, Phys. Rev. Lett. 61, 1396 (1988).
  * (60) G. Y. Panasyuk, G. A. Levin, and K. L. Yerkes, Phys. Rev. E 86, 021116 (2012).
  * (61) J.-T. Lu, M. Brandbyge, P. Hedegard, T. N. Todorov, and D. Dundas, Phys. Rev. B 85, 245444 (2012).
  * (62) J.-S. Wang, J. Wang, and J. T. Lu, Eur. Phys. J. B 62, 381 (2008).