Let $\Omega$ be a bounded Lipschitz domain (smoother if needed), and consider the heat equation $$u_t - \Delta u = 0$$ $$\frac{\partial u(t,x)}{\partial \nu(x)} = a(t,x) - b(t,x)u(t,x)$$ $$u(0) = u_0$$ where $u \in H^1(0,T;L^2) \cap L^2(0,T;H^1)$ is a weak solution for given $u_0 \in L^\infty(\Omega) \cap H^1(\Omega)$, and $a, b \in H^1(0,T;L^2(\partial\Omega)) \cap L^2(0,T;H^1(\partial\Omega)) \cap L^\infty(0,T;L^\infty(\partial\Omega)).$

Also, all data are non-negative hence as is the solution, which is also bounded.

A regularity result is that in fact $u \in C^0([0,T]\times \overline\Omega)$.

Is it true that $u \in C^{1,2}((0,T)\times \Omega)$?

That is, once differentiable in time and twice differentiable in space.

I would like this result in order to apply the strong maximum principle (and doing so would enable me to solve my previous question).

I expect it to be true due to parabolic smoothing but I haven't come across a proof or explanation for this. The closest I could find was this thread, where some commenters suggested that it should hold but I don't know what kind of BC they had in mind, and a little more detail would be useful.