Put $r_d(x)=\|x\|=(\sum_{i=1}^dx_i^2)^{1/2}$. The Laplacian of $1/r_d$ in $\mathbb{R}^d$ is given by $$\Delta(\frac{1}{r_d})=-\frac{d-3}{r_d^3}$$ as a direct calculation shows. Thus, $1/\|x\|$ is harmonic in $\mathbb{R}^d\backslash\{0\}$, if and only if $d=3$.
The Laplacian of a function measures how much does the value of the function at some point differ from the average of its values on the boundary of a small ball around the point. Is there a geometric argument that shows why the average of $1/r_d$ on the boundary of a small ball around a point in $3$-space equals the value at the center, as opposed to all the other-dimensional cases? In the special cases of $d=2,d=3$, is there a physical argument? I can integrate $1/r$ on a boundary of a ball, but I would like to be able to explain this in other ways, e.g, geometrical, or physical.
Here is a geometric (maybe physical?) reason why $d=3$ is the unique possible dimension for which $\frac{1}{r_d}$ can be harmonic. It does not prove that it is harmonic though.
The gradient of $\frac{1}{r_d}$ is $$ \nabla\left(\frac{1}{r_d}\right) = -\frac{1}{r_d^2}\nabla r_d. $$ Moreover, $\nabla r_d$ is the unit normal to concentric spheres. Let us fix $0<r<R$. By the divergence Theorem, \begin{align} \int_{r\leqslant\|x\|\leqslant R} \Delta\left(\frac{1}{r_d}\right) &= \int_{\|x\|=R} \nabla\left(\frac{1}{r_d} \right)\cdot \nabla r_d - \int_{\|x\|=r} \nabla\left(\frac{1}{r_d} \right)\cdot \nabla r_d\\ &= \int_{\|x\|=R} -\frac{1}{r_d^2} \underbrace{\nabla r_d \cdot \nabla r_d}_{=1} + \int_{\|x\|=r} \frac{1}{r_d^2} \underbrace{\nabla r_d \cdot \nabla r_d}_{=1}\\ &= -\frac{1}{R^2} \mathrm{area}(S(0,R)) + \frac{1}{r^2}\mathrm{area}(S(0,r)). \end{align} It can be easily shown that there exist a constant $\omega_{d-1}>0$ depending only on the dimension such that $$ \mathrm{area}(S(0,t)) = \omega_{d-1}t^{d-1}. $$ Hence, $$ \int_{r\leqslant\|x\|\leqslant R} \Delta\left(\frac{1}{r_d}\right) = \omega_{d-1}\left(r^{d-3} - R^{d-3} \right). $$ From there, if $d\neq 3$, $\int_{r\leqslant \|x\|\leqslant R} \Delta\left(\frac{1}{r_d}\right) \neq 0$, and $\Delta\left(\frac{1}{r_d}\right)$ cannot identically vanish. The function $\frac{1}{r_d}$ can then not be harmonic in dimension $d\neq 3$.
In physical terms, we have shown that the flux of $\nabla\left(\frac{1}{r_d}\right)$ is not preserved along concentric spheres. Therefore, $\nabla\left(\frac{1}{r_d}\right)$ cannot be divergence-free, and $\frac{1}{r_d}$ cannot be harmonic.
The same argument shows that $\frac{1}{r_d^{k-2}}$ cannot be harmonic if $d\neq k$.