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.
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
To please the PDE people out there, here is an alternative proof I will share for the sake of variety that uses that weakly harmonic functions are strongly harmonic. I first learnt this argument from a close friend from graduate school, Joel.
We are mathematicians, so new proofs are always welcome, but in this case the proof is even more remarkable because this argument works with even milder requirements over $u$, we require it only to satisfy the mean value property and to be in $L^1_{\text{loc}}(\Omega)$. With this in mind, it suffices to prove that: $$\int_{\Omega} \Delta \varphi \:u\:dx\quad \forall \varphi \in C_0^\infty(\Omega)$$
We need to make a mean over the sphere appear. To this end we notice that: $$\Delta \varphi(x)=\text{trace}(\text{Hess}(\varphi(x)))=\nabla_y\cdot (\text{Hess}(\varphi(x))(y))$$ Now, upon taking a mean over $y$ we obtain that: $$\Delta \varphi(x)=\frac{1}{|B_\varepsilon(0)|}\int_{B_\varepsilon(0)}\nabla_y\cdot (\text{Hess}(\varphi(x))(y))dy=\frac{1}{\varepsilon|B_\varepsilon(0)|}\int_{\partial_{\varepsilon} B(0)}\langle \text{Hess}(\varphi(x))(y),y\rangle dy$$ Now a finite difference trick (just taylor expand the terms over $x$) yields that: $$\varphi(x+y)+2\varphi(x)+\varphi(x-y)=$$ $$\langle \text{Hess}(\varphi(x))(y),y\rangle+\frac{1}{6}(D^3 \varphi(\theta x+(1-\theta)y)(y,y,y)+D^3 \varphi(\tilde{\theta} x+(1-\tilde{\theta})y)(y,y,y))$$ Where $0\leq \theta,\tilde{\theta}\leq 1$ possibly depends on $x$ and $y$ and we will denote by $$r_x(y)=\frac{1}{3}D^3 \varphi(\theta x+(1-\theta)y)(y,y,y)+D^3 \varphi(\tilde{\theta} x+(1-\tilde{\theta})y)(y,y,y))$$ . Substituting this above and making a change of variables yields that: $$\Delta \varphi(x)=\frac{1}{\varepsilon| B_\varepsilon(0)|}\int_{\partial B_\varepsilon(0)}\varphi(x+y)+2\varphi(x)+\varphi(x-y)+r_x(y) dy=$$ $$\frac{2}{\varepsilon| B_\varepsilon(0)|}\left(\int_{\partial B_\varepsilon(x)}\varphi(y)dy-\varphi(x)|\partial B_\varepsilon(0)|-\int_{\partial B_\varepsilon(0)}\frac{r(y)}{2} dy\right)$$ All in all, one obtains that: $$\int_{\Omega} \Delta \varphi \:u\:dx=I_1+I_2+I_3$$ Where for $\varepsilon$ small enough so that $\text{supp}(\varphi)+B_\varepsilon(0)\subseteq \Omega$: $$I_1=\int_{\Omega} \frac{-2}{\varepsilon| B_\varepsilon(0)|}\int_{B_\varepsilon(0)}\frac{r_x(y)}{2} dy \:u(x)\: dx\quad;\quad I_2=\int_{\Omega} \frac{-2|\partial B_\varepsilon(0)|}{\varepsilon| B_\varepsilon(0)|}\varphi(x) u(x) dx$$ $$I_3=\int_{\Omega} \frac{2}{\varepsilon| B_\varepsilon(0)|}\int_{\partial B_\varepsilon(0)}\varphi(x+y)dy\: u(x) dx=\frac{2}{\varepsilon| B_\varepsilon(0)|}\int_{\partial B_\varepsilon(0)}\int_{\Omega}\varphi(x)u(x-y) dxdy=$$ $$\frac{2}{\varepsilon| B_\varepsilon(0)|}\int_{\Omega}|\partial B_\varepsilon(0)|\varphi(x)u(x) dydx$$ Therefore, magically, $I_2+I_3=0$ so that one has: $$\left|\int_{\Omega} \Delta \varphi \:u\:dx\right|=\left|\int_{\text{supp}(\varphi)} \frac{-2}{\varepsilon| B_\varepsilon(0)|}\int_{\partial B_\varepsilon(0)}\frac{r_x(y)}{2} dy u(x) dx\right|\leq$$ $$ c_n\varepsilon\|u\chi_{\text{supp}(\varphi)}\|_{L^1}\frac{1}{\varepsilon^{n-1}}\int_{\partial B_\varepsilon(0)}\frac{\sup_{x\in \text{supp}(\varphi)+B_\varepsilon(0)}\|D^3 \varphi(x)\|\|y\|^3}{\varepsilon^3}d\sigma(y)=$$ $$ c_n\varepsilon\|u\chi_{\text{supp}(\varphi)}\|_{L^1}\sup_{x\in \text{supp}(\varphi)+B_\varepsilon(0)}\|D^3 \varphi(x)\|\xrightarrow{\varepsilon\rightarrow 0^+}0$$ whence we conclude that $u$ is weakly harmonic as we wished to prove.