**3. 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 *VG*.

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:

$$T(E) = \frac{4\Gamma^2}{(E-\epsilon)^2 + 4\Gamma^2} \tag{18}$$

where the molecular level energy  is taken as

$$
\epsilon = E\_0 - \mathfrak{a}V\_{\mathbb{G}\prime} \tag{19}
$$

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

$$G = \frac{\partial I}{\partial V\_{bias}}(V\_{bias} = 0, V\_G) = G\_0 \ T(E = \mu), \tag{20}$$

where *G*0 = <sup>2</sup>*e*2/*h* is the quantum of conductance, with *h* Planck constant. Moreover, in the same limit, the Seebeck coefficient *S* is

$$S = -\frac{\pi^2}{3} \frac{k\_B}{|\mathfrak{e}|} k\_B T \frac{\partial \ln T(E = \mu)}{\partial E} = \frac{\pi^2}{3} \frac{k\_B}{|\mathfrak{e}|} k\_B T \frac{2[\mu - \mathfrak{e}]}{[(\mu - \mathfrak{e})^2 + 4\Gamma^2]} \tag{21}$$

where *kB* is the Boltzmann constant. We remark that *kB*/|*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 *VG* 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 *VG* and are still negative for zero *VG*, 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 *VG* 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 − *μ* = *<sup>α</sup>VG*, hence for *VG* 9 V. This is in contrast with the peak of *G* which occurs at *VG* 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 ge<sup>t</sup> the red line reported in the left panel of Figure 1. It is evident that the agreemen<sup>t</sup> between theory and data is poor, and, in particular, the maximum observed for *VG* around 5 V is not recovered. We remark that *kBT* 0.0086 eV represents the smallest energy scale apart from values of *VG* 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 energy *h*¯ *ω*0 has been estimated to be of the order of 5 meV.

**Figure 1.** (**Left**) Conductance *G* (in units of conductance quantum *G*0) as a function of the gate voltage *VG* (in units of V) at *T* = 100 K from experimental data (black circles, see [17] for fullerene *C*60 junctions) and from a curve (red solid line) obtained by using the parameters of the fit to the Seebeck coefficient. (**Right**) Seebeck coefficient *S* (in units of μK/V) as a function of the gate voltage *VG* (in units of V) at *T* = 100 K from experimental data (black circles), and from a fit (red solid line). 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 energy *h*¯ *ω*0 is still smaller than the thermal energy *kBT* corresponding to the temperature *T* = 100 K fixed for the measurements made in [17]. For *kBT* ≥ *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 *<sup>A</sup>*(*<sup>E</sup>*, *x*)

$$A(E, \mathbf{x}) = \frac{4\Gamma}{(E - \epsilon - \lambda\mathbf{x})^2 + 4\Gamma^2}. \tag{22}$$

within these assumptions, in Equation (2), the interaction Hamiltonian *H* ˆ *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 *EP*

$$E\_P = \frac{\lambda^2}{2m\omega\_0^2},\tag{23}$$

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 *EP* corresponding to values comparable with Γ are considered relevant for fullerene molecular junctions. Taking the parameters extracted from the experimental data discussed above, *EP* 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

$$A(E, \mathbf{x}) = [1 - \rho(\mathbf{x})] \frac{4\Gamma}{(E - \epsilon - \lambda \mathbf{x})^2 + 4\Gamma^2} + \rho(\mathbf{x}) \frac{4\Gamma}{(E - \epsilon - \lambda \mathbf{x} - \mathcal{U})^2 + 4\Gamma^2} \tag{24}$$

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 *ρ*(*x*) = +∞ −∞ *dE*2*πi <sup>G</sup>*<sup>&</sup>lt;(*<sup>E</sup>*, *<sup>x</sup>*), with the lesser Green function *<sup>G</sup>*<sup>&</sup>lt;(*<sup>E</sup>*, *x*) = *i*2 [ *fL*(*E*) + *fR*(*E*)]*A*(*<sup>E</sup>*, *<sup>x</sup>*). The above approximation is valid if the electron–oscillator interaction is not too large, such that Γ *EP* << *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<sup>ω</sup>*0 ≤ *kBT* < Γ *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):

$$A(E) = \int\_{-\infty}^{+\infty} dx P(x)A(E, x),\tag{25}$$

where *<sup>P</sup>*(*x*) is the reduced position distribution function of the oscillator. Notice that, in the absence of electron–electron (*U* = 0) and electron–vibration (*EP* = 0) interactions, the spectral function is proportional to the transmission *T*(*E*) given in Equation (18) through the hybridization width Γ: *T*(*E*) = <sup>Γ</sup>*<sup>A</sup>*(*E*).

In the linear response regime (bias voltage *Vbias* → 0<sup>+</sup> and temperature difference Δ*T* → 0<sup>+</sup>), all the electronic transport coefficients can be expressed as integrals of *<sup>A</sup>*(*E*). To this aim, we report the conductance *G* 

$$G = G\_0 \Gamma \int\_{-\infty}^{+\infty} dE A(E) \left[ -\frac{\partial f(E)}{\partial E} \right],\tag{26}$$

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/*kBT*. The Seebeck coefficient is given by *S* = <sup>−</sup>*GS*/*G*, where the charge conductance *G* has been defined in Equation (26), and

$$G\_S = G\_0 \left(\frac{k\_B}{\varepsilon}\right) \Gamma \int\_{-\infty}^{+\infty} dE \frac{(E-\mu)}{k\_B T} A(E) \left[ -\frac{\partial f(E)}{\partial E} \right]. \tag{27}$$

Then, we calculate the electron thermal conductance *Ge<sup>l</sup> K* = *GQ* − *TGS*2, with

$$\mathcal{G}\_{Q} = \mathcal{G}\_{0} \left(\frac{k\_{B}}{\varepsilon}\right)^{2} \Gamma T \int\_{-\infty}^{+\infty} dE \left[\frac{E-\mu}{k\_{B}T}\right]^{2} A(E) \left[-\frac{\partial f(E)}{\partial E}\right]. \tag{28}$$

Therefore, in the linear response regime, one can easily evaluate the electronic thermoelectric figure of merit *ZTe<sup>l</sup>*

$$ZT^{cl} = \frac{GS^2T}{G\_K^{cl}},\tag{29}$$

which characterizes the electronic thermoelectric conversion. We recall that, in this paper, we do not consider the addition contribution coming from phonon thermal conductance *Gp<sup>h</sup> K*.
