On the Role of Local Many-Body Interactions on the Thermoelectric Properties of Fullerene Junctions

The role of local electron–vibration and electron–electron interactions on the thermoelectric properties of molecular junctions is theoretically analyzed focusing on devices based on fullerene molecules. A self-consistent adiabatic approach is used in order to obtain a non-perturbative treatment of the electron coupling to low frequency vibrational modes, such as those of the molecule center of mass between metallic leads. The approach also incorporates the effects of strong electron–electron interactions between molecular degrees of freedom within the Coulomb blockade regime. The analysis is based on a one-level model which takes into account the relevant transport level of fullerene and its alignment to the chemical potential of the leads. We demonstrate that only the combined effect of local electron–vibration and electron–electron interactions is able to predict the correct behavior of both the charge conductance and the Seebeck coefficient in very good agreement with available experimental data.


Introduction
In recent years, the field of molecular thermoelectrics has attracted a lot of attention [1][2][3][4][5][6][7][8][9][10][11][12]. One of the aims is to improve the thermoelectric efficiency of nanoscale devices by controlling the electronic and vibrational degrees of freedom of the molecules. Moreover, useful information on charge and energy transport mechanisms can be extracted by studying the thermoelectric properties of molecular junctions [1,3,4,13,14]. In addition to the charge conductance G, the Seebeck coefficient S is typically measured in these devices. Measurements in junctions with fullerene (C 60 ) have found a high value of thermopower (of the order or even smaller than −30 µV/K) [4]. Understanding the thermopower is also important for helping advances in thermoelectric performance of large-area molecular junctions [15,16]. Moreover, recently, the application of an Al gate voltage at Au-C 60 -Au junction has allowed to achieve the electrostatic control of charge conductance and thermopower with unprecedented control [17]. However, the precise transport mechanisms affecting both G and S remain elusive in these kinds of measurements. Finally, due to experimental challenges [2,[18][19][20], only recently the thermal conductance of single-molecule junctions has been fully characterized [21].
In molecular junctions, relevant contributions to the thermoelectric properties typically result from intramolecular electron-electron and electron-vibration interactions [1,22]. An additional source of coupling between electronic and vibrational degrees of freedom is also provided by the center of mass oscillation of the molecule between the metallic leads [23]. Different theoretical techniques [1,22] have been used to study the effects of local many-body interactions which affect the thermoelectric transport properties [7][8][9][24][25][26][27][28] in a significant way.
In devices with large molecules such as fullerenes or carbon nanotube quantum dots, a non-perturbative treatment of electron-vibration coupling can be obtained within an adiabatic approach which is based on the slowness of the relevant vibrational modes in comparison with the fast electron dynamics [29][30][31][32][33][34][35][36][37][38][39]. The adiabatic approach can also include a strong Coulomb repulsion allowing the self-consistent calculation of thermoelectric properties of massive molecules, such as fullerenes, within the Coulomb blockade regime [40].
In this paper, the thermoelectric properties of a molecular junction are analyzed focusing on the role of electron-electron and electron-vibration interactions. An adiabatic approach developed in the literature 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 [40]. Parameters appropriate for junctions with C 60 molecules are considered in this paper. In particular, a one-level model is taken into account since it describes the relevant transport level of fullerene and its alignment to the chemical potential of the metallic leads.
The aim of this paper is to thoroughly investigate both the charge conductance and the Seebeck coefficient since accurate experimental data are available for Au-C 60 -Au junction in [17] as a function of the voltage gate. We show that an accurate description of the transport properties is obtained in the intermediate regime for the electron-vibration coupling and in the strong coupling regime for the electron-electron interaction. Moreover, we point out that only the combined effect of electron-vibration and electron-electron interactions is able to predict the correct behavior of both the charge conductance and the Seebeck coefficient finding a very good agreement with available experimental data.
The paper is organized as follows. In Section 2, a very general model for many electronic levels and multiple vibrational degrees is considered and the adiabatic approach is exposed. In Section 3, the one-level model is presented. In Section 4, the theoretical results are presented together with the precise comparison with experimental data. Finally, in Section 5, conclusions and final discussions are given.

Model and Method
In this section, we introduce a general Hamiltonian for a multilevel molecule including many-body interactions between molecular degrees of freedom: the local electron-electron interaction and the local electron coupling to molecular vibrational modes. The model simulates also the coupling of the molecule to two leads in the presence of a finite bias voltage and temperature gradient. The total Hamiltonian of the system isĤ =Ĥ mol +Ĥ leads +Ĥ leads−mol , whereĤ mol is the Hamiltonian describing the molecular degrees of freedom,Ĥ leads the leads' degrees of freedom andĤ leads−mol the coupling between molecule and leads. In this paper, we assume, as usual in the field of molecular junctions, that the electronic and vibrational degrees of freedom in the metallic leads are not interacting [1,41]; therefore, the electron-electron and electron-vibration interactions are effective only on the molecule. In Equation (1), the molecule HamiltonianĤ mol iŝ where c m,σ (c † m,σ ) is the standard electron annihilation (creation) operator for electrons on the molecule levels with spin σ = ↑, ↓, where indices m, l can assume positive integer values with a maximum M indicating the total number of electronic levels in the molecule. The matrix ε m,l σ is assumed diagonal in spin space,n l,σ = c † l,σ c l,σ is the electronic occupation operator relative to level l and spin σ, and U represents the Coulomb-Hubbard repulsion between electrons. We assume that only the diagonal part of the matrix ε m,l σ is nonzero and independent of the spin: ε m,m σ = ε m , where ε m are the energies of the molecule levels.
In Equation (2), the molecular vibrational degrees of freedom are described by the Hamiltonian where s = (1, . . . , N), with N being the total number of vibrational modes; m s is the effective mass associated with the sth vibrational mode; andp s is its momentum operator. Moreover, s is the harmonic potential (with k s the spring constants, and the oscillator frequencies ω s 0 = √ k s /m s ),x s is the displacement operator of the vibrational mode s, and X = (x 1 , . . . ,x N ) indicates all the displacement operators.
In Equation (2), the electron-vibration couplingĤ int is assumed linear in the vibrational displacements and proportional to the electron level occupationŝ where s = (1, . . . , N) indicates the vibrational modes of the molecule, l = (1, . . . , M) denotes its electronic levels,n l = ∑ σ n l,σ is the electronic occupation operator of the level l, and λ s,l is a matrix representing the electron-vibrational coupling.
In Equation (1), the Hamiltonian of the electron leads is given bŷ where the operatorsĉ † k,α,σ (ĉ k,α,σ ) create (annihilate) electrons with momentum k, spin σ, and energy ε k,α = E k,α − µ α in the left (α = L) or right (α = R) leads. The left and right electron leads are considered as thermostats in equilibrium at the temperatures T L = T + ∆T/2 and T R = T − ∆T/2, respectively, with T the average temperature and ∆T temperature difference. Therefore, the left and right electron leads are characterized by the free Fermi distribution functions f L (E) and f R (E), respectively, with E the energy. The difference of the electronic chemical potentials in the leads provides the bias voltage V bias applied to the junction: µ L = µ + eV bias /2, µ R = µ − eV bias /2, with µ the average chemical potential and e the electron charge. In this paper, we focus on the regime of linear response that involves very small values of bias voltage V bias and temperature ∆T.
Finally, in Equation (1), the coupling between the molecule and the leads is described bŷ where the tunneling amplitude between the molecule and a state k in the lead α has the amplitude V m k,α . For the sake of simplicity, we suppose that the density of states ρ k,α for the leads is flat within the wide-band approximation: ρ k,α → ρ α , V m k,α → V m α . Therefore, the full hybridization width matrix of the molecular orbitals is Γ m,n = ∑ α Γ m,n α = ∑ α Γ m,n α , with the tunneling rate Γ m,n α = 2πρ α V m * α V n α . In this paper, we consider the symmetric configuration Γ L = Γ R = Γ/2, where, in the following, bold letters indicate matrices.
In this paper, we consider the electronic system coupled to slow vibrational modes: ω s 0 Γ m,n , for each s and all pairs of (m, n). In this limit, we can treat the mechanical degrees of freedom as classical, acting as slow classical fields on the fast electronic dynamics. Therefore, the electronic dynamics is equivalent to a multi-level problem with energy matrix ε m → ε m + λ m x m , where x m are now classical displacements [32,39]. This is called in the literature adiabatic approximation for vibrational degrees of freedom.
Within the adiabatic approximation, one gets Langevin self-consistent equations for the vibrational modes of the molecule [33,39] m sẍs + k s x s = F s (t)+ξ s (t), where the generalized force F s is due to the effect of all electronic degrees of freedom through the electron-vibration coupling [32,39]: with the trace "Tr", taken over the molecule levels, defined in terms of the lesser molecular matrix Green's function G < (t, t ) with matrix elements G < m,l (t, t ) = i c † m,σ (t)c l,σ (t ) . Quantum electronic density fluctuations on the oscillator motion are responsible for the fluctuating force ξ s (t) in Equation (7), which is derived below together with generalized force.
In deriving equations within the adiabatic approximation [39], next, for the sake of simplicity, we do not include explicitly the effect of the Coulomb repulsion on the molecule Hamiltonian. In the next section, we show that, in the case of a single level molecule with large repulsion U, the adiabatic approach works exactly as in the non-interacting case provided that each Green's function pole is treated as a non interacting level [40].
In our notation, G denotes full Green's functions, while G denotes the strictly adiabatic (or frozen) Green's functions, which are calculated at a fixed value of X. Starting from the Dyson equation [32,39,41], the adiabatic expansion for the retarded Green's function G R is given by where G R (E, X) is the strictly adiabatic (frozen) retarded Green's function including the coupling with the leads ε(X) represents the matrix ε m,l σ + ∑ s λ s x s δ l,m and Σ R,leads = ∑ α Σ R,leads α is the total self-energy due to the coupling between the molecule and the leads. For the lesser Green's function G < , the adiabatic approximation involves The electron-vibration induced forces at the zero order of the adiabatic limit (G < G < ) are given by The leading order correction to the lesser Green's function G < provides a term proportional to the vibrational velocity where the tensor θ can be split into symmetric and anti-symmetric contributions [32]: θ = θ sym + θ a , where we have introduced the notation {C s,s } sym,a = 1 2 {C s,s ± C s ,s } sym,a for symmetric and anti-symmetric parts of an arbitrary matrix C. Indeed, there is a dissipative term θ sym and an orbital, effective magnetic field θ a in the space of the vibrational modes.
We can now discuss the stochastic forces ξ s (t) in Equation (7) within the adiabatic approximation. In the absence of electron-electron interactions, the Wick theorem allows writing the noise correlator as where G > (t, t ) is the greater Green's function with matrix elements G > m,l (t, t ) = −i c m,σ (t)c † l,σ (t ) . In the adiabatic approximation, one first substitutes the full Green's function G by the adiabatic zero-order Green's function G and then observes that the electronic fluctuations act on short time scales only. Therefore, the total forces ξ s (t) are locally correlated in time: where Once the forces and the noise terms are calculated, Equation (7) represents a set of nonlinear Langevin equations in the unknown x s . Even for the simple case where only one vibrational degree of freedom is present, the stochastic differential equation should be solved numerically in the general non-equilibrium case [33,37,38]. Actually, one can calculate the oscillator distribution functions , and, therefore, all the properties of the vibrational modes. Using this function, one can determine the average O of an electronic or vibrational observable O(X, V ): The electronic observables, such as charge and heat currents, can be evaluated exploiting the slowness of the vibrational degrees of freedom. In a previous paper [39], we discussed the validity of the adiabatic approximation, stressing that it is based on the separation between the slow vibrational and fast electronic timescales. Actually, physical quantities calculated within the adiabatic approach are very reliable in a large regime of electronic parameters since this self-consistent approach is not perturbative in the electron-vibration coupling. Therefore, the approach is able to overcome the limitations of the perturbative theory typically used in the literature [42,43].

One-Level Model
In the remaining part of the paper, we consider the simple case where the molecule is modeled as a single electronic level (M = 1 in the previous section) locally interacting with a single vibrational mode (N = 1 in the previous section). Therefore, the focus is on a molecular level which is sufficiently separated in energy from other orbitals. In particular, we analyze the C 60 molecule where the lowest unoccupied molecular orbital (LUMO) energy differs from the highest occupied molecular orbital (HOMO) energy for energies of the order of 1 eV [23,44]. Even when the degeneracy of the LUMO is removed by the contact with metal leads, the splitting gives rise to levels which are separated by an energy of the order of a few tenths of eV [44]. Furthermore, the energy of the molecular orbital can be tuned by varying the gate voltage V G .
One-level transport model has been adopted to interpret experimental data of C 60 molecular junctions [17] neglecting altogether the effect of electron-electron and electron-vibrations interactions. This model is clearly valid for energies close to the resonance, therefore it is particularly useful in the case of the experiments in [17] where the molecular energy is tuned around the Fermi energy of the leads. Moreover, the one-level model has to be used in the regime of low temperatures, therefore temperatures up to room temperature can be considered for the interpretation of experimental data. Within this model, the energy-dependent transmission function T(E) is assumed to be well approximated by a Lorentzian function: where the molecular level energy is taken as with E 0 the energetic separation of the dominant transport level with respect to the chemical potential µ, and α the effectiveness of gate coupling. The expression of takes clearly into account the tuning of the molecular level by the gate voltage. By using Equation (18), in the limit of low temperature of the Landauer-Büttiker approach valid in the coherent regime [1,41], the gate voltage-dependent electrical conductance G becomes where G 0 = 2e 2 /h is the quantum of conductance, with h Planck constant. Moreover, in the same limit, the Seebeck coefficient S is where k B is the Boltzmann constant. We remark that k B /|e| 86.17 µV/K sets the order of magnitude (and, typically, the maximum value in modulus) of the thermopower in molecular junctions. In the right panel of Figure 1, we report the experimental data of Seebeck coefficient S as a function of the gate voltage V G taken from [17] for C 60 junctions. The values of S taken at the temperature T = 100 K are quite large in modulus for negative gate. Moreover, the data show a marked change as a function of the gate voltage suggesting that the chemical potential is able to cross a level of the molecule. Since the values of S are negative for small values of V G and are still negative for zero V G , the charge transport is dominated by the LUMO level of C 60 .
Actually, to fit the experimental data shown in the right panel of Figure 1, Equation (21) has been used, getting the positive value E 0 − µ = 0.057 eV [17]. For the optimization of the fit, in the same paper [17], Γ = 0.032 eV and the gate voltage effectiveness α = 0.006 eV/V are also extracted. These three numerical values put in Equation (21) provide the fit curve shown in the right panel of Figure 1. The fit is good, but not excellent.
In the left panel of Figure 1, we report the experimental data of the charge conductance G as a function of the gate voltage V G taken again from experimental data of [17] for C 60 junctions. Even if the temperature is not high (T = 100 K), the values of G are quite smaller than the conductance quantum G 0 . Moreover, if one uses the parameters (E 0 − µ = 0.057 eV, Γ = 0.032 eV, and α = 0.006 eV/V) extracted from the Seebeck data in [17] and reproduced in the right panel of Figure 1, one finds a peak of the conductance for E 0 − µ = αV G , hence for V G 9 V. This is in contrast with the peak of G which occurs at V G 5 V in the experimental data. If we try to describe the experimental data shown in the left panel of Figure 1 by using Equation (20) and the parameters extracted by fitting the Seebeck data, we get the red line reported in the left panel of Figure 1. It is evident that the agreement between theory and data is poor, and, in particular, the maximum observed for V G around 5 V is not recovered. We remark that k B T 0.0086 eV represents the smallest energy scale apart from values of V G very close to the LUMO level. Therefore, the quality of the comparison cannot depend on the low temperature expansion used in Equation (20).
To improve the interpretation of the experimental data, in this paper, we analyze the role of many-body interactions between molecular degrees of freedom. For example, experimental measurements have highlighted that the effects of the electron-vibration interactions are not negligible in junctions with C 60 molecules and gold electrodes [10,23]. In particular, experimental results for C 60 molecules [23] provide compelling evidence for a sizable coupling between the electrons and the center of mass vibrational mode. Indeed, previous studies have shown that a C 60 molecule is held tightly on gold by van der Waals interactions, which can be expressed by the Lennard-Jones form. The C 60 -gold binding near the equilibrium position can be approximated very well by a harmonic potential with angular frequency ω 0 . For C 60 molecules, the center of mass energyhω 0 has been estimated to be of the order of 5 meV. For both, see [17] relative to C 60 molecular junctions.
In this paper, we focus on the center of mass mode as the relevant low frequency vibrational mode for the molecule. The center of mass mode is expected to have the lowest angular frequency ω 0 for large molecules. For fullerene, the energyhω 0 is still smaller than the thermal energy k B T corresponding to the temperature T = 100 K fixed for the measurements made in [17]. For k B T ≥hω 0 , the self-consistent adiabatic approach introduced in the previous section can be used for a non-perturbative treatment of the electron-vibration coupling. Equation (7) reduces in this case to a single Langevin equation [33,36]. We hereby report the expression for the displacement dependent electronic spectral function A(E, x) within these assumptions, in Equation (2), the interaction HamiltonianĤ int reduces to the same interaction term of the single impurity Anderson-Holstein model [1] and the electron-oscillator coupling sets the characteristic polaron energy E P with m mass of the molecule. Actually, an additional electron injected from the leads compresses the C 60 -surface bond shortening the C 60 -surface distance, but not significantly changing the vibrational frequency. Previous studies [10,23] have estimated that the number of vibrational quanta typically excited by the tunnelling electron in fullerene junctions is not large. Therefore, intermediate values of electron-vibration energy E P corresponding to values comparable with Γ are considered relevant for fullerene molecular junctions. Taking the parameters extracted from the experimental data discussed above, E P 0.030 eV sets the order of magnitude.
To improve the analysis of the fullerene molecular junction, in this paper, we study also the role of electron-electron interactions acting onto the molecule. Indeed, the conductance gap observed in the data of C 60 molecules can be interpreted using ideas borrowed from the Coulomb blockade effect [1,23]. Therefore, these features are understood in term of the finite energy required to add (remove) an electron to (from) the molecule. Within the single-level model introduced in the previous section, this charging energy is simulated by fixing the value of the local Hubbard term U in Equation (2). The maximum conductance gap observed in the experimental data [23] indicates that the charging energy of the C 60 molecule can be around 0.27 eV, therefore experiments set the order of magnitude U 0.3 eV.
To include the Coulomb blockade effect within the adiabatic approach discussed previously, we generalize it to the case in which the electronic level can be double occupied and a strong Coulomb repulsion U is added together with the electron-vibration interaction. The starting point is the observation that, in the absence of electron-oscillator interaction, and in the limit where the coupling of the dot to the leads is small Γ << U [41], the single particle electronic spectral function is characterized by two spectral peaks separated by an energy interval equal to U. In the adiabatic regime, one can independently perturb each spectral peak of the molecule [40], obtaining at the zero order of the adiabatic approach where ρ(x) is the electronic level density per spin. In our computational scheme, ρ(x) has to be self-consistently calculated for a fixed displacement x of the oscillator through the following integral The above approximation is valid if the electron-oscillator interaction is not too large, such that Γ E P << U and the two peaks of the spectral function can be still resolved [40]. We remark that, in comparison with our previous work [40], parameters appropriate for junctions with C 60 molecules are considered in this paper focusing on the temperature T = 100 K fixed for the measurements made in [17], smaller than the room temperature, where the adiabatic approach can be still adopted. Therefore, the approach is valid in the following parameter regime:hω 0 ≤ k B T < Γ U [39,40]. Within the adiabatic approach, the actual electronic spectral function A(E) results from the average over the dynamical fluctuations of the oscillator motion, therefore, as a general observable, it is calculated by using Equation (17): where P(x) is the reduced position distribution function of the oscillator. Notice that, in the absence of electron-electron (U = 0) and electron-vibration (E P = 0) interactions, the spectral function is proportional to the transmission T(E) given in Equation (18) through the hybridization width Γ: In the linear response regime (bias voltage V bias → 0 + and temperature difference ∆T → 0 + ), all the electronic transport coefficients can be expressed as integrals of A(E). To this aim, we report the conductance G where A(E) is the spectral function defined in Equation (25), with f (E) = 1/(exp [β(E − µ)] + 1) the free Fermi distribution corresponding to the chemical potential µ and the temperature T, and β = 1/k B T. The Seebeck coefficient is given by S = −G S /G, where the charge conductance G has been defined in Equation (26), and Then, we calculate the electron thermal conductance G el K = G Q − TGS 2 , with Therefore, in the linear response regime, one can easily evaluate the electronic thermoelectric figure of merit ZT el which characterizes the electronic thermoelectric conversion. We recall that, in this paper, we do not consider the addition contribution coming from phonon thermal conductance G ph K .

Results
In this section, we discuss the thermoelectric properties within the single-level model analyzing the role of the electron-electron and electron-vibration interactions between the molecular degrees of freedom. We point out that only the combined effect of these interactions is able to provide a good agreement between experimental data and theoretical calculations.
The level density ρ is shown in the upper left panel of Figure 2, the charge conductance G in the upper right panel, the Seebeck coefficient S in the lower left panel, and the electronic thermoelectric figure of merit ZT el in the lower right panel. All quantities are plotted as a function of level energy at the temperature T = 100 K. For all the quantities, we first analyze the coherent regime (black solid lines in the four panels of Figure 2), which means absence of electron-electron and electron-vibration interactions. Then, we study the effect of a finite electron-vibration coupling E P (red dash lines in the four panels of Figure 2) focusing on the intermediate coupling regime. Finally, we consider the combined effect of electron-vibration and electron-electron interactions for all the quantities (blue dash-dot lines in the four panels of Figure 2) analyzing the experimentally relevant regime of a large Coulomb repulsion U. The level density ρ per spin reported in the upper left panel of Figure 2 shows the expected decreasing behavior with increasing the level energy . The electron-vibration interaction induces a shift of the curve of about E P . In the presence of electron-electron interactions, the behavior is more complex. Actually, in molecular junctions, the strong Coulomb repulsion usually reduces the electronic charge fluctuations and suppresses the double occupation of the electronic levels [1]. For values of smaller than −U, the density is closer to unity, while, for larger than zero, the density vanishes. For between −U and 0, there is a plateau with a value of the density close to 0.5. Indeed, these phenomena are characteristic of Coulomb blockade effects.
The conductance G is shown in upper right panel of Figure 2. At low temperatures, this quantity as a function of the level position provides essentially the spectral function of the molecular level. Indeed, in the coherent low temperature regime, G can be directly related to the transmission with a Lorentzian profile. One of the main effects of an adiabatic oscillator is to shift the conductance peak towards positive energies proportional to the electron-vibration coupling energy E P . Apparently, another expected effect is the reduction of the peak amplitude. In fact, electron-vibration couplings on the molecule tend to reduce the charge conduction. As a consequence, electron-vibration couplings induces somewhat longer tails far from the resonance. These features, such as the peak narrowing, are common to other theoretical approaches treating electron-vibration interactions, among which that related to the Franck-Condon blockade [45]. Actually, in a previous paper [39], we successfully compared the results of the adiabatic approximation with those of the Franck-Condon blockade formalism in the low density limit where this latter approach becomes essentially exact [1].
We note that a finite electron-electron interaction not only suppresses the electronic conduction for small values of , but it is also responsible for a second peak centered at −U. In fact, there is a transfer of spectral weight from the main peak to the interaction-induced secondary peak. We stress that these features are compatible with experimental data since conductance gap ascribed to Coulomb blockade effects have been measured in fullerene junctions [1,23].
We investigate the properties of the Seebeck coefficient S of the junction in the lower left panel of Figure 2. In analogy with the behavior of the conductance, the main effect of the electron-vibration interaction is to reduce the amplitude of the Seebeck coefficient. Moreover, the shift of the zeroes of S is governed by the coupling E P as that of the peaks of G. Therefore, with varying the level energy , if G reduces its amplitude, S increases its amplitude in absolute value, and vice versa. This behavior and the values of S are in agreement with experimental data [4,17]. In the Coulomb blockade regime, S shows a peculiar oscillatory behavior as a function of the energy , with several positive peaks and negative dips. The energy distance between the peaks (or the dips) is governed by the Hubbard term U. Even in this regime, the Seebeck coefficient S is negligible for the level energies where the electronic conductance presented the main peaks, that is at 0 and −U. This property turns out to be a result of the strong electron-electron interaction U [26]. In any case, close to the resonance (zero values of the level energy ), the conductance looks more sensitive to many-body interactions, while the Seebeck coefficient appears to be more robust.
The electronic conductance G, Seebeck coefficient S, and electron thermal conductance G el K combine in giving an electronic figure of merit ZT el . This latter quantity is shown in lower right panel of Figure 2 at the temperature T = 100 K. We stress that, due to the low value of the temperature, the quantity ZT el does not show values comparable with unity. However, it is interesting to analyze the effects of many-body interactions on this quantity. A finite value of the electron-vibration coupling E P leads to a reduction of the height of the figure of merit peaks. It is worth noting that the position of the peaks in ZT el roughly coincides with the position of the peaks and dips of the Seebeck coefficient S. Finally, the electron-electron interactions tend to reduce the amplitude and to further shift the peaks of the figure of merit. After the analysis of the effects of many-body interactions on the charge conductance and Seebeck coefficient, we can make a comparison with the experimental data available in [17] and shown in Figure 1. These data are plotted again in Figure 3 together with the fit discussed in Figure 1. We recall that for fullerene junctions the one-level model discussed in the previous section is characterized by the following parameters: E 0 − µ = 0.057 eV, Γ = 0.032 eV, and α = 0.006 eV/V. We remark that the level energy used in the previous discussion is related to the energy E 0 and the gate voltage V G through Equation (19). Therefore, once the value of E 0 is fixed, one can switch from the energy to the gate voltage V G . Before introducing many-body effects, we consider a slight shift of the level position considering the case E 0 − µ = 0.065 eV reported in Figure 3. This energy shift is introduced to counteract the shifts of the peaks (conductance) or zeroes (Seebeck) introduced by many-body interactions which, in addition, reduce the amplitudes of response functions. The aim of this paper is to provide an optimal description for both charge conductance G and Seebeck coefficient S.  Starting from the level energy E 0 − µ = 0.065 eV, in Figure 3, we analyze the effect of the electron-vibration coupling in the intermediate regime E P = 0.018 eV. The shift induced in the zero of the Seebeck coefficient is still compatible with experimental data. Moreover, the electron-vibration interaction shifts and reduces the peak of the charge conductance in an important way. However, this is still not sufficient to get an accurate description of the charge conductance. One could increase the value of the coupling energy E P , but, this way, the shift of the conductance peak becomes too large with a not marked reduction of the spectral weight.
Another ingredient is necessary to improve the description of both conductance G and Seebeck coefficient S. In our model, the additional Coulomb repulsion plays a concomitant role. Its effects poorly shift the zero of the Seebeck coefficient and slightly modifies the curve far from the zero. Therefore, the description of the Seebeck coefficient remains quite accurate as a function of the gate voltage. On the other hand, it provides a sensible reduction of the conductance amplitude with a not large shift of the peak. Hence, the effects of Hubbard term are able to improve the theoretical interpretation of the experimental data for the conductance G and the Seebeck coefficient S close to the resonance. Far from the resonance, in a wider window of gate voltages, theory predicts the existence of a secondary peak of the conductance and a complex behavior of the Seebeck coefficient due to Coulomb blockade effects. The features are not negligible as a function of the gate voltage.
As far as we know, experimental measurements of the electronic thermal conductance G el K have become only very recently available [21]. Indeed, it is important to characterize this quantity since it allows determining the thermoelectric figure of merit. Therefore, in Figure 4, we provide the theoretical prediction of the electronic thermal conductance G el K as a function of V G starting from the optimized values of the one-level parameters used to describe both charge conductance and Seebeck coefficient in an accurate way. We stress that the plotted thermal conductance is expressed in terms of the thermal conductance quantum g 0 (T) = π 2 k 2 B T/(3h) [46]. The main point is that, in the units chosen in Figure 4, the thermal conductance G el K shows a strong resemblance with the behavior of the charge conductance G in units of the conductance quantum G 0 as a function of the gate voltage V G . We remark that, at T = 100 K, g 0 (T) 9.456 × 10 −11 (W/K) 100 pW/K. The values of the thermal conductance G el K shown in Figure 4 are fractions of g 0 (T), therefore they are fully compatible with those estimated experimentally in hydrocarbon molecules [19] (50 pW/K).  . Electronic thermal conductance G el K in units of thermal conductance quantum g 0 (T) (g 0 (T) = π 2 k 2 B T/(3h)) as a function the voltage gate V G in units of Volt at the temperature T = 100 K: coherent results (black solid line) corresponding to one-level model with energy E 0 − µ = 0.065 eV, effect of the only electron-vibration coupling E P = 0.018 eV (red dash line), and effect of additional electron-electron interaction U = 0.3 eV (blue dash-dot line).

Conclusions
In this paper, we have theoretically analyzed the role of electron-vibration and electron-electron interactions on the thermoelectric properties of molecular junctions focusing on devices based on fullerene. We have used a self-consistent adiabatic approach which allows a non-perturbative treatment of the electron coupling to low frequency vibrational modes, such as those of the molecular center of mass between metallic electrodes. This approach incorporates Coulomb blockade effects due to strong electron-electron interaction between molecular degrees of freedom. We have analyzed a one-level model which takes into account the LUMO level of fullerene and its alignment to the chemical potential. We have stressed that an accurate description of the transport properties is obtained in the intermediate regime for the electron-vibration coupling and in the strong coupling regime for the electron-electron interaction. Moreover, we have demonstrated that only the combined effect of electron-vibration and electron-electron interactions is able to predict the correct behavior of both the charge conductance and the Seebeck coefficient. The theoretical calculations presented in this paper show a very good agreement with available experimental data of both charge conductance and Seebeck coefficient.
In this paper, we have used a one-level transport model as a starting point to address the role of many-body interactions between molecular degrees of freedom. This model is frequently used in all cases where the energy levels can be tuned around the chemical potential and additional spectral features are absent [4,17]. This is the case of the experiments in [17] for the fullerene junctions analyzed in this paper. The one-level model is expected to be valid for energies close to the Fermi level and low temperatures. Actually, a more realistic description of the molecule and its coupling with metallic leads is needed if more complex transport phenomena take place, in particular interference effects [47,48] recently investigated in molecular junctions. Moreover, inclusion of quantum corrections to oscillator dynamics can be important to explore the effects of additional vibrational modes and further electron-vibration regimes [49] (from adiabatic to anti-adiabatic ones) and their relation with strong electron-electron interactions [50,51].