  The aim of this paper is to give a new proof that any very weak $s$-harmonic
function $u$ in the unit ball $B$ is smooth. As a first step, we improve the
local summability properties of $u$. Then, we exploit a suitable version of the
difference quotient method tailored to get rid of the singularity of the
integral kernel and gain Sobolev regularity and local linear estimates of the
$H^{s}_{\rm loc}$ norm of $u$. Finally, by applying more standard methods, such
as elliptic regularity and Schauder estimates, we reach real analyticity of
$u$. Up to the authors' knowledge, the difference quotient techniques are new.
