Propagation of regularity for the heat equation Let $\Omega$ be an open bounded subset of $\mathbb R^N$.
Let $u_0 \in L^\infty(\Omega)$ and $f \in L^\infty((0,T)\times\Omega)$ and consider the following boundary value problem for the heat equation:
$$
\begin{cases} u_t - \Delta u = f \\
u|_{\partial\Omega} = 0\\
u(0) = u_0
\end{cases}
$$
Fix $x_0 \in \Omega$, $U$ a neighborhood of $x_0 \in \Omega$, and assume that $u_0 \in C^k(U)$, and $f \in C^k([0,T) \times U)$. Is it true that there exists a unique (weak) solution of the heat equation that is regular in $V$, that is $u \in C^k([0,T)\times V) \cap L^\infty$, where $V$ is a neighborhood of $x_0$ which is contained in $U$?
Also, do we have regularity up to the boundary of $\Omega$ if we assume $u_0 = \Delta u_0 = 0$?

I'd also like to draw your attention to a more general question that has appeared on MathOverflow.

Note that I'm talking about the behavior of $u(t, \cdot)$ for $t$ in the time interval $[0,T]$, not just in $[\epsilon, T]$. 
That's why I'm talking about the regularity of the initial data. 
 A: Yes, it is true. Consider the two auxiliary problems
$$\left\{
\begin{eqnarray}
&v_t - \Delta v &= 0 \\
&v(0) &= u_0
\end{eqnarray}
\right.$$
and
$$\left\{
\begin{eqnarray}
&w_t - \Delta w &= f \\
&w(0) &= 0
\end{eqnarray}
\right.$$
where we are looking for $v$ and $w$ in $L^\infty _0 (\Omega)$.
The first problem is just the heat equation on $\Omega$: this means that if $h : (0, \infty) \times \Omega \times \Omega \to (0, \infty)$ is the heat kernel of $\Omega$ (which is known to be in $C^\infty (\Omega) \cap L^p (\Omega)$ for all $p>1$), then it is known that
$$v(t,x) = \int _\Omega h(t,x,y) \ u_0 (y) \ \mathrm d y \ .$$
(This mimicks the convolution operator, because there is no convolution on $\Omega$ since $\Omega$ is not a group under the addition of vectors). This operation is "smoothing", i.e. regardless of how ugly $u_0$ is, $v$ will always be smooth (i.e. $C^\infty$).
The second problem is non-homogeneous, so we attack it with the usual tool: Duhamel's principle. This means that we consider yet another problem, namely
$$\left\{
\begin{eqnarray}
&W^s _t - \Delta W^s &= 0 \\
&W^s (0,s) &= f(s,x)
\end{eqnarray}
\right.$$
which again has the (smooth) solution $W^s (t,x) = \int _\Omega h(t,x,y) \ f(s,y) \ \mathrm d y$. From this we obtain (with Fubini's theorem at the end)
$$w(t,x) = \int _0 ^t W^s (t,x) \ \mathrm d s = \int _0 ^t \int _\Omega h(t,x,y) \ f(s,y) \ \mathrm d y \ \mathrm d s = \int _\Omega h(t, x, y) \int _0 ^t f(s,y) \ \mathrm d s \ \mathrm d y$$
which is seen with the naked eye to be smooth.
Finally, the sum $u = v + w$ clearly is a solution to the original problem, and is obviously smooth. This means that you get much more than just $C^k$ regularity, and this is because $h$ is $C^\infty$ and this is all that matters. It is also an elementary application of Hölder's inequality to check that $u$ is bounded (because $v$, $W^s$ and $w$ are and $\Omega$ has finite measure).
To show that $u(t, \cdot) = 0$ on $\partial \Omega$ it is enough to notice that $h(t,x,y) = 0$ as soon as either $x \in \partial \Omega$, or $y \in \partial \Omega$ (again, this is from the general theory concerning the heat kernel). This implies that $v(t, \cdot) = 0$ and $W^s (t, \cdot) = 0$ on $\partial \Omega$, whence $w(t, \cdot) = 0$ on $\partial \Omega$, hence $u(t, \cdot) = 0$ on $\partial \Omega$.
Smoothness in $t=0$ is again easy; I shall prove it for the first derivative, and then a simple induction does the rest:
$$\partial_t u(t,x) = \partial_t (v(t,x) + w(t,x)) = \partial_t \left( \int _\Omega h(t,x,y) \ u_0 (y) \ \mathrm d y + \int _\Omega h(t, x, y) \int _0 ^t f(s,y) \ \mathrm d s \ \mathrm d y \right) = \\
\int _\Omega \partial_t \ h(t,x,y) \ u_0 (y) \ \mathrm d y + \int _\Omega \partial_t \left( h(t, x, y) \int _0 ^t f(s,y) \ \mathrm d s \right) \mathrm d y = \\
\int _\Omega \Delta h(t,x,y) \ u_0 (y) \ \mathrm d y + \int _\Omega \left( \partial_t \ h(t, x, y) \int _0 ^t f(s,y) \ \mathrm d s + h(t, x, y) \ \partial_t \int _0 ^t f(s,y) \ \mathrm d s \right) \mathrm d y = \\
\int _\Omega \Delta h(t,x,y) \ u_0 (y) \ \mathrm d y + \int _\Omega \left( \Delta h(t, x, y) \int _0 ^t f(s,y) \ \mathrm d s + h(t, x, y) \ f(t,y) \right) \mathrm d y$$
which is easily seen to be continuous in $t=0$ (use Lebesgue's dominated convergence theorem and the standard fact that $\lim _{t \to 0} h(t,x,y) = \delta_x$ - Dirac's distribution). (Also, notice that it doesn't matter if you apply $\Delta$ in $x$ or in $y$ - things remain smooth anyway, because $h$ is so.)
For uniqueness, you want to show that the problem
$$\left\{
\begin{eqnarray}
&u_t - \Delta v &= 0 \\
&u(0) &= 0
\end{eqnarray}
\right.$$
has a unique solution in $L^\infty _0 (\Omega)$. But this, again, is standard knowledge; for instance, you could show it with the approach that uses the energy functional. Or you could use the weak maximum principle.
All of the above holds in the context that I am familiar with, namely that $\partial \Omega$ be smooth. I believe, though, that the same general theory holds when $\partial \Omega$ is only Lipschitz, it's just that the proofs become more technical.
