Spherical average over harmonic, decreasing function On page 529, remark 4.9, in Solovej's article on the ionization conjecture in Hartree-Fock theory (https://arxiv.org/abs/math-ph/0012026) the following assertion is made:
Let $V:\mathbb{R}^3\to\mathbb{R}$ be harmonic and continuous for $|x|>R>0$ and $\lim_{|x|\to\infty}V(x)=0$. Then the limit
$$\lim\limits_{r\to\infty}\int_{S^2}rV(r\omega)\mathrm{d}\omega$$
($S^2$ being the unit sphere in three dimensions) exists. In particular the integral does not even depend on $r$.
The last statement is of course trivial for the case where $V$ is spherical symmetric. However I cannot prove it for the general case.
My intuition is that every function that fulfills the above conditions must be (in some sense) comparable to $|x|^{-1}$ or any translates of that, i.e. $|x-y|$, where $y\in B_0(R)$.
Applying the maximum/minimum principle gives me that there are constants $c$ and $c'$, s.t. 
$$\frac{c'}{|x|}\leq V(x)\leq\frac{c}{|x|}.$$
 A: Here is how the proof should go. Set
$$\phi(r) = \int_{S^2} rV(r\omega) \, d\omega,$$
and differentiate in $r$ to obtain
$$\phi'(r) = \int_{S^2} V(r\omega) +r\nabla V(r\omega)\cdot \omega \, d\omega.$$
Make a change of variables to find
$$\phi'(r) = \int_{rS^2} \frac{1}{r^2} V(\omega) + \frac{1}{r}\frac{\partial V}{\partial n}(\omega) \, d\omega,$$
where $rS^2$ is the sphere of radius $r$. Setting $W(z) = |z|^{-1}$ we can rewrite this as
$$\phi'(r) = \int_{rS^2}  -V(\omega)\frac{\partial W}{\partial n}(\omega) + W(\omega)\frac{\partial V}{\partial n}(\omega) \, d\omega.$$
Now since $W$ and $V$ are harmonic, you can use Green's identity and the boundary condition $V\to 0$ as $|x| \to \infty$ to show that the right hand side is zero. You may have to fill in a few details on this last part, but I believe this is the basic idea of the proof.
A: Your proposed lower bound on $|V|$ is not correct. For example, $V(x) = x_1/|x|^3$ is harmonic on $\mathbb R^3\setminus \{0\}.$ I got that example by taking the Kelvin transform of $x_1.$ Many other examples can be obtained this way.
According to Thm 10.1 in http://www.axler.net/HFT.pdf, a harmonic function in an annular domain in $\mathbb R^3$ has an expansion
$$\sum_{m=0}^{\infty} P_m(x) + \sum_{m=0}^{\infty} \frac{Q_m(x)}{|x|^{2m + 1}},$$
where the $P_m, Q_m$ are homogenous harmonic polynomials of degree $m.$ The result you seek can be obtained from this, although there may be simpler ways.
A: Apparently it is enough to assume that the spherical average (let me forget about $4\pi$) over the unit sphere $F(0,r)=\int_{S^2}V(r\omega)d\omega$ of $V$ goes to zero as $r\to\infty$. Here, I defined $$F(x,r):=\int_{|\omega|=1}V(x+r\omega)d\omega.$$
Now, because $V$ is harmonic for $|x|>R$, we can use Darboux's equation which says that
$$\left(\partial_r^2+\frac{2}{r}\partial_r\right)F(x,r)=\Delta_xF(x,r)$$
for $x\in\mathbb{R^3}$, such that $|x+r\omega|>R$ for all $\omega\in\mathbb{R^3}$ with $|\omega|=1$. Because of elliptic regularity, one can pull $\Delta_x$ inside the spherical average, which makes the right hand side vanish. Hence, if we evaluate both sides at $x=0$, we get that
$$\partial_r^2(rF)(0,r)=0.$$ But this tells us that $F(0,r)$ must be of the form $F(0,r)=a+\frac{b}{r}$ for some $a,b\in\mathbb{R}$. But since $F(0,r)$ must vanish at infinity, we have that $a=0$. Therefore $rF(0,r)=b$ which is independent of $r$.
Let me make one remark to Jeff's solution. In the end, where you use Green's second identity, you get that
$$\phi'(r)=\int_{B_0(r)}\frac{1}{|x|}\Delta V(x)dx.$$ But the right hand side only vanishes for the range $r>R$, which is why the right hand side is effectively
$$\phi'(r)=\int_{B_0(R)}\frac{1}{|x|}\Delta V(x)dx,$$
which is not necessarily zero.
