We study the following class of quasilinear degenerate elliptic equations with critical nonlinearity \(\begin{aligned} {\left\{ \begin{array}{ll}-\Delta _{\gamma ,p} u= \lambda |u|^{q-2}u+\left| u\right| ^{p_{\gamma }^{*}-2}u & \text { in } \Omega \subset \mathbb {R}^N,\\ u=0 & \text { on } \partial \Omega , \end{array}\right. } \end{aligned}\) where \(\Delta _{\gamma , p}u:=\sum _{i=1}^N X_i(|\nabla _\gamma u|^{p-2}X_i u)\) is the Grushin p-Laplace operator, \(z:=(x, y) \in \mathbb {R}^N\) , \(N=m+n,\) \(m,n \ge 1\) , where \(\nabla _\gamma =(X_1, \ldots , X_N)\) is the Grushin gradient, defined as the system of vector fields \(X_i=\frac{\partial }{\partial x_i}, i=1, \ldots , m\) , \(X_{m+j}=|x|^\gamma \frac{\partial }{\partial y_j}, j=1, \ldots , n\) , where \(\gamma >0\) . Here, \(\Omega \subset \mathbb {R}^{N}\) is a smooth bounded domain such that \(\Omega \cap \{x=0\}\ne \emptyset \) , \(\lambda >0\) , \(q \in [p,p_\gamma ^*)\) , where \(p_{\gamma }^{*}=\frac{pN_\gamma }{N_\gamma -p}\) and \(N_\gamma =m+(1+\gamma )n\) denotes the homogeneous dimension attached to the Grushin gradient. The results extend to the p-case the Brezis-Nirenberg type results in Alves-Gandal-Loiudice-Tyagi [J. Geom. Anal. 2024, 34(2),52]. The main crucial step is to preliminarily establish the existence of the extremals for the involved Sobolev-type inequality \( \int _{\mathbb {R}^N} |\nabla _{\gamma } u|^p dz \ge S_{\gamma ,p} \left( \int _{\mathbb {R}^N} |u|^{p_\gamma ^*} dz \right) ^{p/p_\gamma ^*} \) and their qualitative behavior as positive entire solutions to the limit problem \(\begin{aligned} -\Delta _{\gamma ,p} u= u^{p_{\gamma }^{*}-1}\quad \hbox {on}\, \mathbb {R}^N, \end{aligned}\) whose study has independent interest.
