We study the existence of nonnegative solutions to the following nonlocal elliptic problem involving singularity \(\begin{aligned} \mathfrak {M}\left( \int _{Q}\frac{|u(x)-u(y)|^p}{|x-y|^{N+ps}}dxdy\right) (-\Delta )_{p}^{s} u&=\frac{\lambda }{u^{\gamma }}+u^{p_s^*-1}~\text {in}~\Omega ,\\ u&>0~\text {in}~\Omega ,\\ u&=0~\text {in}~\mathbb {R}^N\setminus \Omega , \end{aligned}\) where \(\mathfrak {M}\) is the Kirchhoff function, \(Q=\mathbb {R}^{2N}\setminus ((\mathbb {R}^N\setminus \Omega )\times (\mathbb {R}^N\setminus \Omega ))\) , \(\Omega \subset \mathbb {R}^N\) , is a bounded domain with Lipschitz boundary, \(\lambda >0\) , \(N>ps\) , \(0<s,\gamma <1\) , \((-\Delta )_{p}^{s}\) is the fractional p-Laplacian for \(1<p<\infty \) and \(p_s^*=\frac{Np}{N-ps}\) is the critical Sobolev exponent. We employ a cut-off argument to obtain the existence of k (being arbitrarily large integer) solutions. Furthermore, by using the Moser iteration technique, we prove an uniform \(L^{\infty }({\Omega })\) bound for the solutions. The novelty of this work lies in proving the existence of small energy solutions by using symmetric mountain pass theorem in spite of the presence of a critical nonlinear term which, of course, is super-linear.
