This paper is devoted to the study of the existence, qualitative properties, and numerical approximation of normalized solutions for a class of critical Schrödinger–Bopp–Podolsky systems in \(\mathbb {R}^3\) , involving the p-Laplacian operator, logarithmic nonlinearities, critical Sobolev growth, and nonlocal electromagnetic interactions. More precisely, we consider the system \( {\left\{ \begin{array}{ll} -\varDelta _p u + V(\varepsilon x)|u|^{p-2}u + M(x)|u|^{p-2}u + \kappa \phi u \\ \qquad = \lambda |u|^{p-2}u + \vartheta |u|^{p-2}u \log |u|^p + \mu |u|^{q-2}u + |u|^{p^*-2}u, \\ -\varDelta \phi + a^2 \varDelta ^2 \phi = 4\pi u^2, \end{array}\right. } \qquad x\in \mathbb {R}^3, \) subject to the mass constraint \( \int _{\mathbb {R}^3} |u|^p\,dx = \rho ^p, \) where \(\varDelta _p u = \textrm{div}(|\nabla u|^{p-2}\nabla u)\) denotes the p-Laplace operator, \(p\in (1,3)\) , \(p^*=\frac{3p}{3-p}\) is the critical Sobolev exponent, \(\varepsilon >0\) is a small semiclassical parameter, and \(\lambda \in \mathbb {R}\) is a Lagrange multiplier associated with the normalization. The parameters \(\rho , a, \kappa , \vartheta , \mu > 0\) are given constants, \(V:\mathbb {R}^3\rightarrow \mathbb {R}\) is an external potential, and \(M:\mathbb {R}^3\rightarrow \mathbb {R}\) is a variable mass coefficient. The exponent q satisfies \(p<q<p+\frac{p^2}{3}\) . The main novelty of this work lies in the introduction and analysis of a Schrödinger–Bopp–Podolsky system combining variable mass effects and mixed nonlinearities, including logarithmic, subcritical, and critical growth terms, within a normalized framework. This setting significantly extends existing results by incorporating spatial inhomogeneities and competing nonlinear interactions, which fundamentally modify the variational structure of the problem. In particular, the presence of the variable mass term M(x) destroys translation invariance and affects the concentration behavior of solutions, while the logarithmic nonlinearity introduces non-polynomial growth that requires refined functional tools. To overcome these difficulties, we develop a variational approach on a normalized constraint manifold, employ Orlicz space techniques to handle the logarithmic term, use a reduction method for the nonlocal Bopp–Podolsky field, and apply concentration–compactness principles adapted to critical growth. The simultaneous presence of the p-Laplacian, logarithmic nonlinearity, critical exponent, and nonlocal coupling necessitates a nonstandard combination of analytical techniques beyond classical frameworks. We prove the existence and multiplicity of normalized solutions for sufficiently small \(\varepsilon >0\) and show that these solutions concentrate near the global minima of the external potential V. Furthermore, we analyze the asymptotic behavior as the Lagrange multiplier \(\lambda \rightarrow 0\) and as the Bopp–Podolsky parameter \(a\rightarrow 0\) , recovering the corresponding Schrödinger–Poisson system in the limit. Finally, we propose a numerical approximation scheme that preserves the normalization constraint. The numerical results, presented in tabular form, confirm the theoretical findings and illustrate the dependence of solutions on the model parameters, providing a unified analytical and computational framework for the proposed problem.
