Local boundedness of weak solutions of heat equations...? (Please, what is going on?) The following claim is from a paper which was apparently reviewed by László Erdös, Zhongwei Shen, and Bernard Heffler. Someone tell me it's true. Surely it's true. The entire paper depends on it being true.
I am looking at a 2000 article by Kazuhiro Kurata in the J. London Math. Soc.: An Estimate on the Heat Kernel of Magnetic Schrödinger Operators and Uniformly Elliptic Operators with Nonnegative Potentials. The claim is Lemma 1 on page 893. (In the preprint, available online here, the claim appears more casually, inside a proof on the top of page 12.)
Kurata is studying the heat kernels of various Schrödinger operators, the simplest example being $H = -\Delta + V(x)$. His key condition is that $V(x)$ belongs to a reverse Hölder class $RH_{q}$ with $q > n/2$. 
Now suppose $u(x,t)$ is a weak solution of $(\partial_{t} + H)u = 0$ in some cylinder 
$$Q_{2r}(x_{0}, t_{0}) = B(x_{0}, 2r) \times (t_{0}-(2r)^{2}, t_{0})$$ 
According to Kurata, under just these hypotheses, there is a ``standard subsolution estimate'' 
$$\sup_{(x,t) \in Q_{r/2}(x_{0}, t_{0})} |u| \leq \biggl(\frac{C}{r^{n+2}}
\iint_{Q_{2r/3}(x_{0},t_{0})}|u|^2\,dx\,dt\biggr)^{1/2}$$
with $C$ independent of $V$, $t_{0}$, and $r$. He cites Aronson and Serrin's foundational '67 paper ``Local behavior of solutions of quasilinear parabolic equations'', but gives no specifics. 
In every local boundedness estimate I've seen for parabolic equations (including, incidentally, Theorem 2 in Aronson and Serrin '67), the constants seem to depend on $V$, $t_{0}$, and $r$...!?! Is there just something magical about the reverse Hölder hypothesis?
Life, I guess.
 A: I think by "standard subsolution estimate", Kurata means the following:
Let $u\ge0$ be a subsolution in $Q_2:=B_2(0)\times(-4,0)$, i.e. $\partial_t u - \nabla\cdot(A\nabla u)\le 0$. Then $u$ satisfies
$$sup_{Q_{1/2}} u \le C\Big( \int_{Q_1} |u|^2 \;dxdt \Big)^{\frac{1}{2}}$$
for some constant depending only on the ellipticity contrast $\lambda$ and the dimension $d$.
Then the result of Kurata with a constant independent of $r$ follows by scaling space with $r$ and time with $r^2$. Likewise the constant may be chosen independently of $t_0$ and $x_0$ by translation invariance of the problem.
I am guessing that Kurata assumes $u\ge 0$ and $V\ge 0$, so that a solution to his problem is indeed a subsolution: $\partial_t u - \nabla\cdot(A\nabla u)\le 0$.
The above "standard subsolution estimate" may be obtained by parabolic Moser iteration, cf.
J. Moser, A harnack inequality for parabolic differential equations, Comm. Pure Appl. Math., link.
Moser iteration corresponds to testing the equation with clever test functions. The simplest example would be the following: consider an elliptic problem (e.g. $u$ above does not depend on time) $-\nabla\cdot A \nabla u \le 0$ in $B_2$ and test with $\eta^2 u$, where $\eta$ is a test function which is one in $B_{1/2}$ and zero outside $B_1$. Then after integration by parts
$$0\ge-\int \eta^2 u \nabla\cdot A \nabla u \;dx = \int \Big(\eta^2 \nabla u \cdot A \nabla u + 2\eta u\nabla\eta \cdot A \nabla u \Big) \;dx.$$
Ellipticity and Young's inequality yield
$$\int \eta^2 |\nabla u|^2 \;dx\le C(\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$
This is usually called Caccioppoli's inequality. By the product rule, it can be upgraded to
$$\int |\nabla (\eta u)|^2 \;dx\le C(\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$
Now we use the Sobolev embedding $\|v\|_{L^p(R^d)}\le C(d) \|\nabla v\|_{L^2(R^d)}$ with $p=2d/(d-2)$ (let's assume $d>2$) to obtain
$$\Big(\int |\eta|^p |u|^p \;dx\Big)^{\frac{2}{p}}\le C(d,\lambda) \int |\nabla\eta|^2 u^2 \;dx.$$
Finally, by the choice of test function, e.g. $|\nabla \eta| \le C(d)$ and we obtain
$$\Big(\int_{B_{1/2}} |u|^p \;dx\Big)^{\frac{2}{p}}\le C(d,\lambda) \int_{B_1} u^2 \;dx.$$
This was the first step to give you some gain of integrability (giving up a little bit on the domain of integration on the right hand side). For the next step, you test with something like $\eta^2u^{p-1}$...
I hid a lot of technical details and the argument gets more complicated when you iterate higher and if you do parabolic estimates, but the idea is there. Anyway, I've rambled on much longer than I intended to...
PS: The reason why the constant in Aronson and Serrin seems to depend on $V$ and $r$ is that they consider a (much) more general problem which in particular is not scaling invariant in the same way that $\partial_t u - \nabla \cdot A \nabla u \le 0$ is.
