Let$\Omega \subset \mathbb{R}^{N} $, N ≽ 2, be a smooth bounded domain. For s ∈ (1/2, 1), we consider a problem of the form
$$\left\{\begin{array}{@{}ll} (-\Delta)^s u = \mu(x)\, \mathbb{D}_s^{2}(u) + \lambda f(x), & {\rm in}\,\Omega, \\ u= 0, & {\rm in}\,\mathbb{R}^{N} \setminus \Omega,\end{array}\right.$$
where λ > 0 is a real parameter, f belongs to a suitable Lebesgue space, $\mu \in L^{\infty}$ and $\mathbb {D}_s^2$ is a nonlocal ‘gradient square’ term given by
$$\mathbb{D}_s^2 (u) = \frac{a_{N,s}}{2} \int_{\mathbb{R}^{N}} \frac{|u(x)-u(y)|^2}{|x-y|^{N+2s}}\,{\rm d}y.$$
Depending on the real parameter λ > 0, we derive existence and non-existence results. The proof of our existence result relies on sharp Calderón–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.