Interior Calder\'on-Zygmund estimates for solutions to general parabolic equations of p-Laplacian type

We study general parabolic equations of the form $u_t = div A(x,t, u,D u) + div(|F|^{p-2} F)+ f$ whose principal part depends on the solution itself. The vector field $A$ is assumed to have small mean oscillation in $x$, measurable in $t$, Lipschitz continuous in $u$, and its growth in $Du$ is like the $p$-Laplace operator. We establish interior Calder\'on-Zygmund estimates for locally bounded weak solutions to the equations when $p>2n/(n+2)$. This is achieved by employing a perturbation method together with developing a two-parameter technique and a new compactness argument. We also make crucial use of the intrinsic geometry method by DiBenedetto \cite{D2} and the maximal function free approach by Acerbi and Mingione \cite{AM}.