**2. 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

$$
\hat{\mathcal{H}} = \hat{\mathcal{H}}\_{mol} + \hat{\mathcal{H}}\_{lends} + \hat{\mathcal{H}}\_{lends\text{--}mol\text{-}\prime} \tag{1}
$$

where *H*ˆ *mol* is the Hamiltonian describing the molecular degrees of freedom, *<sup>H</sup>*<sup>ˆ</sup>*leads* the leads' degrees of freedom and *<sup>H</sup>*<sup>ˆ</sup>*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 *H*ˆ *mol* is

$$
\hat{H}\_{mol} = \sum\_{m,l,\sigma} \hat{\varepsilon}\_{m,\sigma}^{\dagger} \varepsilon\_{\sigma}^{m,l} \hat{\varepsilon}\_{l,\sigma} + \mathcal{U} \sum\_{m,l} \hat{n}\_{m,\uparrow}^{\dagger} \hat{n}\_{l,\downarrow} + \hat{H}\_{\text{osc}} + \hat{H}\_{\text{int}} \tag{2}
$$

where *cm*,*<sup>σ</sup>* (*c*† *<sup>m</sup>*,*<sup>σ</sup>*) 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 *ε <sup>m</sup>*,*l σ* is assumed diagonal in spin space, *<sup>n</sup>*<sup>ˆ</sup>*l*,*<sup>σ</sup>* = *c*† *<sup>l</sup>*,*<sup>σ</sup>cl*,*<sup>σ</sup>* 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 *ε <sup>m</sup>*,*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

$$
\hat{H}\_{\text{osc}} = \sum\_{s} \frac{\mathfrak{f}\_s^2}{2m\_s} + V(X),
\tag{3}
$$

where *s* = (1, . . . , *<sup>N</sup>*), with *N* being the total number of vibrational modes; *ms* is the effective mass associated with the *s*th vibrational mode; and *p*ˆ*s* is its momentum operator. Moreover, *V*(*X*) = 12 ∑*s ksx*ˆ2*s* is the harmonic potential (with *ks* the spring constants, and the oscillator frequencies *<sup>ω</sup><sup>s</sup>*0 = √*ks*/*ms*), *x*ˆ*s* is the displacement operator of the vibrational mode *s*, and *X* = (*x*ˆ1,..., *<sup>x</sup>*<sup>ˆ</sup>*N*) indicates all the displacement operators.

In Equation (2), the electron–vibration coupling *H* ˆ *int* is assumed linear in the vibrational displacements and proportional to the electron level occupations

$$
\hat{H}\_{\rm int} = \sum\_{s,l} \lambda\_{s,l} \hat{\mathbf{x}}\_s \hat{\mathbf{n}}\_{l,s} \tag{4}
$$

where *s* = (1, . . . , *N*) indicates the vibrational modes of the molecule, *l* = (1, . . . , *M*) denotes its electronic levels, *n*ˆ*l* = ∑*σ nl*,*<sup>σ</sup>* is the electronic occupation operator of the level *l*, and *<sup>λ</sup><sup>s</sup>*,*<sup>l</sup>* is a matrix representing the electron–vibrational coupling.

In Equation (1), the Hamiltonian of the electron leads is given by

$$
\hat{H}\_{\text{lrads}} = \sum\_{k,\mu,\sigma} \varepsilon\_{k,\mu} \hat{\mathfrak{c}}\_{k,\mu,\sigma}^{\dagger} \hat{\mathfrak{c}}\_{k,\mu,\sigma} \tag{5}
$$

where the operators *<sup>c</sup>*ˆ†*k*,*α*,*<sup>σ</sup>*(*c*<sup>ˆ</sup>*k*,*α*,*<sup>σ</sup>*) create (annihilate) electrons with momentum *k*, spin *σ*, and energy *<sup>ε</sup>k*,*<sup>α</sup>* = *Ek*,*<sup>α</sup>* − *μα* in the left (*α* = *L*) or right (*α* = *R*) leads. The left and right electron leads are considered as thermostats in equilibrium at the temperatures *TL* = *T* + Δ*T*/2 and *TR* = *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 *fL*(*E*) and *fR*(*E*), respectively, with *E* the energy. The difference of the electronic chemical potentials in the leads provides the bias voltage *Vbias* applied to the junction: *μL* = *μ* + *eVbias*/2, *μR* = *μ* − *eVbias*/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 *Vbias*and temperature Δ*T*.

 Finally, in Equation (1), the coupling between the molecule and the leads is described by

$$
\hat{H}\_{mol-lcds} = \sum\_{k,\mu,m,\sigma} (V\_{k,\mu}^{m} \mathfrak{E}\_{k,\sigma}^{\dagger} \mathfrak{E}\_{m,\sigma} + h.c.), \tag{6}
$$

where the tunneling amplitude between the molecule and a state *k* in the lead *α* has the amplitude *Vmk*,*α*. For the sake of simplicity, we suppose that the density of states *ρk*,*<sup>α</sup>* for the leads is flat within the wide-band approximation: *ρk*,*<sup>α</sup>* → *ρα*, *Vmk*,*<sup>α</sup>* → *Vmα* . Therefore, the full hybridization width matrix of the molecular orbitals is Γ*<sup>m</sup>*,*<sup>n</sup>* = ∑*α* Γ*<sup>m</sup>*,*<sup>n</sup> α* = ∑*α* Γ*<sup>m</sup>*,*<sup>n</sup> α* , with the tunneling rate Γ*<sup>m</sup>*,*<sup>n</sup> α* = <sup>2</sup>*πραVm*<sup>∗</sup> *α Vnα* . 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: *<sup>ω</sup><sup>s</sup>*0 Γ*<sup>m</sup>*,*n*, for each *s* and all pairs of (*<sup>m</sup>*, *<sup>n</sup>*). 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 *ε<sup>m</sup>* → *ε<sup>m</sup>* + *λmxm*, where *xm* 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 \ddot{\mathbf{x}}\_s + k\_s \mathbf{x}\_s = F\_s(t) + \mathbf{j}\_s(t),\tag{7}$$

where the generalized force *Fs* is due to the effect of all electronic degrees of freedom through the electron–vibration coupling [32,39]:

$$F\_s^{el}(t) = Tr[i\lambda\_s \mathbf{G}^{<}(t, t)],\tag{8}$$

with the trace "Tr", taken over the molecule levels, defined in terms of the lesser molecular matrix Green's function *<sup>G</sup>*<sup>&</sup>lt;(*<sup>t</sup>*, *t* ) with matrix elements *G*<sup>&</sup>lt; *<sup>m</sup>*,*l* (*t*, *t* ) = *<sup>i</sup>c*† *<sup>m</sup>*,*<sup>σ</sup>*(*t*)*cl*,*<sup>σ</sup>*(*<sup>t</sup>* ). 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<sup>R</sup>* is given by

$$\mathcal{G}^{R} \simeq \mathcal{G}^{R} + \frac{i}{2} (\partial\_{E} \mathcal{G}^{R} (\sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s}) \mathcal{G}^{R} - \mathcal{G}^{R} (\sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s}) \partial\_{E} \mathcal{G}^{R}),\tag{9}$$

where G*<sup>R</sup>*(*<sup>E</sup>*, *X*) is the strictly adiabatic (frozen) retarded Green's function including the coupling with the leads

$$\mathcal{G}^{\mathbb{R}}(E, \mathbf{X}) = [E - \varepsilon(\mathbf{X}) - \mathbf{Z}^{R, \text{lends}}]^{-1},\tag{10}$$

*ε*(*X*) represents the matrix *ε <sup>m</sup>*,*l σ* + ∑*s <sup>λ</sup>sxsδl*,*<sup>m</sup>* and **Σ***R*,*leads* = ∑*α* **<sup>Σ</sup>***R*,*leads α* is the total self-energy due to the coupling between the molecule and the leads. For the lesser Green's function *G*<sup>&</sup>lt;, the adiabatic approximation involves

$$\begin{split} \mathcal{G}^{<\cdot} &\simeq \quad \mathcal{G}^{<\cdot} + \frac{i}{2} \Big[ \partial\_{\mathcal{E}} \mathcal{G}^{<} \big( \sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s} \big) \mathcal{G}^{A} - \mathcal{G}^{R} \big( \sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s} \big) \partial\_{\mathcal{E}} \mathcal{G}^{<} \\ &+ \partial\_{\mathcal{E}} \mathcal{G}^{R} \big( \sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s} \big) \mathcal{G}^{<\cdot} - \mathcal{G}^{<} \big( \sum\_{s} \lambda\_{s} \dot{\mathbf{x}}\_{s} \big) \partial\_{\mathcal{E}} \mathcal{G}^{A} \Big], \end{split} \tag{11}$$

with G< = G*<sup>R</sup>***Σ**<sup>&</sup>lt;G *A*.

The electron–vibration induced forces at the zero order of the adiabatic limit ( *G*<sup>&</sup>lt; G<) are given by

$$F\_s^{el(0)}(X) = -k\_s \mathbf{x}\_s - \int \frac{dE}{2\pi i} tr[\lambda\_s \mathcal{G}^{<}].\tag{12}$$

The leading order correction to the lesser Green's function *G*<sup>&</sup>lt; provides a term proportional to the vibrational velocity

$$F\_s^{cl(1)}(X) = -\sum\_{s'} \theta\_{s,s'}(X) \dot{x}\_{s'} \tag{13}$$

where the tensor *θ* can be split into symmetric and anti-symmetric contributions [32]: *θ* = *θsym* + *θa* , where we have introduced the notation {*Cs*,*<sup>s</sup>*}*sym*,*<sup>a</sup>* = 1 2 {*Cs*,*s* ± *Cs*,*<sup>s</sup>*}*sym*,*<sup>a</sup>* 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

$$
\langle \langle \xi\_s^{\rm cl}(t) \xi\_{s'}^{\rm cl}(t') \rangle \rangle = \operatorname{tr} \{ \lambda\_s \mathbf{G}^{\gtrless}(t, t') \lambda\_{s'} \mathbf{G}^{\lessgtr}(t', t) \}, \tag{14}
$$

where *<sup>G</sup>*<sup>&</sup>gt;(*<sup>t</sup>*, *t* ) is the greater Green's function with matrix elements *G*<sup>&</sup>gt; *<sup>m</sup>*,*l* (*t*, *t* ) = <sup>−</sup>*<sup>i</sup>cm*,*<sup>σ</sup>*(*t*)*c*† *<sup>l</sup>*,*<sup>σ</sup>*(*<sup>t</sup>* ). 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:

$$\langle \zeta\_s^{\text{rel}}(t) \zeta\_{s'}^{\text{rel}}(t') \rangle \simeq \text{tr}\{\lambda\_s \mathcal{G}^>(X,t) \lambda\_{s'} \mathcal{G}^<(X,t) \} = \mathcal{D}(X)\delta(t-t'),\tag{15}$$

where

$$D\_{s,s'}(X) = \int \frac{dE}{2\pi} \text{tr}\left\{\lambda\_s \mathcal{G}^{<} \lambda\_{s'} \mathcal{G}^{>}\right\}\_{sym}.\tag{16}$$

Once the forces and the noise terms are calculated, Equation (7) represents a set of nonlinear Langevin equations in the unknown *xs*. 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 *<sup>P</sup>*(*<sup>X</sup>*, *V*) (where *V* = *X*˙ = (*<sup>v</sup>*1,..., *vN*)), and, therefore, all the properties of the vibrational modes. Using this function, one can determine the average *O* of an electronic or vibrational observable *<sup>O</sup>*(*<sup>X</sup>*, *<sup>V</sup>*):

$$\mathcal{O} = \int \int d\mathbf{X} d\mathbf{V} P(\mathbf{X}, \mathbf{V}) \mathcal{O}(\mathbf{X}, \mathbf{V}).\tag{17}$$

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].
