Mean value property implies harmonicity It is fairly easy to show that harmonic functions satisfy the mean value property, but it seems harder to show the converse. I've seen the following theorem without proof:

If $u \in C(\Omega)$ satisfies $$u(z) = \frac{1}{|\partial
 B_r(z)|}\int_{\partial B_r(z)} u\,dS$$ for all $z \in \Omega$ and
  $B_r(z) \subset \Omega$, then $u \in C^\infty$ and $u$ is harmonic on
  $\Omega$.

When I try to prove it myself, I got stuck. Could anyone kindly show me how to prove this? Thanks.
 A: Theorem: Let $\Omega \subset \mathbb{R}^N$, $u \in C(\Omega)$ be such that  $$\frac{1}{|B(x_0,R)|}\int_{B(x_0,R)}u(y)\ dy = u(x_0) = \frac{1}{|\partial B(x_0,R)|}\int_{\partial B(x_0,R)}u\ dS$$ for every ball $\overline{B(x_0,R)} \subset \Omega$. Then $u \in C^{\infty}(\Omega)$ and it is harmonic
Proof: Consider the standard mollifier: $$\rho(x) := \begin{cases}Ce^{-\frac{1}{1 - \|x\|^2}} & \text{if $\|x\|$ < 1} \\0 & \text{otherwise.}
\end{cases}$$ Here $C$ is a constant such that $\|\rho\|_{L^1} = 1.$ Let $\epsilon > 0$ and consider $$\rho_{\epsilon}(x) = \epsilon^{-N}\rho(x\epsilon^{-N}).$$
Set $\Omega_{\epsilon} = \{x \in \Omega : \text{dist}(x,\partial \Omega) > \epsilon\}$ and define for $x \in \Omega_{\epsilon}$$$u_{\epsilon}(x) = \rho_{\epsilon} * u(x) = \int_{\Omega}\rho_{\epsilon}(x - y)u(y)\ dy.$$ The following is a well know theorem in analysis, if it is new to you you can look for a proof Analysis by Lieb and Loss or anywhere else.
**Theorem:***If $u \in C(\Omega)$, then $u_{\epsilon} \to u$ uniformly on compact subsets of $\Omega$, $u_{\epsilon} \in C^{\infty}(\Omega_{\epsilon})$ and for any multindex $\alpha$ we have $$\frac{\partial^{\alpha}u_{\epsilon}}{\partial x^{\alpha}}(x) = \int_{\Omega}\frac{\partial^{\alpha}\rho_{\epsilon}}{\partial x^{\alpha}}(x - y)u(y)\ dy.$$*
Finally we can proceed with the proof!
Fix $x_0 \in \Omega_{\epsilon}$.
$$u_{\epsilon}(x_0) = \int_{B(x_0,\epsilon)}\rho_{\epsilon}(x - y)u(y)\ dy = \int_{B(0,\epsilon)}\rho_{\epsilon}(z)u(x_0 - z)\ dz = $$ $$ = \int_0^{\epsilon}r^{N - 1}\int_{\partial B(0,1)}\rho_{\epsilon}(rw)u(x_0 - rw)\ dS(w)dr = $$ $$ \int_0^{\epsilon}r^{N - 1}\rho(r)\int_{\partial B(0,1)}u(x_0 - rw)\ dS(w)dr = \int_0^{\epsilon}r^{N-1}\rho_{\epsilon}(r)\frac{\alpha_N N}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u(y)\ dS(y)dr $$ $$ = u(x_0)\|\rho\|_{L^1} = u(x_0).$$
This proves that $u = u_{\epsilon}$ and hence $u \in C^{\infty}(\Omega_{\epsilon})$, for every $\epsilon$. We are left to prove that $u$ is harmonic. To this end consider $$f(r) = \frac{1}{|\partial B(x_0,r)|}\int_{\partial B(x_0,r)}u\ dS.$$
By assumption, $f$ is constant, hence $f' \equiv 0$. By the divergence theorem it is easy to show (and I am sure that you have already seen this result since you have proved that if $u$ is harmonic then it satisfies the mean value property) $$0 = f'(r) = \text{constant}\int_{B(x_0,r)}\Delta u(y)\ dy \to \Delta u(x_0),\ \text{if we let}\ r \to 0^+.$$ Thus $u$ is harmonic. QED
