For a Riemannian manifold $M$ with metric $g$ and Laplace-Beltrami operator $-\Delta\_{g}$, what conditions on $M$ guarantee that $-\Delta\_{g} u(x)$ measures the difference between $u(x)$ and the average of $u$ over a geodesic ball (or sphere) centered at $x$? More precisely, what conditions on $M$ guarantee that


$$-\Delta\_{g} u(x) ~ \propto ~ \lim\_{h \to 0}{ \frac{2}{h^2} \left( u(x) - \frac{1}{|B(x,h)|} \int\_{B(x,h)}{ u(y) dy } \right) } \ \ ? $$


Here, $B(x,h)$ is the geodesic ball with center $x$, radius $h$, and measure $|B(x,h)|$. My guess is that this holds on harmonic manifolds (where harmonic functions can be characterized by the mean value property), but I haven't found the result anywhere. 



Laplace-Beltrami and averaging
This holds in complete generality. I sketch the case in which, instead of the ball $B(x,h)$, you consider the metric sphere $S(x,h)$. In particular (here $\dim M = n$ and the Laplacian is $\Delta = \mathrm{div}\circ \mathrm{grad}$):


$$ (\Delta u)(x) = \lim\_{h\to 0} \frac{2n}{h^2}\frac{1}{|S(x,h)|}\int\_{S(x,h)} [u(y)-u(x)] dy .$$


The case of the ball is similar with a different constant. Choose normal coordinates in a neighbhorhood of $x$, such that


$$u(y)-u(x) = \sum\_{i=1}^k (\partial\_i u)(x) y\_i + \frac{1}{2} \sum\_{i,j=1}^n (\partial\_{ij}^2 u)(x) y\_i y\_j + O(|y|^3). $$


The Riemannian metric in these coordinates is $g\_{ij}= \delta\_{ij}+O(|y|^2)$, hence up to higher orders the measure is the Euclidean one. Also the metric sphere, in these coordinates, for small $h$, coincides with the Euclidean sphere $\mathbb{S}^{n-1}$.


When you take the average over the sphere of each term, the linear part averages to zero (as the integral of a linear function over the sphere). The integral of the quadratic part averages to the sum of second derivatives. To see this, use the fact that for any symmetric matrix $M$


$$\int\_{\mathbb{S}^{n-1}} x^\* Q x \,d\mu\_{\mathbb{S}^{n-1}} = \frac{|\mathbb{S}^{n-1}|}{n} \mathrm{Tr}(Q).$$


Hence, at your point $x$, in normal coordinates, you have


$$ \lim\_{h\to 0} \frac{2n}{h^2}\frac{1}{|S(x,h)|}\int\_{S(x,h)} [u(y)-u(x)] dy = \sum\_{i=1}^n (\partial\_i^2 u)(x).$$


At the point $x$ this coincides with $(\Delta u)(x)$ written in normal coordinates centered in $x$. Since your averaging expression is coordinate invariant, the equality holds in general.


