# Interplay between electron-electron and electron-vibration interactions on
the thermoelectric properties of molecular junctions

C. A. Perroni  CNR-SPIN and Dipartimento di Fisica, Universita’ degli Studi di
Napoli Federico II ,
Complesso Universitario Monte S. Angelo, Via Cintia, I-80126 Napoli, Italy  D.
Ninno  CNR-SPIN and Dipartimento di Fisica, Universita’ degli Studi di Napoli
Federico II ,
Complesso Universitario Monte S. Angelo, Via Cintia, I-80126 Napoli, Italy
IMAST S.c.ar.l.-Technological District on Engineering of polymeric and
composite Materials and Structures, Piazza Bovio 22, I-80133 Napoli, Italy  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 linear thermoelectric properties of molecular junctions are theoretically
studied close to room temperature within a model including electron-electron
and electron-vibration interactions on the molecule. A nonequilibrium
adiabatic approach is generalized to include large Coulomb repulsion through a
self-consistent procedure and applied to the investigation of large molecules,
such as fullerenes, within the Coulomb blockade regime. The focus is on the
phonon thermal conductance which is quite sensitive to the effects of strong
electron-electron interactions within the intermediate electron-vibration
coupling regime. The electron-vibration interaction enhances the phonon and
electron thermal conductance, and it reduces the charge conductance and the
thermopower inducing a decrease of the thermoelectric figure of merit. For
realistic values of junction parameters, the peak values of the thermoelectric
figure of merit are still of the order of unity since the phonon thermal
conductance can be even smaller than the electron counterpart.

##  I Introduction

The direct conversion of temperature differences to electric voltage and vice
versa take place in solid state systems. These thermoelectric effects can be
strong enough in some semiconducting materials to allow either the fabrication
of devices converting wasted heat into electrical energy or the realization of
solid-state coolers nolas ; ioffe . A fundamental parameter to quantify the
energy conversion efficiency is the dimensionless figure of merit
$ZT=GS^{2}T/G_{K}$, where $G$ is the electrical conductance, $S$ the
thermopower, $T$ the absolute temperature, and $G_{K}=G_{K}^{el}+G_{K}^{ph}$
is the total thermal conductance, with $G_{K}^{el}$ and $G_{K}^{ph}$ electron
and phonon thermal conductance, respectively. Indeed, in order to improve the
efficiency, mutually contrasting transport properties of the same material
have to be optimized. For instance, in metals, $ZT$ is typically limited by
the Wiedemann-Franz law kittel . Large effort is currently made in material
science to get bulk values of $ZT$ larger than $1$ and to use solid state
systems for actual thermoelectric devices nolas ; shakouri ; biswas .

Recently, the possibility of controlling materials at the nanoscale has been
exploited to optimize the thermoelectric efficiency. shakouri ; koumoto ;
sothmann  For example, a maximum $ZT\simeq 2.4$ has been observed at room
temperature in a thin-film thermoelectric device venka . High values of $ZT$
have been reported in quantum dot superlattices harman  and in semiconductor
nanowires hochbaum , where phonon confinement can lead to a lower phonon
thermal conductance dressel ; murphy1 . Actually, a significant reduction in
lattice thermal conductivity is considered as the main route for having high
$ZT$ in low-dimensional materials venka1 . The improvement of thermoelectric
efficiency can also derive from the discreteness of energy levels in
nanostructures resulting into a violation of the Wiedemann-Franz law Mahan .
Finally, in nanoscopic Coulomb-coupled systems, the thermoelectric properties
can be optimized by exploiting the Coulomb blockade regime and changing the
gate voltage sothmann .

Molecular devices can be efficient for conversion of heat into electric energy
since both phonon and electron properties can contribute to increase the
thermoelectric figure of merit $ZT$ arad ; dubi . Indeed, the emerging field
of molecular thermoelectrics has attracted a lot of attention in recent years
majum1 ; majum2 ; finch ; murphy ; galperin ; koch ; leijnse . The
thermoelectric properties of molecular junctions are also interesting in that
they can provide useful informations on charge and energy transport otherwise
difficult to obtain, such as the type of carriers (electros/holes) dominating
the transport majum1 ; majum2 ; datta ; cuevas ; datta1 . Measurements of
thermoelectric properties have been performed in junctions with fullerene
($C_{60}$) majum2  finding a high value of the molecular thermopower ($S$ of
the order of $-30$ $\mu V$/K). In these experiments, three different metallic
electrodes (platinum, gold, and silver) have been considered achieving a more
controllable alignment between Fermi level and molecular orbitals (whose
energy separation is still of the order of $0.5$ eV). However, the application
of a gate voltage remains elusive in these kinds of measurements. Moreover,
heat transport in molecular devices remain poorly characterized due to
experimental challenges dubi ; wang1 ; wang2 ; meier  or limited to a range
where transport is elastic cuevas1  .

In molecular junctions, intramolecular electron-electron and electron-
vibration interactions typically constitute the largest energy scales
affecting the thermoelectric properties. cuevas ; galperin1 ; costi  Moreover,
the center of mass oscillation of the molecule Park , or thermally induced
acoustic phonons Qin  can be an additional source of coupling between
electronic and vibrational degrees of freedom. The effects of intramolecular
interactions on the transport properties have been studied in the regime of
linear response and fully out-of-equilibrium by different theoretical tools
cuevas ; galperin1 . The thermopower $S$ and the thermoelectric figure of
merit $ZT$ have been found to be sensitive to the strength of intramolecular
interactions koch ; leijnse ; galperin ; yang ; zianni ; aharon ; liu ; ren ;
tagani ; arra . However, the phonon thermal contribution $G_{K}^{ph}$ to the
figure of merit $ZT$ has been calculated only at a perturbative level of the
electron-vibration coupling hsu .

In devices with large molecules or carbon nanotube quantum dots cava1 , a
nonequilibrium adiabatic approach has been introduced for spinless electrons
exploiting the low energy of the relevant vibrational degrees of freedom mozy
; pistol ; hussein ; alberto ; alberto1 . This method is semiclassical for the
vibrational dynamics, but it is valid for arbitrary strength of electron-
vibration coupling. Within this approach, we have recently implemented a self-
consistent calculation for electron and phonon thermal conductance focusing on
the effects of the electron-vibration coupling perroni2 .

In this paper, we have studied the thermoelectric properties of a molecular
junction with electron-electron and electron-vibration interactions within the
linear response regime focusing on a self-consistent calculation of the phonon
thermal conductance $G_{K}^{ph}$ close to room temperature. The nonequilibrium
adiabatic approach is generalized to treat finite strong Coulomb interactions
within a 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 within the Coulomb blockade regime. Parameters appropriate
for junctions with $C_{60}$ molecules are considered. We have found that,
within the intermediate electron-vibration coupling regime, the effects of
electron-electron interactions can enhance $G_{K}^{ph}$, which, as a function
of the gate voltage, acquires a behavior similar to that of electron thermal
conductance. The electron-vibration interaction induces an increase of the
phonon and electron thermal conductance, and, at the same time, a decrease of
both the charge conductance and the thermopower. The overall effect is a
reduction of the thermoelectric figure of merit. Interestingly, for realistic
parameters of the model, the peak values of $ZT$ are still of the order of
unity. This effect is ascribed to the magnitude of the phonon thermal
conductance which can be smaller than the electronic counterpart in a large
range of gate voltages.

The paper is organized as follows. In Sec. II, the model of molecular junction
is proposed. In Sec. III, the adiabatic approach generalized for strong local
Coulomb interactions is explained. In Sec. IV, the results within the
adiabatic approach are discussed. The paper is closed by Appendix A, where the
comparison between different treatments of the large Coulomb repulsion is made
within the Coulomb blockade regime.

##  II Molecular junction model

In this paper, we analyze the Anderson-Holstein model, which is a reference
model for molecular junctions. cuevas ; Haug  The molecule is modeled as a
single electronic level locally interacting with a single vibrational mode. In
junctions with $C_{60}$ molecules, attention can be focused on a molecular
electronic orbital which is sufficiently separated in energy from other
orbitals. perroni2 ; Lu ; natelson ; mra  In this paper, we will consider
model parameters appropriate to $C_{60}$ molecular junctions.

Figure 1: (Color online) Sketch of the molecular junction studied in this
work. The straight 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$). The Hot Left Lead and the Cold Right Lead are kept at
chemical potential $\mu_{L}=eV/2$, temperature $T_{L}=T+\Delta T/2$ and
chemical potential $\mu_{R}=-eV/2$, temperature $T_{R}=T-\Delta T/2$,
respectively, with $e$ modulus of the electron charge, $V$ bias potential,
$\mu=0$ average chemical potential, $T$ average temperature. The term $U$
indicates the presence of electron-electron interactions, while the term
$E_{P}$ electron-vibration interactions on the molecule.

The Hamiltonian $\hat{H}$ is given by

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

where the Hamiltonian ${\hat{H}}_{el}$ takes into account the electronic
degrees of freedom of the leads and the molecule, ${\hat{H}}_{ph}$ the
vibrational degrees of freedom of the leads and the molecule, and
${\hat{H}}_{int}$ the coupling between electronic and vibrational degrees of
freedom (see Fig. 1 for a sketch of the molecular junction model).

The electronic Hamiltonian ${\hat{H}}_{el}$ of Eq. (1) is

$\displaystyle\hat{H}_{el}$ $\displaystyle=$
$\displaystyle\epsilon\sum_{\sigma}{\hat{n}}_{\sigma}+U{\hat{n}}_{\uparrow}{\hat{n}}_{\downarrow}+\sum_{q,\alpha,\sigma}\varepsilon_{q,\alpha}{\hat{n}_{q,\alpha,\sigma}}+$
(2)
$\displaystyle\sum_{q,\alpha,\sigma}\left(V_{q,\alpha}{\hat{c}^{{\dagger}}_{q,\alpha,\sigma}}{\hat{d}}_{\sigma}+h.c.\right),$

where the molecular electronic level has energy $\epsilon$, the $\sigma$ spin
electron density operator is
${\hat{n}}_{\sigma}={\hat{d}^{{\dagger}}}_{\sigma}{\hat{d}}_{\sigma}$, with
${\hat{d}^{{\dagger}}}_{\sigma}({\hat{d}}_{\sigma})$ creation (annihilation)
$\sigma$ spin electron operator 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 Coulomb repulsion on the molecule is simulated with a
Hubbard term $U$, which gives an energy penalty for electron occupations with
spin $\uparrow$ and $\downarrow$. cuevas  The lead density operator is
${\hat{n}_{q,\alpha,\sigma}}={\hat{c}^{{\dagger}}_{q,\alpha,\sigma}}{\hat{c}}_{q,\alpha,\sigma}$,
where the operators
${\hat{c}^{{\dagger}}_{q,\alpha,\sigma}}({\hat{c}}_{q,\alpha,\sigma})$ create
(annihilate) electrons with momentum $q$, spin $\sigma$, and energy
$\varepsilon_{q,\alpha}=\xi_{q,\alpha}-\mu_{\alpha}$ in the left ($\alpha=L$)
or right ($\alpha=R$) free metallic leads, with $\mu_{\alpha}$ chemical
potential of the lead $\alpha$ in equilibrium at the temperature $T_{\alpha}$.
We consider the temperatures $T_{L}=T+\Delta T/2$ and $T_{R}=T-\Delta T/2$,
with $T$ average temperature. Moreover, we fix the chemical potentials
$\mu_{L}=eV/2$ and $\mu_{R}=-eV/2$, with $e$ modulus of the electron charge,
$V$ bias potential, and average chemical potential $\mu=0$. 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 around 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 the following,
we consider the symmetric configuration: $\Gamma_{L}=\Gamma_{R}=\Gamma/2$. In
junctions with $C_{60}$ molecules, $\hbar\Gamma$ has been estimated to be of
the order of $20$ meV. natelson ; mra  Even if the local Coulomb repulsion is
reduced by the screening of the electrodes, the energy $U$ is expected to be
at least one order of magnitude larger than $\hbar\Gamma$. natelson ; mra

The center of mass mode can be considered as the relevant vibrational mode of
the molecule. perroni2  Indeed, experiments have evidenced a coupling between
the center of mass mode and the electron dynamics in junctions with $C_{60}$
molecules. Park  In Eq.(1), the Hamiltonian ${\hat{H}}_{ph}$ describes the
vibrations of the slow 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
operators, respectively, $M$ is the total large mass, $k$ is the effective
spring constant, with frequency $\omega_{0}=\sqrt{k/M}$. 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$. The left and right phonon leads will be considered as
thermostats in equilibrium at the same temperatures $T_{L}$ and $T_{R}$,
respectively, of the electron leads. Finally, 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}$. For large molecules, the center of mass
mode has a low frequency $\omega_{0}$ which is typically smaller than the
Debye frequency $\omega_{D}$ of the metallic leads ($\hbar\omega_{D}\simeq
15-20$ meV for metals like silver, gold, and platinumkittel ). For example,
$\hbar\omega_{0}\simeq 5$ meV in $C_{60}$ junctions Park , hence
$\omega_{0}\simeq 0.25\Gamma$. Therefore, for large molecules, the adiabatic
regime is valid for the center of mass oscillator: $\omega_{0}<<\Gamma$ and
$\omega_{0}<<\omega_{D}$. Within this regime, the effect of the $\alpha$
phonon lead on the center of mass mode provides a constant damping rate
$\gamma_{\alpha}$. weiss  In analogy with the electronic model, we consider
the symmetric configuration: $\gamma_{L}=\gamma_{R}=\gamma/2$. For junctions
with $C_{60}$ molecules and leads of Ag, Au, and Pt, $\hbar\gamma\simeq 3-8$
meV, therefore $\gamma$ is of the same order of $\omega_{0}$ ($\gamma\simeq
0.15-0.40\Gamma$) perroni2 .

Finally, the interaction term ${\hat{H}}_{int}$ in the Anderson-Holstein model
of Eq.(1) is provided by a linear coupling between the total electron density
on the molecule, ${\hat{n}}=\sum_{\sigma}{\hat{n}}_{\sigma}$, and the
$\hat{x}$ operator of the center of mass:

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

where $\lambda$ is the electron-vibration coupling constant. In the following,
the electron-vibration interaction will be described in terms of the coupling
energy $E_{P}=\lambda^{2}/(2k)$.

In this paper, $\hbar\Gamma\simeq 20$ meV will be the energy unit ($\Gamma$
the frequency unit, $1/\Gamma$ the time unit). We will measure lengths in
units of $2\lambda/k$, and temperatures in units of $\hbar\Gamma/k_{B}$, with
$k_{B}$ Boltzmann constant (the room temperature is of the order of $1.25$ in
these units).

##  III Adiabatic approach within the Coulomb blockade regime

The focus of this paper is on charge and heat transport properties close to
room temperature, therefore for parameters appropriate to the Coulomb blockade
regime: $\hbar\omega_{0}\ll\hbar\omega_{D}\simeq\hbar\Gamma\leq k_{B}T\ll U$,
with $U>10\hbar\Gamma$. Besides, the electron-vibration coupling is not weak,
but it is estimated to be in the intermediate regime: $\hbar\omega_{0}\leq
E_{P}\simeq\hbar\Gamma$. Since $\hbar\omega_{0}$ is the lowest energy scale,
the dynamics of the slow center of mass can be treated as classical. In the
following, the position and the momentum of the oscillator will be indicated
by the c-numbers $x$ and $p$, respectively. The parameter regime appropriate
to these junctions requires a generalization of the adiabatic approach to the
physical situation where the Coulomb interaction is finite and large.
Recently, the adiabatic approach has been combined with a treatment of
electron-electron interactions within a slave-boson approach schiro  which is
valid only in the limit of infinite local Coulomb repulsion for energies close
to the chemical potential and low temperatures perroni3 .

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

The electronic dynamics turns out to be equivalent to that of an adiabatically
slow level with energy $E_{0}(t)=\epsilon+\lambda x(t)$ within the Coulomb
blockade regime. fazio ; mucciolo

At the zero order of the adiabatic expansion, the electronic quantities can be
calculated considering an energy level with a fixed oscillator position $x$.
The effects of the strong Coulomb repulsion are treated inserting the first
self-energy correction upon the atomic limit. Haug  Therefore, for the
paramagnetic solution, the level spectral function $A_{0}(\omega,x)$ at zero
order of the adiabatic expansion becomes

$\displaystyle A_{0}(\omega,x)$ $\displaystyle=$
$\displaystyle[1-\rho(x)]\frac{\hbar\Gamma}{(\hbar\omega-\epsilon-\lambda
x)^{2}+(\hbar\Gamma)^{2}/4}+$ (6)
$\displaystyle\rho(x)\frac{\hbar\Gamma}{(\hbar\omega-\epsilon-\lambda
x-U)^{2}+(\hbar\Gamma)^{2}/4},$

where $\rho(x)$ is the level density per spin self-consistently calculated at
fixed position $x$ through the following integral

$\rho(x)=\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi
i}G^{<}_{0}(\omega,x),$ (7)

with the lesser Green function $G^{<}_{0}(\omega,x)$

$G^{<}_{0}(\omega,x)=\frac{i}{2}[f_{L}(\omega)+f_{R}(\omega)]A_{0}(\omega,x),$
(8)

and
$f_{\alpha}(\hbar\omega)=1/(\exp{[\beta_{\alpha}(\hbar\omega-\mu_{\alpha})]}+1)$
Fermi distribution of the lead $\alpha$ ($\beta_{\alpha}=1/k_{B}T_{\alpha}$).
Actually, the spectral function is characterized by a double peak structure
that, for large $U$, is robust against the effects of electron-vibration
coupling which tend to shift and enlarge the single peaks (the single peak
width increases by a factor of the order of $E_{P}$).

In Appendix A, we compare the spectral function of this treatment for strong
Coulomb repulsion with that of another approach which retains additional self-
energy corrections upon the atomic limit in the absence of electron-vibration
coupling. Haug  For large $U$ and room temperature, the approach considered
here is very accurate, therefore, it represents an optimal starting point for
the adiabatic expansion. In this paper, we will study different properties
varying the electronic level occupation. In our model, these variations can be
controlled changing the molecule level energy $\epsilon$ with respect to the
leads chemical potential (average chemical potential $\mu=0$ in this work). In
Appendix A, we report the molecular electron occupation $N$ as a function of
level energy $\epsilon$ showing the typical profiles of the Coulomb blockade.
In particular, the following energies are relevant: $\epsilon=-U/2$ (close to
half-filling $N=1$), $\epsilon=-U$ (transition from level occupation $N=1$ to
$N=2$), $\epsilon=0$ (from level occupation $N=1$ to $N=0$).

Within the adiabatic approach, one can determine the electronic Green
functions and generic electronic quantities making an expansion on the small
oscillator velocity $v=p/m$. In the absence of electron-electron interactions,
the adiabatic expansion can be determined for any strength of electron-
vibration coupling. alberto ; alberto1 ; alberto2 ; perroni ; perroni1  In
this paper, an approach is devised for the case of strong Coulomb repulsion in
order to include the effects of electron-vibration interaction within the
realistic intermediate coupling regime. Actually, the approach used in this
paper is valid as long as the two peaks characteristic of Coulomb blockade can
be resolved, therefore for the physical regime $E_{P}\ll U$. In the next
subsection, we will use the adiabatic expansion of the level occupation to
derive the motion equation of the slow center of mass oscillator in a self-
consistent way.

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

The effect of the molecule electron degrees of freedom and of the phonon baths
in the leads gives rise to the following generalized Langevin equation for the
slow center of mass

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

which 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,$ (10)

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

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

with $F_{\lambda}(x)=-2\lambda\rho(x)$ induced by the electron-vibration
coupling, and, as a result of the adiabatic expansion, a dissipative force
with an effective position dependent positive definite term $A_{eff}(x)$

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

with $A_{\lambda}(x)$

$A_{\lambda}(x)=2\hbar\lambda^{2}\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi
i}G^{<}_{0}(\omega,x)\left[\frac{\partial
A_{0}(\omega,x)}{\partial(\hbar\omega)}\right]$ (13)

due to the electron-vibration interaction. The fluctuating force $\xi(x,t)$ in
Eq.(9) 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,$ (14)

with $D_{\lambda}(x)$

$D_{\lambda}(x)=2\hbar\lambda^{2}\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi
i}G^{<}_{0}(\omega,x)G^{>}_{0}(\omega,x)$ (15)

determined by the electron-vibration coupling and the greater Green function
$G^{>}_{0}(\omega,x)$

$G^{>}_{0}(\omega,x)=-\frac{i}{2}[2-f_{L}(\omega)-f_{R}(\omega)]A_{0}(\omega,x).$
(16)

It is worthwhile pointing out that, in equilibrium conditions at temperature
$T=T_{\alpha}$ and chemical potential $\mu=\mu_{\alpha}=0$, the adiabatic
procedure gives rise to a generalized fluctuation-dissipation relation
$D_{eff}(x)=2k_{B}TA_{eff}(x)$ valid for each fixed position $x$.

The solution of the Langevin equation (9) represents a central step for this
work. This equation has been numerically solved 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)$ are
determined allowing to evaluate static quantities relative to the center of
mass oscillator. Moreover, these distribution functions will allow to make the
average of an electronic observable $O(x,v)$ dependent on oscillator
parameters.

Figure 2: (Color online) Electron-vibration induced damping rate
$\gamma_{\lambda}$ in units of $\Gamma$ (Upper Panel) and reduced position
distribution function $P$ in units of $k/2\lambda$ (Lower Panel) as a function
of oscillator position x (in units of $2\lambda/k$) for different values of
level energy $\epsilon$ (in units of $\hbar\Gamma$). In all the plots,
$U=20\hbar\Gamma$, $E_{P}=\hbar\Gamma$, and temperature
$T=1.25\hbar\Gamma/k_{B}$ (close to room temperature).

Before going to the section about results, we discuss the features of the
electron-vibration induced damping rate
$\gamma_{\lambda}(x)=A_{\lambda}(x)/m$, with $A_{\lambda}(x)$ given in
Eq.(13). The magnitude of $\gamma_{\lambda}(x)$ always gets enhanced with
increasing the electron-vibration coupling $E_{P}$. However, as reported in
the upper panel of Fig. 2, even for the intermediate coupling $E_{P}=1$, the
peak values of $\gamma_{\lambda}(x)$ are always smaller than realistic values
of the lead induced damping rate $\gamma$ ($\gamma=0.15$ will be considered in
this paper). This implies 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 behavior of $\gamma_{\lambda}(x)$
strongly depends on the occupation of the electronic level. We point out that,
in contrast to the spinless case analyzed in a recent paper, perroni2
$\gamma_{\lambda}(x)$ shows a double-peak behavior due to the effect of the
strong Hubbard interaction. Moreover, as reported in the upper panel of Fig.
2, the peaks of $\gamma_{\lambda}(x)$ largely shift passing from the quasi
half-filling case (close to $\epsilon=-10=-U/2$, state with flat occupation)
to conditions out of half-filling (close to $\epsilon=-20=-U$ and
$\epsilon=0=\mu$, state with strong density fluctuations). The self-consistent
calculation of $\gamma_{\lambda}(x)$ provides a direct signature of the strong
local interaction since it is determined by the adiabatic expansion of the
electron occupation.

A comparison of the $x$ dependence between $\gamma_{\lambda}(x)$ and the
calculated oscillator position distribution $P(x)$ will clarify the conditions
under which the electron-vibration interaction can affect the dynamics of the
center of mass oscillator. Therefore, in the lower panel of Fig. 2, we report
the distribution $P(x)$ with varying the level energy $\epsilon$. We notice
that, apart from the shift of the peaks, close to room temperature, the
distribution $P(x)$ is practically the Gaussian of the free harmonic
oscillator at temperature T for any value of the level energy $\epsilon$. In
the quasi half-filled case ($\epsilon=-10$), the peak positions of
$\gamma_{\lambda}(x)$ and $P(x)$ are well separated. Therefore, one expects
that, in this regime, the effects of the electron-vibration coupling on the
oscillator dynamics are weak. We stress that, within the self-consistent
procedure used in this work, the peak of the $P(x)$ directly signs the level
occupation being close to $-N/2$ within the units used in this paper.
Actually, for $\epsilon=-10$, the value close to $-0.5$ of the peak of $P(x)$
is fully compatible with the half-filled case $N=1$. On the other hand, for
$\epsilon=-20$, the peak position of $P(x)$ shifts towards lower values close
to $-0.75$ ($N\simeq 1.5$), and, for $\epsilon=0$, to $0.25$ ($N\simeq 0.5$).
We point out that, for $\epsilon=-20$, the first peak of $\gamma_{\lambda}(x)$
is close to $x=0$, while, for $\epsilon=0$, the second peak of
$\gamma_{\lambda}(x)$ strongly overlaps with the position distribution $P(x)$.
Therefore, out of half-filling, the effects of the electron-vibration coupling
can affect the oscillator dynamics. In contrast with the spinless case,
perroni2  these effects are present not only close to $\epsilon=\mu=0$, but
also to $\epsilon=-U=-20$, as a result of the strong Coulomb interaction.
Therefore, as discussed in detail in the next section, the complex interplay
between electron-electron and electron-vibration interactions opens an entire
energy region where the phonon heat transport can be enhanced.

##  IV Results

In this paper, we will discuss linear response transport properties trying to
clarify the role of the electron-electron and electron-vibration interactions.
In the next subsections, we will analyze the phonon heat transport, the
electronic spectral function, the charge and electronic heat transport, and
thermoelectric figure of merit. In the following, we will assume
$\omega_{0}=0.25\Gamma$, and $\gamma=0.15\Gamma$ (larger values of $\gamma$
have been considered in a recent paper perroni2 ).

###  IV.1 Phonon heat transport

In this subsection, we will focus on the phonon thermal conductance
$G_{K}^{ph}$ calculated within the linear response regime around temperature T
as

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

with $J_{\alpha}^{ph}$ current from the $\alpha$ phonon lead. wang3 ; perroni2

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

The conductance $G_{K}^{ph}$ is expected to be mostly sensitive to the
coupling of the center of mass mode to the phonons of metallic leads through
the damping rate $\gamma$ ($\gamma=0.15\Gamma$ in this work) which is
typically larger than the peak values of electron-vibration induced damping
rate $\gamma_{\lambda}(x)$. As shown in Fig.3, in the regime of weak electron-
vibration coupling $E_{P}$, low level occupation ($\epsilon\gg 0$), and double
level occupation ($\epsilon\ll-U$), $G_{K}^{ph}$ is close to $0.04$
$k_{B}\Gamma$ ($k_{B}\Gamma$ is about $419.8$ pW/K for $\hbar\Gamma\simeq 20$
meV), a numerical value coincident with an analytical estimate of $G_{K}^{ph}$
given in a recent paper perroni2 . This asymptotic value corresponds to the
contribution given by the only phonon leads neglecting the effects of
electron-electron and electron-vibration interactions on the molecule.

In Fig. 3, we show that $G_{K}^{ph}$ always gets larger with increasing the
electron-vibration coupling $E_{P}$. Moreover, this increase of $G_{K}^{ph}$
strongly depends on the value of level energy $\epsilon$. In contrast with the
spinless case (reported for comparison in Fig. 3 at $E_{P}=1$), we stress that
the enhancement of $G_{K}^{ph}$ takes place not only close to $\epsilon\simeq
0$, but also to $\epsilon\simeq-U$. Therefore, the distance between the peaks
of the phonon thermal conductance is controlled by the energy scale $U$. The
peak values are almost coincident (although slightly smaller than the peak
value of the spinless case), and, at $E_{P}=1$, they are of the order of
$0.05k_{B}\Gamma\simeq 20$ pW/K. Therefore, the calculated $G_{K}^{ph}$ is in
very good agreement with the thermal conductance of the order of a few $10$
pW/K measured for molecules anchored to gold wang2 ; meier . In any case, due
to the strong electron-electron interactions, $G_{K}^{ph}$ can be enhanced in
a new large energy region. On the other hand, for $\epsilon\simeq-U/2$,
$G_{K}^{ph}$ is poorly influenced by the electron-vibration effects even if
$E_{P}$ is not small getting a value close to the asymptotic one. From this
analysis emerges that the complex enhancement of the phonon thermal
conductance $G_{K}^{ph}$ as a function of the electron-electron and electron-
vibration interactions can be mostly ascribed to the properties of additional
electron-vibration induced damping rate $\gamma_{\lambda}(x)$ discussed in the
previous section.

###  IV.2 Electronic spectral function

From the solution of the Langevin equation, one can make the average of an
electronic observable $O(x,v)$ over the oscillator distribution function.
First, we discuss the features of the electronic spectral function which is at
the basis of the thermoelectric properties analyzed in the next subsection.

The electronic spectral function $A(\omega)$ is evaluated making the average
of the function $A_{0}(\omega,x)$ in Eq.(6) over $P(x)$:

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

In this section, the spectral function will be discussed in equilibrium
conditions at temperature $T$ ($V=0$ and $\Delta T=0$). We recall that, in
Appendix A, the features of the spectral function are discussed in the absence
of electron-vibration coupling. Actually, the spectral function is
characterized by a structure with two peaks separated by an energy of the
order of $U$, and it is strongly dependent on the value of the level energy
$\epsilon$.

Figure 4: (Color online) Spectral function (in units of $1/\hbar\Gamma$) as a
function of frequency $\omega$ (in units of $\Gamma$) at level energy
$\epsilon=-8\hbar\Gamma$ (Upper Panel) and $\epsilon=8\hbar\Gamma$ (Lower
Panel) for different values of $E_{P}$ (in units of $\hbar\Gamma$). In all the
plots, $U=16\hbar\Gamma$, $T=1.25\hbar\Gamma/k_{B}$ (close to room
temperature).

In this subsection, we analyze the behavior of the spectral function with
varying the electron-vibration coupling $E_{P}$ at a fixed value of Hubbard
energy $U$. In the upper panel of Fig. 4, we show the spectral function for
different values of the electron-vibration coupling in the half-filled case
$\epsilon=-8=-U/2$ (level occupation $N=1$). For comparison, we report the
spectral function relative to the case where electron-electron and electron-
vibration interactions are neglected (indicated as Free in the figure). We
point out that there is a strong transfer of spectral weight for the double
peak structure toward low frequencies with increasing $E_{P}$. In addition to
the shifts of the peaks, the electron-vibration coupling tends to reduce the
height of the peaks and to enlarge them. Actually, the single peaks increase
their width by a factor of order of $E_{P}$. We stress that, for realistic
values of the coupling $E_{P}$, the two Hubbard peaks do not overlap,
therefore the double peak structure due to the large $U$ is quite robust to
the effects of electron-vibration coupling. Finally, we notice that, in the
spinless case (reported for comparison in Fig. 4 at $E_{P}=0.5$), the spectral
function has a single peak, and it is quite sensitive to the effects of the
electron-vibration coupling.

As shown in the lower panel of Fig. 4, a different behavior takes place in the
regime of low level occupation ($\epsilon=8$ in the figure). For the
considered values of $E_{P}$, the spectral function gets enlarged, but its
peak position is quite rigid. Moreover, the differences with the spinless case
are completely negligible. Even in the presence of electron-vibration coupling
$E_{P}$, the behavior of the spectral function is different in the regime of
half-filling and of low or high level occupation.

###  IV.3 Charge and electronic heat transport, and thermoelectric figure of
merit

In this subsection, the focus will be on the regime of linear response around
the average chemical potential $\mu=0$ and temperature T ($\Delta T\rightarrow
0^{+}$, $V\rightarrow 0^{+}$). We will evaluate the electronic conductance $G$

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

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

$G_{S}=\left(\frac{2e}{\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(\hbar\omega)}{\partial(\hbar\omega)}\right].$ (20)

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

$G_{Q}=\left(\frac{2}{\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(\hbar\omega)}{\partial(\hbar\omega)}\right].$ (21)

The total thermal conductance $G_{K}=G_{K}^{el}+G_{K}^{ph}$ makes feasible the
evaluation of the figure of merit $ZT=GS^{2}T/G_{K}$. 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 efficiency.

Figure 5:  (Color online) Electron conductance G in units of $2e^{2}/h$ (Upper
Left Panel), Seebeck coefficient $S$ in units of $k_{B}/e$ (Upper Right
Panel), electron thermal conductance $G_{K}^{el}$ in units of $2k_{B}\Gamma$
(Lower Left Panel), and figure of merit $ZT^{el}$ (Lower Right Panel) as a
function of the level energy $\epsilon$ (in units of $\hbar\Gamma$) for
different values of electron-vibration $E_{P}$ (in units of $\hbar\Gamma$). In
all the plots, $U=20\hbar\Gamma$, $\gamma=0$, and $T=1.25\hbar\Gamma/k_{B}$
(close to room temperature).

As reported in Fig. 5, we analyze the effects of the electron-vibration
coupling on the electronic response functions as a function of the level
energy $\epsilon$ at a fixed value of Hubbard interaction $U$ ($U=20$) in the
absence of coupling to phonon leads ($\gamma=0$) close to room temperature
($T=1.25$). For comparison, we report the transport properties relative to the
case when electron-electron and electron-vibration interactions are neglected
(indicated as Free in the figure).

The charge conductance $G$ is expected to be smaller than the free one due to
the effects of interactions. As shown in the upper left panel of Fig. 5, close
to room temperature, $G$ has peak values of the order of $10^{-1}$ $e^{2}/h$
($e^{2}/h$ is about $3.87\times 10^{-5}$ S). In particular, for
$\epsilon\simeq 20$, we have checked that $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 . As expected, the conductance as a function of the level
energy $\epsilon$ follows a behavior similar to the double-peak structure of
the spectral function as a function of the frequency. Therefore, $G$ has
maxima for $\epsilon\simeq 0=\mu$ and $\epsilon\simeq-U$, and a minimum at
$\epsilon\simeq-U/2$.

As shown in the upper right panel of Fig. 5, the Seebeck coefficient $S$ shows
large variations with changing $\epsilon$. Indeed, $S$ shows two maxima and
two minima whose magnitude is very large at room temperature being of the
order of $2$ $k_{B}/e$ ($k_{B}/e$ is about $86$ $\mu$eV/K). This complex
behavior is due to the role played by the strong electron correlations liu .
Actually, the structure close to $\epsilon=0$ (where $S$ vanishes) is nearly
translated by $-U$ (for $\epsilon\simeq-20$, $S$ goes again to zero).
Therefore, even at $\epsilon\simeq-U/2$, $S$ gets very small values.
Obviously, for large positive values of $\epsilon$, $S$ is small and negative
(n-type behavior). 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 .

As shown in the upper panels of Fig. 5, the most relevant effect of the
coupling $E_{P}$ on the conductance $G$ and the Seebeck coefficient $S$ is to
shift the curves and reduce the magnitude of the response function. The shift
of the conductance peaks and of the zeroes of the Seebeck coefficient is of
the order of $E_{P}$. At fixed level energy, unlike the conductance $G$, the
Seebeck coefficient is more sensitive to the changes of the coupling $E_{P}$.
For example, this occurs for energies close to the minima and the maxima. By
changing the values of $\epsilon$, there is an inversion in the behavior of
$S$ with increasing the electron-vibration coupling $E_{P}$.

As shown in the lower left panel of Fig.5, with varying the level energy
$\epsilon$, the electron thermal conductance $G_{K}^{el}$ shows the
characteristic double peak structure due to correlation effects liu . The
peaks values of $G_{K}^{el}$ are of the order of a few $0.01$ $k_{B}\Gamma$
($k_{B}\Gamma$ is about $4.198\times 10^{-10}$ W/K for $\hbar\Gamma\simeq 20$
meV). Therefore, the peak values are smaller than the thermal conductance
quantum $g_{0}(T)=\pi^{2}k_{B}^{2}T/(3h)$ at the room temperature
$T=1.25\hbar\Gamma\simeq 300$K ($g_{0}(T)\simeq 9.456\times
10^{-13}(W/K^{2})T)$ jezouin . We point out that electron-vibration
interactions affect the thermal conductance $G_{K}^{el}$ in a way completely
different from the charge conductance $G$ (compare left upper and left lower
panel of Fig. 5). Indeed, $G_{K}^{el}$ gets enhanced with increasing the
electron-oscillator coupling $E_{P}$. As discussed in the previous section,
within the adiabatic approach, the molecular effective level is renormalized
by the position variable $x$ which has a larger spreading upon increasing the
electron-vibration coupling.

We stress that the behavior of the electron thermal conductance $G_{K}^{el}$
shown in the lower left panel of Fig. 5 bears a strong resemblance with that
of the phonon thermal conductance $G_{K}^{ph}$ reported in Fig. 3. Both have a
double peak structure, and both are enhanced by the electron-vibration
coupling. Moreover, $G_{K}^{el}$ acquires values larger than those of
$G_{K}^{ph}$ in the energy region $-U\leq\epsilon\leq 0$. Obviously, the
values of these quantities are comparable for the chosen value of phonon
induced damping rate $\gamma=0.15\Gamma$. If one consider larger values of
$\gamma$ (for example $\gamma\simeq 0.4\Gamma$), then $G_{K}^{ph}$ would play
a major role in the total thermal conductance $G_{K}$. In any case, the values
of $G_{K}^{ph}$ and $G_{K}^{el}$ differ for $\epsilon\gg 0$ and
$\epsilon\ll-U$ since $G_{K}^{ph}$ acquires a finite asymptotic value
(obtained even in the absence of interactions on the molecule), while
$G_{K}^{el}$ goes rapidly to zero.

As shown in the lower right panel of Fig. 5, we analyze the behavior of the
electronic thermoelectric figure of merit $ZT^{el}$ neglecting the
contribution from $G_{K}^{ph}$. The quantity $ZT^{el}$ shows four peaks whose
values are larger than $1$, but smaller than the peak value around $3$
obtained in the absence of interactions. We stress that the peak values of
$ZT^{el}$ at room temperature are almost coincident with maxima and minima of
the Seebeck coefficient $S$. Actually, close to room temperature, the small
values of the conductance $G$ are fully compensated by the large values of the
Seebeck coefficient $S$. With increasing the electron-vibration coupling
$E_{P}$, the reduction of $G$ and $S$ combines with the enhancement of
$G_{K}^{el}$ leading to a sensible reduction of the figure of merit $ZT^{el}$.
Therefore, even if one neglects the role of phonon thermal conductance, the
effect of electron-electron and electron-vibration interactions is able to
induce a reduction of the figure of merit.

Figure 6: (Color online) Dimensionless thermoelectric 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$). In the plot, $U=20\hbar\Gamma$, $T=1.25\hbar\Gamma/k_{B}$,
$\gamma=0.15\Gamma$ and $\omega_{0}=0.25\Gamma$.

In Fig. 6, we focus on the total figure of merit $ZT$ as a function of the
level energy $\epsilon$ for different values of electron-vibration coupling
$E_{P}$ at $U=20\hbar\Gamma$ including the effects of the phonon thermal
conductance ($\gamma=0.15\Gamma$). From the comparison with the results
discussed in the previous paragraph, it emerges that the phonon thermal
conductance $G_{K}^{ph}$ induces an additional suppression of $ZT$. For the
realistic value of $E_{P}=0.5$ (intermediate coupling regime), the peak values
of $ZT$ are decreased by a factor of $2$ in comparison with $ZT^{el}$,
therefore the reduction of $ZT$ is not strong. Only for unrealistically large
electron-vibration couplings ($E_{P}$ larger than $1$), $ZT$ acquires peak
values less than unity. Summarizing, the cooperative effects of phonon leads,
electron-electron and electron-vibration interactions on the molecule are able
to weaken the thermoelectric performance of this kind of device. However,
within a realistic regime of parameters, the thermoelectric figure of merit
$ZT$ is still of the order of unity, making these devices a valid choice for
thermoelectric applications.

##  V 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 realistic
electron-electron and electron-vibration interactions. 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 generalized for including the large electron-electron
Coulomb repulsion. Parameters appropriate to $C_{60}$ molecules have been
considered. Within the intermediate electron-vibration coupling regime, the
phonon thermal conductance $G_{K}^{ph}$ is quite sensitive to the changes in
the occupation of electron level. Moreover, apart from an important asymptotic
value, $G_{K}^{ph}$ resembles the electron thermal conductance $G_{K}^{el}$.
With increasing the electron-vibration coupling, the phonon and the electron
thermal conductance get larger, while the charge conductance $G$ and the
thermopower $S$ get smaller. The figure of merit $ZT$ depends appreciably on
the behavior of $G_{K}^{ph}$ and intramolecular interactions. Indeed, for
realistic parameters of the model, $ZT$ can be substantially reduced, but its
peak values can be still of the order of unity indicating that the emerging
field of molecular thermoelectrics can be very interesting for applications.

The parameters of the junction are determined by the coupling between molecule
and metallic leads in the electronic and vibrational channels. For instance,
the strength of the intramolecular couplings depends on the choice of the
leads which screen the electron-electron and electron-vibration interactions.
In order to improve the thermoelectric efficiency, molecules and metallic
leads forming the junction have to ensure a weak phonon-center of mass
coupling (small $\gamma$) and a small strength of the electron-center of mass
interaction (small $E_{P}$). For realistic values of these couplings, the
values of the phonon thermal conductance $G_{K}$ are small compared to bulk
conductances. Therefore, the values of $ZT$ of the order of unity can be found
in molecular junctions since these systems provide a mechanism to keep the
phonon thermal conduction lower than that of bulks and other low-dimensional
structures. Finally, in this paper, we have shown that, for realistic values
of junction parameters, the phonon thermal conductance can be even smaller
than the electron counterpart in a large range of gate voltages.

The electron-vibration interaction of the Anderson-Holstein model analyzed in
this paper is related to the charge density injected by the external leads
onto the molecule. The renormalization of the lead-molecule hopping integral
induced by the center of mass movement could represent another possible source
of electron-vibration coupling koch  and it can be studied within the
adiabatic approach. However, we expect that the coupling through electron
level density plays a major role due to the large mass of the molecules
considered in this work. Finally, we stress that the approach proposed in this
paper can be generalized to the study of more realistic multi-level molecular
models and to cases where the number of atomic units within the molecule can
be varied.

###### Acknowledgements.

This work has been performed in the frame of the project GREEN (PON02 00029
2791179) granted to IMAST S.c.a.r.l. and funded by the MIUR (Ministero
dell’Istruzione, dell’Università e della Ricerca.).

##  Appendix A Comparison between different approaches within the Coulomb
blockade regime

In this Appendix, we compare the approach used in the main text for strong
Coulomb repulsion with that of Lacroix Lacroix  which retains additional self-
energy corrections upon the atomic limit Haug . We will consider the
electronic properties in the absence of electron-vibration coupling since we
are interested only on the effects induced by the electron-electron
interaction in equilibrium conditions at temperature $T=T_{\alpha}$ and
chemical potential $\mu=0=\mu_{\alpha}$. In this Appendix, we will use the
same units of the main text.

In contrast with the main text, in this Appendix, we will use a slightly
different kind of wide-band approximation for the electron leads. Actually, we
will consider an energy dependent tunneling rate $\Gamma_{0}(E)=\Gamma$, for
$-E_{c}\leq E\leq E_{C}$, and zero elsewhere, with $E_{C}$ cutoff energy much
larger than $U$. Therefore, the retarded self-energy of the electron level
$\Sigma_{0}(E)$ due to the effects of the electron leads is

$\Sigma_{0}(E)=\Lambda_{0}(E)-\frac{i}{2}\Gamma_{0}(E),$ (22)

where $\Lambda_{0}(E)$ is the real part of the retarded self-energy

$\Lambda_{0}(E)=\int_{-\infty}^{+\infty}\frac{dE^{\prime}}{2\pi}\frac{\Gamma_{0}(E^{\prime})}{E-E^{\prime}+\mu}=\frac{\Gamma}{2\pi}\ln{\left|\frac{E-E_{C}+\mu}{E+E_{C}+\mu}\right|}.$
(23)

In the limit where $E_{C}\rightarrow\infty$, one recovers the wide band
approximation used in the main text corresponding to a zero real part
$\Lambda_{0}(E)$.

We focus on the retarded Green function $G^{R}_{L}(\omega)$ relative to the
paramagnetic solution in order to calculate the spectral function
$A_{L}(\omega)=-2\Im{G^{R}_{L}(\omega)}$. The retarded Green function within
the Lacroix approximation for large $U$ Lacroix ; Haug  is

$\displaystyle G^{R}_{L}(\omega)$ $\displaystyle=$
$\displaystyle\frac{1-\rho}{\hbar\omega-\epsilon-\Sigma_{0}(\hbar\omega)-\Sigma_{h}(\hbar\omega)}+$
(24)
$\displaystyle\frac{\rho}{\hbar\omega-\epsilon-U-\Sigma_{0}(\hbar\omega)-\Sigma_{p}(\hbar\omega)},$

where $\rho$ is the level density per spin self-consistently calculated
through the following integral

$\rho=\int_{-\infty}^{+\infty}\frac{d(\hbar\omega)}{2\pi i}G^{<}_{L}(\omega),$
(25)

with the equilibrium lesser Green function $G^{<}_{L}(\omega)$

$G^{<}_{L}(\omega)=if(\hbar\omega)A_{L}(\omega),$ (26)

and $f(\hbar\omega)=1/(\exp{[\beta(\hbar\omega-\mu)]}+1)$ the free Fermi
distribution corresponding to the average chemical potential $\mu=0$. In Eq.
(24), the self-energy $\Sigma_{h}(\hbar\omega)$ is

$\Sigma_{h}(\hbar\omega)=-\frac{U\Sigma_{1}(\hbar\omega)}{\hbar\omega-\epsilon-U-\Sigma_{0}(\hbar\omega)-\Sigma_{3}(\hbar\omega)},$
(27)

while the self-energy $\Sigma_{p}(\hbar\omega)$ is

$\Sigma_{h}(\hbar\omega)=\frac{U\Sigma_{2}(\hbar\omega)}{\hbar\omega-\epsilon-\Sigma_{0}(\hbar\omega)-\Sigma_{3}(\hbar\omega)},$
(28)

where the self-energy $\Sigma_{i}(\hbar\omega)$, with $i=1,2,3$, is given by

$\displaystyle\Sigma_{i}(\hbar\omega)=\int_{-\infty}^{+\infty}\frac{dE}{2\pi}\Gamma_{i}(E)\times$
$\displaystyle\left[\frac{1}{\hbar\omega+E-\mu-2\epsilon-U+i\eta}+\frac{1}{\hbar\omega-E+\mu+i\eta}\right],$
(29)

with $\Gamma_{1}(E)=\Gamma_{0}(E)f(E)$, $\Gamma_{2}(E)=\Gamma_{0}(E)[1-f(E)]$,
$\Gamma_{3}(E)=\Gamma_{0}(E)$, and $\eta\rightarrow 0^{+}$. We notice that,
for large $U$, the weights of the poles of the Green function in Eq. (24) are
the same of the Green function examined in the main text. The Green function
within the Lacroix approach has the additional self-energy terms
$\Sigma_{i}(\hbar\omega)$, which take into account tunneling processes back
and forth to the leads.

Figure 7: (Color online) Spectral function $A$ (in units of $1/\hbar\Gamma$)
as a function of frequency $\omega$ (in units of $\Gamma$) for Hubbard
interaction $U=40\hbar\Gamma$ (Upper panel) and $U=16\hbar\Gamma$ (Lower
panel) at level energy $\epsilon=0$ and $T=1.25\hbar\Gamma/k_{B}$ (close to
room temperature) in the absence of electron-vibration coupling. Solid line:
first correction upon the atomic limit (used in the main text); dashed line:
additional correction upon the atomic limit (Lacroix approach).

As shown in Fig.7, we compare the spectral function obtained within the
approach used in the main text and $A_{L}$ within the Lacroix approximation
Lacroix  close to room temperature for two values of $U$ ($U=40$ upper panel,
$U=16$ lower panel). Both spectral functions exhibit a bimodal structure whose
peaks are separated by the energy $U$. The positions of the peaks within the
two approaches are very close, while the heights of the peaks are slightly
different. However, the ratio of the spectral weights of the two peaks does
not significantly depend on the approach. Obviously, the modification of the
isolated resonances is slightly more complicated within the Lacroix approach
than that due to the self-energy $\Sigma_{0}(\hbar\omega)$ alone. Actually,
the peaks within the Lacroix approach tend to be a little bit asymmetric.
Summarizing, the differences between the two approaches are minimal supporting
the use of the Green function method adopted in the present work. Finally, the
small differences between the two approaches are quantitatively similar with
decreasing $U$ from $40$ to $16$.

Figure 8: (Color online) Level density $N$ as a function of level energy
$\epsilon$ (in units of $\hbar\Gamma$) for different values of the Hubbard
interaction $U$ (in units of $\hbar\Gamma$) at $T=1.25\hbar\Gamma/k_{B}$
(close to room temperature) in the absence of electron-vibration coupling.

In this Appendix, we analyze also the total level occupation $N=2\rho$ (within
the paramagnetic solution). This quantity has been calculated by the two
approaches discussed in this Appendix finding minimal differences. In Fig.8,
we report the occupation determined by the approach used in the main text as a
function of level energy $\epsilon$ for different values of $U$. It shows the
typical profiles of the Coulomb blockade. Actually, for level energy
$\epsilon$ around $-U/2$, $N$ is $1$. The energy region with occupation close
to $1$ gets enhanced with increasing the value of $U$. Moreover, for
$\epsilon$ around $\mu=0$, $N$ goes from $1$ to $0$, while, for $\epsilon$
around $-U$, there is the transition from $N=2$ to $N=1$. These particular
values of $\epsilon$ are carefully analyzed in the main-text when the effects
of the electron-vibration coupling are included.

## References

  * (1) G. S. Nolas, J. Sharp and J. Goldsmid, Thermoelectrics: Basic Principles and New Materials Developments (Springer-Verlag, Berlin, Heidelberg, 2001).
  * (2) A. F. Ioffe, Semiconductor Thermoelements and Thermoelectric Cooling (Infosearch Limited, London, 1957).
  * (3) C. Kittel, Introduction to Solid State Physics (John Wiley & Sons, New York, 2004, 8th ed.).
  * (4) A. Shakouri, Annu. Rev. Mater. Res. 41, 399 (2011).
  * (5) K. Biswas, J. He, I. D. Blum, C.-I Wu, T. P. Hogan, D. N. Seidman, V. P. Dravid , and M. G. Kanatzidis, Nature 489, 414 (2012).
  * (6) K. Koumoto and T. Mori, Thermoelectric Nanomaterials - Materials Design and Applications (Springer-Verlag, Berlin, Heidelberg, 2013).
  * (7) B. Sothmann, R. Sanchez, and A. N. Jordan, arXiv:1406.5329.
  * (8) R. Venkatasubramanian, E. Siivola, T. Colpitts, and B. O’Quinn, Nature (London) 413, 597 (2001).
  * (9) T. C. Harman, P. J. Taylor, M. P. Walsh, and B. E. LaForge, Science 297, 2229 (2002).
  * (10) A. I. Hochbaum, R. Chen, R. D. Delgado, W. Liang, E. C. Garnett, M. Najarian, A. Majumdar, and P. Yang, Nature (London) 451, 163 (2008).
  * (11) L. D. Hicks and M. S. Dresselhaus, Phys. Rev. B 47, 16631 (1993).
  * (12) P. G. Murphy, and J. E. Moore, Phys. Rev. B 76, 155313 (2007).
  * (13) R. Venkatasubramanian, Recent Trends in Thermoelectric Materials Research III, Semiconductors and Semimetals Vol. 71 (Academic Press, New York, 2001), pp. 175-201; G. Chen, Recent Trends in Thermoelectric Materials Research III, Semiconductors and Semimetals Vol. 71 (Academic Press, New York, 2001), pp. 203-259.
  * (14) G. D. Mahan and J. Sofo, Proc. Natl. Acad. Sci. USA 93, 7436 (1996).
  * (15) S. V. Aradhya and L. Venkataraman, Nat. Nanotechnol. 8, 399 (2013).
  * (16) Y. Dubi and M. Di Ventra, Rev. Mod. Phys. 83, 131 (2011).
  * (17) P. Reddy, S.-Y. Jang, R. A. Segalman, and A. Majumdar, Science 315, 1568 (2007).
  * (18) S. K. Yee, J. A. Malen, A. Majumdar, and R. A. Segalman, Nano Lett. 11, 4089 (2011).
  * (19) C. M. Finch, V. M. Garcia-Suarez, and C. J. Lambert, Phys. Rev. B 79, 033405 (2009).
  * (20) P. Murphy, S. Mukerjee, and J. Moore, Phys. Rev. B 78, 161406(R) (2008).
  * (21) M. Galperin, A. Nitzan, and M. A. Ratner, Mol. Phys. 106, 397 (2008).
  * (22) J. Koch, F. von Oppen, Y. Oreg, and E. Sela, Phys. Rev. B 70, 195107 (2004).
  * (23) M. Leijnse, M. R. Wegewijs, and K. Flensberg, Phys. Rev. B 82, 045412 (2010).
  * (24) M. Paulsson and S. Datta, Phys. Rev. B 67, 241403(R) (2003).
  * (25) J. C. Cuevas and E. Scheer, Molecular Electronics: An Introduction to Theory and Experiment (World Scientific Publishing Company, Singapore, 1st edition, 2010).
  * (26) S. Datta, Lessons from Nanoelectronics: A New Perspective on Transport (World Scientific Publishing Company, Singapore, 2012).
  * (27) R. Y. Wang, R. A. Segalman, and A. Majumdar, Appl. Phys. Lett. 89, 173113 (2006).
  * (28) Z. Wang, J. A. Carter, A. Lagutchev, Y. K. Koh, N.-H. Seong, D. G. Cahill, and D. D. Dlott, Science 317, 787 (2007).
  * (29) T. Meier, F. Menges, P. Nirmalraj, H. H lscher, H. Riel, and B. Gotsmann, Phys. Rev. Lett. 113, 060801.
  * (30) W. Lee, K. Kim, W. Jeong, L. A. Zotti, F. Pauly, J. C. Cuevas, and P. Reddy, Nature 498, 209 (2013).
  * (31) V. Schoeps, V. Zlatic, and T. A. Costi, J. Phys.: Conf. Ser. 273, 012155 (2011).
  * (32) M. Galperin, A. Nitzan, and M. A. Ratner, J. Phys.: Condens. Matter 19, 103201 (2007).
  * (33) H. Park, J. Park, A. K. L. Lim, E. H. Anderson, A. P. Alivisatos, and P. L. McEuen, Nature 407, 57 (2000).
  * (34) H. Qin, A. W. Holleitner, K. Eberl, and R. H. Blick, Phys. Rev. B 64, 241302(R) (2001).
  * (35) K.-H. Yang, Y.-L. Zhao, Y.-J. Wu, Y.-P. Wu, Phys. Lett. A 374, 2874 (2010).
  * (36) X. Zianni, Phys. Rev. B 82, 165302 (2010).
  * (37) O. Entin-Wohlman, Y. Imry, and A. Aharony, Phys. Rev. B 82, 115314 (2010).
  * (38) J. Liu, Q.-F. Sun, and X. C. Xie, Phys. Rev. B 81, 245323 (2010).
  * (39) J. Ren, J.-X. Zhu, J. E. Gubernatis, C. Wang, and B. Li, Phys. Rev. B 85, 155443 (2012).
  * (40) M. Bagheri Tagani and H. Rahimpour Soleimani, Physica B 413, 86 (2013).
  * (41) L. Arrachea, N. Bode, F. von Oppen, arXiv:1407.3127.
  * (42) B. C. Hsu, C.-W. Chiang, and Y.-C. Chen, Nanotechnology 23, 275401 (2012).
  * (43) F. Remaggi, N. Traverso Ziani, G. Dolcetto, F. Cavaliere, M. Sassetti, New J. Phys. 15, 083016 (2013); N. Traverso Ziani, G. Piovano, F. Cavaliere, M. Sassetti, Rev. B 83, 245311 (2011).
  * (44) D. Mozyrsky, M. B. Hastings, and I. Martin, Phys. Rev. B 73, 035104 (2006).
  * (45) F. Pistolesi, Ya. M. Blanter, and I. Martin, Phys. Rev. B 78, 085127 (2008).
  * (46) R. Hussein, A. Metelmann, P. Zedler, and T. Brandes, Phys. Rev. B 82, 165406 (2010).
  * (47) A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, and V. Cataudella, Phys. Rev. B 83, 115420 (2011).
  * (48) A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, and V. Cataudella, Phys. Rev. B 86, 035420 (2012).
  * (49) C. A. Perroni, D. Ninno, and V. Cataudella, Phys. Rev. B 90, 125421 (2014).
  * (50) H. Haug and A.-P. Jauho, Quantum Kinetics in Transport and Optics of Semiconductors, (Springer, Berlin, 2008).
  * (51) X. Lu, M. Grobis, K. H. Khoo, S. G. Louie, and M. F. Crommie, Phys. Rev. Lett. 90, 096802 (2003).
  * (52) L. H. Yu and D. Natelson, Nano Lett. 4, 79 (2005).
  * (53) J. Mravlje, A. Ramsak, and T. Rejec, Phys. Rev. B 74, 205320 (2006).
  * (54) U. Weiss, Quantum Dissipative Systems, (World Scientific Publishing Company, Singapore, 2008, 3rd edition).
  * (55) M. Schiró and K. Le Hur, Phys. Rev. B 89, 195127 (2014).
  * (56) C. A. Perroni, V. Cataudella, G. De Filippis, and V. Marigliano Ramaglia, Phys. Rev. B 71, 113107 (2005).
  * (57) J. Splettstoesser, M. Governale, J. König, and R. Fazio, Phys. Rev. Lett. 95, 246803 (2005).
  * (58) A. R. Hernández, F. A. Pinheiro, C. H. Lewenkopf, and E. R. Mucciolo, Phys. Rev. B 80, 115311 (2009).
  * (59) A. Nocera, C. A. Perroni, V. Marigliano Ramaglia, G. Cantele, and V. Cataudella, Phys. Rev. B 87, 155435 (2013).
  * (60) C. A. Perroni, A. Nocera, and V. Cataudella, Europhys. Lett. 103, 58001 (2013).
  * (61) C. A. Perroni, F. Romeo, A. Nocera, V. Marigliano Ramaglia, R. Citro, and V. Cataudella, J. Phys.: Condens. Matter 26, 365301 (2014).
  * (62) R. L. Honeycutt, Phys. Rev. A 45, 600 (1992).
  * (63) R. L. Honeycutt, Phys. Rev. A 45, 604 (1992).
  * (64) J.-S. Wang, Phys. Rev. Lett. 99, 160601 (2007).
  * (65) S. Jezouin, F. D. Parmentier, A. Anthore, U. Gennser, A. Cavanna, Y. Jin, and F. Pierre, Science 342, 601 (2013).
  * (66) C. Lacroix, J. Phys. F 11, 2389 (1981).