We show the capabilities of the operatorial approach, based on the equations of motion and the Green’s function formalisms [1-2], to very efficiently handle strongly correlated electron systems, reaching a deep and rich understanding of their unconventional properties. In particular, we propose a solution of the 2D Hubbard model in the framework of the Composite Operator Method (COM) within a novel four-pole (4p) approximation: the COM is highly tunable and expressly devised for the exploration of emergent phases in strongly correlated systems. First we consider the 2-site system, the smallest cluster where all Hamiltonian terms result active (the minimal model), characterize its exact solution and discuss how such an analysis effectively performed on this minimal model can suggest very efficient and controlled approximation schemes for the corresponding bulk system, which remains the final target of the investigation. We identify the contributions of the operators exactly embodying charge, spin, pair and double occupancy degrees of freedom to the relevant physical quantities of the system and clarify the crucial role played by spin fluctuations. Correspondingly, according to the given general prescription, we devise a novel 4p approximation for the 2D Hubbard model with a basis of fields given by two Hubbard operators plus two fields describing electronic transitions dressed by nearest-neighbor spin fluctuations. As well as this 4p approximate solution is in remarkable agreement with the exact one on the 2-site minimal model, the corresponding one proposed for the 2D system performs very well once compared to advanced (semi-)numerical methods, being by far less computational-resource demanding and more accurate in frequency and momentum resolutions. Moreover, this solution, by treating the spin fluctuations with extreme care and by improving the momentum selectivity of the spectral properties, opens up the possibility to directly address the underdoped-cuprate puzzle on a level before possible only including explicitly a residual self-energy in the calculations . Therefore, motivated by this long-standing experimental challenge posed by the single-particle properties of the underdoped cuprates (Fermi arcs, pseudogap, non-Fermi liquid behavior, extreme momentum-dependence of spectral properties, …), we adopt this 4p approximation to study the single-particle properties of the 2D Hubbard model in the strong coupling regime, where the effects of the spin fluctuations, accurately treated in our approach, are more relevant and induce unconventional features in all analyzed spectral properties that can be put in connection with those of the underdoped cuprates.
 F. Mancini, and A. Avella, Adv. Phys. 53, 537 (2004).
 A. Avella, Eur. Phys. J. B, 87, 45 (2014).
 A. Avella, F. Mancini, Phys. Rev. B 75, 134518 (2007).