Laplace-Beltrami and averaging 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. 
 A: 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.
