I want to prove that the mean value property $$u(\textbf{x}_0) = \frac{1}{\pi r^2} \int \int _{\left \{ \left | x_0-x \right < r| \right \}} u(\textbf{x})d\textbf{x}$$ for non-negative harmonic function $u$ so I can eventually prove Harnack inequality.
I know how to prove the mean value property for a regular harmonic function using the Poisson kernel.
I'm confused how a proof for the the regular haramonic function would differ from proving the mean value property for non-negative harmonic function.
Here is a simple proof of the mean-value theorem for harmonic functions of two variables. We start with Green's Third Identity
$$u(\vec x_0)=\oint_C \left(u(\vec x')\frac{\partial G(\vec x_0,\vec x')}{\partial n'}-G(\vec x_0,\vec x')\frac{\partial u(\vec x')}{\partial n'}\right)\,d \ell' \tag 1$$
where $u(\vec x)$ is harmonic in a region $S$ bounded by the smooth contour $C$ and $G(\vec x_0,\vec x')-\frac{1}{2\pi}\log|\vec x_0-\vec x'|$ is the Green's Function for the Laplacian.
If we choose $S$ to be the region $|\vec x_0-\vec x'|\le r$, and thus $C$ to be the circle of radius $r$ and center $\vec x_0$, then on $C$ we have
$$\begin{align} G(\vec x_0,\vec x')&=\frac{1}{2\pi}\log(r) \tag 2\\\\ \frac{\partial G(\vec x_0,\vec x')}{\partial n'}&=\frac{1}{2\pi r} \tag 3 \end{align}$$
Using $(2)$ and $(3)$ in $(1)$ yields
$$\begin{align} u(\vec x_0)&=\frac{1}{2\pi r}\int_0^{2\pi} u(\vec x')\,d\ell'-\frac{1}{2\pi}\log(r)\oint_C \frac{\partial u(\vec x')}{\partial n'}d \ell' \tag 4\\\\ &=\frac{1}{2\pi r}\int_0^{2\pi} u(\vec x')\,d\ell'-\frac{1}{2\pi}\log(r)\int_{S} \nabla '\cdot \nabla 'u(\vec x')\,dS' \tag 5\\\\ &=\frac{1}{2\pi r}\int_0^{2\pi} u(\vec x')\,d\ell'-\frac{1}{2\pi}\log(r)\int_{S} \nabla '^2 u(\vec x')\,dS' \\\\ &=\frac{1}{2\pi r}\int_0^{2\pi} u(\vec x')\,d\ell' \tag 6 \end{align}$$
In going from $(4)$ to $(5)$, we invoked the Divergence Theorem, while in arriving at $(6)$ we used the fact that $u$ is harmonic.
Finally, if we multiply both sides of $(6)$ by $2\pi r$ and integrate with respect to $r$ we find
$$\pi r^2\,u(\vec x_0)=\int_S u(\vec x')\,dS'$$
whereupon dividing by $\pi r^2$ yields the coveted mean-value theorem
$$u(\vec x_0)=\frac{1}{\pi r^2}\int_{|\vec x_0-\vec x'|\le r} u(\vec x')\,dS'$$