  Let $\Omega \subset \mathbb{R}^N$, $N \geq 2$, be a smooth bounded domain.
For $s \in (1/2,1)$, we consider a problem of the form \[
\left\{\begin{aligned} (-\Delta)^s u & = \mu(x)\, \mathbb{D}_s^{2}(u) + \lambda
f(x)\,, & \quad \mbox{in} \Omega,\\ u & = 0\,, & \quad \mbox{in} \mathbb{R}^N
\setminus \Omega, \end{aligned} \right. \] where $\lambda > 0$ is a real
parameter, $f$ belongs to a suitable Lebesgue space, $\mu \in
L^{\infty}(\Omega)$ and $\mathbb{D}_s^2$ is a nonlocal "gradient square" term
given by \[ \mathbb{D}_s^2 (u) = \frac{a_{N,s}}{2}\mbox{p.v.}
\int_{\mathbb{R}^N} \frac{|u(x)-u(y)|^2}{|x-y|^{N+2s}} dy \,. \] Depending on
the real parameter $\lambda > 0$, we derive existence and non-existence
results. The proof of our existence result relies on sharp Calder\'on-Zygmund
type regularity results for the fractional Poisson equation with low
integrability data. We also obtain existence results for related problems
involving different nonlocal diffusion terms.
