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.