$L^\infty-L^2$ smoothing for heat equation on manifold using Nash-Moser-De-Giorgi technique Let $M$ be a compact and closed smooth Riemannian manifold, and consider weak solution $u$ of the equation
$$u_t - \Delta u = f$$
given $f \in L^2(Q)$ and $u(0)=u_0 \in L^\infty(M)$.
I'm looking for a reference or sketch proof using Nash-Moser-De-Giorgi iterations (and not heat kernels or semigroup stuff, since I wish to adapt it to a different setting) of an $L^\infty-L^2$ smoothing effect for $u$, i.e., 
$$\lVert u \rVert_{L^\infty(Q)} \leq C(\lVert f \rVert_{L^2(Q)}, \lVert u_0 \rVert_{L^\infty(M)})$$
where $Q:=[0,T]\times M$.
Note that $f$ is only in $L^2$. 
I tried books by Davies etc. but for the proofs involving bounded  domains it seems the proofs really depend on the Sobolev inequalities, which varies on manifolds.
 A: Here's a rough calculation that might give you what you want (but the constant depends on the Sobolev constant!): The basic iterative step in Moser iteration for this particular equation is:
\begin{align*}
\frac{1}{p}\frac{d}{dt}\int |u|^p &\le \int |u|^{p-2}u\Delta u + |u|^{p-1}f\\
&\le -\int |\nabla |u|^{p/2}|^2 + |u|^{p-1}f\\
&\le -C^{-1}\left(\int |u|^{pn/(n-2)}\right)^{(n-2)/n}\\
& + \left(\int |u|^{pn/n-2)}\right)^{(n-2)(p-1)/(pn)}\left(\int |f|^{(n+2p-2)/(np)}\right)^{np/(n+2p-2)}\\
&\le (2C)^{p-1}\left(\int |f|^{(n+2p-2)/(np)}\right)^{n/(n+2p-2)}.
\end{align*}
This implies an inequality for $\|u\|_p$, which you can take the limit $p\rightarrow 0$. It appears to me that you don't need to iterate. 
Actually, there's a good chance the Sobolev constant factor converges to $1$ in the limit $p\rightarrow\infty$, so the final $L^\infty$ inequality does not depend on the Sobolev constant. In fact, I used that observation once.
A: To do De Giorgi-Nash-Moser it is important that $f$ is in $L^q$ for $q$ larger than half the dimension. A heuristic is scaling: if $v$ solves $-\Delta v = f$, then the right side for the rescaling $v(\epsilon x)$ is $\epsilon^2 f(\epsilon x)$, whose $L^q$ norm is like $\epsilon^{2-n/q}$, so zooming in helps when $q > n/2$.
When in dimension $4$ or higher, one can take e.g. $v = \sum h_k\varphi(2^k x)$ for some smooth $\varphi$ that is $1$ in $B_1$ and vanishes outside $B_2$. Taking $h_k = 1/k$ makes $v$ unbounded ($\log\log$ growth near the origin) while $f := -\Delta v$ is $L^2$ (the integral of $f^2$ is like $\sum h_k^2 2^{(4-n)k}$).
It seems to me that for such a choice of $f$, the solution $u$ will immediately become unbounded since $u - v$ is unbounded initially and solves the heat equation. 
Note also that indeed, in Deane's computation, the (reciprocal) exponent on $f$ looks like $n/2$ for large $p$. It is important to use that $f \in L^q$ for $q > n/2$, for instance as follows (in the elliptic case for simplicity): Let $A_p = \left(\int u^{p\chi}\right)^{\frac{1}{p \chi}}$, with $\chi = n/(n-2)$. Applying Holder to the last term in the first line of Deane's computation gives an inequality like
$$A_p \leq (pC(f,S))^{1/p} A_{\gamma p},$$
where $\gamma < 1$ depends on $q/(q-1) < \chi$, and $C(f,S)$ depends on $f$ and the Sobolev constant. Iteration gives
$$A_{2\gamma^{-k}} \leq C(f,S)^{\sum j\gamma^j} \|u\|_{L^2(Q)}.$$
Taking $k \rightarrow \infty$ gives an $L^{\infty}$ bound in terms of the desired quantities and the Sobolev constant. I'm not sure how to remove dependence on the Sobolev constant. 
