Summary: We study a nonlinear Schrödinger-Poisson system which reduces to the nonlinear and nonlocal PDE \[ - \Delta u+ u + \lambda^2 \left( \frac{1}{\omega |x|^{N-2}}\star \rho u^2\right) \rho (x) u = |u|^{q-1} u \quad x \in\mathbb{R}^N, \] where \(\omega = (N-2)|\mathbb{S}^{N-1} |\), \(\lambda >0\), \(q\in (1,2^* -1)\), \(\rho :\mathbb{R}^N \rightarrow\mathbb{R}\) is nonnegative, locally bounded, and possibly non-radial, \(N=3,4,5\) and \(2^*=2N/(N-2)\) is the critical Sobolev exponent. In our setting \(\rho\) is allowed as particular scenarios, to either (1) vanish on a region and be finite at infinity, or (2) be large at infinity. We find least energy solutions in both cases, studying the vanishing case by means of a priori integral bounds on the Palais-Smale sequences and highlighting the role of certain positive universal constants for these bounds to hold. Within the Ljusternik-Schnirelman theory we show the existence of infinitely many distinct pairs of high energy solutions, having a min-max characterisation given by means of the Krasnoselskii genus. Our results cover a range of cases where major loss of compactness phenomena may occur, due to the possible unboundedness of the Palais-Smale sequences, and to the action of the group of translations.