Assume that $ (M,g) $ is a smooth closed manifold and $ \mathbb{S}^{L-1} $ is a unit sphere with dimension $ L-1 $ in $ \mathbb{R}^{L} $. Consider the equation of harmonic heat flow $$ \partial_tu-\Delta_g u=|\nabla u|^2u, $$ where $ u=(u^1,u^2,...,u^{L}) $. For the initial problem $$ \left\{\begin{matrix} \partial_tu-\Delta_g u=|\nabla u|^2u&\text{in}&M\times(0,T),\\ u(\cdot,0)=u_0&\text{on}&M. \end{matrix}\right.\quad(*) $$ A well-known result is the if $ u_0\in C^{\infty}(M) $ and $ T $ is sufficiently small, then the above initial problem has a unique smooth solution. Such results can be find in the book "The analysis of harmonic maps and their heat flows" by F. Lin and C. Wang, Page 111. The methods in that book is by using contract map theorem. I want to go another way. Firstly, standards arguments of Galerkin methods implies that $ (*) $ has weak solution for any $ T>0 $, i.e. $ \partial_tu\in L^2(0,T;L^2(M)) $, $ \nabla u\in L^{\infty}(0,T;H^1(M)) $, $ $ for any $ \phi\in C_c^{\infty}(M\times(0,T)) $, we have $$ \int_0^T\int_{M}\partial_tu\phi \mathrm{d}\operatorname{Vol}_g+\int_{0}^T\int_M\langle\nabla u,\nabla\phi\rangle_g\mathrm{d}\operatorname{Vol}_g=\int_0^T\int_{M}|\nabla u|^2u\phi \mathrm{d}\operatorname{Vol}_g. $$ I want to enhance the regularity of this weak solution for sufficiently small $ T $ but I do not know how to start. Can you give me some hints or references?