  We prove the existence of ground state solution to the following problem.
\begin{align*}
(-\Delta)^{s}u+u&=\lambda|u|^{-\gamma-1}u+P(x)|u|^{p-1}u,~\text{in}~\mathbb{R}^N\setminus\Omega\\
N_su(x)&=0,~\text{in}~\Omega \end{align*} where $N\geq2$, $\lambda>0$,
$0<s,\gamma<1$, $p\in(1,2_s^*-1)$ with $2_s^*=\frac{2N}{N-2s}$. %
$0<s^-=\underset{{(x,y)\in\Omega\times\Omega}}{\inf}\{s(x,y)\}\leq s(x,y)\leq
s^+=\underset{{(x,y)\in\Omega\times\Omega}}{\sup}\{s(x,y)\}<1$,
$0<\gamma^-=\underset{{x\in\Omega}}{\inf}\{\gamma(x)\}\leq \gamma(x)\leq
\gamma^+=\underset{{x\in\Omega}}{\sup}\{\gamma(x)\}<1$,
$1-\gamma^-<1<p^-=\underset{x\in\Omega}{\inf}\{p(x)\}\leq p(x)\leq
p^+=\underset{x\in\Omega}{\sup}\{p(x)\}<2_{s^-}^*=\underset{{x\in\Omega}}{\inf}\{2_s^*(x)\}$
with $2_s^*(x)=\frac{2N}{N-2\tilde{s}(s)}$ where $\tilde{s}(x)=s(x,x)$.
  Moreover, $\Omega\subset\mathbb{R}^N$ is a smooth bounded domain,
$(-\Delta)^s$ denotes the $s$-fractional Laplacian and finally $N_s$ denotes
the nonlocal operator that describes the Neumann boundary condition which is
given as follows. \begin{align*}
N_{s}u(x)&=C_{N,s}\int_{\mathbb{R}^N\setminus\Omega}\frac{u(x)-u(y)}{|x-y|^{N+2s}}dy,~x\in\Omega.
\end{align*} We further establish the existence of infinitely many bounded
solutions to the problem.
