Quasilinear problems with critical Sobolev exponent for the Grushin p-Laplace operator

We study the following class of quasilinear degenerate elliptic equations with critical nonlinearity \begin{align*} \begin{cases}-\Delta_{\gamma,p} u= \lambda |u|^{q-2}u+|u|^{p_{\gamma}^{*}-2}u & \text{ in } \Omega\subset \mathbb{R}^N, \\ u=0 & \text{ on } \partial \Omega, \end{cases} \end{align*} where $\Delta_{\gamma, p}v:=\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 \geq 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\}\neq \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 extends 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 \begin{equation*} \int_{\mathbb{R}^N} |\nabla_{\gamma} u|^p dz \geq S_{\gamma,p} \left ( \int_{\mathbb{R}^N} |u|^{p_\gamma^*} dz \right )^{p/p_\gamma^*} \end{equation*} and their qualitative behavior as positive entire solutions to the limit problem \begin{equation*} -\Delta_{\gamma,p} u= u^{p_{\gamma}^{*}-1}\quad \mbox{on}\, \mathbb{R}^N, \end{equation*} whose study has independent interest.