16 the strong non-linearity involved in solving the problem self-consistently. while containing diagrams higher than second order, will probably break important physical properties of the sys- tem at T = 0, because they lack important electronic processes, to second-order the DX vertex correction-like diagram or DPH diagram respectivelly. In the resonant (or quasi-resonant) regime, polarisa- tion effects and vertex corrections play an equally impor- tant role in the electronic structure and transport prop- erties of electron-vibron interacting nanojunctions which electron-hole symmetry. Furthermore, the effects of the second-order diagrams also depend on the transport regime: at and near equilib- rium, the DX diagram dominates over the DPH diagram in the off-resonant transport regime (essentially because there are not many electron-hole excitations available in this regime). In the resonant case, however, both DX and DPH play an equally important role. For large non- equilibrium conditions ( V > ω ), both second-order di- agrams play an important role, since the corresponding spectral functions look qualitatively similar to those ob- tained the resonant regime. We have thus shown that it is indeed necessary to go beyond SCBA to obtain the correct results for a wide range of parameters. This has also been confirmed by comparing our NEGF results to an exact calcula- tion in terms of electron-vibron interacting (the multi- channel inelastic scattering technique MCIST $^{64}$), though it should be noted that comparison between NEGF and MCIST calculations are only valid for the off-resonant transport regime. In this paper, we have presented a method based on NEGF to calculate the equilibrium and non-equilibrium electronic structure of electron-vibron interacting single- molecule junctions. We have applied the method to a model system which consists of a single electronic level coupled to a single vibration mode in the central region, the latter in contact with two non-equilibrium electron reservoirs. We have also studied in detail the effects of self- consistency on the calculations, especially for the cases including the second-order diagrams (see Section III D 1 and Figure 11). We have found a solution to an old prob- lem, well known in many-body perturbation theory: in order to avoid negative spectral densities when including higher order diagrams in the calculations, it seems more appropriate to start the self-consistent loop with SCBA (Hartree-Fock like) electron Green’s functions . In comparison to previous studies performed within a similar approach $^{20,25,26,28,30–34,38,44–46}$, the novelty of our method lies in the fact that it goes beyond the conventionally-used self-consistent Born approximation (SCBA). Higher-order diagrams for the electron-vibron interacting have been implemented in our calculations. Finally, we have studied the linear response (linear con- ductance G$_{lin}$ ) of the nanojunctions, and have found that in the off-resonant regime, the value of G$_{lin}$ is governed by the behaviour of the tail of the main peak at the Fermi level, hence it depends on both the position and the width of this peak. In the off-resonant regime, the DPH dia- gram contribution is negligible and the DX diagram gives a better polaron shift (better position of the main peak) and hence the second-order (mostly DX) calculation gives a better agreement for G$_{lin}$ with the exact result. For the near-resonant regime, the contribution of DPH will be- come more important and both DX and DPH will play an important role in determining the value of G$_{lin}$ by changing both the position and the width of the main peak. In this paper we have considered the second-order dia- grams which contain two-vibron processes: the so-called double exchange DX diagram which is part of the vertex corrections to SCBA, and the partially dressed DPH dia- gram in which the vibron propagator is renormalized by one electron-hole bubble. We have studied the effects of the first and second- order diagrams on the spectral functions for a large set of parameters and for different transport regimes (reso- nant and off-resonant cases) at equilibrium and in the presence of a finite applied bias driving the system out of equilibrium. We have shown the important role played by the Hartree diagram in calculations based only on first-order diagrams. Such a diagram should not be neglected un- less one works within an applied bias range for which the correspond Hartree potential is constant. For calculations including both first and second-order diagram we have found that the effects of the individ- ual second-order diagram are as follows: the DX dia- gram reduces the width of the peaks in the spectral func- tions, with an increase of their height, while the DPH has the opposite effect: it increases the width and decreases the height. Furthermore, the DX diagram moves the vi- bron side-band peak position towards the main peak, and hence gives a better peak separation as should be ob- tained from exact calculations. Calculations performed with both DX and DPH diagram give intermediate re- sults for the spectral functions, which are not simply an average/superposition of the individual effects because of For the resonant regime with electron-hole symmetry, it is necessary to perform the calculations by includ- ing both DX and DPH diagrams in order to conserve the Fermi-liquid properties of the system. We anticipate that this will be also true for higher-order diagrams and calculations need to be performed by including all dia- grams of the same order. Calculations performed with partial resummation of a subset of diagrams will break the expected Fermi-liquid properties of the system and will probably lead to an incorrect value of the linear con- ductance. Finally, we expect that the effects of the second-order IV. CONCLUSION