Recently, I learn a theorem of Bocher from Axler-Bourdon-Ramey's book (Theorem 3.9 and 3.14) on harmonic functions which states that
Theorem
A positive harmonic functions on $\mathbb{R}^n\setminus\{ 0\}$ must be the form $a_0 + a_1 |x|^{2-n}$, where $a_i$ are constants.
My question is that : Can the above theorem be generalized to the multi-singularities case?
Conjecture
A positive harmonic functions on $\mathbb{R}^n\setminus\{ 0,1\}$ must be the form $a_0 + a_1 |x|^{2-n} + a_2 |x-1|^{2-n}$, where $a_i$ are constants.
My apptempt is to prove that such harmonic function can be decomposed to two functions, one is harmonic on $\mathbb{R}^n\setminus\{ 0\}$, the other one is $\mathbb{R}^n\setminus\{ 1\}$. But I have no idea to prove this.
I also see the authors provide the following exercise
Exercise 3.17
Let $P=\{ p_1,p_2,\cdots\}$ denote a discrete set of $\mathbb{R}^n$. Characterize the positive harmonic functions on $\mathbb{R}^n \setminus P$.
I guess its answer is a moditication of my conjecture. But I have no idea to prove it too.
Many thanks for any discussion.