Suppose $f \; : \; \mathbb{R}^n \to \mathbb{R}$ is an $\mathbb{L}^{1}_{loc}$ function (defined everywhere!) $\mathcal{R}\subset(0,+\infty)$ is a set such that (*) holds, meaning : $$f(x) = avg\int_{B(x,r)}{f(y)dy} \; \forall r \in \mathcal{R} \; , \; \forall x \in \mathbb{R}^n (*)$$
I would like to study the property of functions satisfying $(*)$ for a general $\mathcal{R}$, in particular when does $(*)$ imply that $f$ is an harmonic function?
I showed that it does it when $\mathcal{R}$ accumulates to $0$.
Question : Is it true that $\Delta f = 0$ when $\mathcal{R} = \{1\}$ and when $\mathcal{R} = [1,\infty)$ ? If it's true prove it, if not find a counterexample.
You may assume, if you wish, that $f \in C^\infty(\mathbb{R}^n)$
Now I show you what I managed to show so far
First of all if (*) hold for $f$ then it holds also for $f^{\epsilon} := f * \eta_\epsilon$ since $$\int_{B(x,r)}{f^\epsilon(y) dy} = \int_{B(x,r)}{ \int_{B(0,\epsilon)}{ f(y - z)\eta_{\epsilon}(z) dz} dy} = \int_{B(0,\epsilon)}{ \int_{B(x,r)}{ f(y-z)\eta_{\epsilon}(z) dy} dz} = \int_{B(0,\epsilon)}{ |B(x,r)| f(x-z)\eta_{\epsilon}(z) dz} = |B(x,r)| f^\epsilon(x)$$
So it should be fair to assume $f \in C^\infty$
Now assume $f \in C^{1}$, I have
$$f(x) = avg\int_{B(0,r)}{f(x+y)dy}$$
Therefore
$$\nabla f(x) = avg\int_{B(0,r)}{\nabla f(x+y) dy} = nr\cdot avg\int_{\partial B(x,r)}{f(y) n_e dS(y)}$$
So from this I know that $f \in C^{k} \implies \nabla f\in C^{k}$, so $f \in C^{0} \implies f \in C^{\infty}$ (for this case you may also use the previous result to show $f \in C^{0} \implies f \in C^{1}$ )
Now I show that $(*) \implies \Delta f = 0$ when $\mathcal{R}$ accumulates to $0$ (and when $f$ is $C^3$ )
Let $x \in \mathbb{R}^n$ fixed but arbitrary, let $g(y) := f(x+y) - f(x) - \nabla f(x) \cdot y $ clearly $g$ satisfies (*) (and $g(0) = 0$ which will be very useful ), so it's enough to show $\Delta g(0) = 0$
I have
$$g(y) = \frac{1}{2} \nabla^2 g(0) \cdot y^2 + O(|y|^3) = \frac{1}{2}\sum_{i,j = 1}^{n}{ \frac{\partial^2 g(0)}{\partial y_i \partial y_j} y_iy_j } + O(|y|^3)$$
Let $H := \nabla^2 g(0)$, $H$ is a matrix ( a $2$-form) and $tr(H) = \Delta g(0)$.
Now I integrate over $B(0,r)$ (with $r \in \mathcal{R}$ )
And I find
$$0 = \int_{B(0,r)}{ H(y,y) dy} + \int_{B(0,r)}{ O(|y|^3) dy}$$ The second term will be negligible in the end so I can ignore it, the first term is invariant under trasformations $y \leftarrow Ry$ with $R \in O(n)$, applying these transformations properly (and multiplying both sides by $n$ ) I find
$$0 = tr(H)\int_{B(0,r)}{ |y|^2 dy} + \int_{ B(0,r)}{O(|y|^3) dy}$$
finally I divide bot sides by $\int_{B(0,r)}{|y|^2 dy}$ and let $r \to 0$ (which I can do since $\mathcal{R}$ accumulates to $0$) to get $$0 = tr(H)$$