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.
This theorem should answer the question above.
Theorem: Suppose $n > 2$ and $p_1, \dots, p_m$ are distinct points in $\mathbf{R}^n$. Suppose $u \colon \mathbf{R}^n \setminus \{p_1, \dots, p_m\}\to \mathbf{R}$ is harmonic. Then $u$ is bounded below if and only if there exist $a_0, a_1, \dots, a_m \in \mathbf{R}$ such that $a_j \ge 0$ for $j = 1, \dots, m$ and $$ u(x) = a_0 + a_1 |x-p_1|^{2-n} + \cdots + a_m |x - p_m|^{2-n} $$ for all $x \in \mathbf{R}^n \setminus \{p_1, \dots, p_m\}$.
Proof: Use induction on $m$. For $m = 1$, the desired result follows from the theorem quoted in the question.
Now suppose $m > 1$ and the desired result holds for $m-1$. By Theorem 3.9 in Harmonic Function Theory (second edition), there is a constant $a_m \ge 0$ and a harmonic function $v$ defined on some neighborhood of $p_m$ such that $$ u(x) = v(x) + a_m |x - p_m|^{2-n} $$ for all $x \ne p_m$ in a neighborhood of $p_m$.
Extend $v$ to be defined on $\mathbf{R}^n \setminus \{p_1, \dots, p_{m-1}\}$ by defining $$ v(x) = u(x) - a_m |x - p_m|^{2-n}. $$ Then $v$ is harmonic on $\mathbf{R}^n \setminus \{p_1, \dots, p_{m-1}\}$ (because $|x-p_m|^{2-n}$ is harmonic). By our induction hypothesis, $v$ has the desired form, which can then be plugged into the equation above, showing that $u$ has the desired form. QED