Consider the heat equation on a (smooth) domain in $\mathbb{R}^n$ with homogeneous Neumann BCs:

$$u_t - \Delta u = f$$ $$\partial_\nu u = 0$$ $$u|_{t=0} = u_0$$ where $f \in L^p(0,T;L^r(\Omega))$ and $u_0 \in L^q(\Omega)$ can be smoother if required.

Under which conditions on $p,r,q$ and $n$ do we get that the weak solution is in $L^\infty((0,T)\times \Omega)$?

I can find something in the LSU book -- but I'm sure there have been improvements. What's the maximal regularity result for this?