Yukawa potential in higher dimensions I am trying to calculate the integral
\begin{align}
E_n(\mathbf{r}) = \int \frac{d^n \mathbf{k}}{(2\pi)^n} \frac{ e^{i\mathbf{k}\cdot\mathbf{r}} }{ \mathbf{k}^2 + m^2 }
\end{align}
for $n > 2$ (the case $n=3$ is well-known and calculated here).
Related is this post where it is stated that for $n \neq 2$ this should be proportional to $e^{-m r}/{r^{n-2}}$. This intuitively makes sense (as a screened version of the $n$ dimensional Coulomb potential). However, when I perform the integrations, I don't seem to be able to produce this result. One thing to note is that, similar to the Fourier transform of a Coulomb potential (post), this integral may formally diverge, but one can still obtain finite results by regularizing it.
Trying to generalize the steps outlined in the Wiki page, let's use the $n$-dimensional spherical coordinates, where from the $n-1$ angular variables we only need to perform the integration over one because of the symmetry
\begin{align}
E_n (\mathbf{r} ) &= \frac{S_{n-1}}{(2\pi)^n} \int_0^\infty k^{n-1} dk \int_0^\pi \sin^{n-2}(\theta) \, d\theta \, \frac{ e^{i k r \cos(\theta) }}{k^2 + m^2} 
\end{align}
Wolfram says the $\theta$ integral gives a bunch of Bessel functions, but I don't know how to proceed with that.
Update
A different approach is to represent the denominator as an exponential integral
\begin{align}
E_n &= \int \frac{d^n k}{(2\pi)^n} e^{i\mathbf{k}\cdot\mathbf{r}} 
\int_0^\infty d\alpha \, e^{-\alpha ( \mathbf{k}^2 + m^2 )}\\
&= \int_0^\infty d\alpha \, e^{-\alpha m^2} \, \frac{ e^{-\frac{\mathbf{r}^2}{4\alpha}}}{(4\pi\alpha)^{n/2}} \\
&= \frac{1}{(2\pi)^{n/2}} \left( \frac{m}{r} \right)^{\frac{n-2}{2}} \,
K_{\frac{n-2}{2}}(m r),
\end{align}
where from the first to the second line we perform the Gaussian integration and $K$ in the third line is the modified Bessel function of the second kind.
One can now expand the Bessel function for $mr\ll 1$ or $mr\gg 1$; for small distances, substituting the expansion gives a power-law decaying as $r^{-n+2}$, and for large distances, it gives an exponential decay $e^{-mr}$. Although this is very close, but I still can't seem to get the claimed behavior $e^{-mr}/{r^{n-2}}$.
 A: Let us start from the integral
$$I=\int_0^{\infty}d\alpha\,\alpha^{-n/2}\exp\left(-\frac{\alpha m^2}{2}-\frac{r^2}{2\alpha}\right).$$
We are interested in the limit of large $r$, so the integral saturates around one point. In order to find this point, we write
$$\alpha^{-n/2}=\exp\left(-\frac{n}{2}\ln \alpha\right),$$
so we have the integral
$$\int_{0}^{\infty}d\alpha\,e^{+f(\alpha)},\,\,f(\alpha)=-\frac{n}{2}\ln\alpha-\frac{\alpha m^2}{2}-\frac{r^2}{2\alpha}.$$
As was mentioned, the integral saturates near the point $f'(\alpha_0)=0$, which is (in the limit of large $r$) $\alpha_0=r/m$. Then,
$$f(\alpha_0)=-mr-\frac{n}{2}\ln(r/m),\quad f''(\alpha_0)=-m^3/r.$$
The value of the integral $I$ is
$$I\approx \sqrt{\frac{2\pi}{|f''(\alpha_0)|}}e^{f(\alpha_0)},$$
which is
$$\boxed{I\sim \frac{e^{-mr}}{r^{(n-1)/2}}.}$$
It seems that the final expression correctly reproduces screened Coulomb law in 3D,
$$n=3\rightarrow I\sim\frac{e^{-mr}}{r}$$
Additional question is what about case of $m=0$? In case of $m=0$, we have
$$\alpha_0=\frac{r^2}{n},\quad f''(\alpha_0)=-\frac{n^3}{2r^4},$$
so
$$\left.I\right|_{m=0}\sim \frac{1}{r^{n-2}}$$
and coincides with usual Coulomb law.
