Asymptotic behaviour of $I_{\alpha}(x):=\int_{\mathbb{R}^{n}}\frac{e^{-|x-y|^2}}{|y|^{\alpha}}dy$ with $0<\alpha<1$ as $\alpha\rightarrow 0$ I am trying to obtain the asymptotic behaviour of the integral
$$I_{\alpha}(x):= \int_{\mathbb{R}^{n}}\frac{e^{-|x-y|^2}}{|y|^{\alpha}}dy$$ explicitly as $\alpha\rightarrow 0^{+}$.
Clearly, by the dominated convergence theorem
$$I_{\alpha}(x)\rightarrow
\int_{\mathbb{R}^{n}}e^{-|x-y|^2}dy=\int_{\mathbb{R}^{n}}e^{-|y|^2}dy=c.$$
My naive attempt is  to calculate $I_{x}$ explicitly in the hope to get an asymptotic. Here is what I got:
Since $|x-y|^2=|x|^2-2 x\cdot y+|y|^2$ then
$$I_{\alpha}(x)= \int_{\mathbb{R}^{n}}\frac{e^{-|x-y|^2}}{|y|^{\alpha}}dy=e^{-|x|^2}\int_{\mathbb{R}^{n}}\frac{e^{-2 x\cdot y+|y|^2}}{|y|^{\alpha}}dy.$$ Using spherical coordinates, we get
$$I_{\alpha}(x)= e^{-|x|^2}\int_{\mathbb{S}^{n-1}}
\int_{0}^{\infty}
e^{-2 x\cdot \omega r+r^2}r^{n-1-\alpha}dr d\sigma(\omega).$$
By Fubini's theorem we have
$$I_{\alpha}(x)= e^{-|x|^2}
\int_{0}^{\infty}\int_{\mathbb{S}^{n-1}}
e^{-2 x\cdot \omega r} d\sigma(\omega)
e^{-r^2}r^{n-1-\alpha}dr.$$
Using the formula (Grafakos, Classical Fourier Analysis, Appendix D)
$$C\int_{\mathbb{S}^{n-1}} F(x\cdot \omega)d\sigma(\omega)= \int_{-1}^{1}F(s|x|)(\sqrt{1-s^2})^{n-3} ds$$
we have
$$I_{\alpha}(x)= e^{-|x|^2}
\int_{0}^{\infty}\int_{-1}^{1}
e^{-2 s|x| r}(\sqrt{1-s^2})^{n-3} ds 
e^{-r^2}r^{n-1-\alpha}dr.$$
Applying Fubini's theorem one more time
$$I_{\alpha}(x)= e^{-|x|^2}
\int_{-1}^{1}\int_{0}^{\infty}
e^{-2 s|x| r-r^2}r^{n-1-\alpha}dr(\sqrt{1-s^2})^{n-3} ds.$$
Is there a way to calculate
$$\int_{0}^{\infty}
e^{-2 s|x| r-r^2}r^{n-1-\alpha}dr$$
Mathematica gives some kind of the hypergeometric function. It is difficult to understand how Mathematica's answer behaves in terms of $\alpha$.
 A: Expanding the factor $|y|^{-\alpha} = \exp(-\alpha\log|y|)$ in $\alpha$ using the Taylor series of t$\exp(\cdot)$ and invoking the Fubini's Theorem, we may expand
$$ I_{\alpha}(x) = \sum_{k=0}^{\infty} \frac{(-\alpha)^k}{k!} \underbrace{\int_{\mathbb{R}^n} e^{-|x-y|^2} (\log|y|)^k \, \mathrm{d}y}_{=:J_k(x)}. $$
Now applying the computation similarly as in OP, we get
\begin{align*}
J_k(x)
&= \int_{\mathbb{R}^n} e^{-|x-y|^2} (\log|y|)^k \, \mathrm{d}y \\
&= e^{-|x|^2}\int_{0}^{\infty} \biggl( \frac{2\pi^{(n-1)/2}}{\Gamma(\frac{n-1}{2})} \int_{-1}^{1} (1-s^2)^{\frac{n-3}{2}} e^{-2|x|rs} \, \mathrm{d}s \biggr) e^{-r^2} r^{n-1}(\log r)^k \, \mathrm{d}r \\
&= e^{-|x|^2}\int_{0}^{\infty} \biggl( \frac{2\pi^{(n-1)/2}}{\Gamma(\frac{n-1}{2})} \int_{0}^{1} t^{-\frac{1}{2}}(1-t)^{\frac{n-3}{2}} \cosh(2|x|r\sqrt{t}) \, \mathrm{d}t \biggr) e^{-r^2} r^{n-1}(\log r)^k \, \mathrm{d}r \\
&= e^{-|x|^2} \sum_{j=0}^{\infty} \frac{(2|x|)^{2j}}{(2j)!} \biggl( \frac{2\pi^{(n-1)/2}}{\Gamma(\frac{n-1}{2})} \int_{0}^{1} t^{j-\frac{1}{2}}(1-t)^{\frac{n-3}{2}} \, \mathrm{d}t \biggr) \int_{0}^{\infty}  e^{-r^2} r^{2j+n-1}(\log r)^k \, \mathrm{d}r \\
&= e^{-|x|^2} 2\pi^{(n-1)/2} \sum_{j=0}^{\infty} \frac{(2|x|)^{2j}}{(2j)!} \frac{\Gamma(j+\frac{1}{2})}{\Gamma(j+\frac{n}{2})} \int_{0}^{\infty} e^{-r^2} r^{2j+n-1}(\log r)^k \, \mathrm{d}r \\
&= \frac{\pi^{n/2}}{2^k} e^{-|x|^2}  \sum_{j=0}^{\infty} \biggl( \frac{1}{\Gamma(j+\frac{n}{2})} \int_{0}^{\infty} u^{j+\frac{n}{2}-1}e^{-u} (\log u)^k \, \mathrm{d}u \biggr) \frac{|x|^{2j}}{j!} \\
&= \frac{\pi^{n/2}}{2^k} e^{-|x|^2}  \sum_{j=0}^{\infty} \frac{\Gamma^{(k)}(j+\frac{n}{2})}{\Gamma(j+\frac{n}{2})} \frac{|x|^{2j}}{j!} .
\end{align*}
In the intermediate steps, we utilized the substitutions $t = s^2$ and $u = r^2$. Also, we applied the beta function identity to compute the inner integral. Also, $\Gamma^{(k)}(\cdot)$ denotes the $k$-th derivative of the gamam function.
As a sanity check, let us plug $k=0$. Then
\begin{align*}
J_0(x)
= \pi^{n/2} e^{-|x|^2}  \sum_{j=0}^{\infty} \frac{|x|^{2j}}{j!}
= \pi^{n/2},
\end{align*}
which is the same as the $I_{0}(x)$ as expected. Then $J_1(x)$ is
\begin{align*}
J_1(x)
= \frac{\pi^{n/2}}{2} e^{-|x|^2} \sum_{j=0}^{\infty} \psi^{(0)}\left(j+\frac{n}{2}\right) \frac{|x|^{2j}}{j!},
\end{align*}
where $\psi(s) = \frac{\mathrm{d}}{\mathrm{d} s} \log\Gamma(s)$ is the digamma function. In general, each $J_k(x)$ can be written in a similar form involving polygamma functions.
Summarizing, the above formula shows that $I_{\alpha}(x)$ can be expanded as a power series in $\alpha$ and then each coefficient $\frac{(-1)^k}{k!} J_k(x)$ can be further expanded as a power series in $|x|$.
