  The purpose of this paper is to study and classify singular solutions of the
Poisson problem $$ %\begin{equation}\label{eq 0.1} \left \{ \begin{aligned}
{\mathcal L}^s_\mu u = f \quad\ {\rm in}\ \, \Omega\setminus \{0\},\\ u =0
\quad\ {\rm in}\ \, {\mathbb R}^N \setminus \Omega\ %\\ %\liminf_{x \to
0}\:|u(x)| /\Phi_\mu(x) = k. \end{aligned} \right. $$ for the fractional Hardy
operator ${\mathcal L}_\mu^s u= (-\Delta)^s u +\frac{\mu}{|x|^{2s}}u$ in a
bounded domain $\Omega \subset {\mathbb R}^N$ ($N \ge 2$) containing the
origin. Here $(-\Delta)^s$, $s\in(0,1)$, is the fractional Laplacian of order
$2s$, and $\mu \ge \mu_0$, where $\mu_0 =
-2^{2s}\frac{\Gamma^2(\frac{N+2s}4)}{\Gamma^2(\frac{N-2s}{4})}<0$ is the best
constant in the fractional Hardy inequality. The analysis requires a thorough
study of fundamental solutions and associated distributional identities.
Special attention will be given to the critical case $\mu= \mu_0$ which
requires more subtle estimates than the case $\mu>\mu_0$.
