Mean value property of harmonic functions on manifolds A well-known feature of harmonic functions on (domains of) $\mathbb{R}^n$ is the mean-value property: that is, if $\Delta u = 0$, then
$$ u(x_0) = \frac{1}{\text{Vol}(\partial B_r(x_0))}\int_{\partial B_r(x_0)}{u\,dS} = \frac{1}{\text{Vol}(B_r(x_0))}\int_{B_r(x_0)}{u\,dV}. $$
Is the same true on manifolds in general? That is, given a manifold $(M,g)$ and a smooth function $u:M\rightarrow\mathbb{R}$ satisfying
$$\Delta_gu = 0\quad\text{where}\quad\Delta_g = \frac{1}{\sqrt{\det g}}\frac{\partial}{\partial x^i}g^{ij}\sqrt{\det g}\frac{\partial}{\partial x^j},$$
is it true that
$$ u(x_0) = \frac{1}{\text{Vol}(\partial B_r(x_0))}\int_{\partial B_r(x_0)}{u\,dS}\quad\text{or}\quad u(x_0) = \frac{1}{\text{Vol}(B_r(x_0))}\int_{B_r(x_0)}{u\,dV}?$$
Here, $B_r(x_0)$ is the geodesic ball of radius $r$ around $x_0$. If the equalities do not hold, how would the average value (either over a sphere or over a ball) change with $r$? For example, on a surface the Gaussian curvature appears in the higher-order terms in the Taylor expansion of the length/area of a circle/ball of radius $r$--does a similar phenomenon occur for the mean value?
 A: After doing a quick Google search I came across this blog post: https://cuhkmath.wordpress.com/2015/08/14/mean-value-theorems-for-harmonic-functions-on-riemannian-manifolds/
The main points of the article are as follows:


*

*The mean value property of harmonic functions holds on an arbitrary manifold $M$ only when for every point $p\in M$ every geodesic sphere near $p$ has constant mean curvature. Such manifolds are called harmonic manifolds. This coincides with the intuition I had in my comments above.

*For any smooth function $u$ on an $n$-manifold $(M,g)$ it holds that $$\frac{1}{\mathrm{vol}(\partial B_r(p)}\int_{\partial B_r(p)}u(y)~d\sigma(y)=u(p)+\frac{\Delta u(p)}{2n}r^2 + \mathscr{O}(r^4).$$
The $\mathscr{O}(r^4)$ term is given by $B(n)r^4 + \mathscr{O}(r^6)$, where $$B(n):=\frac{1}{24(n+2)}\left(3\Delta^2u - 2\langle\nabla^2 u,\rho\rangle-3\langle\nabla u,\nabla\tau\rangle + \frac{4}{n}\tau\Delta u\right).$$
In the above $\rho(x,y)=\mathrm{tr}\left(R(\cdot,x,y,\cdot)\right)$, $R$ is the Riemann curvature tensor, and $\tau=\mathrm{tr}(\rho)$.

*We have a similar mean value property for geodesic balls. For any smooth function $u$ it holds that $$\frac{1}{\mathrm{vol}(B_r(p))}\int_{B_r(p)}u(y)~d\mu(y)=u(p) + \frac{\Delta u(p)}{2(n+2)}r^2 + B(n+2)r^4 + \mathscr{O}(r^6).$$

*You have similar sub-mean value properties for subharmonic functions given bounds on the Ricci and Riemann curvatures.
The blog provides plenty of references and proofs (they seem to be correct, but I haven't had the time to work through them yet).
