  The paper is concerned with higher order Calderon-Zygmund estimates for the
$p$-Laplace equation $$
  -\textrm{div}(A(\nabla u))
  := -\textrm{div}{(|\nabla
  u|^{p-2}\nabla u)}=-\textrm{div} F, \qquad 1<p<\infty. $$ We are able to
transfer local interior Besov and Triebel-Lizorkin regularity up to first order
derivatives from the force term $F$ to the flux $A(\nabla u)$. For $p\geq 2$ we
show that $F \in B^s_{\rho,q}$ implies $A(\nabla u) \in B^s_{\rho,q}$ for any
$s \in (0,1)$ and all reasonable $\rho,q \in (0,\infty]$ in the planar case.
The result fails for $p<2$. In case of higher dimensions and systems we have a
smallness restriction on $s$. The quasi-Banach case $0<\min\{\rho,q\} < 1$ is
included, since it has important applications in the adaptive finite element
analysis. As an intermediate step we prove new linear decay estimates for
$p$-harmonic functions in the plane for the full range $1<p<\infty$.
