Let $\Omega =\mathbb{R}^d \setminus \overline{B(0,1)}$. Suppose $u \in C^2(\Omega;\mathbb{R})$ ) is a harmonic function. Prove that, for every $x \in \mathbb{R}^d$ , the function $\phi: (|x| + 1, \infty) \to\mathbb{R}$ given by,
$$\phi(r) = \frac{1}{|\partial B(x,r)|} \int_{\partial B(x,r)} u(y) \mathrm{d}S(y)$$
is linear.
I was told to mimic the proof of the Mean Value Property. But I am stuck.
You may have a difficulty proving this statement because it's not true. Consider $u(x) = |x|^{2-d}$ (assuming $d\ge 3$); the averages of $u$ over concentric spheres are strictly positive and tend to zero as $r\to\infty$. It's not linear.
Let's see what is true. Since the gradient $\nabla u$ has divergence zero, its flux out of any domain in $\Omega$ is zero. Consider the spherical shell $A=\{y: r<|y-x|<R\}$; then $$ \int_{\partial A} \nabla u\cdot dS = 0 $$ which can be written as $$ R^{d-1}\phi'(R) - r^{d-1}\phi'(r) = 0 $$ the factors like $r^{d-1}$ coming from the surface area. So the conclusion is that $r^{d-1}\phi'(r)$ is constant, hence $$ \phi(r) = A + Br^{2-d} \qquad (A+B\log r \ \text{ if }d=2) $$ for some constants $A,B$.