So, it is a famous fact that if $u:\mathbb{R}^n \to \mathbb{R}$ is an harmonic function, then its Kelvin transform $$ (Ku)(x) := \frac{1}{|x|^{n-2}} u\left(\frac{x}{|x|^2} \right) $$ is harmonic too. Apparently, all the proofs of this fact that I have seen are some variation of "do the computations and all the terms simplify". Moreover, I already did the computation once and so I am not interested in seeing the computations done in some other equivalent form.
My point is the following: why, really, is this function harmonic? I don't buy the fact that the first person that discovered this just randomly did the computations with the correct power of the norm in front. For $n=2$ this can be proved very easily by properties of holomorphic functions and indeed there's a geometric motivation, in this case $n=2$, on why it sends harmonic functions to harmonic functions. But what for general $n$? Does anyone have at least some geometric/physical intuition on why $Ku$ is harmonic if $u$ is?
I tried many characterizations of harmonic functions like the mean value property or integrating against test functions, and all seem not to bring the result in a clean way. Indeed, at some point, there's always a big computation (being that some tangential Jacobian or Laplacian of composition) that I find basically equivalent to doing the computation from the beginning.
I am searching for a proof or at least some motivation on why would someone think that this transformation is the natural candidate to send harmonic functions to harmonic functions.