Proving the mean value property of harmonic functions using distributions?

A professor I talked to showed me a proof of the mean value property. (He actually showed it for functions solving the heat equation instead of Laplace's equation, but it seems like the argument is the same.) The proof involves distributions, which I am not very familiar with, so there is a step that I am confused about.

Here is how the proof goes:

Let $\Phi(\vec x)$ be the fundamental solution of Laplace's equation.

For simplicity, suppose we are in two dimensions. Then $\Phi(\vec x) = - \frac{1}{2\pi} \log |\vec x|$, so $-\Delta \Phi(\vec x) = \delta_0(\vec x)$, where $\delta_0$ is the delta distribution.

Fix $k \in \mathbb R$ and let $B = \{ \vec x \in \mathbb R^n : \Phi(\vec x) > k \} = \{ \vec x \in \mathbb R^n : |\vec x| < e^{-2\pi k} \}$. Let $\chi_B(\vec x)$ be the characteristic function of $B$, and consider the function $\Theta(\vec x) = \chi_B(\vec x)(\Phi(\vec x) - k)$. Note in particular that $\Theta$ is continuous at the boundary $\partial B$.

If we take the Laplacian of $\Theta$, we will have a $\delta_0$ due to $\Phi$, but also, since the first derivative of $\Theta$ has a jump discontinuity at $\partial B$, there will be another distribution that only involves $\partial B$. Thus, $-\Delta\Theta(\vec x) = \delta_0(\vec x) + f(\vec x)$, where $f$ is a distribution that only involves $\partial B$. (Here is where I start to get confused.)

If $u$ is a harmonic function, then $(\Delta\Theta, u) = (\Theta, \Delta u) = (\Theta, 0) = 0$. But $-\Delta\Theta = \delta_0 + f$, so $(\delta_0, u) = - (f, u)$. Since $(\delta_0, u) = u(0)$, we need to show that $(f, u) = \frac{1}{2\pi e^{-2\pi k}} \int_{\partial_B} u$ to complete the proof.

My question is, how do we show that $f$ satisfies $(f, u) = \frac{1}{2\pi e^{-2\pi k}} \int_{\partial_B} u$? I know this comes from differentiating $\Theta$ at the boundary, but I am not that familiar with distributions. (I can see why in one dimension, the distributional derivative of the Heaviside function is the delta function, but I am not sure how to proceed in higher dimensions.)

-

By definition, the distributional Laplacian satisfies $\int (\Delta \Theta)\varphi =\int \Theta (\Delta \varphi)$ for any smooth compactly supported $\varphi$. This means we want to prove that $\int \Theta (\Delta \varphi)=-\varphi(0)+\frac{1}{2\pi R}\int_{\partial B}\varphi$ where $R$ is the radius of $B$.
(Update) Here is a streamlined proof for all dimensions, based on Green's second identity $$\int_{\partial \Omega}\left( u\frac{\partial v}{\partial n}-v\frac{\partial u}{\partial n}\right)=-\int_\Omega (u\Delta v-v\Delta u)$$ where $n$ is the inward normal and $u,v$ are smooth in the closure of domain ${\Omega}$, which itself has a smooth boundary.
Let $\Omega=\{x:\epsilon<|x|<R\}$, $u=\log |x|-\log R$, and $v=\varphi$ (smooth test function, as above). We get $$\int_{|x|=\epsilon} \left((\log \epsilon-\log R)\varphi_r -\epsilon^{-1}\varphi\right) + \int_{|x|=R} R^{-1}\varphi = -\int_\Omega (\log|x|-\log R)\,\Delta \varphi$$ As $\epsilon\to 0$, this becomes $$-c_n \varphi(0) + \frac{1}{R}\int_{|x|=R} \varphi = \int_{|x|<R} (\log|x|-\log R)\,\Delta \varphi$$ where $c_n$ is the area of unit sphere. QED
(Previous version) In 2 dimensions, write the Laplacian in polar coordinates $(r,t)$: $$\Delta \varphi=\varphi_{rr}+\frac{1}{r}\varphi_r+\frac{1}{r^2}\varphi_{tt}$$ and split the integral into two parts accordingly: radial and tangential derivatives. Integrate the tangential part with respect to $t$ first: it gives $0$. Integrate the radial part with respect to $r$ first, not forgetting to include the Jacobian $r$: $$\int_0^R (\log r-\log R) (r\varphi_{rr}+\varphi_{r})\,dr = \int_0^R (\log r-\log R) (r\varphi_{r})_r\,dr$$ and now by parts: $$\dots =(\log r-\log R) (r\varphi_{r}) \bigg|_{r=0+}^{r=R} -\int_0^R (\log r-\log R)_r (r\varphi_{r})\,dr \\= 0-0-\int_0^R \varphi_{r}\,dr \\ = -\varphi(Re^{it})+\varphi(0)$$ It remains to integrate over $t$ and recall the factor of $-1/(2\pi)$ in $\Theta$ that was ignored in the calculation so far. We get $$-\varphi(0)+\frac{1}{2\pi}\int_0^{2\pi}\varphi(Re^{it})\,dt$$ as desired.
Thank you! This was really easy to follow! One question I have is: how do we generalize this argument to higher dimensions? I know the radial component of $\Delta \phi$ will be $\phi_{rr}+\frac{n-1}{r}\phi_r$, but why is the integral of the remaining part zero? – Alan C Aug 9 '12 at 0:13
@AlanC Maybe this way... Laplacian is the divergence of $\nabla \phi$. By the divergence theorem, the integral of $\Delta\phi$ over a thin spherical shell $\{\rho<|x|<\rho+h\}$ is equal to the difference of integrals of $\phi_r$ over two spheres. In polar form this difference looks like $(\rho+h)^{n-1}\phi_r(\rho+h)-\rho^{n-1}\phi_r(\rho)$. Dividing by $h$ and letting $h\to 0$ we get $\rho^{n-1}\phi_{rr}+(n-1)\rho^{n-2}\phi_r$. So this is the integral of $\Delta \phi$ over the sphere $\{|x|=\rho\}$. – user31373 Aug 9 '12 at 2:22