Kelvin inversion transforms harmonic functions

Put $y:=x-a$, $r:=\mathrm{\mid }y\mathrm{\mid }$ and $z:=I_R(x)=a+R^{2}y\mathrm{/}{r}^2$. On the open set $V:=I_R^{-1}(U)$ we have $x\mathrm{\ne }a$, $r>0$ and $z\in U$, and $x\mapsto z$ is smooth there, being built from the smooth coordinate functions $y_i$ and $r^2$ and the smooth factor $R^2\mathrm{/}{r}^2$; hence $K_{R}u=R^{n-2}{r}^{2-n}(u\circ z)$ is $C^2$ on $V$, and no value is taken at $r=0$.

Coordinate differentiation of $z$ gives, for all $1\le i,k\le n$, ${\mathrm{\partial }}_i{z}_k=R^2({\delta }_{ik}\mathrm{/}{r}^2-2y_k{y}_i\mathrm{/}{r}^4)$ and $\mathrm{\Delta }{z}_k=-2(n-2)R^2{y}_k\mathrm{/}{r}^4$, together with the auxiliary identities ${\sum }_i{\mathrm{\partial }}_i{z}_k{\mathrm{\partial }}_i{z}_l=R^4{\delta }_{kl}\mathrm{/}{r}^4$ and ${\sum }_i{y}_i{\mathrm{\partial }}_i{z}_l=-R^2{y}_l\mathrm{/}{r}^2$; every occurrence of $r$ is positive on $V$.

Put $w_0(x):=r^{2-n}u(z)$, so that $K_{R}u=R^{n-2}{w}_0$. The product and chain rules give ${\mathrm{\partial }}_i{w}_0=(2-n)r^{-n}{y}_{i}u(z)+r^{2-n}{\sum }_k{\mathrm{\partial }}_{k}u(z){\mathrm{\partial }}_i{z}_k$.

The chain and power rules give ${\mathrm{\partial }}_i(r^{2-n})=(2-n)r^{-n}{y}_i$ and ${\mathrm{\partial }}_i^2(r^{2-n})=(2-n)(r^{-n}-n{r}^{-n-2}{y}_i^2)$; summing and using $r^2=\mathrm{\mid }y{\mathrm{\mid }}^2$ gives $\mathrm{\Delta }(r^{2-n})=(2-n)(n{r}^{-n}-n{r}^{-n-2}\mathrm{\mid }y{\mathrm{\mid }}^2)=0$. The Laplacian product rule $\mathrm{\Delta }(fg)=f\mathrm{\Delta }g+2\mathrm{\nabla }f\cdot \mathrm{\nabla }g+g\mathrm{\Delta }f$ applied to $f=r^{2-n}$ and $g=u\circ z$ therefore gives $\mathrm{\Delta }{w}_0=r^{2-n}\mathrm{\Delta }(u\circ z)+2\mathrm{\nabla }(r^{2-n})\cdot \mathrm{\nabla }(u\circ z)$.

Two chain-rule evaluations. First, $\mathrm{\Delta }(u\circ z)={\sum }_{k,l}{\mathrm{\partial }}_k{\mathrm{\partial }}_{l}u(z){\sum }_i{\mathrm{\partial }}_i{z}_k{\mathrm{\partial }}_i{z}_l+{\sum }_k{\mathrm{\partial }}_{k}u(z)\mathrm{\Delta }{z}_k=R^4\mathrm{\Delta }u(z)\mathrm{/}{r}^4-2(n-2)R^{2}y\cdot \mathrm{\nabla }u(z)\mathrm{/}{r}^4$ by step 2.1. Second, since $\mathrm{\nabla }(r^{2-n})=(2-n)r^{-n}y$, step 2.1 gives $\mathrm{\nabla }(r^{2-n})\cdot \mathrm{\nabla }(u\circ z)=(2-n)r^{-n}{\sum }_l{\mathrm{\partial }}_{l}u(z){\sum }_i{y}_i{\mathrm{\partial }}_i{z}_l=(n-2)R^{2}y\cdot \mathrm{\nabla }u(z)\mathrm{/}{r}^{n+2}$.

Substituting step 3.2 into step 3.1, the first-order terms $-2(n-2)R^2{r}^{-n-2}y\cdot \mathrm{\nabla }u(z)$ and $+2(n-2)R^2{r}^{-n-2}y\cdot \mathrm{\nabla }u(z)$ cancel, leaving $\mathrm{\Delta }{w}_0=R^4\mathrm{\Delta }u(z)\mathrm{/}{r}^{n+2}$.

Restoring the factor of step 2.2 gives $\mathrm{\Delta }(K_{R}u)(x)=R^{n-2}\mathrm{\Delta }{w}_0(x)=R^{n+2}{r}^{-n-2}(\mathrm{\Delta }u)(I_R(x))=(R\mathrm{/}\mathrm{\mid }x-a\mathrm{\mid }{)}^{n+2}(\mathrm{\Delta }u)(I_R(x))$ for every $x\in I_R^{-1}(U)$.

Since $R^{n+2}{r}^{-n-2}>0$ on $V$, step 5.1 shows that $\mathrm{\Delta }u$ vanishes at $z=I_R(x)$ exactly when $\mathrm{\Delta }(K_{R}u)$ vanishes at $x$; both sides are evaluated only at points with $x\mathrm{\ne }a$, and $I_R$ is an involution exchanging the two punctured domains, so inversion transfers harmonicity in both directions [F1]. No choice principle and no measure-theoretic input is used.