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

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.

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 1,2 . 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 el K + G ph K is the total thermal conductance, with G el K and G ph K 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 3 . 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 1,4,5 .
Recently, the possibility of controlling materials at the nanoscale has been exploited to optimize the thermoelectric efficiency. 4,6,7 For example, a maximum ZT ≃ 2.4 has been observed at room temperature in a thin-film thermoelectric device 8 . High values of ZT have been reported in quantum dot superlattices 9 and in semiconductor nanowires 10 , where phonon confinement can lead to a lower phonon thermal conductance 11,12 . Actually, a significant reduction in lattice thermal conductivity is considered as the main route for having high ZT in lowdimensional materials 13 . 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 14 . Finally, in nanoscopic Coulomb-coupled systems, the thermoelectric properties can be optimized by exploiting the Coulomb blockade regime and changing the gate voltage 7 .
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 15,16 . Indeed, the emerging field of molecular thermoelectrics has attracted a lot of attention in recent years [17][18][19][20][21][22][23] . 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 17,18,[24][25][26] . Measurements of thermoelectric properties have been performed in junctions with fullerene (C 60 ) 18 finding a high value of the molecular thermopower (S of the order of −30 µ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 16,[27][28][29] or limited to a range where transport is elastic 30 .
In molecular junctions, intramolecular electronelectron and electron-vibration interactions typically constitute the largest energy scales affecting the thermoelectric properties. 25,31,32 Moreover, the center of mass oscillation of the molecule 33 , or thermally induced acoustic phonons 34 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 25,32 . The thermopower S and the thermoelectric figure of merit ZT have been found to be sensitive to the strength of intramolecular interactions [21][22][23][35][36][37][38][39][40][41] . However, the phonon thermal contribution G ph K to the figure of merit ZT has been calculated only at a perturbative level of the electron-vibration coupling 42 .
In devices with large molecules or carbon nanotube quantum dots 43 , a nonequilibrium adiabatic approach has been introduced for spinless electrons exploiting the low energy of the relevant vibrational degrees of freedom [44][45][46][47][48] . 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 49 .
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 ph K 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 ph K , 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. 25,50 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. 49 this paper, we will consider model parameters appropriate to C 60 molecular junctions. The HamiltonianĤ is given bŷ where the HamiltonianĤ el takes into account the electronic degrees of freedom of the leads and the molecule, H ph the vibrational degrees of freedom of the leads and the molecule, andĤ int the coupling between electronic and vibrational degrees of freedom (see Fig. 1 for a sketch of the molecular junction model).
The electronic HamiltonianĤ el of Eq. (1) iŝ where the molecular electronic level has energy ǫ, the σ spin electron density operator isn σ =d † σd σ , withd † σ (d σ ) creation (annihilation) σ 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 ǫ. 25 The Coulomb repulsion on the molecule is simulated with a Hubbard term U , which gives an energy penalty for electron occupations with spin ↑ and ↓. 25 The lead density operator isn q,α,σ =ĉ † q,α,σĉq,α,σ , where the operatorsĉ † q,α,σ (ĉ q,α,σ ) create (annihilate) electrons with momentum q, spin σ, and energy ε q,α = ξ q,α − µ α in the left (α = L) or right (α = R) free metallic leads, with µ α chemical potential of the lead α in equilibrium at the temperature T α . We consider the temperatures T L = T + ∆T /2 and T R = T − ∆T /2, with T average temperature. Moreover, we fix the chemical potentials µ L = eV /2 and µ R = −eV /2, with e modulus of the electron charge, V bias potential, and average chemical potential µ = 0. The electronic tunneling between the molecular dot and a state q in the lead α has the amplitude V q,α . As usual for metallic leads, the density of states ρ q,α is assumed flat around the small energy range relevant for the molecular orbital, making valid the wide-band limit: ρ q,α → ρ α , V q,α → V α . Therefore, the full hybridization width of the molecular orbital is Γ = α Γ α , with Planck constant and the tunneling rate Γ α = 2πρ α |V α | 2 / . In the following, we consider the symmetric configuration: Γ L = Γ R = Γ/2. In junctions with C 60 molecules, Γ has been estimated to be of the order of 20 meV. 52,53 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 Γ. 52,53 The center of mass mode can be considered as the relevant vibrational mode of the molecule. 49 Indeed, experiments have evidenced a coupling between the center of mass mode and the electron dynamics in junctions with C 60 molecules. 33 In Eq.(1), the HamiltonianĤ ph describes the vibrations of the slow center of mass mode, the free phonon modes of the leads, and the coupling between them: (3) The center of mass hamiltonianĤ cm iŝ wherep andx are the center of mass momentum and position operators, respectively, M is the total large mass, k is the effective spring constant, with frequency ω 0 = k/M . In Eq.(3), the operatorsâ † q,α (â q,α ) create (annihilate) phonons with momentum q and frequency ω q,α in the lead α. 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 α is given by the elastic constant C q,α . For large molecules, the center of mass mode has a low frequency ω 0 which is typically smaller than the Debye frequency ω D of the metallic leads ( ω D ≃ 15 − 20 meV for metals like silver, gold, and platinum 3 ). For example, ω 0 ≃ 5 meV in C 60 junctions 33 , hence ω 0 ≃ 0.25Γ. Therefore, for large molecules, the adiabatic regime is valid for the center of mass oscillator: ω 0 << Γ and ω 0 << ω D . Within this regime, the effect of the α phonon lead on the center of mass mode provides a constant damping rate γ α . 54 In analogy with the electronic model, we consider the symmetric configuration: γ L = γ R = γ/2. For junctions with C 60 molecules and leads of Ag, Au, and Pt, γ ≃ 3 − 8 meV, therefore γ is of the same order of ω 0 (γ ≃ 0.15 − 0.40Γ) 49 .
Finally, the interaction termĤ int in the Anderson-Holstein model of Eq.(1) is provided by a linear coupling between the total electron density on the molecule,n = σn σ , and thex operator of the center of mass: where λ is the electron-vibration coupling constant. In the following, the electron-vibration interaction will be described in terms of the coupling energy E P = λ 2 /(2k).
In this paper, Γ ≃ 20 meV will be the energy unit (Γ the frequency unit, 1/Γ the time unit). We will measure lengths in units of 2λ/k, and temperatures in units of Γ/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: Besides, the electron-vibration coupling is not weak, but it is estimated to be in the intermediate regime: Since ω 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 55 which is valid only in the limit of infinite local Coulomb repulsion for energies close to the chemical potential and low temperatures 56 .

A. 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) = ǫ + λx(t) within the Coulomb blockade regime. 57,58 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. 50 Therefore, for the paramagnetic solution, the level spectral function A 0 (ω, x) at zero order of the adiabatic expansion becomes where ρ(x) is the level density per spin self-consistently calculated at fixed position x through the following integral with the lesser Green function G < 0 (ω, x) and . Actually, the spectral function is characterized by a double peak structure that, for large U , is robust against the effects of electronvibration 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. 50 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 ǫ with respect to the leads chemical potential (average chemical potential µ = 0 in this work). In Appendix A, we report the molecular electron occupation N as a function of level energy ǫ showing the typical profiles of the Coulomb blockade. In particular, the following energies are relevant: ǫ = −U/2 (close to half-filling N = 1), ǫ = −U (transition from level occupation N = 1 to N = 2), ǫ = 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. 47,48,[59][60][61] 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 ≪ 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.

B. 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 fol-lowing generalized Langevin equation for the slow center of mass which has the deterministic force F det (x, v) and the position dependent fluctuating force ξ(x, t). The deterministic force can be decomposed into a generalized force F gen (x) with F λ (x) = −2λρ(x) induced by the electron-vibration coupling, and, as a result of the adiabatic expansion, a dissipative force with an effective position dependent pos- due to the electron-vibration interaction. The fluctuating force ξ(x, t) in Eq. (9) is such that where the effective position dependent noise term with D λ (x) determined by the electron-vibration coupling and the greater Green function G > 0 (ω, x) It is worthwhile pointing out that, in equilibrium conditions at temperature T = T α and chemical potential µ = µ α = 0, the adiabatic procedure gives rise to a generalized fluctuation-dissipation relation 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. 47,62,63 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. Before going to the section about results, we discuss the features of the electron-vibration induced damping rate γ λ (x) = A λ (x)/m, with A λ (x) given in Eq. (13). The magnitude of γ λ (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 γ λ (x) are always smaller than realistic values of the lead induced damping rate γ (γ = 0.15 will be considered in this paper). This implies that the effects due to the electronvibration 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 γ λ (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, 49 γ λ (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 γ λ (x) largely shift passing from the quasi half-filling case (close to ǫ = −10 = −U/2, state with flat occupation) to conditions out of half-filling (close to ǫ = −20 = −U and ǫ = 0 = µ, state with strong density fluctuations). The self-consistent calculation of γ λ (x) provides a direct signature of the strong local interaction since it is deter-mined by the adiabatic expansion of the electron occupation.
A comparison of the x dependence between γ λ (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 ǫ. 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 ǫ. In the quasi half-filled case (ǫ = −10), the peak positions of γ λ (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 ǫ = −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 ǫ = −20, the peak position of P (x) shifts towards lower values close to −0.75 (N ≃ 1.5), and, for ǫ = 0, to 0.25 (N ≃ 0.5). We point out that, for ǫ = −20, the first peak of γ λ (x) is close to x = 0, while, for ǫ = 0, the second peak of γ λ (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, 49 these effects are present not only close to ǫ = µ = 0, but also to ǫ = −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 electronelectron 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 ω 0 = 0.25Γ, and γ = 0.15Γ (larger values of γ have been considered in a recent paper 49 ).

A. Phonon heat transport
In this subsection, we will focus on the phonon thermal conductance G ph K calculated within the linear response regime around temperature T as with J ph α current from the α phonon lead. 49,64 The conductance G ph K 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 γ (γ = 0.15Γ in this work) which is typically larger than the peak values of electron-vibration induced damping rate γ λ (x). As shown in Fig.3, in the regime of weak electronvibration coupling E P , low level occupation (ǫ ≫ 0), and double level occupation (ǫ ≪ −U ), G ph K is close to 0.04 k B Γ (k B Γ is about 419.8 pW/K for Γ ≃ 20 meV), a numerical value coincident with an analytical estimate of G ph K given in a recent paper 49 . 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 ph K always gets larger with increasing the electron-vibration coupling E P . Moreover, this increase of G ph K strongly depends on the value of level energy ǫ. In contrast with the spinless case (reported for comparison in Fig. 3 at E P = 1), we stress that the enhancement of G ph K takes place not only close to ǫ ≃ 0, but also to ǫ ≃ −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 Γ ≃ 20 pW/K. Therefore, the calculated G ph K is in very good agreement with the thermal conductance of the order of a few 10 pW/K measured for molecules anchored to gold 28,29 . In any case, due to the strong electronelectron interactions, G ph K can be enhanced in a new large energy region. On the other hand, for ǫ ≃ −U/2, G ph K 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 ph K 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 γ λ (x) discussed in the previous section.

B. 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(ω) is evaluated making the average of the function A 0 (ω, x) in Eq.(6) over P (x): In this section, the spectral function will be discussed in equilibrium conditions at temperature T (V = 0 and ∆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 ǫ.
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 ǫ = −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 (ǫ = 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.

C. 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 µ = 0 and temperature T (∆T → 0 + , V → 0 + ). We will evaluate the electronic conductance G where f ( ω) = 1/(exp [β( ω − µ)] + 1) is the free Fermi distribution corresponding to the average chemical potential µ = 0. Then, we will calculate the Seebeck coef- (20) Finally, we will determine the electron thermal conduc- The total thermal conductance G K = G el K + G ph K 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 (γ = 0), G K = G el K , so that ZT = ZT el , which can be used to characterize the electronic thermoelectric efficiency.
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 ǫ at a fixed value of Hubbard interaction U (U = 20) in the absence of coupling to phonon leads (γ = 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 × 10 −5 S). In particular, for ǫ ≃ 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 18 . As expected, the conductance as a function of the level energy ǫ follows a behavior similar to the double-peak structure of the spectral function as a function of the frequency. Therefore, G has maxima for ǫ ≃ 0 = µ and ǫ ≃ −U , and a minimum at ǫ ≃ −U/2.
As shown in the upper right panel of Fig. 5, the Seebeck coefficient S shows large variations with changing ǫ. 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 µeV/K). This complex behavior is due to the role played by the strong electron correlations 38 . Actually, the structure close to ǫ = 0 (where S vanishes) is nearly translated by −U (for ǫ ≃ −20, S goes again to zero). Therefore, even at ǫ ≃ −U/2, S gets very small values. Obviously, for large positive values of ǫ, S is small and negative (ntype behavior). In particular, for ǫ = 20, S is about −0.45k B /e ≃ −38.5µV /K in agreement with the magnitude of experimental data in C 60 18 . 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 ǫ, 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 ǫ, the electron thermal conductance G el K shows the characteristic double peak structure due to correlation effects 38 . The peaks values of G el K are of the order of a few 0.01 k B Γ (k B Γ is about 4.198 × 10 −10 W/K for Γ ≃ 20 meV). Therefore, the peak values are smaller than the thermal conductance quantum g 0 (T ) = π 2 k 2 B T /(3h) at the room temperature T = 1.25 Γ ≃ 300K (g 0 (T ) ≃ 9.456 × 10 −13 (W/K 2 )T ) 65 . We point out that electron-vibration interactions affect the thermal conductance G el K in a way completely different from the charge conductance G (compare left upper and left lower panel of Fig. 5). Indeed, G el K 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 el K shown in the lower left panel of Fig.  5 bears a strong resemblance with that of the phonon thermal conductance G ph K reported in Fig. 3. Both have a double peak structure, and both are enhanced by the electron-vibration coupling. Moreover, G el K acquires values larger than those of G ph K in the energy region −U ≤ ǫ ≤ 0. Obviously, the values of these quantities are comparable for the chosen value of phonon induced damping rate γ = 0.15Γ. If one consider larger values of γ (for example γ ≃ 0.4Γ), then G ph K would play a major role in the total thermal conductance G K . In any case, the values of G ph K and G el K differ for ǫ ≫ 0 and ǫ ≪ −U since G ph K acquires a finite asymptotic value (obtained even in the absence of interactions on the molecule), while G el K 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 ph K . 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 el K 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 electronvibration interactions is able to induce a reduction of the figure of merit.
In Fig. 6, we focus on the total figure of merit ZT as a function of the level energy ǫ for different values of electron-vibration coupling E P at U = 20 Γ including the effects of the phonon thermal conductance (γ = 0.15Γ). From the comparison with the results discussed in the previous paragraph, it emerges that the phonon thermal conductance G ph K 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 ph K 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 ph K is quite sensitive to the changes in the occupation of electron level. Moreover, apart from an important asymptotic value, G ph K resembles the electron thermal conductance G el K . 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 ph K 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 phononcenter of mass coupling (small γ) 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 electronvibration coupling 22 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.

ACKNOWLEDGMENTS
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 66 which retains additional self-energy corrections upon the atomic limit 50 . 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 α and chemical potential µ = 0 = µ α . 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 Γ 0 (E) = Γ, for −E c ≤ E ≤ E C , and zero elsewhere, with E C cutoff energy much larger than U . Therefore, the retarded selfenergy of the electron level Σ 0 (E) due to the effects of the electron leads is where Λ 0 (E) is the real part of the retarded self-energy (A2) In the limit where E C → ∞, one recovers the wide band approximation used in the main text corresponding to a zero real part Λ 0 (E).
We focus on the retarded Green function G R L (ω) relative to the paramagnetic solution in order to calculate the spectral function A L (ω) = −2ℑG R L (ω). The retarded Green function within the Lacroix approximation for large U 50,66 is where ρ is the level density per spin self-consistently calculated through the following integral with the equilibrium lesser Green function G < L (ω) and f ( ω) = 1/(exp [β( ω − µ)] + 1) the free Fermi distribution corresponding to the average chemical potential µ = 0. In Eq. (A3), the self-energy Σ h ( ω) is while the self-energy Σ p ( ω) is where the self-energy Σ i ( ω), with i = 1, 2, 3, is given by has the additional self-energy terms Σ i ( ω), which take into account tunneling processes back and forth to the leads. 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 66 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 Σ 0 ( ω) 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. In this Appendix, we analyze also the total level occupation N = 2ρ (within the paramagnetic solution). This quantity has been calculated by the two approaches dis- cussed 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 ǫ for different values of U . It shows the typical profiles of the Coulomb blockade. Actually, for level energy ǫ around −U/2, N is 1. The energy region with occupation close to 1 gets enhanced with increasing the value of U . Moreover, for ǫ around µ = 0, N goes from 1 to 0, while, for ǫ around −U , there is the transition from N = 2 to N = 1. These particular values of ǫ are carefully analyzed in the main-text when the effects of the electron-vibration coupling are included.