Why is $\frac{1}{|x|^{n-2}}u(\frac{x}{|x|^2})$ harmonic if $u$ is harmonic? I found myself trying to prove the following, but I had to compute everything explicitly.
It is well known that if $u:\mathbb{R}^n\to\mathbb{R}$ is an harmonic function on $\mathbb{R}^n$, then the so-called Kelvin transform of $u$
$$(Ku)(x):=\frac{1}{|x|^{n-2}}u\biggl(\frac{x}{|x|^2}\biggr)$$
is harmonic in $\mathbb{R}^n\setminus\{0\}$.
Question: Is there a way to prove that $Ku$ is harmonic in $\mathbb{R}^n\setminus\{0\}$ WITHOUT making the explicit computations?
I know that there are tricks to soften the calculations (for example noticing that $1/|x|^{n-2}$ is harmonic in $\mathbb{R}^n\setminus\{0\}$), however I wasn't able to find a "clever" way to do it. I noticed that $x\mapsto x/|x|^2$ is a special transformation but I wasn't able to go anywhere, not even using the mean value property.
 A: An explanation can be the following. Take a harmonic function of the form $u(x)=r^\alpha P(\omega)$, with $r=|x|$, $|\omega|=1$ and $P$ a polynomial. Then $0=\Delta u=r^{\alpha-2}\left (\alpha(N-2+\alpha)P(\omega)+\Delta_S P(\omega)\right)$ with $\Delta_S$ the Laplace-Beltrami on the unit sphere $S$. Then $\Delta_S P=-\alpha(N-2+\alpha)P$. However we know that the eigenvalues of $\Delta_S$ are given by $-k(N-2+k)$ with $k$ a nonegative integer but the equation $\alpha(N-2+\alpha)=k(N-2+k)$ has another root $\alpha=2-N-k$ which corresponds to the Kelvin transform $Ku(x)=|x|^{2-N}u(x/|x|^2)$. This shows that $Ku$ is harmonic whenever $u$ is an harmonic polynomial and then (approximating locally) for any harmonic function $u$.
A: Here is another explanation which has various generalizations.
Denote by $L = \Delta + \frac{n-2}{4(n-1)}R$ the conformal Laplacian of a Riemannian manifold $(M^n,g)$ (my convention is that $\Delta$ is a nonnegative operator).  If $\Phi \colon M \to M$ is a conformal diffeomorphism — i.e. $\Phi^\ast g = e^{2\phi}g$ for some $\phi \in C^\infty(M)$ — then
$$ (Lu) \circ \Phi = e^{-\frac{n+2}{2}\phi} L\left( e^{\frac{n-2}{2}\phi} (u \circ \Phi) \right) . $$
Let's now specialize to $\mathbb{R}^n$ with $g$ the Euclidean metric. Then $\Phi(x) := x/\lvert x \rvert^2$ is a conformal diffeomorphism of $\mathbb{R}^n \setminus \{0\}$ with $\Phi^\ast g = \lvert x \rvert^{-4} g$. Therefore
$$ (\Delta u)\left( \frac{x}{\lvert x\rvert^2}\right) = \lvert x \rvert^{n+2} \Delta \left( \lvert x \rvert^{2-n} u \bigl( \frac{x}{\lvert x\rvert^2} \bigr) \right) $$
for any $u \in C^2(\mathbb{R}^n \setminus \{0\})$.
Specializing to the case when $u$ is harmonic recovers your result.
One generalization of this is the following:
There is a family of operators $L_{2\gamma} = \Delta^\gamma + \mathrm{l.o.t.}$, $\gamma \in (0 , n/2]$, on $(M^n,g)$ such that if $\Phi^\ast g = e^{2\phi}g$, then
$$ \tag{$\ast$} \label{eqn} (L_{2\gamma}u) \circ \Phi = e^{-\frac{n+2\gamma}{2}\phi} L_{2\gamma}\left( e^{\frac{n-2\gamma}{2}\phi} (u \circ \Phi) \right) . $$
On Euclidean space, $L_{2\gamma} = \Delta^\gamma$, so one concludes that if $u \in \ker \Delta^\gamma$, then $\lvert x\rvert^{2\gamma-n}u\bigl( x/\lvert x\rvert^2 \bigr) \in \ker \Delta^\gamma$.
While there are relatively simple proofs of \eqref{eqn}, they all require developing a fair bit of machinery besides doing a direct computation in Euclidean space by some other mechanism.
