Let $u$ be harmonic in $\Omega$, and let $B_R(x_0) \subset \Omega$ be a ball of radius $R$ centered at $x_0$.

Claim: $u(x_0) = \frac{1}{|\partial B_R|}\int_{\partial B_R} u(y)\,dS_y$

Step 1: Define the Average Function

Define: $\phi(r) = \frac{1}{|\partial B_r|}\int_{\partial B_r(x_0)} u(y)\,dS_y$ for $0 < r \leq R$.

We will show that $\phi$ is constant, which implies $\phi(r) = \lim_{r \to 0^+} \phi(r) = u(x_0)$.

Step 2: Compute the Derivative

In spherical coordinates centered at $x_0$, using $dS = r^{n-1}\omega_{n-1}d\omega$ where $\omega$ represents angular variables: $\phi(r) = \frac{1}{\omega_{n-1}r^{n-1}}\int_{\partial B_r} u(x_0 + r\omega)\,r^{n-1}d\omega = \frac{1}{\omega_{n-1}}\int_{S^{n-1}} u(x_0 + r\omega)\,d\omega$

Taking the derivative: $\phi'(r) = \frac{1}{\omega_{n-1}}\int_{S^{n-1}} \nabla u(x_0 + r\omega) \cdot \omega\,d\omega$

Step 3: Apply Divergence Theorem

By the divergence theorem applied to the region between spheres of radius $\epsilon$ and $r$: $\int_{\partial B_r} \frac{\partial u}{\partial n}\,dS = \int_{B_r \setminus B_\epsilon} \nabla^2 u\,dx + \int_{\partial B_\epsilon} \frac{\partial u}{\partial n}\,dS$

Since $\nabla^2 u = 0$ in $B_r$: $\int_{\partial B_r} \frac{\partial u}{\partial n}\,dS = \int_{\partial B_\epsilon} \frac{\partial u}{\partial n}\,dS$

As $\epsilon \to 0$, the right side vanishes (assuming $u$ is $C^1$), so: $\int_{\partial B_r} \frac{\partial u}{\partial n}\,dS = 0$

Step 4: Relate to $\phi'(r)$

Note that $\frac{\partial u}{\partial n} = \nabla u \cdot \frac{x - x_0}{|x - x_0|} = \frac{1}{r}\nabla u \cdot \omega$ on $\partial B_r$.

Therefore: $0 = \int_{\partial B_r} \frac{\partial u}{\partial n}\,dS = \int_{\partial B_r} \frac{1}{r}\nabla u \cdot \omega\,r^{n-1}d\omega = r^{n-2}\int_{S^{n-1}} \nabla u \cdot \omega\,d\omega$

This gives: $\phi'(r) = \frac{1}{\omega_{n-1}}\int_{S^{n-1}} \nabla u \cdot \omega\,d\omega = 0$

Step 5: Conclude

Since $\phi'(r) = 0$ for all $0 < r \leq R$, the function $\phi(r)$ is constant. By continuity: $u(x_0) = \lim_{r \to 0^+} \phi(r) = \phi(R) = \frac{1}{|\partial B_R|}\int_{\partial B_R} u(y)\,dS_y$

This completes the proof of the spherical mean value property.