Neumann Li-Yau gradient estimate under integral Ricci curvature bounds

We prove a Li-Yau gradient estimate for positive solutions to the heat equation, with Neumann boundary conditions, on a compact Riemannian submanifold with boundary ${\mathbf{M}}^n\subseteq {\mathbf{N}}^n$, satisfying the integral Ricci curvature assumption: \begin{equation} D^2 \sup_{x\in {\bf N}} \left( \oint_{B(x,D)} |Ric^-|^p dy \right)^{\frac{1}{p}} < K \end{equation} for $K(n,p)$ small enough, $p>n\mathrm{/}2$, where $diam(\mathbf{M})\le D$. The boundary of $\mathbf{M}$ is not necessarily convex, but it needs to satisfy the interior rolling $R-$ball condition.