Laplacian in 4-dimensions I know for 3-D $$\nabla^2 \left(\frac1r\right)=-4\pi\, \delta(\vec{r})\,.$$
I would like to know, what is $$\text{Div}\cdot\text{Grad}\left(\frac{1}{r^2}\right)$$ in 4-Dimensions ($r^2=x_1^2+x_2^2+x_3^2+x_4^2$)?
 A: Let $\sigma_k$ denote the hypersurface area measure on the $k$-sphere.  For example,
$$\text{d}\sigma_1(\varphi_1)=\text{d}\phi_1\,,\,\,\text{d}\sigma_2(\varphi_1,\varphi_2)=\sin(\varphi_2)\,\text{d}\varphi_1\,\text{d}\varphi_2\,,$$
and 
$$\text{d}\sigma_3(\varphi_1,\varphi_2,\varphi_3)=\sin(\varphi_2)\,\sin^2(\varphi_3)\,\text{d}\varphi_1\,\text{d}\varphi_2\,\text{d}\varphi_3\,.$$
In general,
$$\text{d}\sigma_k(\varphi_1,\varphi_2,\ldots,\varphi_k)=\prod_{j=1}^k\,\left(\sin^{j-1}(\varphi_j)\,\text{d}\varphi_j\right)\text{ for all }k=1,2,3,\ldots\,.$$
For simplicity, write $\Phi_k$ for the angular tuple $\left(\varphi_1,\varphi_2,\ldots,\varphi_k\right)$. 
The volume element in the polar coordinates of $\mathbb{R}^n$ is given by $$\text{d}\lambda_n(\mathbf{x})=r^{n-1}\,\text{d}r\,\text{d}\sigma_{n-1}\left(\Phi_{n-1}\right)\,,$$
if $\mathbf{x}\in\mathbb{R}^n$ is represented by the polar coordinates $(r,\varphi_1,\varphi_2,\ldots,\varphi_{n-1})$.  Define
$$\Psi_n(\mathbf{x}):=\frac{1}{\|\mathbf{x}\|_2^{n-2}}\text{ for all }\mathbf{x}\in\mathbb{R}^n\text{ for }n>2\,,$$
where $\|\_\|_2$ is the usual Euclidean norm on $\mathbb{R}^n$.
For a differentiable function $f:\Omega\to\mathbb{R}$, where $\Omega$ is an open region in $\mathbb{R}^n$, the gradient $\boldsymbol{\nabla}f:\Omega\to\mathbb{R}^n$ is given by
$$(\boldsymbol{\nabla}f)(\mathbf{x}):=\sum_{j=1}^n\,\left(\frac{\partial f}{\partial x_j}(\mathbf{x})\right)\,\mathbf{e}_j\text{ for all }\mathbf{x}\in\Omega\,,$$
where $\mathbf{x}=\left(x_1,x_2,\ldots,x_n\right)$ and $\mathbf{e}_1,\mathbf{e}_2,\ldots,\mathbf{e}_n$ are the usual standard basis vectors of $\mathbb{R}^n$.  For a differentiable vector function $\mathbf{v}:\Omega\to\mathbb{R}^n$, the divergence $\boldsymbol{\nabla}\cdot\mathbf{v}:\Omega\to\mathbb{R}$ is the function
$$(\boldsymbol{\nabla}\cdot\mathbf{v})(\mathbf{x}):=\sum_{j=1}^n\,\left(\frac{\partial v_j}{\partial x_j}(\mathbf{x})\right)\,,$$
where $\mathbf{v}=\left(v_1,v_2,\ldots,v_n\right)$.  The Laplacian operator $\nabla^2$ is defined by
$$\nabla^2f:=\boldsymbol{\nabla}\cdot(\boldsymbol{\nabla}f)\,.$$
We want to evaluate
$$L(f):=\int_{\mathbb{R}^n}\,\Psi_n(\mathbf{x})\,(\nabla^2f)(\mathbf{x})\,\text{d}\lambda_n(\mathbf{x})$$
for a well behaved function $f$ (say, $f$ is sufficiently fast decaying at large distances from the origin $\boldsymbol{0}_n$ of $\mathbb{R}^n$, at least equipped with second weak derivatives, and with bounded first weak derivatives).  We can see that the distribution $\nabla^2\Psi_n$ satisfies
$$\int_{\mathbb{R}^n}\,f(\mathbf{x})\,(\nabla^2\Psi_n)(\mathbf{x})\,\text{d}\lambda_n(\mathbf{x})=(-1)^2\,\int_{\mathbb{R}^n}\,\Psi_n(\mathbf{x})\,(\nabla^2f)(\mathbf{x})\,\text{d}\lambda_n(\mathbf{x})=L(f)\,.$$
(The equality above is where the assumption that $f$ is fast decaying at large distances comes into play.)  Note that
$$\left(\boldsymbol{\nabla}\Psi_n\right)(\mathbf{x}) = -(n-2)\,\left(\frac{\mathbf{x}}{\|\mathbf{x}\|_2^{n}}\right)\text{ for all }\mathbf{x}\neq \boldsymbol{0}_n\,,$$
and so
$$(\nabla^2\Psi_n)(\mathbf{x})=-(n-2)\,\left(\frac{n}{\|\mathbf{x}\|_2^n}-n\,\sum_{j=1}^n\,\frac{x_j^2}{\|\mathbf{x}\|_2^{n+2}}\right)=0\text{ for every }\mathbf{x}\neq \boldsymbol{0}_n\,.$$ 
That is, 
$$L(f)=\lim_{R\to0^+}\,\int_{\mathbb{R}\setminus B_R(\boldsymbol{0}_n)}\,\Psi_n(\mathbf{x})\,\big(\nabla^2 f\big)(\mathbf{x})\,\text{d}\lambda_n(\mathbf{x})\,,\tag{*}$$
where $B_\rho(\mathbf{y})$ is the open ball centered at $\mathbf{y}\in\mathbb{R}^n$ with radius $\rho>0$.
We can write
$$f(\mathbf{x})\,(\nabla^2\Psi_n)(\mathbf{x})=\big(\boldsymbol{\nabla}\cdot(f\,\boldsymbol{\nabla}\Psi_n)\big)(\mathbf{x})-\big((\boldsymbol{\nabla}f)(\mathbf{x})\big)\cdot\big((\boldsymbol{\nabla}\Psi_n)(\mathbf{x})\big)$$
and
$$\Psi_n(\mathbf{x})\,(\nabla^2f)(\mathbf{x})=\big(\boldsymbol{\nabla}\cdot(\Psi_n\,\boldsymbol{\nabla}f)\big)(\mathbf{x})-\big((\boldsymbol{\nabla}\Psi_n)(\mathbf{x})\big)\cdot\big((\boldsymbol{\nabla}f)(\mathbf{x})\big)\,.$$
Thus,
$$f(\mathbf{x})\,(\nabla^2\Psi_n)(\mathbf{x})=\Psi_n(\mathbf{x})\,(\nabla^2f)(\mathbf{x})+(\boldsymbol{\nabla}\cdot \mathbf{v})(\mathbf{x})\,,$$
where
$$\mathbf{v}:=f\,(\boldsymbol{\nabla}\Psi_n)-\Psi_n\,(\boldsymbol{\nabla}f)\,.$$
For $\mathbf{x}\neq \boldsymbol{0}_n$, we obtain
$$\Psi_n(\mathbf{x})\,(\nabla^2f)=-(\boldsymbol{\nabla}\cdot \mathbf{v})(\mathbf{x})\,.$$  From (*), we get
$$L(f)=-\lim_{R\to0^+}\,\int_{\mathbb{R}^n\setminus B_R(\boldsymbol{0}_n)}\,(\boldsymbol{\nabla}\cdot\mathbf{v})(\mathbf{x})\,\text{d}\lambda_n(\mathbf{x})\,.$$
Using the Divergence Theorem, we obtain
$$L(f)=\lim_{R\to 0^+}\,\int_{\partial B_R(\boldsymbol{0}_n)}\,\mathbf{v}(\mathbf{x})\cdot\mathbf{x}\,\|\mathbf{x}\|_2^{n-2}\,\text{d}\sigma_{n-1}(\Phi_{n-1})\,.$$
That is,
$$\begin{align}L(f)&=\lim_{R\to 0^+}\,\int_{\partial B_R(\boldsymbol{0}_n)}\,\Big(f(\mathbf{x})\,(\boldsymbol{\nabla}\Psi_n)(\mathbf{x})-\Psi_n(\mathbf{x})\,(\boldsymbol{\nabla}f)(\mathbf{x})\Big)\cdot\mathbf{x}\,\|\mathbf{x}\|_2^{n-2}\,\text{d}\sigma_{n-1}(\Phi_{n-1})\\
&=\lim_{R\to0^+}\,\left(\int_{\partial B_R(\boldsymbol{0}_n)}\,(-n+2)\,f(\mathbf{x})\,\text{d}\sigma_{n-1}(\Phi_{n-1})-\int_{\partial B_R(\boldsymbol{0}_n)}\,\mathbf{x}\cdot(\boldsymbol{\nabla}f)(\mathbf{x})\,\text{d}\sigma_{n-1}(\Phi_{n-1})\right)
\\
&=-(n-2)\,f(\boldsymbol{0}_n)\,\int_{\partial B_1(\boldsymbol{0}_n)}\,\text{d}\sigma_{n-1}(\Phi_{n-1})-0=-(n-2)\,\Sigma_{n-1}\,f(\boldsymbol{0}_n)\,,\end{align}$$
where $\Sigma_k$ denotes the hypersurface measure of the unit $k$-sphere.  (For example, $\Sigma_1=2\pi$, $\Sigma_2=4\pi$, $\Sigma_3=2\pi^2$, and for every $k=1,2,3,\ldots$, $\Sigma_k=\dfrac{2\pi^{\frac{k+1}{2}}}{\Gamma\left(\frac{k+1}{2}\right)}$, where $\Gamma$ is the gamma function.)
In conclusion,
$$\left(\nabla^2\Psi_n\right)(\mathbf{x})=-\frac{2(n-2)\pi^{\frac{n}{2}}}{\Gamma\left(\frac{n}{2}\right)}\,\delta_n(\mathbf{x})\,,\tag{#}$$
where $\delta_n$ is the $n$-dimensional Dirac delta distribution.  In particular,
$$\left(\nabla^2\Psi_3\right)(\mathbf{x})=-4\pi\,\delta_3(\mathbf{x})\text{ and }\left(\nabla^2\Psi_4\right)(\mathbf{x})=-4\pi^2\,\delta_4(\mathbf{x})\,.$$
In fact, it also makes sense to consider $\Psi_1$, where $\Psi_1(x)=|x|$ for all $x\in\mathbb{R}$.  The same formula (#) works and we can readily check that 
$$\left(\nabla^2\Psi_1\right)(x)=2\,\delta_1(x)\text{ for all }x\in\mathbb{R}\,.$$  For $\Psi_2$, (#) also works trivially, noting that $\Psi_2\equiv 1$ almost everywhere, whence $\nabla^2\Psi_2\equiv 0$ almost everywhere.
If you want to obtain a similar result in the $2$-dimensional case, then you can take $$\Xi(\mathbf{x}):=\ln\big(\|\mathbf{x}\|_2\big)\text{ for all }\mathbf{x}\in\mathbb{R}^2\setminus\{\boldsymbol{0}_2\}\,.$$
Then,
$$L(f):=\int_{\mathbb{R}^2}\,\Xi(\mathbf{x})\,(\nabla^2f)(\mathbf{x})\,\text{d}\lambda_2(\mathbf{x})=\int_{\mathbb{R}^2}\,f(\mathbf{x})\,(\nabla^2\Xi)(\mathbf{x})\,\text{d}\lambda_2(\mathbf{x})$$
for all well behaved functions $f:\mathbb{R}^2\to\mathbb{R}$.  Note that
$$(\boldsymbol{\nabla}\Xi)(\mathbf{x})=\frac{\mathbf{x}}{\|\mathbf{x}\|_2^2}\text{ and }(\nabla^2\Xi)(\mathbf{x})=0\text{ for all }\mathbf{x}\neq \boldsymbol{0}_2\,.$$
We perform the same trick as before to get
$$\begin{align}L(f)&=\lim_{R\to0^+}\,\int_{\partial B_R(\boldsymbol{0}_2)}\,\Big(f(\mathbf{x})\,(\boldsymbol{\nabla}\Xi)(\mathbf{x})-\Xi(\mathbf{x})\,(\boldsymbol{\nabla}f)(\mathbf{x})\Big)\cdot\mathbf{x}\,\text{d}\sigma_1(\Phi_1)
\\&=f(\boldsymbol{0}_2)\,\int_{\partial B_1(\boldsymbol{0}_2)}\,\text{d}\sigma_1(\Phi_1)-0=\Sigma_1\,f(\boldsymbol{0}_2)=2\pi\,f(\boldsymbol{0}_2)\,.
\end{align}$$
That is,
$$(\nabla^2\Xi)(\mathbf{x})=2\pi\,\delta_2(\mathbf{x})\,.$$
A: Since the Fourier transform of a radial function is also a radial function and the transform of $r^\lambda$ is a homogeneous function of degree $-\lambda - n$ for $(-\lambda - n)/2 \notin \mathbb N^0$,
$$\mathcal F[r^\lambda] =
(r^\lambda, e^{i \boldsymbol x \cdot \boldsymbol \xi}) =
C_{\lambda, n} \rho^{-\lambda - n}, \\
\mathcal F[r^{2 - n}] = C_n \rho^{-2}, \\
\mathcal F[\nabla^2 r^{2 - n}] = -\rho^2 \mathcal F[r^{2 - n}] = -C_n, \\
\nabla^2 r^{2 - n} = -C_n \delta(\boldsymbol x).$$
The constant $C_n = 4 \pi^{n/2}/\Gamma(n/2 - 1)$ can be found by taking a simple test function and using the identity $(r^\lambda, \mathcal F[\phi]) = (\mathcal F[r^\lambda], \phi)$.
