We consider the following nonlocal Brézis-Nirenberg type critical Choquard problem involving the Grushin operator \(\begin{aligned} \left\{ \begin{aligned} -\Delta _\gamma&u =\lambda u + \left( \displaystyle \int _\Omega \frac{|u(w)|^{2^{*}_{\gamma ,\mu }}}{d(z-w)^\mu }dw\right) |u|^{2^{*}_{\gamma ,\mu }-2}u \quad & \text {in} \ \Omega ,\\ u&= 0 \quad & \text {on} \, \partial \Omega , \end{aligned} \right. \end{aligned}\) where \(\Omega \) is an open bounded domain in \(\mathbb {R}^N\) , \(N \ge 3\) , with \(\Omega \cap \{ x=0\} \ne \emptyset \) , and \(\lambda >0\) is a parameter. Here, \(\Delta _\gamma \) represents the Grushin operator, defined as \( \Delta _\gamma u(z) = \Delta _x u(z) +(1+\gamma )^2 |x|^{2\gamma } \Delta _y u(z), \quad \gamma \ge 0, \) where \(z=(x,y)\in \Omega \subset \mathbb {R}^m\times \mathbb {R}^n\) , \(m+n=N \ge 3\) and \(2^{*}_{\gamma ,\mu }= \frac{2N_\gamma -\mu }{N_\gamma -2}\) is the Sobolev critical exponent in the Hardy-Littlewood context with \(N_\gamma = m+(1+\gamma )n\) is the homogeneous dimension associated to the Grushin operator and \(0<\mu <N_\gamma \) . The homogeneous norm related to the Grushin operator is denoted by \(d(\cdot )\) . In this article, we prove the existence of bifurcation from any eigenvalue \(\lambda ^*\) of \(-\Delta _\gamma \) under Dirichlet boundary conditions. Furthermore, we show that in a suitable left neighborhood of \(\lambda ^*\) , the number of nontrivial solutions to the problem is at least twice the multiplicity of \(\lambda ^*\) .
