3 B. Non-equilibrium electron Green’s functions and electron-vibron self-energies where the interacting electrons in the central region are created (annihilated) on electronic level n by the opera- tors d † n ( d$_{n}$ ). In principle, our Hamiltonian for the central region may contain electron-electron interactions, and is built from a complete set of single-electron creation and anni- hilation operators. In our current work we do not include any electron-electron interactions, and so the electronic part of the total Hamiltonian for the central region H$_{C}$ becomes Non-equilibrium Green’s functions (NEGF) within the Keldysh formalism 66–69 represent an extremely useful tool for studying the non-equilibrium properties of many- particle systems. The Green’s functions are functions of two space-time coordinates, and are obviously more com- plicated than the one-particle density which is the main ingredient of density-functional-based theories. One of the great advantage of NEGF techniques is that one can improve the calculations in a systematic way by taking into account specific physical processes (represented by Feynman diagrams) which is what we do in this paper for the electron-vibron interaction. The Green’s functions provide us directly with all expectation values of one- body operators (such as the density and the current), and also the total energy, the response functions, spec- tral functions, etc. Meanwhile, the vibron part of the Hamiltonian is rep- resented by where a † λ creates (annihilates) a vibron in vibron mode λ with frequency ω$_{λ}$ . The electron-vibron coupling term is taken to be linear in the vibron displacement and its most general expression is then given by 16 In Appendix A, we provide more details about NEGF and how to obtained the electron-vibron self-energies from a Feynman diagrammatic expansion of the electron- vibron interaction. We now briefly describe how we apply the NEGF formalism to the SSSM model. H e C = ∑ n ε$_{n}$d † $_{n}$d$_{n}$ . (6) H$_{vib}$ = ∑ λ ℏ ω$_{λ}$a † $_{λ}$a$_{λ}$ , (7) H$_{e}$$_{-vib}$ = ∑ λ,n,m γ$_{λnm}$ ( a † λ + a$_{λ}$ ) d † $_{n}$d$_{m}$ , (8) G r,a = g r,a C + g r,a C Σ $^{r,a}$G$^{r,a}$, (11) G >,< = (1 + G $^{r}$Σ$^{r}$) g >,< C (1 + Σ $^{a}$G$^{a}$) + G $^{r}$Σ$^{>,<}$G$^{a}$. (12) Σ x = Σ x L + Σ x R + Σ x $_{int}$. (13) The Green’s functions are calculated via Dyson-like equations. For the retarded and advanced Green’s func- tions G r,a these are where g r,a, C is the non-interacting Green’s function for the isolated central region. For the greater G > and lesser G < Green’s functions, we use a quantum kinetic equation of the form where γ$_{λnm}$ is the coupling constant for exciting the vi- bron mode λ by electronic transition between the elec- tronic levels n and m . We concentrate on the simplest version of the Hamil- tonian of the central part: the single-site single-mode SSSM model, in which one considers just one electron level coupled to one vibration mode. Despite the sim- plicity of this model, remarkably it not only contains all the physics we require but also ensures that we isolate the properties we are interested in without the complication of the added electronic levels. The reasons for this are as follows. Firstly, when the Fermi energy of the leads is pinned around the midgap, then at low and intermedi- ate biases one of the frontier orbitals (either HOMO or LUMO) dominates the transport properties. Secondly in conjugated organic molecules (mostly used in single- molecule junction experiments), it is known that the op- tically relevant vibration modes are mostly coupled to either the HOMO or LUMO levels $^{65}$. We can therefore be confident that our model contains the relevant physics. Here Σ $^{x}$, ( x = r, a, >, < ) is a total self-energy consisting of a sum of the self-energies from the constituent parts of the system: Σ x L,R are the self-energies arising from the non- interacting leads α = L, R and as such are simple to calculate: The total Hamiltonian for the central region thus be- comes where we now have just one electronic level ε$_{n}$ → ε$_{0}$ and just one vibron mode a$_{λ}$ → a , coupled via the electron- vibron coupling constant γ$_{0}$ . The expression for the lead- central-region coupling (equation (5)) also simplifies to become where f$_{α}$ is the Fermi-Dirac distribution for lead α , with Fermi level µ$_{α}$ = µ eq + η$_{α}$eV and temperature T$_{α}$ . The fraction of potential drop at the left contact is η$_{L}$ = ± η$_{V}$ and η$_{R}$ = ∓ (1 − η$_{V}$ ) at the right contact $^{61}$, hence η$_{L}$ − η$_{R}$ = eV is indeed the applied bias, and η$_{V}$ ∈ [0 , 1]. The component of the retarded Green’s function for the isolated (non-interacting) lead α corresponding to the where we have replaced the coupling potentials V$_{αn}$ with hopping integrals t$_{0}$$_{α}$ . Σ r α = t 2 0 $_{α}$g r 0 $_{α}$( ω ) , (14) Σ a α = (Σ r $_{α}$) $^{∗}$, (15) Σ > α = 2 i ℑ m [Σ r $_{α}$( ω )] (1 − f$_{α}$ ( ω )) , (16) Σ < α = − 2 i ℑ m [Σ r $_{α}$( ω )] f$_{α}$ ( ω ) , (17) H$_{C}$ = ε$_{0}$d $^{†}$d + ℏ ω$_{0}$a $^{†}$a + γ$_{0}$ ( a † + a ) d $^{†}$d, (9) V$_{LC}$ + V$_{CR}$ = ∑ α = L,R t$_{0}$$_{α}$ ( c † $_{α}$d + d $^{†}$c$_{α}$ ) , (10)