In this paper, we introduce a new class of Schrödinger–Poisson type systems characterized by multi-scale logarithmic nonlinearities and nonlocal interactions. The model is governed by the coupled system \( \left\{ \begin{array}{ll} -\Delta u + V(x)u + \phi u = \lambda u + g(x,u,\phi ), & x \in \mathbb {R}^3, \\ -\Delta \phi = u^2, & x \in \mathbb {R}^3, \end{array}\right. \) where \(u:\mathbb {R}^3 \rightarrow \mathbb {R}\) represents the wave function, \(\phi :\mathbb {R}^3 \rightarrow \mathbb {R}\) denotes the electrostatic potential generated by the particle density, \(V:\mathbb {R}^3 \rightarrow \mathbb {R}\) is a continuous external potential satisfying \(V(x)\ge V_0>0\) , and \(\lambda \in \mathbb {R}\) is a real parameter. The nonlinear term \(g:\mathbb {R}^3 \times \mathbb {R}\times \mathbb {R}\rightarrow \mathbb {R}\) is defined by \( g(x,u,\phi ) = u \log (1+|u|^\alpha ) + \beta (x)\,u \log \log (e+|u|^\gamma ) + \eta \,\phi u \log (1+|u|^2), \) where \(\alpha ,\gamma >0\) , \(\beta :\mathbb {R}^3 \rightarrow \mathbb {R}\) is a bounded function, and \(\eta >0\) is a parameter describing the strength of the nonlocal coupling. This formulation introduces a new interaction mechanism between the wave function and the self-consistent electrostatic field, significantly extending classical Schrödinger–Poisson and Schrödinger–Bopp–Podolsky models studied in the literature. In contrast to recent works that mainly consider power-type or single logarithmic nonlinearities, the present framework incorporates double-logarithmic growth together with nonlocal logarithmic coupling, leading to a richer and more intricate variational structure. From a mathematical perspective, the presence of iterated logarithmic terms requires the development of refined analytical tools beyond standard Sobolev settings. By combining variational methods, generalized growth inequalities, and concentration–compactness techniques, we establish the existence of nontrivial weak solutions and prove the existence of infinitely many geometrically distinct solutions via genus theory. Moreover, we analyze the concentration behavior of solutions, showing that they localize near the global minima of the potential V(x). The results obtained in this work significantly extend existing theories for Schrödinger–Poisson type systems by capturing multi-scale nonlinear effects and strong nonlocal interactions. Finally, numerical simulations presented in tabular form confirm the theoretical findings, illustrating convergence properties, multiplicity of energy levels, and the influence of the parameters \(\eta \) and \(\gamma \) on the qualitative behavior of solutions.
