Consider the inhomogeneous linear heat equation
$$\partial_tu-\Delta u=F$$
on $\mathbb R^n\times [0,1]$ (say) with zero initial data. Assume $F$ is very nice (say Schwarz), so that we have a nice solution $u$. It is quite standard that
$$\|\nabla^2 u\|_{L_x^2L_t^2}\ll \|F\|_{L_x^2L_t^2}. $$
I'm now asking if there is a similar estimate
$$\|\nabla^2 u\|_{L_x^2L_t^p}\ll_p \|F\|_{L_x^2L_t^p} $$
for $1<p<\infty$.
Remark: It is false for $p=\infty$, because the $L^2\to L^2$ norm of second derivatives of the heat kernel $K_t$ (i.e. the $L^\infty$ norm of its spatial Fourier transform) blows up like $1/t$ as $t\to 0$, which is not integrable in $t$. Also by duality it is also false for $p=1$, so I'm interested in what happens in between.