  In this article, we discuss the non-trivial collective charge excitations
(plasmons) of the extended square-lattice Hubbard model. Using a fully
non-perturbative approach, we employ the hybrid Monte Carlo algorithm to
simulate the system at half-filling. A modified Backus-Gilbert method is
introduced to obtain the spectral functions via numerical analytic
continuation. We directly compute the single-particle density of states which
demonstrates the formation of Hubbard bands in the strongly-correlated phase.
The momentum-resolved charge susceptibility is also computed on the basis of
the Euclidean charge density-density correlator. In agreement with previous
EDMFT studies, we find that at large strength of the electron-electron
interaction, the plasmon dispersion develops two branches.
