How to get the idea of the formula for the mean value property for the heat equation From the mean-value property of the Laplace's equation, we have the following mean-value property:
$$
u(x)=\frac{1}{a(n)r^n}\int_{B(x,r)}u\,dy.
$$
But for the mean-value property of the Heat equation, Evans' book defines a heat ball:
$$
E(x,t,r)=\left\{(y,s)\in R^{n+1}\bigg|s\leq t, \Phi(x-y,t-s)\geq \frac{1}{r^n}\right\}.
$$
Then, the theorem claims that if $u\in C^2_1(U_T)$ solves the heat equation. Then,
$$
u(x,t)=\frac{1}{4r^n}\iint\limits_{E(x,t,r)}u(y,s)\frac{|x-y|^2}{(t-s)^2}\,dy\,ds.
$$
My question is: is there any explanation (or a guessed one) about the discover of this theorem? The mean-value property is intuitive. But how can we know that we can achieve the goal by making the integrand as the multiplication of $u(y,s)$ with such a strange factor, $\frac{|x-y|^2}{(t-s)^2}$, and a nonintuitive heat ball? 
 A: I highly recommend you skim through this original paper,

*

*Fulks, W., A mean value theorem for the heat equation, Proc. Am. Math. Soc. 17, 6-11 (1966). ZBL0152.10503.

The idea is to find an analogy of the Green's 2nd identity (with $\Delta$ replaced by the heat operator $H = \Delta - \partial_t$), so by choosing a suitable region of integration (not too surprisingly we will choose the heat ball $E$ as it is the symmetry inherent in the heat equation) and some suitable test functions we arrive at
\begin{equation}
u(x, t) = \int_{\partial E(x,t,1/c)} Q(x-y,t-s)u(y,s)\;d\mathcal{H}^1\end{equation}
where $Q(y,s)=cx^2[4x^2t^2 + (2t-x^2)^2]^{-1/2}$ (see the paper for the details)
This is not yet the volume integral we want. A remarkable observation is that the formula above does not depend on the choice of the constant $c$, so by the coarea formula we have
\begin{multline}
u(x,t) = \int^\infty_1\frac{1}{c^2}\left(\int_{\partial E(x,t,1/c)}  Q(x-y,t-s)u(y,s)\;d\mathcal{H}^1\right) dc \\= \int_{E(x,t,1)}u(y,s)\frac{\vert x - y\vert^2}{4(t-s)^2}\;dyds\end{multline}
The weight $1/c^2$ I append in this weighted average is to kill the additional $c^2$ in this calculation.
Personally I prefer this derivation than that presented by Evans, as it seems we may apply this method (Green's 2nd identity) to derive mean value formulas for a wide class of PDE.
