Mean value property for harmonic functions Consider a bounded harmonic function $u:\mathbb{R}^p \to \mathbb{R}$ (i.e. $u$ is a $C^2$ function such that the Laplacian $\Delta u=0$).
Prove, without using Liouville's theorem, the following version of the mean value property:
$$\forall x \in \mathbb{R}^p,\; u(x)=\frac{1}{2^p}\int\limits_{[-1,1]^p}u(y+x)dy$$
How can we prove it?
 A: I believe no proof of this statement will be shorter than the following classical proof of a strong form of Liouville's theorem. (This proof is well known, but should be even better known.) Given a bounded function, adding a constant yields a positive function, so the next statement implies that bounded harmonic functions on ${\mathbb R}^d$ are constant.   Let $B_R(x)$ be the ball of radius
$R$ around $x$ in ${\mathbb R}^d$, and let $|B_R|=R^d |B_1|$ be its volume.
Claim: If $u:{\mathbb R}^d \to [0,\infty)$ is harmonic, then it is constant.
Proof:  Given $x,y \in {\mathbb R}^d$ with $|x-y|=\delta$, we have
$$u(x)=\frac{1}{|B_R|} \int_{B_R(x)} u \,dz \le
 \frac{1}{|B_R|} \int_{B_{R+\delta}(y)} u \,dz =  \frac{|B_{R+\delta}|}{|B_R|} u(y) \,.$$
Taking  $R \to \infty$ yields $u(x) \le u(y)$. The same argument also gives $u(y) \le u(x)$.
A: I don't think there's a way to avoid Liouville completely. At most we can disguise it's proof to give a direct argument. This is what I do.
By translation invariance, we only need to prove this for $x=0$.
Call $C_r u$ the average on cubes, i.e.
$$
C_ru:= (2r)^{-p} \int_{[-r,r]^p} u(y)\, dy.
$$
If $|u(x)|\leq M$, then $|\nabla u(x)|\leq c(p)M$ (this follows from the usual estimates for harmonic functions). With this and the mean value theorem you can bound
$$
|u(0)-C_ru|\leq c(p)Mr.
$$
Set $H_M$ to be the set of harmonic functions bounded by $M$. Then the above says that for any $\varepsilon>0$ and $u\in H_M$ we have
$$
|u(0)-C_{\varepsilon/c(p)M} u|\leq \varepsilon.
$$
Notice that the scaling transformation $u_r(x)=u(x/r)$ leaves the set $H_M$ invariant for any $r>0$, and so we get, for any $u\in H_M$,
$$
|u(0)-C_1u|\leq \varepsilon.
$$
Since $\varepsilon$ was arbitrary the result follows.
This works with $C_r$ replaced by any nice enough averaging operator. The point is that the translation and scale invariance of $H_M$ is way too good to be true; only constant functions will be there (e.g. a common proof of Liouville exploits this at the level of the gradient estimate $|\nabla u(x)|\leq c(p)M$, where instead we get $|\nabla u(x)|\leq c(p)M/r$ for any $r>0$). So even though we "avoided" Liouville, the basic principle behind it is very much the key to getting the formula.
