How to get the idea of the formula for the mean value property for the heat equation
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.