Existence and multiplicity of solutions for a critical Grushin problem with a singular nonlinearity

We investigate the existence and multiplicity of positive solutions to the problem \begin{equation} \begin{cases} \begin{aligned} - \Delta_{\gamma} u &= \lambda u^{p} + u^{-\delta} &\quad \text{in } \Omega, \quad u &= 0 &\quad \text{on } \partial \Omega, \end{aligned} \end{cases} \end{equation} where $\Delta_{\gamma}$ denotes the Grushin operator defined by \begin{equation} \Delta_{\gamma} := \Delta_x + (1+\gamma)^2 |x|^{2\gamma}\Delta_y, \end{equation} with $\gamma>0$, $z=(x,y)\in \mathbb{R}^N$, $N=n+m$, $n \geq 1$, $m\geq 1$, $\Omega \subset \mathbb{R}^N$ a smooth bounded domain, $\lambda>0$, $1<p<\infty$, and $\delta>0$. The analysis depends on the exponent $p$, which may be subcritical, critical, or supercritical, that is, $p<2_\gamma^*-1$, $p=2_\gamma^*-1$, or $p>2_\gamma^*-1$, respectively, where $2_\gamma^*=\frac{2Q}{Q-2}$ is the critical Sobolev exponent associated with the Grushin operator, and $Q=m+(1+\gamma)n$ is the corresponding homogeneous dimension.