  We consider the following fractional semilinear Neumann problem on a smooth
bounded domain $\Omega\subset\mathbb{R}^n$, $n\geq2$, $$\begin{cases}
(-\varepsilon\Delta)^{1/2}u+u=u^{p},&\hbox{in}~\Omega,\\ \partial_\nu
u=0,&\hbox{on}~\partial\Omega,\\ u>0,&\hbox{in}~\Omega, \end{cases}$$ where
$\varepsilon>0$ and $1<p<(n+1)/(n-1)$. This is the fractional version of the
semilinear Neumann problem studied by Lin--Ni--Takagi in the late 80's. The
problem arises by considering steady states of the Keller--Segel model with
nonlocal chemical concentration diffusion. Using the semigroup language for the
extension method and variational techniques, we prove existence of nonconstant
smooth solutions for small $\varepsilon$, which are obtained by minimizing a
suitable energy functional. In the case of large $\varepsilon$ we obtain
nonexistence of nonconstant solutions. It is also shown that as
$\varepsilon\to0$ the solutions $u_\varepsilon$ tend to zero in measure on
$\Omega$, while they form spikes in $\overline{\Omega}$. The regularity
estimates of the fractional Neumann Laplacian that we develop here are
essential for the analysis. The latter results are of independent interest.
